Steady laminar flow simulation through a backward facing step

Introduction:

In the geometry shown in the figure below, the area of intrest is the step which is 0.08 metres from the inlet. As there is sudden change in cross section, fluid particles rotate and eddy formation happens.

When the fluid has low enough viscosity, eddies form right after the fluid crosses the step. This occurs because of seperation of boundary layer, where negative velocity gradient happens at the step and eddies forms.

Boundary layer separation occurs due to adverse pressure gradient, i.e. when there is a sudden increase in pressure, the velocity head comes down and this leads to negative velocity gradient, which causes eddies.

In this simulation we use 'icoFoam' solver to run the simulation. icoFoam is a transient solver, meaning that it time marches the Navier Stokes equation. It is used in case of incompressible, transient states and for Newtonian fluids.

Process:

1. Study the geometry and fix the location of the origin, to define the geomtry
2. Copy the predefined 'cavity' folder from Openfoam tutorials to '$FOAM_RUN'
3. In the copied cavity folder, open 'blockMeshDict' file using 'gedit' command, which is located inside the 'system' folder
4. Define the origin point and define all other vertices with respect to the origin (which here acts as a global coordinate system)
5. Now define each 'block' and divide the blocks by entering the desired number of cells in x, y, z directions, also a desired grading factor is given to refine the mesh
6. After defining the blocks, define the 'type' and enter the vertices of each face in a way that the normal vector of the plane is oriented outwards
7. Save the blockMeshDict file, come back one step from the directory, and enter the command 'blockMesh' to generate the mesh
8. Go back to the system folder, open 'controlDict' file, give the time span of simulation and the elemental size of time. Save the file and close
9. Come back a step out and go to the '0' folder (intial conditions - time = 0), open the pressure file (p) and enter the inlet and outlet pressure values at time zero, save and close the file. In a similar way, open velocity file (U) and enter the inlet and outlet values at time 0.
10. Come out of the 0 folder and enter the command 'icoFoam' to start the simulation
11. After the simulation is done, post process it to get the desired data by using Paraview. To open and view the results using Paraview, enter the command 'ParaFoam'

NOTE:

• The values of grading factor, time step and velocity are to be chosen by considering the CFL number, such that the CFL number should be less than one, for a stable solution
• In defining the blocks, note that the direction of vector joining first and second points acts as x axis, direction of vector joining second and third points acts as y axis and the direction of vector joining first and fifth points acts as z axis
• In this simulation the value of kinematic viscosity, which also plays a prominent role, the value of which is 0.01 (here), is not altered

BlockMeshDict file:

In FVM as the elements are volumes, we have finite volumes and hence blocks. The file which contains information and through which we can create, control the number, meshing of these blocks is the blockMeshDict file.

In this case, the geometry is divided in to six blocks.

Note that the analysis is done for different values of grading factor, 0.2, 0.5, 0.8. The file remains same for all other parameters.

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  7
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 0.01;

vertices
(
    (0 0 0)
    (8 0 0)
    (8 0.5 0)
    (0 0.5 0)
    (8 1 0)
    (0 1 0)
    (20 1 0)
    (20 0.5 0)
    (20 0 0)
    (20 -0.5 0)
    (8 -0.5 0)
    (8 -1 0)
    (20 -1 0)
    (0 0 0.1)
    (8 0 0.1)
    (8 0.5 0.1)
    (0 0.5 0.1)
    (8 1 0.1)
    (0 1 0.1)
    (20 1 0.1)
    (20 0.5 0.1)
    (20 0 0.1)
    (20 -0.5 0.1)
    (8 -0.5 0.1)
    (8 -1 0.1)
    (20 -1 0.1)   
);

blocks
(
    hex (0 1 2 3 13 14 15 16) (200 10 1) simpleGrading (0.2 5 1)
    hex (3 2 4 5 16 15 17 18) (200 10 1) simpleGrading (0.2 0.2 1)
    hex (2 7 6 4 15 20 19 17) (200 10 1) simpleGrading (5 0.2 1)
    hex (1 8 7 2 14 21 20 15) (200 10 1) simpleGrading (5 5 1)
    hex (10 9 8 1 23 22 21 14) (200 10 1) simpleGrading (5 1 1)
    hex (11 12 9 10 24 25 22 23) (200 10 1) simpleGrading (5 5 1)
);

