2D Backward Facing Step (BFS) Flow Simulation

OBJECTIVE

The objective of this project is to setup a 2D backward facing step flow problem with three different mesh size cases as mentioned below and solve using Converge Flow Solver.

Comparison of the results between all the three simulations.

Case 1: dx = 2.0e-3 m; dy = 2.0e-3 m; dz = 2.0e-3 m

Case 2: dx = 1.5e-3 m; dy = 1.5e-3 m; dz = 1.5e-3 m

Case 3: dx = 1.0e-3 m; dy = 1.0e-3 m; dz = 1.0e-3 m

NOTE: During the simulation setup, learnt to restart simulation. Initially, tried to run simulations for just 7000 cycles. Using the feature restart, continued to run simulation from 7000 cycles to 15000 cycles. Then, merged the ‘two’ set of output files into ‘single’ output file, to plot the results for mass flow rate, grid cells count, etc.

Understanding of Problem:

Backward facing step channel (considered for simulation) dimensions of the bounding box, i.e., 0.272 × 0.02 × 0.025 (all dimensions considered in metres).

Pressure and Temperature considered at different zones are:

• At inlet, 110325 Pa and 300 K
• At outlet, 101325 Pa and 300 K
• In volumetric region, 101325 Pa and 300K

The flow direction would be from inlet to outlet. Because, pressure difference between inlet and outlet is 9000 Pa.

Here, the flow encounters a sudden backward facing step, i.e. increase of cross-sectional area causing the separation of flow which we are going to observe through the simulation.

Mechanical Features of Backward Facing Step:

Reference: A review of Backward-Facing Step (BFS) flow mechanisms, heat transfer and Control, Len Chen et.al.

While drawing inferences, we will try to capture these features through our simulation outcomes.

Case Setup

Three standard steps for performing any simulation:

1. Pre-Processing
2. Solver
3. Post-Processing

Pre-Processing

Importing Geometry

Just drag and drop

Every entity in converge is represented by triangles (or discretised into minimum possible triangles), which are further used to define boundaries.

Boundary Creation

Geometry > Boundary > Flag > Add

For the channel flow simulation, the complete geometry is defined through five boundary names:

1. Inlet
2. Outlet
3. Top-Bottom-Walls
4. Front2D
5. Back2D

Each boundary is assigned with their respective boundary ID number and defined using selection or grouping of triangles.

For simulating in 2D, when the defined geometry is three dimensional, we need to define the symmetry conditions or define axis for infinite length which would replicate the same exact results along the z axis for any cross section or slice.

Case Setup (Some of the settings are only being mentioned/explained here)

Application:

• General Flow

Material:

• Air (predefined mixture)
• Gas simulation
• Reaction mechanism

Simulation Parameters:

• Pressure-based steady solver
• Start-End time (0-7000 cycles)
• Initial time step (1.0e-6 s)
• Minimum time step (1.0e-6 s)
• Maximum time step (1.0 s)

Boundary:

• Inlet – INFLOW – 110325 Pa – Air
• Outlet – OUTFLOW – 101325 Pa – Air
• Top-Bottom-Walls – WALL – Law of Wall
• Front2D – TWO_D
• Back2D – TWO_D

Inlet and Outlet pressure differs by 9000 Pa, by which we can say that flow would be turbulence initially. Thus, requires turbulence models for simulation.

Regions and initialization:

Based on defined boundaries, the boundaries are grouped to form volumetric zones called as regions.

This region is pre-set (default) to values for pressure and temperature, i.e., 101325 Pa and 300 K respectively.

This region is defined or contained with species that is air. Air comprises of 23% oxygen and 77% nitrogen by mass fraction.

Base Grid:

It gives the control to define base grid size, i.e., dx, dy, and dz which is considered to be the input for discretization of geometry into cells or to create mesh. Cell count estimation helps to predict cell counts after discretization

• (2.0e-3, 2.0e-3, 2.0e-3) >>>>>>> 12557.7 (estimated cell counts)
• (1.5e-3, 1.5e-3, 1.5e-3) >>>>>>> 29766.4 (estimated cell counts)
• (1.0e-4, 1.0e-4, 1.0e-4) >>>>>>> 100462 (estimated cell counts)

Fixed Embedding:

In the grid control, it provides the control to refine the grid with reference to base grid

• Entity Type – Boundary
• Boundary ID – top-bottom-walls
• Mode – Permanent
• Scale – 2 (defines level of refinement)
• Layer – 6 (defines number of layers)

Post Variable Selection: (Default)

Output files:

Time interval to writing 3D output data files – 100.0 cycles

Time interval for writing text output – 1.0 cycles

Time interval for writing restarting output – 1000.0 cycles

Solver

