Week 9 - Mini project - Flow over an Airfoil

In this project, a trasient simulation of a NACA 2412 airfoil will be performed in converge. The maximum camber is 2% of chord length at 40% of the chord length. The maximum thickness is 12% of chord length at 30% of chord length. I will also be looking at the effect of angle of attack on Co-efficient of lift and drag.

An airfoil is the cross-sectional shape of an object whose motion through a gas is capable of generating significant lift. It is used in tubines, aircraft winds, propellors, etc.

Geometry Creation

The geometry can either be created manually using converge or it can also be imported from another CAD software, in my case I created the geometry. The geometry was sourced from: http://airfoiltools.com/airfoil/details?airfoil=naca2412-il

The dat file can either be imported into a CAD software, or a vertex can be created for each point in CONVERGE. Then, the vertices can be joined, and triangles can be created.

This is the resultant geometry in converge.

This airfoil's chord length was 1m long, the entire geometry was scaled down by 0.1. Hence, final length is 0.1m.

After the airfoil is created, a wind tunnel has to be created. The tunnel should be atleast 50 chord lengths (5m) and 30 chord length (3m). In this case, the wind tunnel is a cubioid of length 5.1m, breath 3m and width 0.1m. (width can be anything, as simulation will be 2D) The airfoil is a hole in the cubiod.

This is the geometry with the wind tunnel in converge.

To modify the geometry for different angles of attack, the center of the airfoil was roughly determined using the measure tool. The vertex points of only the airfoil were selected and rotated for the required angle with the center earlier determined being the center of rotation. A full diagnosis was conducted and errors if any were rectified.

Boundary

This is the boundary setup as seen in converge:

Case Setup

Now, begin the case setup:

Application Type

Materials: Select Air as predefined mixture.

 

Gas Simulation, Global Transport Parameters and Reaction Mechanisms were set to default.

In species, O2 and N2 were added.

Simulation Parameters

Default values were used for Run Parameters

Simulation Time Parameters

The inlet velocity is 31.5 m/s and the length of the domain is 5.1m, this means the time taken for each cycle would be around 0.15s, 3 cycles would be more than enough: 3*0.15 = 0.45 ~= 0.5s.

Note: a maximum convection CFL limit is required else the solution will never converge.

Default values for Solver parameters were used.

Boundary Conditions

inlet: INFLOW type with 0.315m/s velocity in the x axis. Temperature of 300K. Other values were set to zero neumann condition.

Calculation of Inlet Velocity

We know, `Re = (rho * v * L) / mu`

where, `L = 0.1m, rho = 1.177 kg//m^3, mu = 1.853 xx 10^-5  N*s//m^2 and Re = 200000`

`V = (Re * mu)/(rho*L)`

`V = (200000 * 1.853 xx 10^-5) / (1.177*0.1)`

`V = 31.48683092608326 m//s ~= 31.5m//s`

outlet: OUTFLOW type with 101325 Pa of Pressure. Other values were set to zero neumann condition.

front: TWO_D

back: TWO_D

topAndBottom: SYMMETRY

airfoil:

Initial Conditions:

The inlet velocity was used as the initial velocity

Physical Models

Turbulence model used was `k-omega SST` as it the best model for this type of flow.

The values were not changed.

Grid Control: Fixed Embededding was enabled to ensure a finer mesh is used in regions near the airfoil.

Note: AMR was not used as the area of interest of this simulation is properties on and near the surface of the airfoil only. AMR used considerable resources to precisely calculate regions of the domain I was not interested in.

Base Grid: This is the step where sizes of each element is provided.

Output/Post Processing

All the default variables were selected for post variable selection.

Under Output Generation, Wall Output: Boundaries Only was chosen to get values of flow variables near the airfoil boundary

Around 100 timesteps were required to ensure the post processed result is smooth enough for a good animation, hence the time interval was chosen accordingly.

Now, our case setup is complete. The files will be exported into a folder using the Files Export tool (File -> Export->Export input files)

In total 12 files were exported, these are:

These files contain all the necessary information for the simulation.

Running the Simulation

1. Open cygwin

2. Navigate to directory where case files were exported

3. Run the following command

mpiexec.exe -n 4 "C:\Program Files\Convergent_Science\CONVERGE\3.0.16\bin\intelmpi\converge.exe" restricted </dev/null> logfile.txt &

This will take a while, you can view the progress in task manager or by opening the logfile. CPU usage is usually maxed out.

Once CPU usage drops from 100%, the output files are generated. This simulation took around 45 minutes. To view them in paraview, we must export them to a format which is supported by paraview.

Go to 3D-post processing in converge,

These files can be read by paraview.

Post-Processing

In Paraview

