Simulation of flow through a pipe in OpenFoam using symmetry boundary condition
OBJECTIVE:
1) To develop a MatLab program that can generate the computational mesh automatically (in OpenFoam) for given angles using a symmetry boundary condition.
2) To simulate the flow through the pipe in OpenFOAM using icoFoam solver and post-process the result in paraFoam.
3) Compare the results obtained from the Hagen Poisuelle\'s equation, using wedge boundary condition, using a symmetry boundary condition.
Introduction: Take a cylindrical pipe and assume that the liquid exhibits laminar flow in the pipe. When this fluid enters into the pipe so there is a bunch of circular layers called lamina of liquid and each layer having a velocity which determined by the radial distance from the center of the tube.
We completed the simulation through the pipe in OpenFOAM using wedge boundary conditions for angle 4, so we will use the direct result here.
for Symmetry BC we use symmetry for front and back.
Flow properties:
1) Geometry: Axi-symmetric
2) Reynolds number Re: 2100 (Laminar Flow)
3) Working fluid: water (Density`rho`=997 kg /m^3
4) Fluid temperature: 25 degree Celsius
5) Diameter of pipe = 0.02 m (R=0.01 m)
5) Dynamic viscosity(`mu`) = 0.890*10^-3 pa.s
6) Kinematic viscosity (`nu`) = `frac{mu}{rho}`
= `0.8927*10^-6 m^2s^-1`
7) theta = 10,25,45
and calculate the length of the pipe using entry length formula
L = 0.06*Re*d
= 0.06*2100*0.02
=2.52 meters.
Hagen Poiseuille Flow:
HP flow applies to the Laminar flow of an incompressible fluid of viscosity through a tube of length and given radius.
According to this,
Flow rate: `Q= V_(avg)A_e=frac{DeltaPpiD^4}{128muL}`
and Average Velocity `V_(avg)=frac{DeltaPD^2}{32muL}`
Analytical Solution
1) Average Velocity: `V_(avg) = frac{R_emu}{rhoD}=0.09373 ms^-1`
2) Max Velocity: `V_(max)=2.V_(avg)=0.1875 ms^-1`
3) Pressure Drop: `DeltaP =frac{32muLV_(avg)}{D^2}=16.817 pa`
4) Kinematic Pressure drop: `frac{DeltaP}{rho}=0.01687 m`
Geometry and Mesh Generation
We use the icoFoam solver for incompressible flow and geometry will make with the use of MatLab code which further generates the mesh for the geometry.
MatLab Script:
This MatLab script generates a blockMesh txt file which can be further converted into a c++ file and use for generating a mesh for the geometry.
% Matlab scripts for generating the blockMeshDict file for a wedge angle
clear all
close all
clc
% Given Inputs
% wedge angle
theta = input(\'Enter the value of wedge angle = \');
Re = 2100;
D=0.02; %Diameter of the pipe
L = 0.06*Re*D; % Length of the pipe
R=D/2;
mu = 0.89e-3; %Viscosity of water at 25 deg C
rho = 997.0; % Densiy of water at 25 deg C
nu = mu/rho; % Kinematic viscosity of water at 25 deg C
% Calculation of pressure ans velocity
v_avg = Re*mu/(D*rho); %Average Velocity
v_max = 2*v_avg ; %Max Velocity
del_P = (32*mu*v_avg*L)/(D^2); % Pressure drop
kin_P = del_P/rho; % Kinematic Pressure drop
n1 = 50; % Number of cells along the pipe
n2 = 1; % Number of cells radially
n3 = 300; % Number of cells on the curve
gfactor=0.5;
% Vertices of the wedge
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)];
% Printing the file
text1=\'/*-----------------------------*- C++ -*----------------------------*\\\';
text2=\' ========= | \';
text3=\' \\\\ / F ield | OPENFOAM: The Open Source CFD Toolbox \';
text4=\' \\\\ / O peration | Version: 7 \';
text5=\' \\\\ / A nd | Web: www.OpenFOAM.ORG \';
text6=\' \\\\/ M anipulation | \';
text7=\'\\*------------------------------------------------------------------*/\';
text8=\' FoamFile\';
text9=\'{\';
text10=\' version 2.0;\';
text11=\' format ascii;\';
text12=\' class dictionary;\';
text13=\' object blockMeshDict;\';
text14=\'}\';
text15=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
text16=\'//*****************************************************************//\';
% Creating file
f1 = fopen(\'blockMeshDict.txt\',\'wt\');
fprintf(f1,\'%s\\n\',text1);
fprintf(f1,\'%s\\n\',text2);
fprintf(f1,\'%s\\n\',text3);
fprintf(f1,\'%s\\n\',text4);
fprintf(f1,\'%s\\n\',text5);
fprintf(f1,\'%s\\n\',text6);
fprintf(f1,\'%s\\n\',text7);
fprintf(f1,\'%s\\n\',text8);
fprintf(f1,\'%s\\n\',text9);
fprintf(f1,\'%s\\n\',text10);
fprintf(f1,\'%s\\n\',text11);
fprintf(f1,\'%s\\n\',text12);
fprintf(f1,\'%s\\n\',text13);
fprintf(f1,\'%s\\n\',text14);
fprintf(f1,\'%s\\n\',text15);
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',\'convertToMeters 1;\');
fprintf(f1,\'%s\\n\',\'vertices\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'%s (%f %f %f) \\n\',blanks(5),V0(1),V0(2),V0(3));
fprintf(f1,\'%s (%f %f %f) \\n\',blanks(5),V1(1),V1(2),V1(3));
fprintf(f1,\'%s (%f %f %f) \\n\',blanks(5),V2(1),V2(2),V2(3));
fprintf(f1,\'%s (%f %f %f) \\n\',blanks(5),V3(1),V3(2),V3(3));
fprintf(f1,\'%s (%f %f %f) \\n\',blanks(5),V4(1),V4(2),V4(3));
fprintf(f1,\'%s (%f %f %f) \\n\',blanks(5),V5(1),V5(2),V5(3));
fprintf(f1,\'%s\\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',\'blocks\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'%s %s\\n\',blanks(5),\'hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)\');
fprintf(f1,\'%s\\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',\'edges\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'%s arc 1 2 (%f %f %f) \\n\',blanks(5),0,R,0);
fprintf(f1,\'%s arc 4 5 (%f %f %f) \\n\',blanks(5),L,R,0);
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',\'boundary\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'inlet\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'{\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'type patch;\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'faces\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'(\');
fprintf(f1,\'%s %s \\n\',blanks(15),\'(0 1 2 0)\');
fprintf(f1,\'%s %s \\n\',blanks(10),\');\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'}\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'outlet\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'{\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'type patch;\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'faces\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'(\');
fprintf(f1,\'%s %s \\n\',blanks(15),\'(3 4 5 3)\');
fprintf(f1,\'%s %s \\n\',blanks(10),\');\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'}\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'front\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'{\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'type symmetry;\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'faces\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'(\');
fprintf(f1,\'%s %s \\n\',blanks(15),\'(0 3 5 2)\');
fprintf(f1,\'%s %s \\n\',blanks(10),\');\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'}\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'back\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'{\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'type symmetry;\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'faces\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'(\');
fprintf(f1,\'%s %s \\n\',blanks(15),\'(0 1 4 3)\');
fprintf(f1,\'%s %s \\n\',blanks(10),\');\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'}\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'pipe\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'{\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'type wall;\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'faces\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'(\');
fprintf(f1,\'%s %s \\n\',blanks(15),\'(1 4 5 2)\');
fprintf(f1,\'%s %s \\n\',blanks(10),\');\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'}\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'axis\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'{\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'type empty;\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'faces\');
fprintf(f1,\'%s %s \\n\',blanks(10),\'(\');
fprintf(f1,\'%s %s \\n\',blanks(15),\'(0 3 3 0)\');
fprintf(f1,\'%s %s \\n\',blanks(10),\');\');
fprintf(f1,\'%s %s \\n\',blanks(5),\'}\');
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'mergePatchPairs\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',text16);
fclose(f1);
Generate BlockMeshDict file for angle 10 degree
/*-----------------------------*- C++ -*----------------------------*\\
========= |
\\\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Version: 7
\\\\ / A nd | Web: www.OpenFOAM.ORG
\\\\/ M anipulation |
\\*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009962 -0.000872)
(0.000000 0.009962 0.000872)
(2.520000 0.000000 0.000000)
(2.520000 0.009962 -0.000872)
(2.520000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
//*****************************************************************//
Generated BlockMeshDict file for angle 25 degree
/*-----------------------------*- C++ -*----------------------------*\\
========= |
\\\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Version: 7
\\\\ / A nd | Web: www.OpenFOAM.ORG
\\\\/ M anipulation |
\\*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009763 -0.002164)
(0.000000 0.009763 0.002164)
(2.520000 0.000000 0.000000)
(2.520000 0.009763 -0.002164)
(2.520000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
//*****************************************************************//
Generated BlockMeshDict file for Angle 45 degree
/*-----------------------------*- C++ -*----------------------------*\\
========= |
\\\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Version: 7
\\\\ / A nd | Web: www.OpenFOAM.ORG
\\\\/ M anipulation |
\\*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009239 -0.003827)
(0.000000 0.009239 0.003827)
(2.520000 0.000000 0.000000)
(2.520000 0.009239 -0.003827)
(2.520000 0.009239 0.003827)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
//*****************************************************************//
Generated Mesh for angle 10 degree
Generated Mesh for angle 25 degree
Generated Mesh for 45 Degree
Initial conditions and Boundary Conditions
Pressure
/*--------------------------------*- 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
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //
Velocity
/*--------------------------------*- 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);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
pipe
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //Set up conditions for time step dt, start time and end time
/*--------------------------------*- 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 10;
deltaT 0.005;
writeControl timeStep;
writeInterval 25;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Also, use the
Kinematic viscosity (`nu`) = `frac{mu}{rho}`
= `0.8927*10^-6 m^2s^-1` condition
Simulation Result:
For angle 10 degree
For angle 25 degree
For angle 45 degree
Velocity variation at the inlet
Velocity variation at the outlet (Fully developed region)
Velocity contour at different positions for wedge angle 4 deg
1) at 0.1m from the Inlet
2) at 0.5m from the Inlet
3) At 1 m from the Inlet
4) At 2 m from the inlet
4) At 2.52 m from the Inlet
Shear stress profile:
We can define shear stress analytically from the given formula
`tau=mufrac{du}{dy}`
where,
`tau` = shear stress(F/A)
`mu` = dynamic viscosity
du/dy = rate of shear deformation
Numerically we can find the shear stress by multiplying velocity gradient with the dynamic viscosity of fluid under the study. Here we use the water as the fluid and his dynamic viscosity is 0.000890 ps-s and plot the graph between shear stress and radial distance.
From the graph, we can say that at the end of the pipe as we go towards the center of the pipe then shear force decreases while when we go towards the wall of pipe it gradually increases due to the fully developed velocity region.
Result and Comparision:
|
Wedge angle (Degree) |
Max Velocity obtained from the simulation |
Max Velocity obtained from Hagen Poiseuille’s Equation |
Courant number Mean |
Courant Number Max |
Max shear stress |
|
4 (Wedge BC) |
0.175932 |
0.1875 |
0.0384577 |
0.143645 |
0.0339956 |
|
10 (symmetry BC) |
0.176095 |
0.1875 |
0.0384587 |
0.143895 |
0.034081 |
|
25 (symmetry BC) |
0.177121 |
0.1875 |
0.0384646 |
0.145835 |
0.0346265 |
|
45 (symmetry BC) |
0.179638 |
0.1875 |
0.0384791 |
0.150125 |
0.0360188 |
We can see in the table that at an angle 45 degrees (symmetry BC) we get the highest value for maximum velocity so we can that as we increase the wedge angle the maximum velocity closes to the theoretical value. We also get the maximum shear stress value for the angle 45 degree which uses the symmetrical BC.
Courant number also increases as we increase the angle of wedge.
Conclusions:
1) Matlab program is successfully developed and from that, we generate the BlockMesh in OpenFoam for different angles (10,25,45) using symmetry BC that is one kind of data processing fro Matlab to OpenFoam.
2) Hagen Poiseullie\'s method solves this flow problem analytically and In the simulation, we use the Wedge BC method and symmetry BC method.
3) Symmetry BC is slightly better than the wedge BC as we get the result closer to the analytical results in symmetry BC.
4) Symmetry BC took some more time to complete the simulation as compare to wedge BC but give the closest result to analytical results.
5) As we increase the wedge angle, the Courant number also increases.