Using Cygwin, a dos-based solver is used to solve the problem with the help of numerical technique.

The input files directory for backward facing step is to be set as present working directory, through which solver is executed using the below mentioned command line.

mpiexec.exe -n 4 converge.exe </dev/null> logfile &

This command line helps to solve using parallel processing with better usage of system cores and also saves ‘logfile’, easy to understand and helpful for future reference.

Post-Processing

1. Grid (Group of Cells)

Based on three cases, we can see that cell count increases as the base grid size decreases shown in below figure and data has also been tabulated below.

Case No.	Base Grid Size (m)	Total Cell Count
1	2.0e-3	9884
2	1.5e-3	12524
3	1.0e-3	21885

The below figure shows the refinement of mesh based on grid size definition and fixed embedding at boundary with two level of refinement and six layers. For ease in visual inspection and comparison, all three cases have been magnified to same level.

Case 1: 2.0e-3

Case 2: 1.5e-3

Case 3: 1.0e-3

2. Mass Flow Rate

Based on the observation from mass flow rate plot for each case shown below, we can say that as the total cell count is increased for any case setup, the cycles required to reach the steady state also increases. Here, in all the three cases, we are able to achieve steady state much before the completion of cycles, i.e. 15000.

During the initial cycles, the occurrence of occasional disturbances can be observed. These disturbance gets damped over a period of cycle as reached the steady state depending upon the consideration of base grid size (cells density).

While comparison, its good practise to compare with variation of only one variable. Thus, here only base grid size is varied while cycles are fixed to 15000 irrespective of achieving the steady state for all cases that here is not a problem.

3. Conservation of Mass Flow Rate

Previous section shows the plot for mass flow rate at both inlet and outlet. It can be seen that inlet mass flow rate is negative while at outlet is positive. It must be a preliminary check after getting solutions for any simulation. As, we know that mass flow rate must be conserved (i.e., what goes in…. must come out….). Since, inlet and outlet have opposite signs, then on their addition it must be zero. Since, the solver employs numerical techniques to reach towards the result, there would be some errors but are in acceptable or considerable range, i.e. can be seen from 9000 to 15000 cycles in the shown below plot that errors are well within considerable limits of `±0.0004`.

4. Computational Time

The computational time taken by solver to compute for all three cases has been tabulated below:

Case No.	Total Time (s)	Cycles Run
1	703.082	15000
2	1091.178	15000
3	2534.128	15000

We can infer from the plot (obtained from table shown above) that as the base grid size is decreased the computational time increases. But, able to achieve more accurate results with refined grid. Here, all the three cases were run for same cycles, i.e. 15000 cycles and are able to achieve steady state for all three cases. Based on the variation of computational time w.r.t base grid size, we could observe the exponential rise in computational time with decrease in grid size (which means increase in cell count).

Normalized plot for computational time variation helps to understand much more easily. The computational time and base grid size of case 1 has been considered as a reference value. Thus, just through observation as base grid reduced to half, the computational time increased by the factor of ‘3.60’.

NOTE: In further study, results has been considered for case 3, as it is most approximate results due to finer cells especially at boundaries.

5. Flow Field: Vector and Streamline Plot

In backward facing step, sudden expansion of area takes place due to which flow filed of channel gets effected. At the point of backward facing step, the phenomena of flow separation take place and as well as initialization of bubble formation. As the flow field is fully developed along the span of the channel, the mean position could be located where reattachment occurs and further, the boundary layer is recovered. Through the vector plot, the location could be identified at the wall where the vector points towards both the positive as well as negative direction. At the reattachment location, negative velocity vector implies that some flow field gets reversed along the wall. Due to which the wall experiences drag force, thus resulting in wall shear. As this flow field again reaches backward facing step, the kinetic energy gets completely converted into force called as pressure force, i.e. velocity becomes zero. Thereafter, it again continues to move along the main flow field. This phenomenon takes place layer-by-layer. Thus, the recirculation zone could be observed through the streamline plot. The recirculation zone and main flow filed is separated by another layer termed as separated shear layer. This layer acts as a transition layer being turbulent, the layer experiences drag force called shear stress, through this layer momentum diffusion takes place from higher to lower velocity flow field.

Thus, recirculation region is bad for mixing and transportation problem. The recirculation region can be avoided in channel flow by eliminating sudden step, instead considering steep or curved step which guides boundary layer while dodging boundary separation.

6. Velocity Contour and Profile

The velocity contour has been shown above, representing the variation of velocity along channel’s span. The velocity profile has been plotted at 5 different section as per the velocity contour, to understand the change in velocity profile due backward facing step.

