Week 11 - Simulation of Flow through a pipe in OpenFoam
1. Aim:
To write a program in MATLAB that can generate the computational mesh automatically for any wedge angle and grading schemes and post-process Velocity profile and shear stress. Also, to compare the obtained velocity profile with velocity obtained from Hagen poiseuille's equation.
2. Equations solved:
a. Continuity Equation:
`(delu)/(delx)+(delv)/(dely)+(delw)/(delz)=0`
b. Momentum Equations:
`(delu)/(delt)+u⋅(∇u)=-(delp)/(delx)+μ*((del^2u)/(delx^2)+(del^2u)/(dely^2)+(del^2u)/(delz^2))`
`(delv)/(delt)+u⋅(∇v)=-(delp)/(dely)+μ*((del^2v)/(delx^2)+(del^2v)/(dely^2)+(del^2v)/(delz^2))`
`(delw)/(delt)+u⋅(∇w)=-(delp)/(delz)+μ*((del^2w)/(delx^2)+(del^2w)/(dely^2)+(del^2w)/(delz^2))`
c. Shear stress equation (solved in post-processing for shear stress plot):
`tau=mu(delu)/(delr)`
3. Procedure:
- First of all decide the geometry dimensions (radius of pipe, length of pipe, angle of wedge), flow properties (Temperature, Reynolds number), meshing properties (grading). The length of pipe should be chosen such that it is sufficiently greater than entrance length and so the flow is fully developed. Find out the density, dynamics viscosity and kinematic viscosity at the decided temperature. Also, find out the average flow velocity from Reynolds number equation.
- Write a MATLAB code using file parsing that takes angle of wedge, length of pipe, radius of pipe and grading in x,y,z directions as input and in return generates the blockMesh file for these geometry features. The block should be of the shape of a wedge. The curved surface is defined using 'arc' keyword inside edges part of the code.
- Define the initial pressure and velocity distribution at all boundaries and in internal region at inital time step (t=0) in the p and U file.
- In the transportProperties file, assign the appropriate value of kinematic viscosity to the variable nu.
- In controlDict file, specify time step size and endTime for simulation. Time step size should be such that the Courant number stays less than 1.
- Use blockMesh command in Terminal to generate mesh. Use icoFoam command to run the simulation and use paraFoam to view results in ParaView.
- Plot velocity along radial direction in fully developed region. Export this data and also plot the theoretical velocity using Hagen–Poiseuille equation and compare the two.
- Plot pressure along the length of pipe to find the pressure drop across the length of pipe. Compare it with pressure drop obtained from Hagen–Poiseuille equation.
- Finally plot shear stress variation in radial direction by calculating gradient of velocity and multiplying it by dynamic viscosity.
4. Calculations:
For the given problem, we assume the flowing fluid is water and simulation is performed at 25 degree Celsius and flow is laminar. So following are the properties of fluid (water):
Density: `rho = 997 kg//m^3`
Dynamic Viscosity: `mu=0.0008891Ns//m^2`
Kinematic Viscosity: `nu=mu/rho=0.00000089177 m^2//s`
Reynolds Number: `Re=2100`
Diameter of Pipe: `D=0.02m`
Entrance Length: `L_e=0.05*Re*D=2.1m`
Length of Pipe: `L=3m`
Velocity at Entrance: `V_(avg)=(Re*mu)/(rho*D)=0.936364m//s`
Max Theoretical Velocity: `V_(max)=2*V_(avg)=0.1872728m//s`
Velocity distribution in radial direction according to Hagen-Poiseuille equation: `u=V_(max)(1-r^2/R^2)`
Pressure drop according to Hagen-Poiseuille equation: `DeltaP=(32muLv_(avg))/D^2`
Kinematic pressure drop according to Hagen-Poiseuille equation: `DeltaP_k=(DeltaP)/rho`
5. MATLAB File for generating blockMeshfile:
clc
clear all
close all
theta=2; %Angle of Wedge cut out of cylinder (in degrees)
L=3; %Length of pipe (in m)
r=0.01; %Radius of Wedge (in m)
g_r=0.1; %Grading in radial direction
g_t=1; %Grading perpendicular to radial direction
g_a=1; %Grading in axial direction
f=fopen('blockMeshDict.txt','w'); %Open file in write mode
fprintf(f,'%s\n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f,'%s\n','| ========= | |');
fprintf(f,'%s\n','| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |');
fprintf(f,'%s\n','| \\ / O peration | Website: https://openfoam.org |');
fprintf(f,'%s\n','| \\ / A nd | Version: 8 |');
fprintf(f,'%s\n','| \\/ M anipulation | |');
fprintf(f,'%s\n\n','\*---------------------------------------------------------------------------*/');
fprintf(f,'%s\n','FoamFile');
fprintf(f,'%s\n','{');
fprintf(f,'\t%s\n','version 2.0;');
fprintf(f,'\t%s\n','format ascii;');
fprintf(f,'\t%s\n','class dictionary;');
fprintf(f,'\t%s\n','object blockMeshDict;');
fprintf(f,'%s\n','}');
fprintf(f,'%s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f,'%s\n','vertices');
fprintf(f,'%s\n','(');
fprintf(f,'\t(%d %d %d)\n',0,0,0);
fprintf(f,'\t(%d %d %d)\n',r*cosd(theta/2),-r*sind(theta/2),0);
fprintf(f,'\t(%d %d %d)\n',r*cosd(theta/2),r*sind(theta/2),0);
fprintf(f,'\t(%d %d %d)\n',0,0,L);
fprintf(f,'\t(%d %d %d)\n',r*cosd(theta/2),-r*sind(theta/2),L);
fprintf(f,'\t(%d %d %d)\n',r*cosd(theta/2),r*sind(theta/2),L);
fprintf(f,'%s\n\n',');');
fprintf(f,'%s\n','blocks');
fprintf(f,'%s\n','(');
fprintf(f,'\thex (0 1 2 0 3 4 5 3) (50 1 300) simpleGrading (%d %d %d)\n',g_r,g_t,g_a);
fprintf(f,'%s\n\n',');');
fprintf(f,'%s\n','edges');
fprintf(f,'%s\n','(');
fprintf(f,'\tarc 1 2 (%d %d %d)\n',r,0,0);
fprintf(f,'\tarc 4 5 (%d %d %d)\n',r,0,L);
fprintf(f,'%s\n\n',');');
fprintf(f,'%s\n','boundary');
fprintf(f,'%s\n','(');
fprintf(f,'\t%s\n','inlet');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type patch;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 2 1 0)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','curvedSurface');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type wall;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(1 2 5 4)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','front');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type wedge;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 1 4 3)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','back');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type wedge;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 3 5 2)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','outlet');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type patch;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(3 4 5 3)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','axis');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type empty;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 3 3 0)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'%s\n\n',');');
fprintf(f,'%s\n','mergePatchPairs');
fprintf(f,'%s\n','(');
fprintf(f,'%s\n\n',');');
fprintf(f,'%s\n','// ************************************************************************* //');6. Codes Files for Mesh, Time Control and Initial Conditions:
a. blockMesh 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
vertices
(
(0 0 0)
(0.00999848 -0.000174524 0)
(0.00999848 0.000174524 0)
(0 0 3)
(0.00999848 -0.000174524 3)
(0.00999848 0.000174524 3)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 300) simpleGrading (0.1 1 1)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
curvedSurface
{
type wall;
faces
(
(1 2 5 4)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(0 3 5 2)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //b. 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 startTime;
startTime 0;
stopAt endTime;
endTime 40.0;
deltaT 0.02;
writeControl timeStep;
writeInterval 25;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //c. U File (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
{
inlet
{
type fixedValue;
value uniform (0 0 0.0936364);
}
curvedSurface
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //d. p File (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
{
inlet
{
type zeroGradient;
}
curvedSurface
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
outlet
{
type fixedValue;
value uniform 0;
}
axis
{
type empty;
}
}
// ************************************************************************* //e. transportProperties 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.00000089177;
// ************************************************************************* //7. Output:
a. Mesh (x,y are radial direction and z is flow direction):
b. Velocity Contour near entrance of pipe:
c. Velocity contour near exit region of pipe:
d. Velocity distribution contour in radial direction at exit (flow is fully developed):
e. Velocity variation plot in radial direction in fully developed region (z=2.7m):
f. Comparison of velocity obatined in simulation with the theoretical velocity (Hagen–Poiseuille equation) at z=2.7m:
g. Pressure variation in flow direction (relative to pressure at exit):
h. Shear stress distribution in radial direction in fully developed region (z=2.7m)
8. Inference and Conclusion:
The MATLAB code has successfully generated the desired mesh. From the first figure we can see that the mesh is graded in radial direction and is closely spaced near the outer surface. The second and third figure shows velocity contour near the entrance and exit. Here we can see that the value of velocity near the entrance is close to 0.1 m/s (0.0936364 precisely). And further ahead in the flow direction, the velocity close to the outer wall starts becoming zero due to no-slip condition. At the exit, we see that velocity near the central axis is close to 0.19 m/s which is almost double of the entrance (average velocity). This is expected as `V_(max)=2*V_(avg)`. The fourth figure shows the velocity contour in radial direction at exit. In this contour we can see the effect of boundary layer phenomenon and how it affects the velocity profile in fully developed region. The velocity is maximum near the central axis and drops to zero near the boundary.
After this we have plotted velocity variation in radial direction in the fully developed flow region and in the next figure compared it with the theoretical velocity. We can see that the two almost overlap and hence the velocity profile obtained is correct. Next figure shows the pressure variation (relative to pressure at exit) along the flow direction. The drop in kinematic pressure from 0.24m to 3m (fully developed region) as observed from the data is 0.004048 `m^2//s^2`. From the Hagen–Poiseuille equation pressure drop between 2.4m and 3m is calculated to be 0.004008 `m^2//s^2`. Both these values are very close. Also, the shear stress is dependent on the gradient of velocity and is expected to be higher near the boundary of pipe which is evident from the shear stress plot. From all these observations, we can say that the flow is fully developed and simulation is accurate.