Simulation of Flow through a pipe in OpenFoam
Axisymmetric Flow & Wedge boundary condition:
Axisymmetric Flow:
The flow which possesses an axis of symmetry along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles with respect to this axis are called as axisymmetric flow. The cylindrical coordinate system `(r,theta,z)` is found to be convenient for analyzing such problems.
Wedge boundary condition:
This boundary condition is specified by two planes that must be selected on separate sides of the domain (referred to as front and back) on surfaces running along the axis.
Given condition:
1.Flow is laminar
2.flow through constant crossectional area
3. Reynold's number at inlet(Re)=2100
Hagen- Poiseuille's flow:
There is no exact solution for Navier Stroke equation.But inorder to derive exact solution,assume a special cases(Parallel flow).
Special cases:
Parallel Flow:
A flow which contains only one velocity component(`v_z`) and other two(`v_r=0,v_theta=0`)a re negligible,such flow are called as Parallel flow.
Examples of parallel flow:
- Couette Flow
- Poiseuille Flow
Poiseuille flow:
Poiseuille flow is pressure-induced axisymmetric flow in an long plate of Length(L) .If Poiseuille flow is occur in circular pipe of length(L) and Radius(R) is called as Hagen- Poiseuille's flow.
Hagen-Poiseuille flow is analyzed using cylindrical polar coordinates (r, θ, z) with origin on the center-line of the pipe entrance and z-direction aligned with the center-line.
Derivation of Hagen-Poiseuille flow:
Assumptions
- Parallel flow(`v_z!=0,v_r=v_theta=0`)
- Incompressible Flow(density(`rho`)=constant)
- Viscous and steady flow
- laminar flow(Re<2300)
- Fully developed flow(ie the velocity profile does not change in the fluid flow direction )
Continuity Equation(Polar coordinates):
for steady state and compressible flow
`(delv_r)/(delr)+v_r/r+1/r(delv_theta)/(deltheta)+(delv_z)/(delz)=0`
here `v_r=v_theta=0`
`(delv_z)/(delv_z)=0`
hence `v_z!=v(z)`(it is not a function of "Z")
we assumed Axisymmetric flow .hence all parameter will be constant in all angles
`del/(deltheta)(anyparameter)=0`
hence `v_z!=v_z(theta)`(It is not a function of `theta`)
then,`v_z`is only function to radius
`v_z=v_z(r)`
N-S Equation(Polar coordinates):
for incompressible flow
where `nabla.vecv=0`
Apply steady state condition`(del(any parameter))/(delt)=0`
apply parallel flow `(v_theta=v_r=0)`
assume body force terms are negligible
Along "r" direction
`(delp)/(delr)=0`
Along "`theta`" direction
`(delp)/(deltheta)=0`
then pressure(p) is function of only z direction.
`p=p(z)`
Along "z" direction:
`-(dp)/dz+mu((d^2v_z)/(dr^2)+1/r(dv_z)/(dr))=0`
`((d^2v_z)/(dr^2)+1/r(dv_z)/(dr))=1/mu(dp)/dz`
`(d)/(dr)((dv_z)/(dr))+1/r(dv_z)/(dr)=1/mu(dp)/dz`
multiply "r" on both sides
`(d)/(dr)(r(dv_z)/(dr))=1/mu(dp)/dzr`
Integrate above equation wrt "dr"
`int(d)/(dr)(r(dv_z)/(dr))dr=int1/mu(dp)/dzrdr`
`(dv_z)/(dr)=1/(2mu)(dp)/dzr+c_1/r`
Integrate above equation wrt "dr"
`int(dv_z)/(dr)dr=int1/(2mu)(dp)/dzrdr+intc_1/rdr`
`v_z=1/(4mu)(dp)/dzr^2+c_1lnr+c_2`
Apply Boundary conditions
At r=0,`(dv_z)/(dr)=0`
`c_2=0`
At r=R,`v_z=0`
`c_1=(-1/(4mu)(dp)/dzR^2)/lnR`
substitute `c_1&c_2`values in above equation
`v_z=R^2/(4mu)(-(dp)/dz)(1-r^2/R^2)`
this is Hagen -Poiseuille's equation for velocity profile distribution
Hence At r=0, the velocity is maximum
`(v_z)_max=R^2/(4mu)(-(dp)/dz)`
-ve sign indicates as pressure gradient decreases,its velocity will increases.
Flow Rate(Q):
`Q=(v_z)_(avg)*`Area of pipe
`Q=intv_z*`cross sectional area
`Q=intv_z(r)*2pir*dr`
`Q=intR^2/(4mu)(-(dp)/dz)(1-r^2/R^2)*2pir*dr`
`Q=(piR^4)/(8mu)(-(dp)/dz)`
`(v_z)_(avg)=Q/(piR^2)`
`(v_z)_(avg)=((piR^4)/(8mu)(-(dp)/dz))/(piR^2)`
`(v_z)_(avg)=(R^2)/(8mu)(-(dp)/dz)`
`(v_z)_(avg)=(v_z)_(max)/2`
SHEAR STRESS at wall:
Newton's law of viscosity
`T_(wall)=-mu((dv_z)/(dr))_(r=R)`
-ve sign indicates,as velocity gradient increases,Shear stress decreases.
`v_z=R^2/(4mu)(-(dp)/dz)(1-r^2/R^2)`
`(dv_z)/(dr)=R^2/(4mu)(-(dp)/dz)(-(2r)/R^2)`
`((dv_z)/(dr))_(r=R)=R/(2mu)(-(dp)/dz)`
`T_(wall)=-mu*R/(2mu)(-(dp)/dz)`
`T_(wall)=R/(2)(-(dp)/dz)`
we know
`(v_z)_(avg)=(R^2)/(8mu)(-(dp)/dz)`
`(-(dp)/dz)=((v_z)_(avg)*8mu)/(R^2)`
substitute `(-(dp)/dz)`value in shear stress equation
`T_(wall)=(4*mu*(v_z)_(avg))/R`
Skin friction coefficient(`c_f`):
`c_f=T_(wall)/(1/2*rho*(v_z)_(avg)^2)`
`c_f=((4*mu*(v_z)_(avg))/R)/(1/2*rho*(v_z)_(avg)^2)`
`c_f=(8*mu)/(rho*(v_z)_(avg)*R)`
multiply and divide by 2 on RHS side
`c_f=(16)/((rho*(v_z)_(avg)*2R)/mu)`
`c_f=(16)/((rho*(v_z)_(avg)*D)/mu)`
`c_f=(16)/(Re)`
where
`Re=(rho*(v_z)_(avg)*D)/mu`
CONCEPT:
- write Matlab code to calculate theoritical values of Hagen-Poiseuille flow at fully developed region and to generate automatically an fully functional blockMeshDict file.
- Update constant properties and initial condition for flow
- Update ControlDict file
- Run blockMesh command to check mesh was generated as per requirement.(ensure boundaryfile is empty in blockMeshDict file)
- Run icoFoam command
- Run paraFoam command to navigate to paraview.
MATLAB CODE:
clear all
close all
clc
%constant properties
Re=2100; %Reynolds number
D=0.02; %Diameter of pipe
R=D/2; %Radius of pipe
L=3; %length of pipe
mu=1.002*10^(-3); %Dynamic viscosity
rho = 998; %Density
theta=3; %angle of wedge
nu=mu/rho; %Kinematic viscosity
%vertices
v0=[0 0 0];
v1=[R*cosd((theta/2)) -R*sind((theta/2)) 0];
v2=[R*cosd((theta/2)) R*sind((theta/2)) 0];
v3=[0 0 L];
v4=[R*cosd((theta/2)) -R*sind((theta/2)) L];
v5=[R*cosd((theta/2)) R*sind((theta/2)) L];
% Mesh details
nx=10;ny=1;nz=100;%number of grid points along respective axis
gx=1;gy=1;gz=1; % grading factor along respective axis
%analytical solution
L_e=0.05*Re*D
vz_avg=(Re*mu)/(rho*D) % average velocity
vz_max=2*vz_avg %maximum velocity
dp= (-8*mu*vz_avg/(R^2))*L;
r=linspace(0,R,nx);
v_z=vz_max*(1-(r.^2./R^2))
T_wall=4*mu*vz_avg/R
plot(r,v_z,'linewidth',2)
xlabel('radius')
ylabel('velocity')
title('velocity profile at fully developed region')
f1=fopen('blockMeshDict.txt','w');
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,'%s \n\n','convertToMeters 1 ;');
fprintf(f1,'%s \n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,' (%d %d %d) \n', [v0(1) v0(2) v0(3)]);
fprintf(f1,' (%d %d %d) \n', [v1(1) v1(2) v1(3)]);
fprintf(f1,' (%d %d %d) \n', [v2(1) v2(2) v2(3)]);
fprintf(f1,' (%d %d %d) \n', [v3(1) v3(2) v3(3)]);
fprintf(f1,' (%d %d %d) \n', [v4(1) v4(2) v4(3)]);
fprintf(f1,' (%d %d %d) \n', [v5(1) v5(2) v5(3)]);
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s',' hex (0 1 2 0 3 4 5 3)');
fprintf(f1,' (%d %d %d)',[nx ny nz] );
fprintf(f1,' %s (%d %d %d)','simpleGrading', [gx gy gz]);
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',' inlet');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type patch;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 1 2 0)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' outlet');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type patch;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (3 5 4 3)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' wedgeside1');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type wedge;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 3 4 1)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' wedgeside2');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type wedge;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 2 5 3)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' pipewall');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type wall;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (2 1 4 5)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' axis');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type empty;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 0 3 3)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');
fclose(f1);
THEORITICAL CALCULATION OF POISEUILLE'S FLOW:
Reynold's Number(Re)=2100;
Diameter of Pipe(D)=0.02m;
Radius of Pipe(R)=D/2=0.01m;
Dynamic viscosity(`mu`)=1.002*10^(-3) kg/m.s
Density of Fluid(`rho`)=998 kg/m^3
kinematic viscosity(`nu`)=`mu/rho`=1.0040*10^(-6)
Above code results,blockMeshDict file as below:
/*--------------------------------*- 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)
(9.996573e-03 -2.617695e-04 0)
(9.996573e-03 2.617695e-04 0)
(0 0 3)
(9.996573e-03 -2.617695e-04 3)
(9.996573e-03 2.617695e-04 3)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (10 1 100) simpleGrading (1 1 1));
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wedgeside1
{
type wedge;
faces
(
(0 3 4 1)
);
}
wedgeside2
{
type wedge;
faces
(
(0 2 5 3)
);
}
pipewall
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 0 3 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //INITIAL CONDITION:
Velocity Updation(Initial condition):
/*--------------------------------*- 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 0 0.11);
}
outlet
{
type zeroGradient;
}
wedgeside1
{
type wedge;
}
wedgeside2
{
type wedge;
}
pipewall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //Pressure Updation(Initial condition):
/*--------------------------------*- 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 1;
}
wedgeside1
{
type wedge;
}
wedgeside2
{
type wedge;
}
pipewall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //ControlDict Updation:
/*--------------------------------*- 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 30;
deltaT 0.005;
writeControl timeStep;
writeInterval 30;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Results:
After updation of blockMesh:
Mesh was created as per number of Nodes and Grading factor.
Velocity profile:
Velocity Profile at different length of pipe:
Observed velocity profile was unchanged after z=2.4.so Hydro-dynamic length (Le) =2.4m.After Hydro-dynamic length (Le=z=2.4),the velocity profile was similar.
As we move away from center line,intensity of velocity was decreasing.Finally it attains NoSlip condition at walls.
Theoritical Velocity profile(Fully developed)
Here theoritical velocity profile was calculated by formulas in matlab code.
Shear stress(at wall):
As we move away from center axis,Shear stress is increased.Hence near wall of pipe,shear stress is more.
Theoritical stress is calculated from formulae(T_wall =0.423)
Validation:
| parameters | Analytical solution | Numerical solution |
| Hydro-dynamic length(Le) | 2.1m | 2.4m |
| velocity profile(vz) | Parabolic shape | Parabolic shape |
| maximum velocity(`(v_z)_max`) | 0.21m/s | 0.218m/s |
| wall shear stress `T_(wall)` | 0.0423N/m^2 | 0.03N/m^2 |