edges
(
);

boundary
(
	inlet
	{
		type patch;
		faces
		(
			(0 13 16 3)
			(3 16 18 5)
		);
	}
	
	outlet
	{
		type patch;
		faces
		(
			(20 7 6 19)
			(21 8 7 20)
			(22 9 8 21)
			(25 12 9 22)
		);
	}
	
	noslip
	{
		type wall;
		faces
		(
			(13 0 1 14)
			(14 1 10 23)
			(23 10 11 24)
			(24 11 12 25)
			(17 4 5 18)
			(19 6 4 17)
		);
	}
	
	frontandback
	{
		type empty;
		faces
		(
			(13 14 15 16)
			(16 15 17 18)
			(15 20 19 17)
			(14 21 20 15)
			(23 22 21 14)
			(24 25 22 23)
			(0 3 2 1)
			(3 5 4 2)
			(2 4 6 7)
			(1 2 7 8)
			(10 1 8 9)
			(11 10 9 12)
		);
	}
   
);

mergePatchPairs
(
);

// ************************************************************************* //

1. Faces after defining their type, note that red colour denotes patch type, blue denotes wall type and the grey denotes empty type.

Controldict file:

The controldict file contains the information about the run time, time step, type of application. In this case, the run time is taken as 0.6 seconds and the time step value if 1e-5 (so that Courant number safely lies less than one).

The Courant number will be same for all the grading factors 0.2, 0.5, 0.8. Although the element size changes (dx and dy) with grading factor, as the value of dt is very less, for all of the three values of grading factors, CFL number lies lesser than one. The control dict file will be as follows,

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  7
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.6;

deltaT          0.00001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //

Initial conditions:

Initial conditions define the state of the system right before the start of simulation. The parameters involved in this case are velocity and pressure. The initial conditions here are such that, at the inlet one need to be 'fixed' (Dirichlet) and the other 'float' (Neumann). And it is same for the outlet too.

It is because, at the inlet and at the outlet, the characteristic lines are in two directions, i.e they are coming out and going in. The lines which come out tend to resist the flow and hence a gradient exists, and so the intial condition also need to be floating i.e Neumann condition. For the 'going in' condition, a Dirichilet condition is given.

The practice here is that, at inlet velocity is fixed (as we are intrested in how the velocity behaves along the considered geometry) at the inlet and pressure is given a gradient value. And the vice-versa is at the outlet. That is velocity is given a gradient value and pressure, a fixed value.

The intial conditions are also same regardless of the grading factor. It is to be noted that the value of velocity is to be considered carefully, so as to maintain the CFL number less than one. Also, for the simplicity of the problem, the value of presure is considered zero, at the inlet and at the outlet.

The intial conditions in this case are as follows,

A. Pressure:

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  7
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
        
    }

    noslip
    {
        type            zeroGradient;
    }

    frontandback
    {
        type            empty;
    }
    outlet
    {
    	type		fixedValue;
    	value 		uniform 0;
    }
}

// ************************************************************************* //

B. Velocity:

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  7
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (2 0 0);
    }

    noslip
    {
        type            noSlip;
        value 		 uniform (1 0 0);
    }

    frontandback
    {
        type            empty;
    }
    
    outlet
    {
    	type		 zeroGradient;    	    	
    }
}

// ************************************************************************* //

Results:

After running the simulation using 'icoFoam' command, result data is analysed using post processor, Paraview, by entering the command 'ParaFoam'.

1. Velocity profile after the system has reached steady state:

2. Pressure profile after the system has reached steady state:

Inferences:

• The eddies formation can be seen clearly right after the step in the velocity figure. Later the velocity has reached a constant value after that, with the constant development in boundary layer
• In the pressure figure too, the value is high at the corners and and slowly started decreasing and reached zero later along the longitudinal length. The decreasing pressure gradient is because to compensate for the loss in momentum due to wall shear stress

Notice that in the blockMeshDict file, a reciprocal of grading factor is given in case when a bigger cell is needed. For example, 0.2 and 1/0.2 i.e 5 is maintained.

1. For a grading factor of 0.2

