Simulation of a Backward-Facing Step in OpenFOAM

SIMULATION OF A BACKWARD-FACING-STEP IN OpenFOAM

 

AIM:

1. Simulate an incompressible-laminar-viscious flow through the backward facing step geometry.
2. Perform transient simulation; choose the solver based on the described physics of flow in problem statement. 
3. Explain simulation process and way of setting up the problem statement in openFOAM and also, explain every changes thats made inside the reference case files obtained from OpenFOAM tutorials. 
4. Next, perform two case study, firstly, form simulation of the flow without using any grading factor(i.e. GF=1). Secondly, simulate the flow with grading factor of 0.2. The cells should be finer near the walls. 

 

OBJECTIVE: 

1. Open OpenFOAM tutorials & find the solver to solve incompressible-laminar-viscious flow probelm. 
2. Next step is to copy the required file in solver and paste it in OpenFOAM's run folder so that the official file of OpenFOAM remains untouched. 
3. Find out the block-mesh script file in order to create domain space and provide mesh specification.
4. Find 'controlDict' file to edit start and end time for the simulation. 
5. In the initial folder '0' edit initial parameter for velocity and pressure according to faces of the given domain.
6. Then, use the 'blockMesh' command to generate mesh and chcekMesh to verify it.
7. After that, Open the simulation file in 'paraView' to plot the velocity profile at point 0.085m from the inlet.
8. Also, Repeat the above procedure for the grading factor 0.2 to analyse the differences in results. 
9. Write a detailed report explaining all the command and codes used in the project. 

 

PROBLEM STATEMENT: 

Simulate an incompressible-laminar-viscious flow through the backward-facing step geometry & perform two case studies as described below and compare the results. 

CASE 1: Simulate the flow without using any grading factor(i.e., GF=1). 

CASE 2: Simulate the flow with grading factor of 0.2. The cells should be finer near the walls (including the step wall). 

 

Mesh Specificaitons:

1. Number of cells along the X-direction(longer dimension) = 200. 
2. Number of cells along the Y-direction = 20

 

THEORY: 

A flow is classified as being compressible or incompressible, depending on the level of variation of density during the flow. Incompressibility is an approximation, in which the flow is said to be incompressible if the density remains nearly constant over the course of its motion. For a density to remain constant, the control volume `V = m/rho`has to remain constant. Even though the pressure changes, the density will be constant for an incompressible flow as the variation of density due to pressure change is negligible.  

                                                 

 

Some flows are smooth and orderly while others are rather chaotic. The highly ordered fluid motion characterized by smooth layers of fluid is called laminar. The flow of high-viscosity fluids such as oils at low velocities is typically laminar. 

                     

 This picture shows the difference between turbulent and laminar flow from a faucet. The left side of the picture shows the turbulent flow, and the right side shows the laminar flow. 

Now, When two fluid layers move relative to each other, a friction force develops between them and the slower layers tries to slow down the faster layer. This internal resistance to flow is quantified by the fluid property viscosity, which is a measure of internal stickiness of the fluid. Viscosity is caused by cohesive forces between the molecules in liquids and by molecular collisions in gases. There is no fluid with zero viscosity, and thus all flui flows involve  viscous effects to some degree. Flows in which the frictional effects are significant are called viscous flows. 

                     

Across any length of pipe there will be a drop in pressure due to frictional shear acting within the fluids . 

`Deltap = Punderset(IN)`-`Punderset(OUT)`

The pressure drop in turbulent flow in much larger than in laminar flow. 

Here, we essentially have to model the fluid domain around the blockmesh as a mesh of discrete elements, define boundary conditions and fluid properties and apply an appropirate assessment technique to find a solution. 

 

Firstly, to begin with the simulation process we're opening the terminal window in Linux OS by pressing 'Ctrl+Alt+t' buttons or by simply right clicking on the desktop window and selecting the option 'Open Terminal'. 

We'll see the terminal window shown as follows: 

Now, in order to simulate our problem statement in OpenFOAM, we're going to use the 'tutorials' present in OpenFOAM by typing 'tut' to open OpenFOAM tutorial directory. Then, we can type 'ls' to view all the tutorials present inside this directory.  

 We have to solve incompressible flow simulation for this problem hence we're going to select the 'incompressible' directory by simply typing 'cd incompressible' and then type 'ls' again to view various folders inside the incompressible directory.  

