Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim:
To simulate an axi-symmetric flow through a pipe by applying the wedge boundary condition
Objective:
- To validate Hydro-dynamic length with the numerical result
- To validate the fully developed flow velocity profile with its analytical profile
- To validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
Assumptions:
- The flow is steady and laminar
- Pipe cross section is constant
- Flow is incompressible
Theory:
Hydro-dynamic length
In fluid dynamics, the entrance length is the distance a flow travels after entering a pipe before the flow becomes fully developed
For laminar flow entry length is given by:
`L_e=0.04*Re*D`
Hagen- Poiseuille's equations
`V_(avg)=(Re*mu)/(rho*D)`
`V_max=2*V_(avg)`
`triangleP=(8*mu*L*V_(avg))/r^2`
`triangleP_k=(triangleP)/rho`
`tau=(2*mu*V_max*R)/r^2`
where `Re=2100` is reynold's number
`mu=0.00089` is dynamic viscosity
D=0.02 is diameter of the pipe
R=0.01 is radius of the pipe
`rho=1000` is density of the fluid
`L_e=1.68` is length of the pipe
Simulation:
The simulation is carried in openFoam. For the given flow the icoFoam solver is chosen from the openFoam solvers and is then copied to the Foam_run folder for simulation.
For creating the geometry using File parsing in matlab the blockMeshDict file is generated.The geometry is as shown below
Code:
clear all
close all
clc
%Fluid properties
Re=2100;
rho=1000;
mu=0.00089;
nu=mu/rho;
%Pipe dimensions
D=0.02;
r=D/2;
L_e=0.06*Re*D;
L=L_e+0.7;
theta=3;
%Hagen- Poiseuille's equations
v_avg=(Re*mu)/(rho*D);
v_max=2*v_avg;
del_p=(8*mu*L*v_avg)/r^2;
p_k=del_p/rho;
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2=' ========= |';
h3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4=' \\ / O peration | Website: https://openfoam.org';
h5=' \\ / A nd | Version: 8';
h6=' \\/ M anipulation |';
h7='\*---------------------------------------------------------------------------*/';
h8=' version 2.0;';
h9=' format ascii;';
h10=' class dictionary;';
h11=' object blockMeshDict;';
h12='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%coordinates of mesh
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 0 0];
v4=[L r*cosd(theta/2) -r*sind(theta/2)];
v5=[L r*cosd(theta/2) r*sind(theta/2)];
f1=fopen('blockMeshDict.txt','w')
fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'FoamFile\n{\n');
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'}\n');
fprintf(f1,'%s\n\n',h12);
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n(\n');
fprintf(f1,'\t(%d %d %d)\n',v0);
fprintf(f1,'\t(%d %d %d)\n',v1);
fprintf(f1,'\t(%d %d %d)\n',v2);
fprintf(f1,'\t(%d %d %d)\n',v3);
fprintf(f1,'\t(%d %d %d)\n',v4);
fprintf(f1,'\t(%d %d %d)\n',v5);
fprintf(f1,');\n\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (%d %d %d) ',300,30,1);
fprintf(f1,'simpleGrading (1 0.1 1)\n');
fprintf(f1,');\n\n');
fprintf(f1,'edges\n');
fprintf(f1,'(\n');
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,L,r,0);
fprintf(f1,');\n\n');
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 1 2 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\ttop\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(1 4 5 2)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 4 1)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 2 5 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
h13='// ************************************************************************* //';
fprintf(f1,'%s\n',h13);
fclose(f1)blockMeshDict file:
The following file is generated and used for the simulation
/*--------------------------------*- 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 9.996573e-03 -2.617695e-04)
(0 9.996573e-03 2.617695e-04)
(3.220000e+00 0 0)
(3.220000e+00 9.996573e-03 -2.617695e-04)
(3.220000e+00 9.996573e-03 2.617695e-04)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (300 30 1) simpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.220000e+00 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
controlDict file:
The controlDict file is used to give simulation time and the time step values
/*--------------------------------*- 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 10;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Initial conditions:
Pressure(p)
For the pressure its given as kinematic pressure value calculated from the above equations
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0214247208;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Velocity(U)
/*--------------------------------*- 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.09345 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.09345 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //transportProperties file:
In this file the value for kinematic viscosity is given which is calculated using
`nu=mu/rho`
where `mu` is dynamic viscosity and `rho` is density of the fluid
/*--------------------------------*- 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] 8.9e-07;
// ************************************************************************* //Post-processing:
Mesh
Velocity contour
At inlet
At outlet
Results:
At x=0.01
At x=0.2
At x=0.5
At x=1
At x=1.5
At x=2
At x=2.5
As we can see that after x=1.5 there is not much variation in the velocity profile which indicates us that the flow is fully developed.
Pressure
Pressure difference across the length of the pipe. since there is pressure difference the flow is possible
Shear stress
Shear stress is maximum at the pipe wall and minimum at the center of the pipe
Conclusion:
Velocity profile of fully developed flow
- Analytical
- Simulation
Note: This profile is for the half cross section of the pipe
Validation
| Variables | Analytical value | Simulation value |
| Hydro-dynamic entry length | 1.68 | 1.5 |
| Maximum velocity | 0.186 | 0.18 |
| Pressure difference | 0.021 | 0.026 |
| Shear stress | 0.033 | 0.034 |