Grading factor relates the starting cell and ending cell by equating the length of starting cell to the grading factor times the length of last cell. We can roughly calculate the value of CFL number considering the least value of the size of the element after applying grading factor, which comes out 8e-5 (dx value), and as time stpe is 1e-5 and velocity is 2, in x - direction, the CFL number will be ~0.025, which is safe.

The meshed geometry with a grading factor of 0.2 will be as follows,

• The mesh here is more concentrated near walls and at the step. This is to have more data from closely packed cells is used near walls and at the step so that more 'finer' results can be obtained near them

NOTE: The velocity profile is plotted after the system has reached steady state and the 'length' in abscissa denotes transverse length of the geometry at 0.085 metres from the inlet.

A. Velocity profile at 0.085 metres from the inlet:

B. Pressure variation along the longitudinal length:

C. Variation of magnitude of velocity variation along the longitudinal length:

D. Velocity in x, y, z directions:

Inferences:

• As there are eddy formation right after the step (at the middle of the height), a gradient is beign developed as the height is increased. As we know that a positive velocity gradient represents the growth of boundary layer, it can be inferred that the boundary layer formation has started from 0.085 metres downstream to the inlet
• The minimum velocity is at the walls, at the bottom and the at the top. (In case of laminar flow, the velocity profile will be parabolic at a transverse section with minimum (~zero) at the walls). This is because of the wall shear stress, which is because of rate of shear strain (Newton's law of viscosity), the velocity will be zero at the walls
• Also there exists a decreasing pressure gradient along the longitudinal length, which is because to compensate for the loss in momentum due to wall shear stress. The increment during the start of graph is because of formation of eddies. Or on the contrary, the eddy formation is because of positive pressure gradient, which causes boundary layer separation
• In the velocity variation along the length, the value is highest near 0.1 metres, which is halfway along the total length, because after the eddy formation at 0.08 metres, the boundary layer starts developing again and the positive velocity gradient exists
• In the plot showing velocity in different directions, it can be clearly seen how the velocity in transverse direction is having negative values and reaching zero at almost 0.12 metres from the beginning. This negative values are because of eddy formation, which is formed for the reasons mentioned in above points

2. For a grading factor of 0.5:

The meshed geometry is as follows,

• The mesh is finer at the walls and also at the step, so that more data is taken from the 'nearest' poins to refine the data. And ofcourse, in case of grading factor 0.2, where the mesh is finer, grading factor 0.5 gives the date that may vary very little

A. Velocity variation at 0.085 metres from the inlet:

Note that length here means transverse length.

B. Pressure variation along the longitudinal length:

C. Variation of magnitude of velocity along the longitudinal length:

D. Velocity variation in x, y, z directions:

Inferences:

• The variation in velocity is for the similar reasons mentioned above, but the values will be a little far from the corresponding parameters when the grading factor is 0.2

3. For a grading factor of 0.8:

The meshed geometry is as follows,

• The mesh concentration is visibly less finer than with grading factors 0.5 and 02.

A. Velocity variation at 0.085 metres from the inlet:

B. Pressure variation along the longitudinal length:

C. Variation of magnitude of velocity along the longitudinal length:

D. Velocity variation in x, y, z directions:

Inferences:

• The inferences are same as in above mentioned points

Conclusion:

• When exported the data of variation of velocity (our parameter of intrest is velocity), there is no point exactly located at 0.085 metres (from the inlet). So the velocity value from the point right next to 0.085 metres, which is at 0.085023 metres, is taken. Similarly, the velocity value at the point right before 0.085 metres, which is at 0.084822 metres, is taken
• Using these data from points before and after, the value at 0.085 metres from the inlet is linearly interpolated for different mesh grading factors, which are as follows

For GF 0.2, magnitude of velocity at 0.085 metres is 1.0425228 metres/second.

For GF 0.5, magnitude of velocity at 0.085 metres is 1.0495545 metres/second.

For GF 0.8, magnitude of velocity at 0.085 metres is 1.0524793 metres/second.

• It can be safely concluded from the values above that, the value refinement happens as the mesh becomes finer, although it may increase the load on the machcine. So, based on our requirement and the accuracy needed, we can choose the grading factor value