Week 10: Modeling and Simulation of flow around an Ahmed Body

Ahmed body is a simplified vehicle body used for capturing basic yet charecteristic features seen in objects used in the automobile industry.

Calculation of Reynolds Number

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 u = 40 m//s`

`Re = (1.177 * 40 *  0.288) / (1.853 xx 10^-5)`

`Re = 7.31 xx 10^5`

Geometry Creation

The geometry can be created in any CAD software. This is the geometry:

I used FreeCAD due to the lack of a better alternative. The process is quite easy. Draw a cubiod of dimensions 1044*288 units with the bend. Extrude it on both sides for 194 units. On the front end, use the fillet tool to create fillets of radius 100 units for all edges. Then, the legs/cylinders of each 30 units diameter can be added. The angle of bend (`varphi`) is 35 degrees for this geometry.

This is the geometry seen in FreeCAD:

Export the geometry as an .stl file and import it in converge. Scale down the geometry if required. Now, a virtual wind tunnel must be created. The dimensions depend on the processing capabilities of the computer used. In this case, the distance in-front of the body to the inlet is 3 times the body length (~3 metres) and distance from the end of the body to the outlet is 8 times the body length (~8 metres). The width and height of the tunnel are 3m and 3m respectively.

Since this simulation will be compared to experimental data, the experimental conditions must be followed as much as possible. One of the conditions is that a boundary layer of 30mm (on the floor/bottom) is present 400mm from the front end of the body.

Calculation of beginning of Free-Slip (Slip) Boundary

To find distance from the front of the body at which the bottom wall needs to be no-slip so that a boundary layer of 30mm (on the floor/bottom) is present 400mm from the body.

Recall flat plate boundary layer thickness, `delta = 0.382 x // Re ^ (1/5)`

Where `delta` is the boundary layer thickness (= 0.03m), `x` is the distance from the beginning of the boundary layer (=0.4m) and `Re` is Reynolds Number (=`7.31 xx 10^5`)

`x = (delta * Re^ (1/5)) /0.382`

`x = 1.80m`

Distance from origin = 2.844m (since origin is at the end of the body)

This is the geometry with the wind tunnel in converge.

(some surfaces hidden to show ahmed body)

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.

Simulation Parameters

Default values were used for Run Parameters. Steady State solver was checked under Misc.

Steady-State Monitor

Simulation Time Parameters

The inlet velocity is 40 m/s and the length of the domain is 12m, this means the time taken for each cycle would be around 0.3s, 3 cycles would be more than enough: 3*0.3 = 0.9 ~= 1s.

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

Default values for Solver parameters were used.

Boundary Conditions

ahmedBody: WALL type with No-slip surface. Wall Model Dissipation BC and Standard wall function for the same.

freeSlipBottom: WALL type with Slip surface.  zero neumann Dissipation BC and Standard wall function for the same.

bottom: WALL type with No-slip surface. Wall Model Dissipation BC and Standard wall function for the same.

Capturing flow variables in the viscous sublayer for this flow (high Re), is impossible with the current mesh used, hence we use the `k- omega SST` turbulence model and wall model BC to ensure we can skip calculating the turbulent boundary layer and move on to finding flow variables in the logarithmic region.

top: SYMMETRY type

front: SYMMETRY type

back: SYMMETRY type

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

outlet: OUTFLOW type with 101325Pa pressure. Temperature of 300K. Other values were set to zero neumann condition.

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 external flows.

The values were not changed.

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

Two fixed embedding were used. One for the cubiodal/box region around the ahmed body and one adjacent to the walls of the ahmed body.

This is the box seen in converge:

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 2h 30 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

This is the mesh seen in paraview

zooming in,

The mesh is finer than the base mesh grid in a big region around the ahmed body, this is due to fixed embedding 1.

The mesh is even more finer in a small region around the ahmed body, it is about 10 layers thick. This is due to fixed embedding 2. The size of the smallest elements is 0.0125m (= 0.8/2^6).

2. Variation of Drag Force with Simulation Time

This was taken in converge's line plotting tool.

