Week 11 - Simulation of Flow through a pipe in OpenFoam
*AIM:-
-To write a program in matlab that can generate the computational mesh automatically for any wedge angle and grading schemes
-To calculate the length of the pipe using the entry length formula for laminar flow through a pipe
-To compare velocity profiles at different positions along the length of the pipe
*Introdution:-
-Many different enterance lengths exist to describe a variety of flow conditions. The entrance length refers to the length of the entry region where effects from the pipe interior wall spread as an expanding boundary layer in the flow.
-The entrance length describes the distance travelled by the flow after entering the pipe before the flow becomes fully developed.
-Hydrodynamic entrance length describes the formation of a velocity profile caused by viscous forces propagating from the pipe wall.
-After the fluid becomes fully developed the flow characteristics no longer changes with the increased distance along the length of pipe.
-Entrance length for the laminar flow is given below:-
Entrance length=Le=0.06*Re*D
*Data:-
Reynold's number(Re):- 2100
Diameter of the pipe(D)= 0.02m
Radious of the pipe(R)= D/2= 0.01m
Density (rho)= 997 kg/m^3
Kinematic viscocity(v):- 0.00089 n/m^2
*Hagen-poisselus equation:-
- Poiseuille equation 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.
-Poiseuille's equation describes the pressure drop due to the viscocity of the fluid
-Average velocity:- `V_(avg)=(R_e∙μ)/(ρ∙D)`
-Maximum velocity:- `V_max=2∙V_(avg)`
-Pressure difference:- `∆p=(8∙μ∙L∙V_(avg))/r^2`
-Kinematic pressure difference:- `∆P_k=(∆p)/ρ`
-Shear stress:-`τ=(2∙μ∙V_max∙R)/r^2`
*Procedure for numerical setup and simulation:-
-First of all our simulation is incompressible-laminar-viscous flow, So we google and search "which solver use for incompresible laminar flow in openfoam"
-Here we seen that this condition is satisfied for the icoFoam solver, so we use this solver for our simulation.
-After that open terminal in the linux operating system.
-Then open the tutorial file, For that we use 'tut' shortcut key.
-Inside the tutorial directory type "ls" for the seen the directory and files inside that directory.
-Then type "cd incompressible/" for go to inside the incompressible directory.
-Again type "ls" and seen the directory and files inside that directory.Here we seen that icoFoam also available as directory.
-Type ''cd icoFoam/" and then use "ls" command, Here we seen that cavity files is available.
-We copy this file and run folder, Because every simulations are run in to only run folder.
-For that use this command:-
"cp -r cavity/ $FOAM_RUN/flowpipe_1
Here, cp= Copy any file or folder
-r= Entire directory
cavity= name of file
$FOAM_RUN= paste inside the openFoam run folder.
flowpipe_1= New name provide of the file
-After that go to the run directory using this command:- "cd $FOAM_RUN" , we seen that our das file present there.
-In this directory we changes setting depends on our requirements.
1) blockmeshDict:- It has the geometric and mesh charasteristics of the domain.
2) transportproperties:- It has information of transport quantity. In our case transport object is air and transport property is kinematic viscocity.
3) controlDict:- It has control of time.It use for start time and end time of the simulation.
4) fvschemes:- It has set the numerical schemes for the term such as derivatives in equations there are calculate the simulation.
5) fvsolutions:- It governs the equations, solver and numerical algorithm.
-After that run the simuation using th"icoFoam" command.
-Then use paraview software for seen the post-processing result. (Shortcut key:-"paraFoam")
1) blockmeshDict:-
/*--------------------------------*- 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 9.993e-03 3.489e-04) //1
(0 9.993e-03 -3.489e-04) //2
(2.82 0 0) //3
(2.82 9.993e-03 3.489e-04) //4
(2.82 9.993e-03 -3.489e-04) //5
);
blocks
(
hex (0 3 5 2 0 3 4 1) (500 25 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 1e-2 0)
arc 4 5 (2.82 1e-2 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
(
);
// ************************************************************************* //2) transportproperties:-
/*--------------------------------*- 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.9256e-07;
// ************************************************************************* //3) 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 0.8;
deltaT 1e-3;
writeControl timeStep;
writeInterval 10;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //4) fvschemes:-
/*--------------------------------*- 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 fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //5) fvsolutions:-
/*--------------------------------*- 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 fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //6) 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;
}
outlet
{
type fixedValue;
value uniform 5.377;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //7) 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.09373 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.09373 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //*Post-processing:-
-Before the simulation we know about the simulation is good or not, For that we use this two commands:-
i) checkMeah :- For mesh generation is good or not
ii) blockMesh:- For geometry is good or not
-Then open paraview software using the "paraFoam" Shortkey.
-We first of all seen the result of meshes using opction "surface with faces'
-Then create pressure and velocity profile.
1) Surface only:-
2) Surface with edge:-
3) Pressure:-
4) Velocity:-
5) Velocity profile at inlet:-
6) Velocity profile at 0.2m from inlet:-
7) Velocity profile at 0.5m from inlet:-
*Matlab code:-
clear all
close all
clc
%Reynold's number
Re=2100;
%Diameter of pipe
D=0.02;
%Entry length of pipe
L=0.06*Re*D;
%Total length of pipe
L_P=L+3;
%Radious of pipe
r=D/2;
%Angle of wedge
theta=4;
%Dynamic viscocity
mu=0.89e-03;
%Density
rho=997;
%Hagen poissuel flow equation
%Average velocity
v_avg=(Re*mu)/(rho*D);
%Maximum velocity
v_max=2*(v_avg);
%Pressure difference
del_p=(8*mu*L_P*v_avg)/(r^2);
%Kinematic pressure
kinematic_p=del_p/rho;
%Calculating shear stress
R=linspace(0,(D/2),1000);
for i=1:length(R)
%Shear stress
tau(i)=2*mu*v_max*(abs(R(i)))/r^2;
end
plot(R,tau)
xlabel('radious');
ylabel('shear stress');
axis([0 0.01 0 0.1]);
grid on
legend('shear stress');
header1='/*--------------------------------*- C++ -*----------------------------------*\'
header2=' ========= | '
header3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox '
header4=' \\ / O peration | Website: https://openfoam.org '
header5=' \\ / A nd | Version: 8 '
header6=' \\/ M anipulation | '
header7='\*---------------------------------------------------------------------------*/'
f1=fopen('blockMwshdict','w')
fprintf(f1, '%s\n', header1)
fprintf(f1, '%s\n', header2)
fprintf(f1, '%s\n', header3)
fprintf(f1, '%s\n', header4)
fprintf(f1, '%s\n', header5)
fprintf(f1, '%s\n', header6)
fprintf(f1, '%s\n', header7)
fprintf(f1, 'FoamFile\n{\n')
fprintf(f1, '%s \t 2.0;\n','version')
fprintf(f1, '%s \t ascii;\n','format')
fprintf(f1, '%s \t dictionary;\n','class')
fprintf(f1, '%s \t blockMeshDict;\n','object')
fprintf(f1, '}\n')
header8='\*---------------------------------------------------------------------------*/'
fprintf(f1, '%s\n\n', header8)
fprintf(f1, 'conertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1, '\t (%d %d %d)\n',0,0,0)
fprintf(f1, '\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',L_P,0,0)
fprintf(f1, '\t (%d %d %d)\n',L_P,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',L_P,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, ');\n\n')
fprintf(f1,'block\n')
fprintf(f1,'(\n')
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',500,25,1)
fprintf(f1, 'simpleGrading (1 0.2 1)\n\n')
fprintf(f1, ');\n\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,L_P,r,0)
fprintf(f1, ');\n\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1, '\t inlet \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type patch;\n')
fprintf(f1, '\t\t faces \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, '\t outlet \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type patch;\n')
fprintf(f1, '\t\t faces \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, '\t top \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wall;\n')
fprintf(f1, '\t\t faces \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, '\t front \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wall;\n')
fprintf(f1, '\t\t faces \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, '\t axis \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type empty;\n')
fprintf(f1, '\t\t faces \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,'mergePatchParis\n')
fprintf(f1, '(\n')
fprintf(f1, ');\n')
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%s\n', header9)
fclose(f1)*Plot of shear stress vs Radious:-
