Simulation and Analysis of flow characteristics of a Backward Facing Step in OPENFOAM.

Aim: To simulate and analyze flow characteristics of a backward facing step in OpenFOAM.

 

Introduction:

The study of fluid motion across various geometries is one of the major aspect of a CFD engineer. The fluid flowing over different geometries results into different flow fields (velocity field, pressure field and more). By using this flow field data product manipulations are done for the higher efficacy. In this case, a backward facing step is used as a geometry through which fluid will be analyzed. The flow over backward facing step is one of the contemporary case through which various fluid flow characteristics can be studied. 

 

Objectives:

• To simulate the flow through a backward facing step (Transient simulation).
• To simulate the same flow by changing the grading factor to 0.2.
• Measure and compare the velocity profile at 0.085m (8.5cm) from the inlet of geometry for both case studies.
• To understand the significance of grading/growth factor.

 

Assumptions for the given problem:

• The fluid flowing is newtonian fluid.
• The flow is incompressible and laminar.
• The kinematic viscosity is taken at air at `T = 20^o c`,  `nu =`1.516e-5 `m^2/s` (Source: Air standard table`
• Zero frictional losses along the domain.
• The backward facing step is to be analyzed in 2-dimensions (x and y) and we will assign negligible value to the z-direction and its respective boundaries will be kept empty.

 

Software used: OpenFOAM and Paraview

 

Solver used: icoFoam

 

Brief course of action:

1. Selecting an appropriate solver for our case i.e icoFoam

2. Select the appropriate solver tutorial for modifying our geometry and case conditions.

3. Modifying initial and boundary conditions.

4. Adjusting time step with respect to CFL criteria.

5. Initiating execution and visualizing in Paraview.

6. Plotting results and analyse the flow field variables.

 

Geometry of backward facing step:

The above geometry has to be discretized into the set of blocks and combined together to create rigid domain. So now we will discretize the domain into set of 5 blocks.

 

 

Procedure for Numerical set-up and simulation:

1. Open a terminal in the linux operating system.

2. Type in "tut" which is a shortcut to the tutorial directory.

3. Inside the "tutorial" directory, type in "ls" which gives list of other directories/files inside that directory, then type in "cd incompressible/" to access to all the tutorials pertaining to the incompressible flows and its respective solver.

4. Again type in "ls", then access the icoFoam directory and copy the cavity directory by typing (cp -r cavity/ $FOAM_RUN/cavity_test). Here cp is syntax used to copy any file or folder, "-r" is used to denote the entire directory with all the files then "$FOAM_RUN" is the run directory inside the OpenFOAM which allows the user to run user modified problems or tutorials. And finally we rename our cavity file as cavity_test to avoid name rewriting.

5. Now, all we have to play with is the FOAM_RUN directory. Here we are having various directories stating different significance.

a) blockMeshDict (cavity_test/system/blockMeshDict) = It has the geometric and mesh characteristics of the domain.

b) transportProperties (cavity_test/constant/transportProperties) = It has the information of the transport quantity. In our case, the transport object is air and transport property is kinematic viscosity (`nu`).

c) controlDict (cavity_test/system/controlDict) = It has the time control (features like start time, end time, deltaT, writeControl, writeInterval and many more). It governs the time adjustments and modifications.

d) fvSchemes (cavity_test/system/fvSchemes) = It sets the numerical schemes for terms such as derivatives in equations that are calculate the simulation. It has the categories such as gradient, divergence, laplacian , time derivatives etc.

e) fvSolution (cavity_test/system/fvSolution) = It governs the equations, solvers and numerical algorithm.

So the stated above are the files which need to be modified for our case.

6. So starting with the "blockMeshDict", where all the geomteric and mesh modifications are done.

blockMeshDict edit:

Vertice edit: First we wil, edit vertices according to our geometry dimensions. The x y and z co-ordinates are defined in centi-meters (cm). As far as our geometry is concerned, we have a total of 21 vertices.

 

Note: As OpenFOAM is written in C++ language, it takes 0 as an index.

Block edit:

Here, in block edit, we will combine all the inidividual blocks as per their ordered vertices to another blocks. The starting letter "hex" defines the "hexahedral" type of mesh for the block meshing. The total number of cells along x-dimesnion is given to be as 200 and for y-dimension it is given as 20. And here simple grading of (1 1 1) is used, which means the cells along the three dimensions are equally spaced (uniform).

Boundary edit:

There are various boundaries in our domain which are inlet, outlet, front, back and no slip wall ( velocity of the immediate fluid layer along the wall surface tends to be zero). Faces are referred to be as "patch" in OpenFOAM. So the inlet and outlet face is patch from where fluid enters and exits. The front and back boundary is kept to be as "empty"since we are not interested in 3D analysis of backward facing step.The top, bottom and step face is referred to be as "no slip wall", where boundary condition of "wall" is assigned. OpenFOAM will treat as wall with no slip phenonemenon at its surface.

7. Modifying the transport properties which in the constant directory. 

transportProperties:

In this directory we will edit the kinematic viscosity of to 1.516e-5 `m^2/s` from the initial value. As we are assuming air to be at T = `20 ^oc`so accordingly took value from the air standard table. 

8. Then similarly modifying controlDict, 

controlDict:

The controlDict directory controls the time functions such time stepping, deltaT , writeControl and other. So time step modification is done by taking time step as 5e-3 (as it affects the CFL criterion). Such time step should be rationally guessed or one may try with some initial guess and according to expected CFL, if not satisfied or error occurs then may decrease it accordingly.   

 

9. The cavity directory copied from tutorial directory has a "0" folder, which contains two files called as "p" (pressure) and "U" (velocity). These are the initial conditions for our problem. We need to modify these initial conditions according to our flow problem.

Initial condition:

For velocity:

We assume uniform velocity at the inlet condition giving it as U = (0.1 0 0) (U = 0.1 m/sec in x-direction and "0" in y-z direction) (Dirichlet boundary condition). Outlet will be kept as zeroGradient at initial state but will be calculated as function of time step advancement (Neumann boundary condition). Front and back patch will be kept as empty (No velocity field is calculated at this patch). Finally the noslipwall (explained in previously) is kept at "noSlip" condition.

 

For pressure:

However, pressure at inlet is stated to be as zeroGradient (Neumann B.C.) and for outlet pressure is explicitly stated as "uniform 0" (at 1 atm) (Dirichlet B.C.). Front and back patches are kept empty same as of velocity. And for the no slip wall boundary, pressure is defined to be as zeroGradient ( Neumann B.C.) 

 

10. Then we will first try to visualize the mesh in Paraview, to execute this, we first have to be in cavity_test directory, then type in blockMesh and "polyMesh" file be created inside the constant directory. Accordingly, type in "paraFoam" and click apply, then mesh can be visulaized in the Paraview, 

Note: More of cells are alloted to the increased area to capture the physics.

 

I. Case 1: Grading factor (1):

In this case, grading factor is kept to be as one (1) for all the x,y and z dimensions.

Now, we will process this problem in Paraview, with all the boundary, initial and time conditions. To process further, we must recall the solver used to solver our problem, which is nothing but "icoFoam".

Steps:

a) Make sure cavity_test directory is opened, here we have three sub-directories named as 0, constant and system.

b) Type in "icoFoam" and results will get printed on the terminal window.

c) Look at the CFL number changing with iterations calculating mean and max CFL, this approach will provide intuition about providing an appropriate time-step.

d) Note that icoFoam is an implicit solver, which means it can handle larger CFL, but it is good practice to limit CFL as low as possible.

e) After plotting all the results on the terminal window, type in "paraFoam", which will take us to Paraview and results can be visualized.

f) Click "apply" in the left hand side corner of the Paraview interface, now the geometry is visible then click on to surface with edges and U to visualize the visualize the geometry ready for the execution.

g) Now start the simulation to visualize results by clicking on the start button.

As a part of result, we will plot velocity contour,

Velocity contour:

The velocity contour plot is resulted showing the velocity distribution, flow separation, boundary layer formation, free-stream velocities, recirculation zone and much more physics.

 

Simulation:

 

 Plot of velocity along y-dimension at x = 0.085 (Grading factor (1 1 1)):

 

Case 2: Using grading factor of 0.2

In case 2, grading factor helps to stretch the mesh along one direction. The mesh grading is given rationally to capture the physics at certain areas. As per the conditions of given problem, grading is required at the step wall, bottom wall, top wall and at the no slip wall. As more of physics is concerned at this locations.

blockMeshedit:

Mesh grading is given as a factor of 0.2 or 20, which same but changes as per the orientation of the domain (i.e position of starting and ending cell).

Graded mesh:

Graded mesh can be visualize in the Paraview interface and we can see that more of finer mesh is provided at the step wall, because it is location at which more of physics must get captured. Similarly, at the top and bottom wall boundary layer formation is captured.

Note: More of finer mesh can be provided to backward facing step at some of locations to study pressure and velocity variation along the domain.

 

 Plot of velocity along y-dimension at x = 0.085 (Grading factor 0.2):

 

Comparison of both the cases on the basis of grading factor:

As we know that, more the finer mesh, more is accurate result (least residuals in the approxiation). The velocity profile after the step is measured in both the cases, as per the previous statement, the case 2 (with grading) produced more accurate results than that of case1 (no grading). Lets have a look at the velocity profile measured along y distance at x = 0.085m. Additionally, if the mesh is not fine enough at location where is fluid behaviour is more of concerned, it will give average results because it has less elements to capture the physics at each node or at each element volume. So in order to get most accurate results, it is required to have sufficiently finer mesh size to capture the flow properties of the fluid.

By comparing we have a intuition that mesh size or cell size will have significant consequence on our results. At x = 0.085, the incoming fluid passed the way of boundary layer separation. So from the bottom of the enlarged s/c (0.085,0) to (0.085,0.02), velocity will increase in the y-dimension. However, after undergoing flow separtion, the velocity increament in the y-direction won't be uniform, it should have patterns of variations along the y-direction. The velocity is recirculating in this section after flow separation, hence it cannot have uniform increament in the y-direction.Moreover, the velocity distribution at the top wall must show variation in velocity because of the presence of boundary layer, and appropriate velocity distribution is shown by case 2. Since the velocity distribution shown resulted by case 2 at x = 0.085m shows more accurate results than that of the case 1.

We can call this study as grid dependency test, as the flow field variables change significantly as a function of the mesh or grid size.

 

Case 3: 

In this case, finer mesh is used to capture more approximate physics for understanding the fluid behaviour over the backward facing step.

blockMeshedit:

jaKJ'

Velocity contour:

Velocity plot:

A] Inlet velocity along Y-distance:

Discussion:

In this plot, velocity is calculated along the Y-distance to understand some fundamental physics. The resulted plot is parabolic in nature, which certainly show the effect of boundary layer at the wall surface. Here 0 and 0.01 are the wall surfaces, where velocity is zero (No slip condition). As we move in positive Y direction, there is increase in velocity in parabolic manner. The maximum velocity is 0.15 which is nothing but the freestream velocity which is unaffected by the wall shear stress. 

 

B] Outlet velocity along Y distance:

 

Discussion:

Similarly as explained above there is also a parabolic variation in velocity profile at the outlet condition. However the maximum velocity reached is 0.075 m/s. There is decrease in the freestream velocity and this can be explained by the mass conservation principle. As mass is conserved in our control volume. We can state that 

Rate of change of mass  = Inflow - Outflow, as the rate of change of mass is equal to zero.

`A_1V_1 = A_2V_2`

Since there is increase in area, the velocity has to be reduced to compensate for rate of change of mass inside the domain. 

 

Velocity contour plot using glyphs:                                                                      C] Velocity distribution along Y distance at x = 0.085m:

        

 

 

Discussion:

The velocity distribution is plotted after fluid overcomes the backward step. At this condition, more fluid physics can be observed because there is sudden disruption in uniform fluid motion. Let's understand the physics at this state.

Initially, as we saw above, the velocity is following nature at the inlet condition. And there is sudden disruption in the geometry (start of the backward step) due to increased cross-section. This sudden enlargement in the cross-section breaks the boundary layer at the step corner. The significance of wall shear stress gets zero. The velocity of fluid layer adjacent to wall surface stays no longer zero. It has to change as its contact with the wall surface is disturbed. Now, the velocity of the bottommost fluid layer has two choices: 1) To move further along with flow in positive x-direction or 2) To flow in the negative x-direction. 

Here at this condition, the bottomost fluid layer will try to flow in negative x-direction . This is because of some unusual phenomenon called as adverse pressure gradient. Before that, let's understand that the pressure gradient variation in a fluid flow. Generally for any fluid to flow from one end to another, it is required to have one end at high pressure and another at low pressure. Then the fluid will flow from high pressure end to low pressure end. In this process of fluid flow, we can state that, as the fluid is flowing from high pressure end to low pressure end, at each and every step the pressure difference (pressure gradient, `nablap`) is getting decreased between the both ends. As the fluid approaches to the end condition, the resulting pressure difference (`nablap`) becomes very less as compared to the pressure difference (`nablap`) at the inlet condition. So at this condition, pressure gradient is favourable pressure gradient.

However, there is another phenomenon called as "Adverse pressure gradient". In an adverse pressure gradient, the static pressure of the fluid increases in the direction of flow. As the pressure increases, there is also an increament in the potential energy which leads to reduction in kinetic energy or deacceleration of the fluid element.

            

Fig: Transformations in boundary layer under the influence of pressure gradient

Annotations for boundary layer transformation:

1-2: Normal boundary layer formation along the fluid flow.

3: Increase in static pressure which change the boundary layer (Flow separation)

4: Reversed flow (Negative velocity gradient) (Recirculation region)

Fig: Flow characteristics in backward facing step

Under such conditions, the flow is struggling and the local velocity gradient approaches zero. At this point the flow detaches from the surface and this detachment is referred to as flow separation. This flow detachment is generally caused by sudden geometric changes which breaks the boundary layer at the step, which also causes flow separation. And as the momentum of the fluid is so large that the flow cannot remain attached to the surface of the body.

Beyond this separation point, a reversed flow occurs near the surface and the velocity gradient becomes negative. Generally the bottommost fluid layer with less velocity will try to move in negative x-direction. And there is formation of recirculation region at step zone. This recirculation makes some disturbances in the fluid motion at step height and corner eddies at sharp angles (corner of the step height). And these disturbances has nothing to do with the turbulence as the inlet velocity is way smaller to create turbulence. However after some time this recirculating fluid join the streamlines to make the fluid stream. The point at which these fluid streams rejoin together is called as reattachment zone. Finally after reattachment the fluid streams flow together in steady conditions towards the outlet.

 

The pressure distribution along the backward facing step in x-direction is resulted as follows:

As per the above discussion, the pressure gradient must continuously decrease as function of flow direction. In this scenario, this flow property is followed just before the step and after the step encounter the pressure gradient significantly increases. This is because of the adverse pressure gradient phenomenon. The increment can be seen to a certain length (x -direction) and after the reattachment again the pressure gradient (favourable) is dropping out till the exit.

 

D] Velocity distribution after reattachment:

 

Discussion:

In this plot, we can see the normalized fluid velocity flowing to the exit of the domain. The fluid velocity is parabolic in nature, whcih shows the steady flow of the fluid by following boundary layer characteristics at the wall surface and unaffected freestream velocity.

 

Simulation:

 

Conclusion:

Initially backward facing step was analyzed using simple grading of 1 and with grading of 0.2. By which an observation was resulted that, using mesh grading factor of 0.2, the control volumes are made smaller to approximate more accurate physics at the step height and boundary walls. Then case 3 study was done to understand the physics of fluid flow when it encounters a sudden disruption in geometry. Mass conservation principle was used to understand the compensation of velocity due to increased in area. Moreover the effect of adverse pressure gradient and boundary layer separation significantly affects the fluid flow resulting in recirculation of fluid and formation of eddies at the step height. These distrubances are normalized along the flow and reattachment makes everything uniform and steady till the fluid approaches trailing end of the domain.

 

References:

1. https://cfd.direct/openfoam/user-guide/

2. https://www.cfd-online.com/Forums/openfoam-solving/

3. https://www.simscale.com/docs/validation-cases/flow-reattachment-flow-over-a-backward-facing-step/#:~:text=The%20phenomenon%20of%20flow%20reattachment,Coefficient%20of%20Pressure

4. https://aip.scitation.org/doi/abs/10.1063/1.5141565

5. https://en.wikipedia.org/wiki/Adverse_pressure_gradient