Formation Flight Computational Fluid Dynamics

     This is a post about computational fluid dynamics analysis of formation flight.

     The results from an analysis of three unmanned combat aerial vehicles (UCAVs) flying in a V-type formation are presented. The chosen UCAV configuration, shown in Figs. 1-2 and available for download here, is named SACCON (Stability And Control Configuration) UCAV. This configuration is used because of the availability of the geometric and aerodynamic data, used in the validation and verification of the numerical analysis. The SACCON UCAV is designed by NATO's (North Atlantic Treaty Organization) RTO (Research and Technology Group) under Applied Vehicle Task Group (AVT-161) to assess the performance of military aircrafts.

Fig. 1, SACCON UCAV

Fig. 2, Technical drawing for the SACCON UCAV

     The simulation is performed using commercially available computational fluid dynamics code.  The details about the solver and the discretization schemes are presented next. The simulation is performed using SIMPLE-R solver for pressure-velocity coupling. The diffusive terms of the Navier-Stokes equations are discretized using central differentiating scheme while the convective terms are discretized using the upwind scheme of second order. The κ-ε turbulence model with damping functions is implemented to model turbulence. The simulation predicts three-dimensional steady–state flow over UCAVs.

     The Reynolds and Mach number of flow are set at 1e6 and 0.15, respectively. The two trailing UCAVs are placed 3 wingspans behind the leading UCAV. The trailing UCAVs are a wingspan apart. The V-type configuration is chosen because it is the most common observation in birds (author's observation). The V-type formation is shown in Fig. 3. All three UCAVs are at 5° angle-of-attack.

Fig. 3, The V-type formation (top view)

     The boundaries of the computational domain are located at a distance equivalent to 10 times the distance between the nose of the leading UCAV and the tail of the trailing UCAV. The mesh is made of 821,315 cells. Mesh controls are used to refine the mesh in the areas of interest i.e. on the surfaces of the UCAVs and in the wake of all three UCAVs. A cartesian mesh with immersed boundary method is used for the present study. The computational domain with mesh is shown in Fig. 4 while a closeup of the mesh is shown in Fig. 5.

Fig. 4, The computational domain and mesh

Fig. 5, Closeup of mesh, notice the refined wakes of the UCAVs

     For validation and verification, the lift and drag forces from the present study are compared with studies [1-2]. The results are in close agreement with [1,2] As a result of flying in a formation, an improvement in the lift-to-drag ratio of 10.05% is noted. The lift-to-drag ratio of the trailing UCAVs is at 11.825 in comparison with a lift-to-drag ratio of a single UCAV, i.e. 10.745. The lift coefficient is increased by 7.43% while the drag coefficient decreased by 2.174%. The reason(s) to why the efficiency increases will be looked upon later, if ever the author has the time and will power ☺.

     The results from post processing of the simulations are presented in Figs. 6-7. The pressure iso-surfaces colored by velocity magnitude are shown in Fig. 6. While the velocity iso-surfaces colored by pressure magnitude are shown in Fig. 7.

Fig. 6, Flight direction towards the reader

Fig. 7, Flight direction away from the reader

     Thank you very much for reading. If you would like to collaborate on research projects, please reach out.

[1] https://doi.org/10.2514/1.C031386

[2] https://doi.org/10.1155/2017/4217217

Hypersonic Flow over a Two Dimensional Heated Cylinder

     This post is about the simulation of hypersonic flow over a heated circular cylinder, in two dimensions.

     Equation 1 is used as a relationship between Mach and the Reynold number.

M= Re*μ*√(R*T) ÷ d*P*√γ     (1)

     w.r.t. equation 1, the parameters represent the following quantities.

     M     Freestream Mach number at 17.6
     Re    Reynolds number at 376,000
     μ     Dynamic viscosity at 1.329045e^-5 Ns.m^-2
     R     Specific gas constant at 286.9 J.(kg.K)^-1
     T     Freestream temperature 200 K
     d     Cylinder diameter at 5.6730225e^-4 m
     P     Freestream pressure at 101325 Pa
     γ     Specific heat ratio at 1.4
     Tw  Wall temperature of cylinder at 500 K
     Pr    Prandtl number at 0.736

     The boundary conditions were taken from [1]. A comparison with [1] is shown in Fig. 1. Inside Fig. 1, the red dotted line with circles represents the data from [1]. The black solid line represents the data from the present simulation. Within Fig. 1, 0° represents the stagnation point. The velocity, pressure, Mach number and temperature contours are shown in Fig. 2.

Fig. 1 A comparison with previous research [1].

Fig. 2, Top Row, L-R: Velocity and pressure contours. Bottom Row, L-R: Mach number and temperature contours.

The computational mesh and the computational domain with boundary conditions visible are shown in Fig. 3-4, respectively. The computational domain had a size of 20D x 20D. The mesh had 836,580 total cells and 944 cells were located at the solid fluid boundary. Several local mesh controls were employed to capture the shockwave properly.

Fig. 3, The computational mesh.

Fig. 4, The computational domain.

     The solution method is Finite Volume method. SIMPLE-R is the solver employed. Implicit central difference scheme for diffusion terms, second-order Upwind scheme for convective terms and first-order implicit for temporal terms are used. The mesh created uses the Cartesian mesh with Immersed Boundary method.

     Reference:

     [1] https://arc.aiaa.org/doi/10.2514/6.2004-2371 

     Thank you for reading. If you would like to collaborate on research projects, please reach out. I am looking for a PhD position, any guidance would be appreciated.

Vertical Axis Wind Turbine Computational Fluid Dynamics Analysis

     This post is be about the validation and verification of the computational fluid dynamics analysis of a three blade vertical axis wind turbine. The turbine had a diameter of 2 m with each blade being 1 m tall. The blades had an NACA-0018 airfoil cross section.

     The computational fluid dynamics analysis employed the κ-ε turbulence model with damping functions as the turbulence model, SIMPLE-R as the numerical algorithm. The spatial discretization schemes for the convective fluxes and diffusive terms used are the second order upwind and central approximations, respectively. An implicit first-order Euler scheme is employed to approximate the time derivatives.

     The Cartesian computational mesh with immersed boundary method had a total of 769,357 cells. Among those 769,357 cells, 166,188 cells were around the turbine blades. Mesh controls were employed to refine the mesh near the turbine blades. A time step of 3e^-3 was employed. The computational domain inlet was 1.5 D away from the turbine and the outlet was 3D away. The computational domain walls on the sides were 1D x 1.5D, where D represents the turbine diameter. The mesh and the computational domain are shown in Fig. 1. The vertical teal arrow represents the force of gravity, the curved teal arrow represents the direction of turbine rotation. The dark blue arrow represents the direction of free stream velocity.

Fig. 1, Mesh and computational domain.

     The simulations ran at a tip-speed ratio of 1.87 at a wind speed of 4.03 m.s^-1. The velocity distribution around the turbine after 4 revolutions is shown in Fig. 2. Validation of the numerical analysis was carried out using [1]. The results of power produced by the turbine were with in 4% of the experimental results [1]. An animation of the numerical analysis is also shown.

Fig. 1, Flow field around the turbine.

     Thank you for reading. If you would like to contribute to the research, both financially and scientifically, please feel free to reach out.





[1] Yi-Xin Peng, You-Lin Xu, Sheng Zhan and Kei-Man Shum, High-solidity straight-bladed vertical axis wind turbine: Aerodynamic force measurements, Journal of Wind Engineering and Industrial Aerodynamics, January 2019.

Improvement of the Volume Flow Rate Through a Blower Fan

     In this post, an improvement in the volume flow rate through the blower fan assembly made are presented. The only thing changed in the blower fan was the cross section of the fan blades. In the previous version, the blade cross section resembled a flat plate with fillets at the leading edge. The trailing edge in the previous design was blunt. In the modified design, there aero-foils were selected, namely NACA 9410, NACA 9420 and the NACA 9430. All the other parameters were kept the same to the previous case. The CAD models of the modified fan blades are shown in Fig. 1.

Fig. 1, Fan blade geometries.

     The velocity contours are shown in Fig. 2 while the pressure contours are shown in Fig. 3, super imposed with velocity vectors and the computational mesh. The volume flow rate was the most for the fan with blades having cross-section of NACA 9410 aero-foil, followed by the fan with blades having cross-section of NACA 9420 and the NACA 9430 cross sections, respectively.

Fig. 2, Pressure contours. Row 1, L-R; fan with the NACA 9430 and NACA 9420 cross sections. Row 2, fan with NACA 9410 cross sections.

Fig. 3, Velocity contours. Row 1, L-R; fan with the NACA 9430 and NACA 9420 cross sections. Row 2, fan with NACA 9410 cross sections.

     Thank you for reading. If you would like to collaborate, both scientifically and financially, on research projects, please reach out.

Computational Fluid Dynamics Analysis of a Blower/Centrifugal Fan: Update 01

     In this post the results from a CFD analysis of a blower fan are presented. The fan had a diameter of 66 mm and a height of 12.57 mm. The fan's rotational velocity was at 10,000 rpm. The CAD model is shown in Fig. 1.

Fig. 1, CAD Assembly of the Blower Fan.

     The simulations were completed in SolidWorks Flow Simulation Premium code. The code employs immersed boundary method to create a Cartesian mesh. The sliding mesh feature was employed to simulate the rotation of the fan at atmospheric conditions. The code employs κ-ε model with Two-Scales Wall Functions approach as the turbulence model. The numerical algorithm implemented is the SIMPLE-R, modified. The second-order upwind discretization scheme is used to approximate the convective fluxes while the diffusive terms are approximated using the central differencing scheme. The time derivatives are approximated with an implicit first-order Euler scheme.

     The numerical model for the fan had 816,994 cells of which 209,421 cells were at the solid-fluid interface. Two mesh controls were employed to refine the mesh near the blades of the fan and at the boundary of the stationery and the rotating domains. The results were indeed, mesh independent. Due to the fact that this was an internal flow problem, domain independence test was not applicable. The mesh and the computational domain is shown in Fig. 2. The curved teal arrow represents the direction of rotation of the fan. The blue arrows represent the pressure boundary conditions at the inlet and at the outlet of the fan assembly. The straight teal arrow represents the force of gravity (the arrow is inverted).

Fig. 2, The mesh and the computational domain.

     The pressure and velocity plots are shown in Fig. 3-4.

Fig. 3, Pressure contours.

Fig. 4, Velocity contour

     Thank you very much for reading. If you would like to collaborate on research projects or want a tutorial for the setup of the numerical simulations such as this one, please reach out.

Update 01

     CAD files are available here.

Computational Fluid Dynamics Analysis of a Symmetrical Wing, Update 01

     This post is about the computational fluid dynamics analysis of a wing. The wing analyzed employed the NACA 0021 section throughout. The wing had a span of 4 m and a chord length of 1 m. The Reynolds number was kept at 3,000,000. The software employed was SolidWorks Flow Simulation Premium.

     The mesh had a total of 385,064 cells of which 84,826 cells were in contact with the wing surface, as shown in Fig. 1. The results are, indeed, mesh independent. Mesh controls were employed to refine the mesh near the wing surface. The computational domain employed was of cylindrical shape.

 

Fig. 1, The computational mesh around the wing.

 

     The velocity variation at various angles of attack around the wing cross-section is shown in Fig. 3 while the pressure variation on the wing surface is shown in Fig. 4. The results were validated against experiments conducted by [1].

 

Fig. 2, Velocity variation around the wing at 0-25 degree AOA, 5 degree increments.

 

Fig. 3, Pressure variation at the wing surface at 0-25 degree AOA, 5 degree increments.

     The purpose of this blog is maintain my online portfolio. I did this analysis because I realized I haven't written anything of this nature before. All of my previous simulations and/or blog entries were from the propulsion, renewable energy and turbo-machinery areas.
 

     Update 01

     CAD files are available here.

 

    
     Thank you for reading. If you would like to collaborate on research projects, please feel free to contact.

     [1] Fernando A. Rocha, Adson A. de Paula, Marcos d. Sousa, André V. Cavalieri, and Vitor G. Kleine, "Lift enhancement by wavy leading edges at Reynolds numbers between 700,000 and 3,000,000," Proceedings of the 2018 Applied Aerodynamics Conference, AIAA AVIATION Forum, Atlanta, GA, 2018.

     https://doi.org/10.2514/6.2018-3816

Comparison of VAWT Blade Designs (Leading-Edge Tubercle, Leading and Trailing-Edge Tubercle, Unmodified) (Update 05)

Numerical simulations were run on SolidWorks Flow Simulation Premium (model files are available here) software to compare the torque characterizes of three distinct vertical-axis wind turbine blade designs shown in Fig. 1. The torque characteristics are shown in Fig. 2.

This publication was used to verify and validate the numerical methodology. The results were within 8% of the publication's results at the design point of TSR of 1.2 at 90 RPM and 7.85 m/s wind speed. The dimensions of the turbine, the  blades and the cross section used are mentioned in the publication.

Fig. 1. Top Row, L-R: VAWT with blades having tubercles at the leading edge (ten tubercles per blade span, configuration name 10T), VAWT with blades having tubercles at both the leading and the trailing edge (ten tubercles per blade span). Bottom Row, VAWT with blades having no modifications.

It is clear from the Fig. 2 that the baseline design provides the most stable torque. On average the turbine with no modifications on the blades produced 5.31 Nm torque in one complete rotation, while the turbine with tubercles at the leading edge only, produced 5.20 Nm torque. The turbine with tubercles added to both the leading and the trailing edge produced 5.09 Nm torque in one complete rotation.

The peak torque was maximum for the turbine with the leading edge tubercles, followed by the turbine with the tubercles added to both the leading and the trailing edge of the turbine blades and the turbine with no modifications on the blades at 21.59 Nm, 21.45 Nm and 20.58 Nm respectively.

Fig. 2. Top Row, L-R: Torque curves for VAWT with blades having tubercles at the leading edge, Torque curves for VAWT with blades having tubercles at both the leading and the trailing edge. Bottom Row, Torque curves for VAWT with blades having no modifications. Three colors denote each of the blades in the turbine.

CFD post processing will be added later (may be next week). The effect of leading edge tubercle geometry will be investigated next. The blade design with tubercles added to both the leading and the trailing edge will not be investigated further because it produced the lowest average torque and second highest peak torque.

Update 01:

Decreased the number of tubercles per unit length of the blade, i.e. made the wavelength of the tubercles longer, kept the sweep angle same. As a result, the average and peak torque decreased to 4.53 Nm, and 19.33 Nm, respectively. The figure is attached.

Fig. 3. T-B: Torque curves for VAWT with blades having large wavelength tubercles at the leading edge (five tubercles per blade span, configuration name 5T45). Three colors denote each of the blades in the turbine. Render of the blades.

Update 02:

Increased the number of tubercles per blade span, i.e. made the wavelength of the tubercles smaller, kept the sweep angle same. As a result, the average and peak torque increased to 5.80 Nm, and 23.36 Nm, respectively. The figure is attached.

Fig. 4. T-B: Torque curves for VAWT with blades having smaller wavelength tubercles at the leading edge (fifteen tubercles per blade span, configuration name 15T45). Three colors denote each of the blades in the turbine. Render of the blades.

Update 03:

Again, increased the number of tubercles per blade span, i.e. made the wavelength of the tubercles smaller, kept the sweep angle same. As a result, the average and peak torque increased to 6.1 Nm, and 24.12 Nm, respectively. The figure is attached.

Fig. 5. T-B: Torque curves for VAWT with blades having smaller wavelength tubercles at the leading edge (twenty tubercles per blade span, configuration name 20T45). Three colors denote each of the blades in the turbine. Render of the blades.

Update 04:

Once more, increased the number of tubercles per blade span, i.e. made the wavelength of the tubercles smaller, kept the sweep angle same. As a result, the average and peak torque increased to 6.42 Nm, and 24.63 Nm, respectively. The figure is attached.

Fig. 6. T-B: Torque curves for VAWT with blades having smaller wavelength tubercles at the leading edge (twenty-five tubercles per blade span, configuration name 25T45). Three colors denote each of the blades in the turbine. Render of the blades.

A table for the tubercle geometry is shown below.

Table 01, Tubercle Geometry

Configuration Name	Amplitude (m)	Wavelength (m)	Sweep Angle (°)
Baseline	0	0	0
5T45	0.12777778	0.25555556	45
10T45	0.06052632	0.12105263	45
15T45	0.03965517	0.07931034	45
20T45	0.02948718	0.05897436	45
25T45	0.02346939	0.04693878	45

It is evident from Table 2 that adding more tubercles to the wind turbine's blade causes an increase in both the peak and the average torque. But it is also clear from the Table 2 that the percentage difference in both the average and the peak torque from the previous configuration (less tubercles per blade span) decreases as the number of tubercles per blade span is increased. It appears to be converging to a value.

Table 02, Tubercle Efficiency

Configuration Name	Peak Torque (Nm)	Average Torque (Nm)	Percentage Difference in the Average Torque from the Previous Configuration	Percentage Difference in the Average Torque from then Baseline Configuration
Baseline	20.58	5.31	N/A	N/A
5T45	19.33	4.53	-17.22	-17.22
10T45	21.59	5.2	12.89	-2.12
15T45	23.36	5.8	10.35	8.45
20T45	24.12	6.1	4.92	12.95
25T45	24.63	6.42	4.98	17.29

I think the difference between both the peak and the average torque produced by 25T45 and 20T45 configuration is comparable, up next, a new sweep angle.

Update 05

Following are my publications relating to the subject of this post.

Butt, F.R., and Talha, T., "A Numerical Investigation of the Effect of Leading-Edge Tubercles on Propeller Performance," Journal of Aircraft. Vol. 56, No. 2 or No. 3, 2019, pp. XX. (Issue/page number(s) to assigned soon. Active DOI: https://arc.aiaa.org/doi/10.2514/1.C034845)

Butt, F.R., and Talha, T., "A Parametric Study of the Effect of the Leading-Edge Tubercles Geometry on the Performance of Aeronautic Propeller using Computational Fluid Dynamics (CFD)," Proceedings of the World Congress on Engineering, Vol. 2, Newswood Limited, Hong Kong, 2018, pp. 586-595, (active link: http://www.iaeng.org/publication/WCE2018/WCE2018_pp586-595.pdf).

Butt, F.R., and Talha, T., "Optimization of the Geometry and the Span-wise Positioning of the Leading-Edge Tubercles on a Helical Vertical-Axis Marine Turbine Blade ," AIAA Science and Technology Forum and Exposition 2019, Turbomachinery and Energy Systems, accepted for publication.

Thank you for reading.

SolidWorks Animation: Transient NREL Phase-VI Wind Turbine CFD Simulation [Validated]

     10 KW wind turbine CFD simulation using Flow Simulation Premium. Design points: 10 m/s wind speed, rotational velocity 7.5 rad/s.

     The rendered volume shows vorticity (curl of the velocity field). It is colored by dynamic pressure. Low pressure in the center of the helix shows very small wind speed.

     Power from the CFD analysis was 9,854.96 W while the experimental power is 10,000 W, a difference of only 1.45 %, that too by using only 693,141 cells in the mesh.

     Do you want me to make a tutorial about the simulation setup with SolidWorks Flow Simulation Premium?