We're now able to see a lot of solvers inside the incompressible folder. For this particular solver we're needed to solve for laminar-viscious flow. Hence we're choosing the 'icoFoam' solver. icoFoam is solver in OpenFoam which is used as a transient solver for incompressible, laminar flow of newtonian fluids. Newtonian fluid is that type of fluid in which the viscious stress arising from its flow, at every point, are linearly correlated to the local strain rate- the rate of change of its deformation over time. Newtonian fluids are the simplest mathematical models of fluids that account for viscosity. 

Hence, we'll select icoFoam folder by typing 'cd icoFoam' and then again we'll type 'ls' to explore that folder in the terminal window. 

 

From the above image, we can see that there's two options available to choose from. 'cavity' & 'elbow'. They are both part of the icoFoam solver. 

Cavity's section is needed to be used to create the required domain and blockmesh. So, now we are again changing the directory towards cavity's folder. 

Inside the cavity's folder, we're again going to see a cavity folder where we have to define our fluid domain and then generate the blockmesh. 

Since OpenFOAM dosen't let us to make any changes to their original tutorial's section, we're going to make a copy of this folder using the command 'cp -r cavity' & then we have to paste this information inside the OpenFOAM's Run folder by typing 'cp -r cavity $FOAM_RUN'/cavity_test'. 

Now, we are free to make any changes to this cavity_test file and it'll not interfere with the original tutorial's files. 

We're now going to open this cavity_test folder in order to provide specifications for our blockmesh in the system folder. 

From the above image, we can observe that there are  four script files inside the system's folder. 

The blockMesh utility creates parametric meshes with grading and curved edges. blockMesh reads this dictionary, generates the mesh and writes out the mesh data to points and faces, cells and boundary files in the same directory.

We have to create the geometry through blockmesh as shown below:

In order to generate the above geometry, we're considering to create the geometry using five blocks as shown in the image below: 

We're then going to provide vertices for each and every block by referring from the sktech shown below:

 

Now, to create a domain and provide required specifications for our blockmesh we're going to open 'blockMeshDict' using the 'gedit' command. Text Editor (gedit) is the default GUI text editor in the Ubuntu operating system. 

 

When we open this file, at first, we're going to see the header information in the form of comment which cannot be executed.

 After that, we're going to see the file identifier using which, OpenFOAM understands the type of file it is reading which is in this case 'blockMeshDict' and it contains the information in a dictionary format which essentially means in key and its value's format. Here, version is designated as the key and it's value is 2.0 which is why it is called as a key-value pair. 

We won't be making any changes here. 

 

Next, we set units in meter, by converting all type of length into meter using the command:

convertTometers 0.01;

After that we're providing all the respective vertices for each & every blocks. 

Mesh genaration without any grading factor: 

Inside the same blockMesh file, we have to goto the lines where we can see the word 'hex'. That is the line where we specify the faces of the blocks and also the number of cells in X,Y and Z direction. 

 Mesh grading factor is given as: 

Grading Factor = Size of the last cell/ Size of the first cell

After that, we're then specifying the faces as the inlet, outlet and walls as shown in the image below:

The blockMeshDict file is edited as shown in the above image. 

 

After that, we going to change our directory and goto inside the 'cavity_test' folder in order to go inside a folder named as '0' which consist of two script files, one for pressure & one for velocity. 

Using 'gedit' command we're going to open the script file for velocity. We're then going to edit our boundary conditions and set the environmental variables for our patch. 

The environmental booundary condition we're using for our inlet is 'fixedValue'. Along the X-axis the 'fixedValue' is stated as 1 in the global-coordinate system which means that the flow will be flowing towards the X-direction. If X's value was stated as -1 then the flow would tend to go in the opposite direction. 

Now, for the outlet, if we give a fixed value then we need to specify the value which we can't obtain because we cannot predict the values at the outlet from beforehand. Hence, we're going to calculate is from the inner cells of the solution. 

So, the value of the outer node can be given as, `Vunderset(n)` = `2Vunderset(n-1)`-`Vunderset(n-2)`. This is known as the zero-gradient boundary condition. Linear interpolation is being used here. 

Now, for the walls, we know that there exists a 'no-slip' condition which is a special type of environmental boundary condition for walls.

We're then going to copy the entire set of boundary conditions from our velocity script, navigate ourself to the pressure script & paste it inside our pressure script like in the image shown below: 

