Simulation of Flow through Pipe: Part-1
A. INTRODUCTION:
This project simulates the flow through a cylindrical pipe. The flow through a cylindrical pipe is based on pressure difference and the pressure difference is given by Hagen Poiseuille equation:
`trianglep =frac{8muLQ}{piR^4}`
where `triangleP`is the pressure difference between the lengths; L is length of pipe; `mu` is dynamic viscosity; Q is volumetric flow rate; A is area of cross section of pipe.
In this project, the pressure difference is obtained numerically and is compared with the analytical value and checked for equality. Also it is shown that the entry level is enough to produce a fully developed flow.
A1. ABOUT THE CHALLENGE:
The distance the fluid has to travel to become fully developed is given by the Entry Length formula.
Given Reynolds number is 2100. This implies flow is laminar.
Hence Entry length formula for Laminar is given by: `frac{L}{d} = 0.06.Re`
where L is effective length; d is dia of pipe and Re is reyonolds number.
Assuming the dia of pipe as 10mm; substituting the values, we get L = 1.26m;
Hence a comparison of the velocity profiles have to be made before and after the critical distance to prove the entry length.
After the entry length, the Hagen Poiseuille equation can be applied and comparitive study between the numerical and analytical values can be made.
A2. ASSUMPTIONS
- Flow along pipe is axi-symmetrical. Hence considering wedge shape of geometry.
- Incompressible, Laminar flow;
- Working fluid: WATER
B. THEORY
B1. ABOUT SOLVER
IcoFoam solver solves the incompressibe Navier Stokes equations using PISO Algorithm. PISO stands for Pressure Implicit with Splitting of Operators. The algorithm involves pressure-velocity calculation for Navier Stokes equation with one predictor step and two corrector step.
B2. GEOMETRY & MESHING
The geometry of the pipe is symmetrical along its axis. Since the simulation of the pipe is a very tedious and resource consuming process, the advantage of the axi-symmetricity. The axi-symmetric element is a wedge and this simulation is based on geometry of wedge.
The geometry is input in the 'BlockMeshDict' file which is present inside the 'System' folder of the tutorial of the respective solver. The code of blockmeshDict file is attached below. The complete is defined through a right handed coordinate system into the solver via the block command. Once the vertices are defined, they are then converted into a block.
The wedge angle considered in this problem is 4 degrees. OpenFOAM recommends choosing angle within 5 degrees for axi-symmetrical cases.
After the geometry is defined, it is necessary to define the mesh points for discrerization. The mesh used here is hexahedron mesh. The geometry is meshed along the length of the pipe and along the radius of the pipe. The grading of the meshes is done near the walls and along the length of the pipe.
After the mesh is created, the boundaries are defined according to their function;
The blockmeshDict code is attached below:
/*--------------------------------*- 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.001; //scales to metres
vertices
(
(0 0 0) //point0
(-0.1744 4.9969 0) //point1
(0.1744 4.9969 0) //point2
(0 0 1700) //point3
(-0.1744 4.9969 1700) //point4
(0.1744 4.9969 1700) //point5
);
blocks
(
hex (0 2 1 0 3 5 4 3)(25 1 800) simpleGrading (0.1 1 1)
);
edges
(
arc 2 1 (0 5 0)
arc 4 5 (0 5 1700)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
Inclined1
{
type wedge;
faces
(
(4 3 0 1) //LHS
);
}
Inclined2
{
type wedge;
faces
(
(3 5 2 0) //RHS
);
}
Top
{
type wall;
faces
(
(2 5 4 1) //Top
);
}
Axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
The blockMeshDict file is also generated by the following MATLAB/OCTAVE code as attached below; The input of the code is wedge angle and the grading factor along the direction of radius of wedge and along the length of the pipe.
MATLAB/ OCTAVE CODE:
clear all
close all
clc
ang = input('Enter the wedge angle: ');
r = 5; %[assuming the radius of wedge as 5mm]
mesh_rad = input('Enter the mesh grade along the radius: ');
x0 = 0;
y0 = 0;
z0 = 0;
x1 = -(r*sind(ang/2));
y1 = (r*cosd(ang/2));
z1 = 0;
x2 = (r*sind(ang/2));
y2 = (r*cosd(ang/2));
z2 = 0;
x3 = 0;
y3 = 0;
z3 = 1700;
x4 = -(r*sind(ang/2));
y4 = (r*cosd(ang/2));
z4 = 1700;
x5 = (r*sind(ang/2));
y5 = (r*cosd(ang/2));
z5 = 1700;
v0 = [x0 y0 z0];
v1 = [x1 y1 z1];
v2 = [x2 y2 z2];
v3 = [x3 y3 z3];
v4 = [x4 y4 z4];
v5 = [x5 y5 z5];
f1 = fopen('E:blockMeshDict.txt','w');
header = "FoamFile n{ nversion 2.0; n format ascii; nclass dictionary; nobject blockMeshDict; n}";
fprintf(f1,header);
conv_m = "nnconvertToMeters 0.001;";
fprintf(f1,conv_m);
vertex_start = "nnvertices n(";
fprintf(f1,vertex_start);
fprintf(f1,"n("); %%vertex 0
fprintf(f1,'%d%c%d%c%d',v0(1),' ',v0(2),' ',v0(3) );
fprintf(f1,")");
fprintf(f1,"n("); %%vertex 1
fprintf(f1,'%d%c%d%c%d',v1(1),' ',v1(2),' ',v1(3) );
fprintf(f1,")");
fprintf(f1,"n("); %%vertex 2
fprintf(f1,'%d%c%d%c%d',v2(1),' ',v2(2),' ',v2(3) );
fprintf(f1,")");
fprintf(f1,"n("); %%vertex 3
fprintf(f1,'%d%c%d%c%d',v3(1),' ',v3(2),' ',v3(3) );
fprintf(f1,")");
fprintf(f1,"n("); %%vertex 4
fprintf(f1,'%d%c%d%c%d',v4(1),' ',v4(2),' ',v4(3) );
fprintf(f1,")");
fprintf(f1,"n("); %%vertex 5
fprintf(f1,'%d%c%d%c%d',v5(1),' ',v5(2),' ',v5(3) );
fprintf(f1,")");
vertex_end = "nn);";
fprintf(f1,vertex_end);
block_start = "nnblocks n(";
fprintf(f1,block_start);
fprintf(f1,'n%s%d%s','hex ( 0 2 1 0 3 5 4 3)(25 1 800) simpleGrading(',mesh_rad,' 1 1)');
block_end = "nn);";
fprintf(f1,block_end);
edge_start = "nnedges n(";
fprintf(f1,edge_start);
fprintf(f1,'n%s%d%c%d%c%d%c','arc 2 1 (',0,' ',r,' ',0,')');
fprintf(f1,'n%s%d%c%d%c%d%c','arc 2 1 (',0,' ',r,' ',1700,')');
edge_end = "nn);";
fprintf(f1,edge_end);
boundary_start = "nnboundary n(";
fprintf(f1,boundary_start);
fprintf(f1,'nInlet ');
fprintf(f1,'n{ ');
fprintf(f1,'n%s','type patch;');fprintf(f1,'nfaces');
fprintf(f1,'n( ');
fprintf(f1,'n%s','( 0 2 1 0 )');
fprintf(f1,'n); ');
fprintf(f1,'n} ');
fprintf(f1,'nOutlet ');
fprintf(f1,'n{ ');
fprintf(f1,'n%s','type patch;');
fprintf(f1,'nfaces');
fprintf(f1,'n( ');
fprintf(f1,'n%s','( 3 4 5 3 )');
fprintf(f1,'n); ');
fprintf(f1,'n} ');
fprintf(f1,'nInclined1 ');
fprintf(f1,'n{ ');
fprintf(f1,'n%s','type wedge;');
fprintf(f1,'nfaces');
fprintf(f1,'n( ');
fprintf(f1,'n%s','( 4 3 0 1 )');
fprintf(f1,'n); ');
fprintf(f1,'n} ');
fprintf(f1,'nInclined2 ');
fprintf(f1,'n{ ');
fprintf(f1,'n%s','type wedge;');
fprintf(f1,'nfaces');
fprintf(f1,'n( ');
fprintf(f1,'n%s','( 3 5 2 0 )');
fprintf(f1,'n); ');
fprintf(f1,'n} ');
fprintf(f1,'nTop ');
fprintf(f1,'n{ ');
fprintf(f1,'n%s','type wall;');
fprintf(f1,'nfaces');
fprintf(f1,'n( ');
fprintf(f1,'n%s','( 2 5 4 1 )');
fprintf(f1,'n); ');
fprintf(f1,'n} ');
fprintf(f1,'nAxis ');
fprintf(f1,'n{ ');
fprintf(f1,'n%s','type empty;');
fprintf(f1,'nfaces');
fprintf(f1,'n( ');
fprintf(f1,'n%s','( 0 3 3 0 )');
fprintf(f1,'n); ');
fprintf(f1,'n} ');
boundary_end = "nn);";
fprintf(f1,boundary_end);
mpp_start = "nnmergePatchPairs n(";
fprintf(f1,mpp_start);
mpp_end = "nn);";
fprintf(f1,mpp_end);
disp('The blockMeshDict file has been generated...');
fclose(f1);
Meshed Geometry:
- The front face of the wedge is shown in the image followed by isometric view.
- The wedge angle is 4 degrees.
- Mesh grading is given only for the walls.
B4. INITIAL AND BOUNDARY FIELDS:
BOUNDARY FIELDS:
Inlet/Outlet condition:
Usually for velocity, the zeroGradient condition which extrapolates the value at the oulet and and a fixedValue of velocity is given at the inlet at the respective axis of flow.
Usually in subsonic conditions, it is assumed that, at the outlet, the total static pressure, P = 0; and pressure is extrapolated at the inlet using Neumann Boundary condition.
Top:
The walls are given noSlip condition meaning the velocity at the walls = 0; which is true for any pipe flow. The flow is maximum at the centre of pipe and 0 at the walls. The pressure is set to zeroGradient which means the normal gradient is 0.
Wedge:
The inclined wedge angled boundaries are given a special boundary type called 'Wedge'. The wedge boundary takes into consideration the axi-symmetrical condition and interpolates the parameters accordingly.
Axis:
The 'empty' boundary condition is applied to the surface of the block which has been collapsed to create the wedge geometry. No information is needed from this plane and hence empty condition is applied here. The solver does not calculate any parameter on this plane.
Initial Conditions:
The pressure and velocity initial conditions are attached below:
VELOCITY INITIAL CONDITION:
/*--------------------------------*- 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.187257);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 0.187257);
}
outlet
{
type zeroGradient;
}
Inclined1
{
type wedge;
}
Inclined2
{
type wedge;
}
Top
{
type noSlip;
}
Axis
{
type empty;
}
}
// ************************************************************************* //
PRESSURE INITIAL CONDITION:
/*--------------------------------*- 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;
}
outlet
{
type fixedValue;
value uniform 0;
}
Inclined1
{
type wedge;
}
Inclined2
{
type wedge;
}
Top
{
type zeroGradient;
}
Axis
{
type empty;
}
}
// ************************************************************************* //
B5. CALCULATION OF TIMESTEP & RUN TIME
The simulation is run until it reaches steady state. As a general rule, the fluid should pass through the domain 10 times to achieve steady state.
The start time is t_start = 0 sec;
The end time is `nXfrac{L}{V}`; where L is length of domain, V is avg velocity of fluid` and n = 10;
The average velocity of the fluid is calculated from the given reynolds number;
Given Re = 2100;
`V = frac{2100 X 8.917 X10^-7}{0.010}`
Hence V = 0.187251 m/s;
And t-end = `10 X frac{1.7}{0.18725}`
t_end~ 91 sec;
Next, the timestep, 'dt' value is to be calculated. For calculating 'dt', the CFL condition is used.
Cell Size must be considered in direction of flow of information; dx is given by L/n where L is length of domain and n is number of meshes required in the direction. Hence considering block , L = 1.7m and n = 800; and hence dx = 0.0021;
Now CFL condition is given by `Co = frac{dt*Umax}{dx}`
where dt is time step, dx is cell size and Co is Courant number whose max value is equal to 1 for Convection cases.
Hence `dt = frac{Co*dx}{Umax}`
`t = (1X0.0021)/(0.18725)`
Hence value of dt ,must be less than this: dt_max = 0.0113
Taking dt = 0.001 for more accuracy and incorporating these values into the 'Controldict' File which is present in the 'system' folder.
The ControlDict code is attached below:
CONTROLDICT:
/*--------------------------------*- 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 91;
deltaT 0.001;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
B6. TRANSPORT PROPERTIES
The transport properties relates the viscosity to the strain rate of the fluid used by the solver. In icoFoam, the Newtonian model is used. The newtonian model assumes kinematic viscosity is constant.
The working fluid used in this project is water. Hence its kinematic viscosity is added into the Transport Property file.
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.917e-7;
// ************************************************************************* //
C. EXECUTION
- First, the geometry is created in the blockmesh Dict file.
- The mesh is generated according to requirement.
- The boundary types are specified.
- The Initial and Boundary conditions are entered as per application.
- The controlDict file is edited for Simulation run time and timestep;
- At any point of time, the status of the geometry can be viewed in paraView software by entering 'ParaFoam' command in the terminal.
- After the geometry is input, the simulation can be started by calling 'icoFoam' command in the terminal.
D. RESULTS
D1. COMPARISON OF VELOCITY PROFILES
Plot1: Velocity Profile at X = 0.1m
Plot2: Velocity Profile at X = 0.5m
Plot3: Velocity Profile at X = 1.1m
Plot4: Velocity Profile at X = 1.3m
Plot5: Velocity Profile at X = 1.4m
Plot6: Velocity Profile at X = 1.7m
INFERENCE:
- It can be observed that the velocity profile is flat near the centre of the pipe in the initial phase of flow as shown in plot 1.
- The flatness of the velocity profile keeps decreasing and the parabolic nature manifests itself in the velocity profile.
- Theoretically, the entry length value for this problem is 1.26m according to the calculation done above.
- It can be observed that the velocity profile does not undergo any big changes after the x = 1.3m mark.
- At x = 1.3m, the velocity profile has attained the complete parabolic shape as well.
- Hence after the entry level distance, the flow is fully developed and is sufficient to produce a fully developed flow.
- Henceforth, Hagen Poiseuille equation is applicable to the pipe.
D2. ANALYTICAL VS NUMERICAL STUDY:
Considering a length of 0.3 meters from X = 1.4m to X = 1.7m, to validate the Hagen Poiseuielle equation.
Numerical Data obtained at X = 1.7m
Numerical Data obtained at X = 1.4m
From the above plot data, the difference in kinematic pressure in the considered length of pipe `triangleP` = 0.161085 Pa (Obtained from the difference of the magnitude of the red line from the previous two plots) and the max velocity, Vmax obtained is 0.371364 m/s numerically.
According to analytical calculations:
- We find that the obtained values from analytical solution matches closely with the numerical solutions post the entry region of the pipe.
- The kinematic pressure drop error percentage obtained is ~0.8% and velocity error is ~0.83%
- Hence, it can be said that the velocity plots matches with the Hagen Poiseuielle equation after the Entry level distance.
D3. SHEAR STRESS CURVE:
INFERENCE:
- The shear stress diagram is obtained by post-processing the velocity results obtained from Openfoam.
- The shear stess of fluid in pipe is given by: `tau = mu.frac{dU}{dy}`
- Shear stress is basically the slope of the velocity profile multiplied by the absolute viscosity. Hence it is calculated by using the formula: `frac{y2-y1}{x2-x1}`
- It is found that the shear stress is zero near the centre of the pipe and it is maximum near the walls of the gometry. This means that the fluid does not expericene any resistant force at the centre of the pipe and hence velocity of the pipe is maximum at the middle.
- The shear stress is a linear plot and the exact shape can be obtained when the best fit line is plotted on the above plot.
E. CONCLUSION:
- The project is solved using axi-symmetrical wedge shaped element since, simulation of the whole pipe would require a lot of computer resource.
- The velocity profiles changes until the flow is fully developed. In this project, the entry length for the flow is 1.24m. Post this distance from inlet, the velocity attains its characteristic parabolic shape and does not change any further.
- The velocity profile obtained is parabolic in nature with max velocity attained at the centre of the pipe.
- Once the flow is fully developed, the Hagen Poiseuielle equation is applicable. Comparison analysis is done between the numerical and analytical solution and it is found that the error percentage is around 1%. Hence, the numerical velocity profiles obtained after the entry length hold good for Hagen Poiseuielle equation.
- The Shear stress plot is obtained after post processing and it is found that the shear stress is maximum at the walls and zero at the centre of the pipe.