In section 1-1’ and 2-2’, laminar flow could be observed. At section 3-3’, wall separation takes place due to which first half height of the channel has velocity profile is similar to laminar flow but the second half have velocity as zero throughout. Intermediate to these profiles, there is a sudden but little drop in velocity with magnitude as negative. This represents the initiation of recirculation zone. While, section 4-4’ and 5-5’ represents flow field being developed and recovering of boundary layer. Here, the velocity profile gets skewed towards the top wall due to the reason of reattachment and recovery of boundary layer at bottom wall.

The velocity contour has been shown above, representing the variation of velocity along channel’s span. The velocity profile has been plotted at 5 different section as per the velocity contour, to understand the change in velocity profile due backward facing step. As per observation from the velocity profile shown below along the channel’s span, no definite inferences could be drawn. So, this plot being irrelevant could be avoided along the longitudinal direction.

7. Pressure Contour and Profile

The pressure contour has been shown above, representing the variation of pressure along the channel’s span. The pressure profile has been plotted at 5 different traverse section as per the pressure contour, to understand the change in pressure profile due to backward facing step. At the point of wall separation due to backward facing step, the significant drop in pressure is observed and the pressure force experienced by the walls can be observed. Between the section 3-3’ and 4-4’, the recirculation region is present where pressure is too low. And, the pressure we are talking about here is static pressure along the channel.

The above shown pressure contour represents the sections along the longitudinal direction based on which pressure profile is plotted. This pressure plot is not required as same thing is quite understandable through the contour and plot shown previously.

8. Input and Computed Parameters at inlet and outlet

As per the case of our simulation, we are able to achieve the steady state which was our prime motive while using steady state solver. Now, lets compare different properties either given as boundary conditions or obtained as output through simulation at inlet and outlet.

	Inlet	Outlet	Δ
Mass Flow Rate (kg/s)	1.37131	1.37128	-0.00003
Average Velocity (m/s)	135.796	58.0168	-77.78
Average Density (kg/m3)	1.1813	1.18189	+0.00059
Total Pressure (Pa)	110325	103616	-6709
Static Pressure (Pa)	9905.1	101325	+91420
Average Temperature (K)	290.861	297.406	+6.54
Average Mach Number (dimensionless)	0.396312	0.167832	-0.23

From the above tabulated data, we can observe that:

1. Mass flow rate is conserved.
2. Due to sudden expansion due BFS model (by continuity equation), velocity decreases at outlet.
3. During simulation parameters definition compressible flow solver is used, but not much change in density is observed at outlet.
4. By total pressure at inlet and outlet, pressure drop of 6709 Pa is observed, which implies flow direction is from inlet to outlet.
5. Initial temperature of 300 K is defined at inlet and outlet. But after steady state, temperature at inlet and outlet dropped by 10 and 3 K respectively. Overall rise in temperature as computed is 6.54 K.
6. At inlet and outlet, Mach number is 0.3963 and 0.1678, i.e. 39.63% and 16.78% of speed of sound. Implies flow is subsonic in nature.
7. Since, Mach number is significantly smaller, i.e., < 0.2-0.3 that means flow is quasi-steady and isothermal. Compressibility effects will be small and simplified incompressible flow equations can be used.

Expansion Ratio:

Continuity equation at inlet and outlet,

Substituting computed results using simulation,

Geometric measurement,

`h=8.5484 mm ;H=1.9975 mm`

So,

`h/H=(8.5484×10^(-3))/(1.9975×10^(-2) )=0.4279`

Reynold’s Number:

`S=H-h=11.427 mm ; μ=1.2×10^(-2) Pa∙s`

`Reynold's Number=(ρvD)/μ=(1.1813×135.796×D)/(1.2×10^(-2) )=13368×D`

As per different literature, characteristic length, D could be defined as either h, or S (step size).

`Re_h=13368×h=114.27`

`Re_S=13368×S=152.75`

Conclusions:

• As the base grid size is decreased, the total cell and the computational time increases.
• Damping of the disturbances depends on the base grid size, smaller the base grid size more cycles taken to damp completely and achieve steady state. Whereas in our case, we are able to achieve steady state for each of them.
• For all the three cases, mass flow rate is conserved through the backward facing step channel flow.
• We were able to observe the key features of backward facing step through simulation results, a classical example of CFD.
  • Able to locate the reattachment point of boundary layer.
  • Able to observe the recovery of boundary layer.
  • Able to depict different forces acting through recirculation region.
  • Abe to understand the concept of momentum diffusion through shear layer which could be either laminar or turbulent based on flow field.
  • At the recirculation zone, the low pressure was observed as well as lower velocity that to opposite to main flow field.
• The total pressure drops at outlet inferred the flow direction from inlet to outlet.
• At outlet, the velocity was lower in comparison to inlet velocity due to BFS, as could also be verified using continuity equation applied to inlet and outlet.