Import these files into paraview

The required plots/animations are generated in paraview.

In converge

Go to Line plotting, select the case folder and plots can be viewed

Outputs and Plots with explanations

1. Mesh

These were taken in paraview.

This is the mesh of the simulation, there are more cells near the airfoil/center due to fixed embeddding.

zooming in,

zooming in,

There are about 10 layers of the smallest cells on the surface of the airfoil. These cells have a size of 0.000625m (0.08 / 2^7). This way, more resources are spent to make the solution more precise only around the airfoil.

2. Pressure, Temperature Contours

These were taken in paraview.

The temperature is slightly higher at forward and backward end, this can be attributed to frictional heating.

Pressure is highest at the forward tip of the airfoil due to accumulated airflow and change in streamlined flow. Pressure is lower on the top of the airfoil and higher on the bottom, this is expected and contrast to velocity. The difference in pressure is what causes the airfoil to generate lift forces, and this is the basis for modern aviation.

3. Pressure and Y+ plots

This graph shows the lift force and drag force (Pres_Y and Pres_X respectively). The drag force increases as the angle of attack increase, this is expected as more eddies and turbulent layers are formed when angle of attack increase (streamlined flow is disturbed). The lift force increases till 10 degrees and drops as the airfoil crosses its stall angle.

All values of y are above 30 indicating cells are not in the transisitonal region, which is a good thing. Except when angle of attack is 15 degrees.

4. Animation

https://youtu.be/d_I54z7Ezdc

The animation shows how velocity flows around the airfoil from the beginning to the end of the simulation. The velocity is higher on the top of the airfoil and low on the bottom. This is in stark contrast to the pressure and is expeceted, according to reynold's law and conservation of momentum.

5. Drag Coefficient v/s Angle of Attack

Drag Coefficient is a dimensionless quantity that is used to quantify the drag or resistance of an object in a fluid environment, such as air or water.

`underset(d)(c) = (2 * underset(d)(F))/(rho*u^2*A)`

where, `underset(d)(c)` is the coefficient of drag, `underset(d)(F)` is the drag force, `rho` is the density of the fluid, `u` is the velocity of the fluid and `A` is the cross-section/reference area.

In this case,

Drag direction is in the x-axis

`underset(d)(c) = (2 * underset(X)(P))/(rho*u^2*A)`

Where, `rho = 1.177 kg//m^3, u = 31.5 m//s and A = 1m^2`

`=> underset(d)(c) = 8.5625 xx 10^-4 * underset(X)(P)`

This expression was used to calculate Drag Coefficient for each angle of attack, the results are shown below:

The coefficient of drag increases as angle of attack increases. This is expected as the angle of attack increases, the airfoil becomes less streamlined to the incoming fluid. This causes the airfoil to generate more resistance to flow than before.

6. Lift Coefficient v/s Angle of Attack

Lift Coefficient is a dimensionless coefficient that relates the lift generated by a lifting body to the fluid density around the body, the fluid velocity and an associated reference area.

The relation is given by:

`underset(L)(C) = (2 * L)/(rho*u^2*A)`

where, `underset(L)(C)` is the coefficient of lift, `L` is the lift force, `rho` is the density of the fluid, `u` is the velocity of the fluid and `A` is the cross-section/reference area.

In this case,

Lift direction is the y-axis

`underset(L)(C) = (2 * underset(Y)(P))/(rho*u^2*A)`

Where, `rho = 1.177 kg//m^3, u = 31.5 m//s and A =1 m^2`

`=> underset(L)(C) = 8.5625 xx 10^-4 * underset(Y)(P)`

This expression was used to calculate Lift Coefficient for each angle of attack, the results are shown below:

The coefficient of lift increases as angle of attack increases and drops after 10 degrees. This is because the airfoil is past its stall angle. A stall angle is the angle at which the coefficient of lift no longer increases and lift generated by the airfoil is surpassed by the opposite force: gravitational force. This causes any aircraft to "lose lift" and the aircraft is losing altitude now. This observation is in line with other studies on the NACA 2412 airfoil. Example: https://ijisrt.com/assets/upload/files/IJISRT20MAR229.pdf.

Effect of other turbulence models on the simulation

Standard k-epsilon model was used as the test case. This model should not be used for this type of simulation because it has poor handling of jets, seperation, etc.

1. Drag Coefficient v/s Angle of Attack

This is the relation between angle of attack and drag coefficient for standard k-epsilon:

The relation has not changed for the new turbulence model.

2. Lift Coefficient v/s Angle of Attack

This is the relation between angle of attack and drag coefficient for standard k-epsilon:

The relation has not changed for the new turbulence model. Lift coefficient values seem to have changed sligthly however.