Now here, for the boundary conditions at inlet, the inlet pressure(`(dP)/(dx)`) can be derived from the inlet velocity by solving through momentum equation. So, as we're getting the value for pressure from the solution itself and hence we are going to use 'zeroGradient' environmental boundary condition for it. 

Now for the outlet, the outside pressure is assumed to be equal to atmospheric pressure. Hence we're using the 'fixedValue' environmental boudary condition for outlet with a vlaue of '0' as we know 0 gauge = 1 bar = 1 atm. 

And finally for the walls we're again using the 'zeroGradient' environmental boundary condition as it can be calculated from the solution. The setup is complete now. 

Now, we're going to change directory and goto our 'cavity_test' folder and run 'blockMesh' command. The reason we're doing this is because when we run our blockmesh command in the terminal, a folder named as 'polymesh' is getting generated inside the 'constant' folder from our folder 'cavity_test' which we are then going to view. We can then enter the 'constant' folder & we'll find that along with 'polyMesh' folder another file is there named as 'transportProperties'. 

We're then going to enter into the 'polyMesh' folder where we'll find files such as 'boundary', 'faces', 'neighbour', 'owner' & 'points'. For every cell, it is going to classify it's own neighbour and owner by itself. 

Now, we're going to use the command 'gedit' and open the script file for boundary. Based on the patches we've created so far, it tells us its type, number of faces and the starting face. We can see the values for inlet, outlet, walls and default values. These value are for OpenFOAM and needs not be edited in any way. 

We're again changing the directory to enter into the 'polyMesh' folder to view the script file for 'transportProperties'. 

All the transport constant can be viewed from here. As 'icoFoam' only solves for pressure and velocity and hence we're supposed to have a contant for momentum diffusivity. 

Here, Kinematic Viscosity is being viewed which can be given as Kinematic Viscosity = Dynamic Viscosity/Density. 

Here the dimensions are in meter-square per second and the value of kinematic viscosity for air can be given as 1e-5 which is a realistic value for room temperature. 

Finally, we're then changing directory to 'cavity_test' folder where we'll then enter into the 'system' folder. Then, we're going to open the script file for 'controlDict'. Then, in order to decide what should be the end time for our solver we're redirecting to our velocity folder(U) where the speed is stated as 1 m/s which is quite high for a geometry like ours. Hence, we're going to modify it by changing it into a value favourable for our geometry. Here, we're stating the value of speed as 0.02 m/s. 

Now, we know that our length is 0.2m and speed is 0.02m/s. So, as Time=Distance/Speed we should get the value of 10 secs for our endtime value but we're then modifying it to 20 secs and we're stating the value of `Deltat` as 1e-2. 

We're then going back to the 'cavity_test' folder and we're going to use the command 'blockMesh' in the terminal to generate our final mesh. It is going to show us the following output in the terminal shown in the images below: 

So here, we can see the faces created by OpenFOAM for inlet, outlet, walls and also OpenFOAM has generated faces for defaultFaces in patch 3. As OpenFOAM is aware of all the local coordinates of the blocks and their orientation in the geometry, it automatially patches it accordingly. 

After that, we're going to run a command called 'checkMesh' which is used in order to check if the mesh which is being generated above is in appropirate condition or not. It shows us the following output from which we can observe that the mesh is deemed to be in correct order. 

 

 Now, we're going to use a command known as 'paraFoam' to open the blockmesh in paraView. So, initially the paraView will open as shown below: 

 We'll click on apply inside the properties toolbar to generate the mesh visually. Then using the toolbar shown below, we'll change the type of surface for the blockMesh to 'Surface with edges' and we'll use a solidcolor to view the mesh. 

 

 

Above image shows the mesh for the grading factor of 0.2. 

Now, we're going to see the profiles for velocity and pressure. 

We can then observe the velocity(U) profile for initial condition from the image given below where we can see that initially the velocity is zero everywhere except at the inlet where the velocity is  m/s. 

  

For pressure at the initial condition, the pressure at outlet is zero and at inlet, the pressure is yet to be found. 

Now, we'll goto the terminal and enter the 'icoFoam' command. 

icoFoam Solver: Used for solving the Navier-Stokes equation for the incompressible laminar flow problems using the PISO algorithm. It can take mesh non-orthogonality into account with successive non-orthogonality iterations. The number of PISO and non-orthogonality corrections are controlled through user input.The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions.

