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.

Flying Wing Design using Computational Fluid Dynamics (Verification and Validation)

     This post is about a transient simulation of the ONERA M-6 flying-wing aircraft with a cross-section of ONERA D airfoil, as shown in Fig. 1.

Fig. 1, The simulated geometry. 

     The mesh had 2,474,614 cells in total with 300,892 cells on the wing surface, as shown in Fig. 2. The computational domain is shown in Fig. 4. The computational domain walls are at a distance equal to ten times the wingspan.

Fig. 2, The computational mesh.

Fig. 3, The computational domain.

     The simulated conditions are taken from [1-4] i.e. a freestream Mach number of 0.8395 at 101,325 Pa and 293.2 K. The angle of attack is 3.06°.

     The lift force from the present numerical simulation is at 19,551.65 N as compared to the lift force of 20,438.53 N as determined by [1-4]. The numerically determined drag force from present simulation is 1,370.88 N as compared to 1,313.24 N, as determined by [1-4].

     The results of lift and drag force from the present simulation are within 4.35% and 4.2% of the results calculated by NASA [4] and [1-3]. Streamlines and pressure surface plot around the aircraft surface are shown in Fig. 4.

Fig. 4, Results.

References

[1] Le Moigne, "A Discrete Navier-Stokes Adjoint Method for Aerodynamic Optimization of Blended Wing-Body, Configurations", PhD thesis, Cranfield University, United Kingdom, 2002.

[2] J. Lee, C. S. Kim, C. Kim, O. H. Rho, and K. D. Lee, "Parallelized Design Optimization for Transonic Wings using Aerodynamic Sensitivity Analysis", AIAA Paper 2002-0264, 2002.

[3] E. J. Neilsen and W. K. Anderson, "Recent Improvements in Aerodynamic Design Optimization on Unstructured Meshes", AIAA Paper 2001-0596, 2001.

[4] 3D ONERA M6 Wing Validation, https://turbmodels.larc.nasa.gov/onerawingnumerics_val_sa.html.

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

Automotive Computational Fluid Dynamics (CFD) Analysis

     This post is about the numerical analysis of an Ahmed body. It was a new experience because of the area between the car's floor and the road which is different from most of the numerical analysis performed in open channel aeronautics and turbo-machinery.

     The numerical analysis was performed using the commercial software, SolidWorks Flow Simulation. The software employs κ − ε turbulence model with damping functions, SIMPLE-R (modified), as the numerical algorithm and second-order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The time derivatives are approximated with an implicit first-order Euler scheme. Flow simulation solves the Navier–Stokes equations, which are formulations of mass, momentum, and energy conservation laws for fluid flows. To predict turbulent flows, the Favre-averaged Navier–Stokes equations are used.

     The software generates Cartesian mesh using immersed boundary method. The mesh had a cell size of 0.035 m in the far field regions within the computational domain. A fine mesh was need between the road the the car's floor to make sure the interaction of the car's floor with the road was captured accurately. Therefore the mesh between the car's floor and the road was refined to have a cell size of 0.00875 m. Another mesh control was applied around the body to refine the mesh with a cells size of 0.0175 m to capture the trailing vortices. The resulting mesh had 209,580 total cells, among those cells, 31,783 cells were at the solid fluid boundary. The computational domain size was ~1L x 1.12L x 3L where L being the vehicle's length. The computational domain along with the computational mesh is shown in Fig. 1.

Fig. 1 Mesh, computational domain and the boundary conditions.

     The red arrows within the Fig. 1 represents the inlet boundary condition of ambient (free-stream) velocity and the blue arrows represent the outlet boundary condition of the ambient pressure. The green arrows represents the co-ordinates axes direction.

     The results from the numerical analysis were compared with [1-3]. The results are within 10% of the experimental results. The velocity (superimposed by the velocity streamlines) and pressure profiles around the car body at various free-stream velocities is shown in Fig. 2.

Fig. 3 Velocity and pressure plots. From the top, Row 1, L-R; ambient velocity of 30 and 40 m/s. Row 2, L-R; ambient velocity of 60 and 80 m/s. Row 3, ambient velocity of 105 m/s.

     It was a good experience learning about automotive CFD after spending a long time in aeronautic/turbo-machinery CFD. Thank you for reading. Please share my work. If you would like to collaborate on a project please reach out.

[1] F.J.Bello-Millán, T.Mäkelä, L.Parras, C.delPino, C.Ferrera, "Experimental study on Ahmed's body drag coefficient for different yaw angles", Journal of Wind Engineering and Industrial Aerodynamics, Volume 157, October 2016, Pages 140-144.

[2] Guilmineau E., Deng G.B., Queutey P., Visonneau M. (2018) Assessment of Hybrid LES Formulations for Flow Simulation Around the Ahmed Body. In: Deville M. et al. (eds) Turbulence and Interactions. TI 2015. Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol 135. Springer, Cham.

