Week 11 - Simulation of Flow through a pipe in OpenFoam
OBJECTIVE.
- To simulate an axis-symmetric laminar flow through the constant cross-sectional of the pipe by applying the wedge boundary condition.
- To validate our results with the Hagen-poiseuille equation.
Hagen-Poiseuille's law:
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross-section. The pressure drop is given by:
`DeltaP=(32muLV)/D^2`
where,
`DeltaP` is the pressure drop.
`mu` is the dynamic viscosity.
L is the length of the pipe.
V is the velocity.
At the entrance, the flow is said to be in hydrodynamic entrance region where the velocity profile gradually takes the shape of a parabola from a straight line. As soon as the parabolic velocity profile becomes steady the flow is said to be in its fully developed state. Normally Hagen Poiseuille's Equation is applicable for a fully developed flow which helps us to calculate the pressure loss.
Assumptions made:
Geometry: Axisymmetric.
Reynolds number: 2100.
Working fluid: Water.
Diameter of the pipe: 0.02 m.
Dynamic viscosity(`mu`)=0.89x10^-3 Pas
Kinematic viscosity(`nu`)=0.8927x10^-6 m^2/s
theta=4
EQUATIONS USED.
Average velocity
`V_(avg)=(R_emu)/(rhoD)` m/s
Maximum velocity
`V_max=2*V_(avg)` m/s
Pressure drop
`DeltaP=(32muLV)/D^2` pa
kinematic pressure drop=`(DeltaP)/rho`
PROCEDURE.
- Will start with copying the cavity file from the tutorial to the FOAM_RUN folder.
- In the run folder only we are going to simulate and whatever we want to do the changes in the blockMeshDict file and pressure and velocity files and controlDict files we will be doing here.
- creating the block in the OpenFoam here we are creating the wedge after making the required changes in the blockMeshDict file to check the geometry run the blockMesh command.
- After creating the wedge assign the boundary conditions, and also sitting the simulation end time and time-step in such a way that the courant number should be less than 0.5 and also in the transportProperties enter the kinematic viscosity.
-
Run the simulation using the icoFoam solver after calculating the time needed for the fluid to travel from inlet to outlet
-
After the simulation, use paraFoam to view the flow inside the wedge and use plot over line command to plot the velocity and shear stress changes at a point after the hydrodynamical entry region so that we analyse the fully developed region of the fluid.
-
Use Matlab to parse the blockMeshDict so that the same file can be used for any wedge values.
MATLAB code to generate the blockMeshDict file.
clear all
close all
clc
% inputs
D= 0.02;
r= D/2;
theta= 4;
Re= 2100;
L_h = 0.05*Re*D;
L_pipe= L_h ;
scale= 1;
%vertices
v0= [0 0 0];
v1= [0 r*cosd(theta/2) r*sind(theta/2)];
v2= [0 r*cosd(theta/2) -r*sind(theta/2)];
v3= [L_pipe 0 0];
v4= [L_pipe r*cosd(theta/2) r*sind(theta/2)];
v5= [L_pipe r*cosd(theta/2) -r*sind(theta/2)];
%number of cells along x,y and z axis
nx= 20;
ny= 1;
nz= 200;
%grading factors
gx= 1;
gy= 1;
gz= 1;
%writing to a file
f1 = fopen('blockMeshDict_45','w');
%header region
fprintf(f1,'%s\n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'%s\n',' ========= |');
fprintf(f1,'%s\n',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'%s\n',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'%s\n',' \\ / A nd | Version: 8');
fprintf(f1,'%s\n',' \\/ M anipulation |');
fprintf(f1,'%s\n','\*---------------------------------------------------------------------------*/');
fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n',' version 2.0;');
fprintf(f1,'%s\n',' format ascii;');
fprintf(f1,'%s\n',' class dictionary;');
fprintf(f1,'%s\n',' object blockMeshDict;');
fprintf(f1,'%s\n','}');
fprintf(f1,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
%initial inputs
fprintf(f1,'%s','convertToMeters ');
fprintf(f1,'%g',scale);
fprintf(f1,'%s\n',';' );
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'(\n');
fprintf(f1,' (%d %d %d)\n',v0);
fprintf(f1,' (%g %g %g)\n',v1);
fprintf(f1,' (%g %g %g)\n',v2);
fprintf(f1,' (%g %g %g)\n',v3);
fprintf(f1,' (%g %g %g)\n',v4);
fprintf(f1,' (%g %g %g)\n',v5);
fprintf(f1,');\n\n');
%Defining the wedge block
fprintf(f1,'%s\n','blocks');
fprintf(f1,'(\n');
fprintf(f1,' %s (%g %g %g %g %g %g %g %g)','hex (0 2 1 0 3 5 4 3)');
fprintf(f1,'%g %g %g)',nx,ny,nz);
fprintf(f1,' %s (%g %g %g)\n','simpleGrading',gx,gy,gz);
fprintf(f1,');');
fprintf(f1,'\n\n');
%edges (if needed to add arc, etc.)
fprintf(f1,'%s\n','edges');
fprintf(f1,'(\n');
fprintf(f1,' %s (%g %g %g)\n','arc 1 2' ,0 ,r, 0);
fprintf(f1,' %s (%g %g %g)\n','arc 1 2' ,L_pipe ,r, 0);
fprintf(f1,');\n\n');
%Boundary Conditions
fprintf(f1,'%s\n','boundary');
fprintf(f1,'(\n');
%Inlet
fprintf(f1,' %s\n','inlet');
fprintf(f1,' {\n');
fprintf(f1,' %s\n','type patch;');
fprintf(f1,' %s\n','faces');
fprintf(f1,' (\n');
fprintf(f1,' %s\n','(0 1 2 0)');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Outlet
fprintf(f1,' %s\n','outlet');
fprintf(f1,' {\n');
fprintf(f1,' %s\n','type patch;');
fprintf(f1,' %s\n','faces');
fprintf(f1,' (\n');
fprintf(f1,' %s\n','(3 4 5 3)');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Front wall
fprintf(f1,' %s\n','front');
fprintf(f1,' {\n');
fprintf(f1,' %s\n','type wedge;');
fprintf(f1,' %s\n','faces');
fprintf(f1,' (\n');
fprintf(f1,' %s\n','(2 5 3 0)');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Back wall
fprintf(f1,' %s\n','Back');
fprintf(f1,' {\n');
fprintf(f1,' %s\n','type wedge;');
fprintf(f1,' %s\n','faces');
fprintf(f1,' (\n');
fprintf(f1,' %s\n','(0 3 4 1)');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Top wall
fprintf(f1,' %s\n','pipe');
fprintf(f1,' {\n');
fprintf(f1,' %s\n','type wall;');
fprintf(f1,' %s\n','faces');
fprintf(f1,' (\n');
fprintf(f1,' %s\n','(1 4 5 2)');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
%Bottom wall
fprintf(f1,' %s\n','axis');
fprintf(f1,' {\n');
fprintf(f1,' %s\n','type empty;');
fprintf(f1,' %s\n','faces');
fprintf(f1,' (\n');
fprintf(f1,' %s\n','(0 3 3 0)');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,');\n\n');
%Merging Patch pairs
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'(\n');
fprintf(f1,');\n\n');
fprintf(f1,'%s','// ************************************************************************* //');
fclose(f1)blockMeshDict 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.0099390827 0.000348994967)
(0 0.0099390827 -0.000348994967)
(2.52 0 0)
(2.52 0.0099390827 0.000348994967)
(2.52 0.0099390827 -0.000348994967)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (20 1 200) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.52 0.01 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(2 5 3 0)
);
}
Back
{
type wedge;
faces
(
(0 3 4 1)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //P 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.01687;
}
front
{
type wedge;
}
Back
{
type wedge;
}
pipe
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
U 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.09373 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
front
{
type wedge;
}
Back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
controlDict
/*--------------------------------*- 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 8;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
RESULTS.
Pressure flow in the wedge for fully developed flow.
Velocity contour diagram.
Shear stress in the wedge.
Velocity profile at the length of x=0.01
Velocity profile at the length of x=2.52 at the fully developed flow region.
Shear stress plot
Pressure drop plot.
Comparison between the analytical and numerical solution.
variables Analytical Numerical
Velocity 0.1874 0.1700
shear stress 0.0334 0.03213