Simulation of a Laminar Flow through a Pipe in OpenFoam
Week 11 Report
- Aim: To simulate an axi-symmetric laminar flow by applying the wedge boundary condition for a flow of Reynolds Number of 2100.
- Theoretical Introduction:
In fluid dynamics, laminar flow is characterized by fluid particles following smooth paths in layers, with each layer moving smoothly past the adjacent layers with little or no mixing. At low velocities, the fluid tends to flow without lateral mixing, and adjacent layers slide past one another like playing cards. There are no cross-currents perpendicular to the direction of flow, nor eddies or swirls of fluids. In laminar flow, the fluid's motion is very orderly with particles close to a solid surface moving in straight lines parallel to that surface.
Difference between Laminar Flow and Turbulent Flow
When a fluid is flowing through a closed channel such as a pipe or between two flat plates, either of two types of flow may occur depending on the velocity and viscosity of the fluid: laminar flow or turbulent flow. Laminar flow occurs at lower velocities, below a threshold at which the flow becomes turbulent. The velocity is determined by a dimensionless parameter characterizing the flow called the Reynolds number, which also depends on the viscosity and density of the fluid and dimensions of the channel.
Reynolds Number:
The Reynolds number (Re) helps predict flow patterns in different fluid flow situations. At low Reynolds numbers, flows tend to be dominated by laminar (sheet-like) flow, while at high Reynolds numbers, flows tend to be turbulent. The turbulence results from differences in the fluid's speed and direction, which may sometimes intersect or even move counter to the flow's overall direction (eddy currents). These eddy currents begin to churn the flow, using up energy in the process, which for liquids increases the chances of cavitation. Reynolds numbers are an important dimensionless quantity in fluid mechanics.
The Reynolds number has wide applications, ranging from liquid flow in a pipe to the passage of air over an aircraft wing. It is used to predict the transition from laminar to turbulent flow and is used in the scaling of similar but different-sized flow situations, such as between an aircraft model in a wind tunnel and the full-size version. The predictions of the onset of turbulence and the ability to calculate scaling effects can be used to help predict fluid behaviour on a larger scale, such as in local or global air or water movement and thereby the associated meteorological and climatological effects.
Axisymmetric Flow:
In hydrodynamics and fluid mechanics, a flow in which the streamlines are symmetrically located around an axis.Every longitudinal plane through the axis would exhibit the same streamline pattern. - Procedure:
1. The first step that we need to do is write a code that will help us automate the creation of the BlockMeshDict file. We must be able to automate it such that we simply give in the basic inputs for the creation of the blockMesh, such as the diameter of the cylindrical pipe or even the length of the pipe.
2. To begin with the coding, we must start by using the commands of clear all, close all, clc to clear any previously stored values and clear the command window as well.
3. Next, we must give in the necessary inputs that are needed to create the block in the desired shape. They are as following:
%Inputs: l= 1.5; %Length of the pipe. theta = 5.0; %Wedge angle. d = 0.01; %Diameter of the pipe. c = 1.0; %Conversion Factor. r = d/2; %Radius of the Cylinder. nx = 500; %number of Mesh points in x-direction on the block. ny= 25; %Number of Mesh points in the y-direction on the block. nz= 1; %number of Mesh Points in the z-direction on the block. gx = 1; %Grading along the x-direction. gy = 1; %Grading along the y-direction. gz = 1; %Grading along the z-direction.
As we can see, we have given the basic inputs such as the number of mesh points and even the grading of the mesh that is being created.
4. Next, we must print all the inputs in the format of a blockMeshDict file line by line. So, to write it in, we must create a file and give the permission as writing and then go on to write the code in the format of a BlockMeshDict file. WE must be careful when we write it line by line and not miss a semi-colon in it.
Since we are creating the block, the vertices and faces must be created appropriately with respect to their normals. We create the vertices with respect to the wedge angle and the Cylindrical pipe radius using basic trigonometry.
So, because we have already conceptualized and made the block on paper, we can therefore create the vertices and assign the faces with their boundary conditions or the type of face that it is.
Like, for example, the front and the back face of the block are assigned as for the wedge boundary conditions. The top wall has to be assigned the wall condition. We must also create the axis and then assign an empty boundary condition to it as well.
The Main-code is given below:
clear all close all clc %Inputs: l= 1.5; %Length of the pipe. theta = 5.0; %Wedge angle. d = 0.01; %Diameter of the pipe. c = 1.0; %Conversion Factor. r = d/2; %Radius of the Cylinder. nx = 500; %number of Mesh points in x-direction on the block. ny= 25; %Number of Mesh points in the y-direction on the block. nz= 1; %number of Mesh Points in the z-direction on the block. gx = 1; %Grading along the x-direction. gy = 1; %Grading along the y-direction. gz = 1; %Grading along the z-direction. f1 = fopen('BlockMeshDict.txt','w'); fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n'); fprintf(f1,' ========= | n'); fprintf(f1,' \ / F ield | OpenFOAM: The Open Source CFD Toolbox n'); fprintf(f1,' \ / O peration | Website: https://openfoam.org n'); fprintf(f1,' \ / A nd | Version: 8 n'); fprintf(f1,' \/ M anipulation | n'); fprintf(f1,'/*------------------------------------------------------------------*/ n'); fprintf(f1, 'FoamFile n'); fprintf(f1,'{ n'); fprintf(f1,' version 2.0; n'); fprintf(f1,' format ascii; n' ); fprintf(f1,' class dictionary; n'); fprintf(f1,' object blockMeshDict; n'); fprintf(f1,'} n'); fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // n'); fprintf(f1,'convertToMeters %0.1f; n',c); fprintf(f1,'n'); fprintf(f1,'vertices n'); fprintf(f1,'( n'); fprintf(f1,'(0 0 0) n'); fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),0); fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),0); fprintf(f1,'(0 0 %0.3f) n',-l); fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),-l); fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),-l); fprintf(f1,'); n'); fprintf(f1,'n'); fprintf(f1,'blocks n'); fprintf(f1,'( n'); fprintf(f1,'hex (0 3 5 2 0 3 4 1) (%0.0f %0.0f %0.0f) simpleGrading (%0.0f %0.0f %0.0f) n',nx,ny,nz,gx,gy,gz); fprintf(f1,'); n'); fprintf(f1,'n'); fprintf(f1,'edges n'); fprintf(f1,'( n'); fprintf(f1,'); n'); fprintf(f1,'n'); fprintf(f1,'boundary n'); fprintf(f1,'( n'); fprintf(f1,'in n'); fprintf(f1,'{ n'); fprintf(f1,'type patch; n'); fprintf(f1,'faces n'); fprintf(f1,'( n'); fprintf(f1,'(0 1 2 0) n'); fprintf(f1,'); n'); fprintf(f1,'} n'); fprintf(f1,'out n'); fprintf(f1,'{ n'); fprintf(f1,'type patch; n'); fprintf(f1,'faces n'); fprintf(f1,'( n'); fprintf(f1,'(3 5 4 3) n'); fprintf(f1,'); n'); fprintf(f1,'} n'); fprintf(f1,'ft n'); fprintf(f1,'{ n'); fprintf(f1,'type wedge; n'); fprintf(f1,'faces n'); fprintf(f1,'( n'); fprintf(f1,'(0 3 4 1) n'); fprintf(f1,'); n'); fprintf(f1,'} n'); fprintf(f1,'bk n'); fprintf(f1,'{ n'); fprintf(f1,'type wedge; n'); fprintf(f1,'faces n'); fprintf(f1,'( n'); fprintf(f1,'(3 0 2 5) n'); fprintf(f1,'); n'); fprintf(f1,'} n'); fprintf(f1,'top n'); fprintf(f1,'{ n'); fprintf(f1,'type wall; n'); fprintf(f1,'faces n'); fprintf(f1,'( n'); fprintf(f1,'(5 2 1 4) n'); fprintf(f1,'); n'); fprintf(f1,'} n'); fprintf(f1,'axis n'); fprintf(f1,'{ n'); fprintf(f1,'type empty; n'); fprintf(f1,'faces n'); fprintf(f1,'( n'); fprintf(f1,'(0 3 3 0) n'); fprintf(f1,'); n'); fprintf(f1, '} n'); fprintf(f1,'); n'); fprintf(f1,'n'); fprintf(f1,'mergePatchPairs n'); fprintf(f1,'( n'); fprintf(f1,'); n'); fprintf(f1,'n'); fprintf(f1,'// ************************************************************************* // n'); %END
5. The next step is to check if the Block Mesh is created properly or not. To check if it is created appropriately, we must copy the blockMeshDict file's output that we have created using MATLAB and then paste it in the blockMeshDict file in OpenFoam. This file is in the system folder. Once this is done, we must then run the blockMesh command. After the command is run, we then check for the mesh in the post-processor, known as Paraview. To access this, we must type in the paraFoam command. Here, we can check if the mesh is okay.
Another command would be of checkMesh. This command checks the mesh for the basic quality criteria, such as minimum mesh length, the aspect ratio etc.
6. After the Mesh is made, we need to alter the Velocity and Pressure quantities' Boundary conditions. To change the Boundary Conditions, we must use the 0 file in the Run folder. In this file, we can change the Boundary Conditions for Velocity and Pressure.
For all the Boundaries that we have created, we must give in some condition for the quantities of Velocity and Pressure. We need to be careful and thoughtful when giving in the Boundary Conditions for these quantities as it is according to the Conditions that we mention, the entire problem will change.
The Boundary Conditions are given as follows:
a. Pressure:
/*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 8 \/ M anipulation | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class volScalarField; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -2 0 0 0 0]; internalField uniform 0; boundaryField { in { type zeroGradient; } out { type fixedValue; value uniform 0; } ft { type wedge; } bk { type wedge; } top { type zeroGradient; } axis { type empty; } } // ************************************************************************* //
b. Velocity:
/*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 8 \/ 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 { in { type fixedValue; value uniform (0 0 -0.21); } out { type zeroGradient; } ft { type wedge; } bk { type wedge; } top { type noSlip; } axis { type empty; } } // ************************************************************************* //
7. After thee Boundary Conditions were specified, the next step is to give in the Timestep and the Start and End time for the simulation. The time step will control the accuracy and the Courant number of the entire simulation. To access the controlDict file, we must use the systems file inside the run folder. Here, we must use the command of gedit to edit the controlDict file. The following were the values entered in the controlDict file:
/*--------------------------------*- C++ -*----------------------------------* ========= | \ / F ield | OpenFOAM: The Open Source CFD Toolbox \ / O peration | Website: https://openfoam.org \ / A nd | Version: 8 \/ M anipulation | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application icoFoam; startFrom 0; startTime 0; stopAt endTime; endTime 50; deltaT 0.001; writeControl timeStep; writeInterval 20; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; // ************************************************************************* //
8. The next step is to run the simulation in the run folder. To run the simulation, we must use the name of the solver in which we want to run it. It is icoFoam in our case as this solver solves for incompressible flow. We need to run this solver in the run folder and then check the results in the post-processor, known as ParaView.
After the simulation is done, we need to open it in the post-processor. We then check the simulation time and the results and also plot certain graphs to study the simulation. - Results:
After the Simulation is run, the next step is to study it and then derive an appropriate understanding from it.
a. To validate the Hydrodynamic length with the Numerical Result.
So, In fluid dynamics, the entrance length is the distance a flow travels after entering a pipe before the flow becomes fully developed. Entrance length refers to the length of the entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer.
Here, numerically, we know that the velocity profile will be different before and after the Hydrodynamic entry length. After the Hydrodynamic entry length, the Velocity profile will be the final parabolic profile that it must have in a fully developed region.
Velocity profile before the Hydrodynamic Length
Here is the velocity profile before the Hydrodynamic Length. We can see that this plot is not properly parabolic.
The location where the velocity was measured for the above plot:
Whereas, if we compare the velocity profile at a location after the hydrodynamic length and in the fully developed region, the velocity profile will be much more defined and proper.
Velocity after the Hydrodynamic length
The location where the velocity was measured for the above plot:
Now, If we wish to calculate the Hydrodynamic length analytically, the formula is as follows:
`L_h = 0.005*R_e* D`
where:
Re = Reynolds number.
D = Diameter.
In our case, Re = 2100 and the Diameter = 0.1m, so the Hydrodynamic Length is:
Lh = 0.005*2100*0.1 = 1.05m.
And from the two graphs above, we can confirm that only after the Hydrodynamic length of approximately 1.0m, the velocity profile is changing into a fully parabolic shape.
b. Fully Developed Velocity Profile:
Velocity after the Hydrodynamic Length
The velocity profile in the fully developed region will be parabolic if we measure across in the radial direction of the Cylinder. We have achieved that and can confirm with the above plot. The velocity will be in this profile because of the friction on the wall and that this friction almost stops the flow of the fluid, and then because of the velocity difference between this layer of fluid and the adjacent layer, the velocity keeps on decreasing. But, the gradient is also decreasing as we move along the radial direction of the cylinder.
c. Velocity and Pressure Drop for the Fully developed region:
Velocity after the Hydrodynamic Length
The above graph is the representation of the Velocity Profile for the fully developed region in the cylinder. The velocity drop is approximately:
0.385-0.009 = 0.376 units.
Pressure drop along the axial direction
The pressure drop is shown in the plot above. The value for the Pressure drop is :
0.46-0.15 = 0.31 units.
Now, analytically, the pressure drop is:
`V_(avg) = triangleP * (D^2) / (32*mu*L)`
Here, `V_(avg) = 0.21`
`D = 0.1`
`mu = 0.00089`
`L = 1.5`
Hence = `triangle P = 0.36 units`
d. Shear Stress Plot:
Shear Stress Plot
Shear Stress is the maximum near the wall because of the friction that the wall has. But, as we move away from the wall, the Shear Force decreases and it can be confirmed with the graph above. - Conclusion: In conclusion, I would like to state that this challenge helped me understand how to simulate an axi-symmetric laminar flow through a cylinder. It also helped me comprehend the other methods, like the plot over line and even the way to measure the shear stress radially along the wall.