The pressure initially flucuates rapidly and slowly converges to a steady form of periodicity. This apparent periodicity is due to vortex shedding, which causes waves or ripples of pressure changes when flow passes against any sharp/rough edges. This periodocity causes the drag force to exhibit a periodic pattern, the pattern initially has a very high amplitude and eventually drops as solution progesses from transient to a steady state. There is some periodicity in the steady state, but it is extremely minor and can be ignored for most purposes.

3. Calculation of Steady State Drag Force and Coefficient of Drag

The data used to generate "Variation of Drag Force with Simulation Time" is exported to excel. The imported data is averaged over the last 0.1s, to get the steady state drag force = 91.881246520202 N.

This is further used to find coefficient of drag for the ahmed body.

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.

Finding cross-sectional area using Paraview

A plane is created which cuts the ahmed body at its thickest point.

A ruler is then used to measure the height and width of the ahmed body.

Measuring height,

Measuring width,

In this case,

Where, `underset(d)(F) = 91.881246520202 N, rho = 1.177 kg//m^3, u = 40 m//s and A = 0.288*0.389 m^2`

`=> (2*91.881246520202)/(1.177*40^2*0.288*0.389)`

`underset(d)(c) = 0.871`

This is the steady state coefficient of drag for the Ahmed Body.

4. Animation of Velocity and Pressure Distribution

These were taken in paraview. A Slice was created on the Y normal (XZ plane), which allowed for a cut through view of the simulation parameters around the ahmed body.

Velocity Distribution:https://youtu.be/DU0Cq5PKl8o

Velocity is low at the front of the body due to presence of accumulated air at the front. Immediately after that, there is a region of very high velocity below the body. There are multiple regions of velocity above the body. There is the free stream region of high velocity which colides with the vortex regions consisting of regions with very low velocities. This type of zone is also seen towards the bend at the back of the body. The back of the body also contains tiny circular pockets of recirculating air. As the simulation approaches steady state, the velocities become more stable.

Pressure Distribution: https://youtu.be/9ifI7YQpRGc

There is a high pressure region at the front of the body, this is due to accumulation of moving air hitting the front of the body. It is immediately followed by low pressure regions on the top and bottom of the body and some faint waves of vortex shedding can be seen. A low pressure region can also be seen at the top of the bend on the top of the body. Less variation is seen in the pressure distrubution as the flow moves from transient to steady-state.

5. Variation of Velocity along the Symmetry axis of the Ahmed Body

These were taken in paraview using the plot over line tool. Multiple lines parallel to the z axis each of varying length on the XZ plane at various positions before, above and beyond the ahmed body. The X component of velocity and length scale were exported at these points to excel. The simulation data was compared to experimental data at the same location for all cases. These were the results:

Note: The vertical bars on the experimental data are the 5% error margins.

At x=-0.263m,

At x=-0.113m,

At x=-0.063,

At x=-0.013,

At x=0.037m,

At x=0.087m,

At x=0.137m,

At x=0.187m,

Upon initial observation, we can see that the simulation data varies wildly from the experimental data. The simulation data seems to captures the general trend of the change in X velocity (such as decreases, increases, etc.) but the accuracy in which the mentioned trends are measured is quite low. In most cases, the experimental velocity quickly increases and plateaus, while the simulation velocity increases slowly but eventually plateaus.

In locations with finer mesh (x = 0.037m), the part of the line with finer mesh is modelled accurately and is in close agreement with the experimental values. The deviation between experimental and simulation data widens as the mesh starts to become more coarse (from fixed embedding 2 to fixed embedding 1).

Therefore it can be concluded that the reason for such divergence is because the mesh is not fine enough.

Conclusion

1. The simulation runs well and desired results are obtained.

2. The results agree with the simulation data to some extent, like validating general trends, expectations and patterns.

3. Due to the complexity of the simulation, the desired data is not of enough significance due to its low precision and accuracy.

4. The reason for this low precision was deduced to be the coarse mesh, and further mesh refining is required which will increase the computational resources needed for the simulation and time needed drastically.