It is used to solve for pressure and velocity in the blockmesh. We're then getting the output at the terminal as follows:

 

From the above image we can observe that since the blockmesh is graded hence for each cell there is a seperate CFL(Courant-Friedrich-Lewy) number. Here, 'mean' tells us the average CFL number from the overall simulation & 'max' tells us the smallest CFL number. If our cell is very small in size than our value for max would be very high. Hence, we have to keep in mind that if our cell is very small than our timesteps has to be relatively lower such that the ratio is maintained. 

IcoFOAM is an implicit solver and hence even if our CFL number becomes higher in size chances are it'll still not crash but once it goes beyond a certain limit then with the increase in value of CFL, our simulation is going to be very much inaccurate and hence knowing that OpenFOAM might crash it since no useful data is being generated from that simulation. For a cell with a CFL value of 1 will bring higher degree of accuracy as compared to a cell with a CFL value of 2. 

We have both convection and diffusion here and hence if we increase our kinematic viscosity, the CFL's value will also increase. Instead of thermal diffusion, we got moment diffusion here. CFL is given by addition of CFL convection and CFL diffusion which has to been less than 1 to 2.

Now, if we type 'ls' command is  the terminal we'll be able to see all the intermediate solution files as shown below:

Now, we're going to refresh the blockmesh in paraView from the properties  so that all the values will get updated for the velocity and pressure generated through the solver. 

 

We're going to now simulate the results for two case studies where the grading factor will vary with one another. 

 

 

CASE 1: HERE THE GRADING FACTOR IS 1

 

 

MESH GRID FOR GRADING FACTOR OF 1:

 

MAGNIFIED VIEW OF THE MESH GRID: 

We can observe from the above image that the mesh structure is fine. 

 

VELOCITY PROFILE FOR OUR DOMAIN: 

 

We can use the magnifying glass from the toolbar to enhance our view of the flow seperation:

 

PRESSURE PROFILE FOR OUR DOMAIN: 

 

MAGNIFIED VIEW OF THE PRESSURE PROFILE: 

The Plot Over Line filter samples the data set attributes of the current data set at the points along a line. The values of the point-centered variables along that line will be displayed in an XY Plot. Now, we are going the use the plot over line command in paraView in order to view our plot at respective points for our velocity profile. 

When we click on the 'plot over line' option it gives us a line with two end-points as shown below: 

   

We're measuring the velocity profiles for before and after the step wall is crossed by the flow: 

Modifications made in Properties toolbar to calculate velocity profile at 0.078:  

Modifications made in Properties toolbar to calculate velocity profile at 0.085:  

 

Next, by clicking 'apply' we'll generate the respective plots at grading factor as 1: 

 

CASE 2: HERE THE GRADING FACTOR IS 0.2 

 

MESH GRID FOR GRADING FACTOR OF 0.2: 

 

MAGNIFIED VIEW OF THE MESH GRID: 

We can observe from the above image that the mesh structure is coarse. 

 

VELOCITY PROFILE FOR OUR DOMAIN: 

 

We can use the magnifying glass from the toolbar to enhance our view of the flow seperation:

 

PRESSURE PROFILE FOR OUR DOMAIN: 

 

MAGNIFIED VIEW OF THE PRESSURE PROFILE: 

We're measuring the velocity profiles for before and after the step wall is crossed by the flow: 

 

RESULTS: 

We have completed the necessary plotting for both case studies, one for without grading factor(GF=1) & one for grading factor of 0.2. From the plots we can see how our velocity profiles changes. 

HERE, 1ST PLOT IS FOR GF = 0.2  & 2ND PLOT IS FOR GF = 1.

1. The plot without the GF(grading factor=1) seems to be more accurate & smoother than the plot with GF=0.2. Hence, we can observe that when the mesh is not uniform, it can disrupt our output results to some extent. Grading factor actually helps us to adjust the sizes of the respective cells from the initial section of the block to the final section of that block. 

2. Mesh-grid needs to be finer in order to capture the flow properties with better accuracy as shown in the plot with GF=1 as the mesh is uniform at all section. 

3. With the grading factor of 0.2 we can notice that the flow quantities were captured well at grid points nearer to the wall where the mesh is coarse and hence resolution of the velocity with mesh-grading have more accuracy compared to when there's no grading. 

4. We can observe from the simulation that when the fluid is going to from inlet towards outlet, the area increases and as the area increases, the velocity profile decreses.