[3] A. Thacker, S.Aubrun, A.Leroy, P.Devinant, "Effects of suppressing the 3D separation on the rear slant on the flow structures around an Ahmed body", Journal of Wind Engineering and Industrial Aerodynamics, Volumes 107–108, August–September 2012, Pages 237-243.

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

Computational Fluid Dynamics (CFD) Analysis of a Horizontal Axis Tidal Turbine (Update 01)

     In this post, the results of numerical simulations carried out on a three-blade horizontal-axis tidal turbine blades are made available. These simulations are a part of larger project relating to a horizontal-axis tidal turbine.

 

     The turbine under investigation had a diameter of 4 m. The aero-foils employed include the NACA 4424 (chord length 0.25 m), NACA 4420 (chord length 0.2312 m), NACA 4418 (chord length 0.2126 m), NACA 4417 (chord length 0.1938 m), NACA 4416 (chord length 0.175 m), NACA 4415 (chord length 0.1562 m), NACA 4414 (chord length 0.1376 m), NACA 4413 (chord length 0.1188 m) and the NACA 4412 (chord length 0.1 m) at a distance of 0.4 m, 0.6 m, 0.8 m, 1 m, 1.2 m, 1.4 m, 1.6 m, 1.8 m, 2 m from the blade root, respectively [1].

 

     The numerical simulations were carried out on SolidWorks Flow Simulation Premium software. The software employs κ-ε turbulence model with damping functions, SIMPLE-R (modified) as the numerical algorithm and second order upwind and central approximations as the spatial discretization schemes for the convective fluxes and diffusive terms. The time derivatives are approximated with an implicit first-order Euler scheme. To predict turbulent flows, the Favre-averaged Navier-Stokes equations are used. The software employs a Cartesian mesh which is created using the immersed boundary method.

 

     The simulations were carried out using the Local rotating region(s) (Sliding) technique. The mesh had a total of 589,763 cells. 132,045 cells were at the turbine blade-rotating region boundary. While the process of mesh generation, a mesh control was employed to refine the mesh near the turbine blades and at the boundary of the rotating region and the stationery computational domain. To refine the mesh for the mesh independent test, the curvature level was increased i.e. the mesh in the areas of interest was refined by a factor of 8. As a result, the refined mesh had a total of 916,577 cells. A total of 221,978 were present around the turbine blade. The computational mesh around the turbine blade is shown in Fig. 1. The refined mesh in the areas of high gradients is also clearly visible. The computational domain had dimensions of 2D x 2D x 2.4D, where D is the diameter of the turbine. The dimensions for the computational domain resemble to those in [2].

     Fig. 2 shows computational domain around the turbine blades. In Fig. 2, the curved teal arrow represents the direction of rotation of the turbine. The green arrows represent the respective co-ordinate directions. The brown arrow represents the direction of the force of gravity. The blue arrow represents the direction of the fluid velocity. The circular disc represents the rotating region.

 

Fig. 1, The computational mesh.

    

Fig. 2, The computational domain and the orientation of the boundary conditions.

 

     The simulations were carried out for a total of 5 tip-speed ratios for the turbine ranging from 2 to 10. The fluid; water, velocity was set at 2 m/s. The results are indeed, mesh independent. The mesh independence test was conducted on the design point of the turbine i.e. at the TSR of 6. The plot between the turbine tip-speed ratio; TSR, and the co-efficient of power is shown in Fig. 3. It can be clearly seen from Fig. 3 that the results are in close agreement with the results from [1, 3]. The CFD results from both studies are lower then the BEM, Blade-Element Momentum, results because the three-dimensional effects are not considered while implementing the BEM method.

 

Fig. 3, Turbine efficiency plot.

 

     Fig. 4 shows velocity streamlines colored by velocity magnitude, both of these features are drawn relative to the rotating reference frame, around the turbine blade cross-section at various tip-speed ratios. It can be seen from Fig. 4 that the turbine stalls at TSR of 2 due to a large positive angle of attack. It can also be seen that as the TSR increases, the angle of attack on the blade decreases, it is because of this reason that the power output from the turbine increases.

 

 

Fig. 4, Row 1, L-R; TSR of 2 and 4. Row 2, L-R; TSR of 6 and 8. Row 3, TSR of 10.

 

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

     [1] Binoe E. Abuan, and Robert J. Howell, "The Influence of Unsteady Flow to the Performance of a Horizontal Axis Tidal Turbine," Proceedings of the World Congress on Engineering. London, 2018.

     [2] Ibrahim, I. H., and T. H. New, "A numerical study on the effects of leading-edge modifications upon propeller flow characteristics," Proceedings of Ninth International Symposium on Turbulence and Shear Flow Phenomena. Melbourne, 2015.
     [3] Bahaj, A., Batten, W., & McCann, J. (2007). Experimental verifications of numerical predictions for the hydrodynamic performance of horizontal axis marine current turbines. Renewable Energy, 2479-2490.

Update 01:

     CFD post processing added. One more TSR simulated.

Steady-State VS Transient Propeller Numerical Simulation Comparison

     This post is about the comparison between steady-state and transient computational fluid dynamics analysis of two different propellers. The propellers under investigation are 11x7 and 11x4.7 propellers. The first number in the propeller nomenclature is the propeller diameter and the second number represents the propeller pitch, both parameters are in inch. The transient analysis was carried out using the sliding mesh technique while the steady-state results were obtained by the local rotating region-averaging method. For details about 11x7 propeller click here, for the details about 11x4.7 propeller, click here.

 

     As expected, the propeller efficiencies of transient and steady-state analysis are within 0.9% of each other, as shown in Fig. 1-2. Therefore, it is advised to simulate propellers and horizontal axis wind turbines using the steady-state technique as long as no time-dependent boundary conditions are employed.

 

Fig. 1, Propeller efficiency plot.

  

 Fig. 2, Propeller efficiency plot.

 

     It can be seen from Fig. 3-4 that time taken by the steady-state simulation to converge is on average 42.37% less that the transient analysis.  The steady-state analysis takes considerably less time to give a solution then a transient analysis.

 

Fig. 3, Solution time.

 

Fig. 4, Solution time.

 

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

1D Transient Diffusion (MATLAB code)

clear; clc;% clear the screen and memory

Nx=500; %number of space nodes

Nt=10000; %number of time nodes

Lx=0.3; %length of space (m)

Lt=10; %length of physical time (s)

dx=Lx/(Nx-1); %grid spacing

dt=Lt/(Nt-1); %time step

c=0.000097; %speed of wave (constant)

a=c*dt/dx.^2;

u=100*ones(Nx,1); %initialization of matrix for the property under investigation

x=zeros(Nx,1); %initialization of space

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

u(1)=-20; %boundary condition

u(Nx)=50; %boundary condition

for t=0:dt:Lt %time loop

    un=u; %u(i,t)

    for i=2:Nx-1 %solution loop, backward in space forward in time

        u(i)=((1-2*a)*un(i))+a*(un(i+1)+un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u) %plotting

    title({t;'Time Elaplsed (s)'},'FontSize',38.5);

    pause(0.001)

end

 

 

1D non-Linear Convection (MATLAB code)

clear; clc;% clear the screen and memory

Nx=41; %number of space nodes

Nt=82; %number of time nodes

Lx=2; %length of space (m)

Lt=1; %length of time (s)

dx=Lx/(Nx-1); %grid spacing

dt=Lt/(Nt-1); %time step

a=dt/dx;

u=ones(Nx,1); %initialization of solution matrix

x=zeros(Nx,1); %initialization of space

hold on

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

for i=0.5/dx:1/dx %initial conditions

    u(i)=2;

end

for t=1:Nt %time loop

    un=u; %u(i,t)

    for i=2:Nx %solution loop, backward in space forward in time

        u(i)=un(i)-u(i)*a*(un(i)-un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u,'k') %plotting

    pause(0.1)

end

1D Linear Convection (MATLAB code) (Updated with CFL condition)

clear; clc;% clear the screen and memory

Nx=41; %number of space nodes

Nt=82; %number of time nodes

Lx=2; %length of space (m)

Lt=1; %length of time (s)

dx=Lx/(Nx-1); %grid spacing

dt=Lt/(Nt-1); %time step

c=2; %speed of wave (constant)

a=c*dt/dx;

u=ones(Nx,1); %initialization of solution matrix

x=zeros(Nx,1); %initialization of space

hold on

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

for i=0.5/dx:1/dx %initial conditions

    u(i)=2;

end

for t=1:Nt %time loop

    un=u; %u(i,t)

    for i=2:Nx %solution loop, backward in space forward in time

        u(i)=un(i)-a*(un(i)-un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u,'k') %plotting

    pause(0.1)

end



New Code, with the CFL condition satisfied.

clear; clc;% clear the screen and memory

Nx=13; %number of space nodes

Lx=2; %length of space (m)

Lt=1; %length of time (s)

dx=Lx/(Nx-1); %grid spacing

c=2; %speed of wave (constant)

dt=dx*0.1/c; %time step, Courant number = 0.1

Nt=(Lt/dt)+1; %number of time nodes

a=c*dt/dx;

u=ones(Nx,1); %initialization of solution matrix

x=zeros(Nx,1); %initialization of space

t=zeros(fix(Nt),1); %initialization of space

hold on

for i=1:Nt-1

    t(i+1)=t(i)+dt;

end

for i=1:Nx-1 %space loop

    x(i+1)=x(i)+dx;

end

for i=0.5/dx:1/dx %initial conditions

    u(i)=2;

end

for j=0:dt:Lt %time loop

    un=u; %u(i,t)

    for i=2:Nx %solution loop, backward in space forward in time

        u(i)=un(i)-a*(un(i)-un(i-1)); %discretized equation, u(i,t+1)

    end

    plot(x,u,'k') %plotting

    pause(0.1)

end