Flow through a pipe computational study in OpenFoam
Objective:
1. To simulate a 2D incompressible laminar flow through a cylinder pipe in OpenFOAM.
2. We can use the MATLAB function to automate the blockMeshDict file.
3. Validate Hydro-dynamic length with the numerical result.
4. Validate the fully developed flow velocity profile with its analytical profile.
5. Validate maximum velocity and pressure drop for fully developed flow.
6. Post-process Shear stress and validate wall shear stress for fully developed flow.
Assumptions:
- Diameter of pipe=0.02m
- wedge angle `theta`=3 deg
- length=2.7m
- fluid=water
- density=`1000(kg)/(m^3)`
- dynamic viscosity=`mu=0.00089Pa.s`
- kinematic viscosity=`nu=0.00000089(m^2)/(s)`
- axisymmetric,laminar,constant cross-section
- Reynolds number `Re=2100`
As per the calculations, the keep the flow laminar the average velocity of the flow is taken to be 0.0935 m/s and thus Reynold's Number, Re = 2100. For a circular pipe, the flow becomes turbulent after when Reynold's number becomes greater than 4000 and in between 2300-4000 the flow remains in the transition state.
At the entrance, the flow is said to be in the 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 fully developed flow which helps us to calculate the pressure loss. The distance required for a flow to develop is called the entrance length(Le), and depends on the shear stresses that act on the pipe wall.
For a circular pipe, we can calculate the average and maximum velocity, maximum shear stress and pressure loss analytically which are as follows.
`underset(avg)(v)=(muRe)/(rho underset(h)(D)`
`underset(max)(v)=2 underset(avg)(v)`
`underset(max)(tau)=(2muv)/(underset(h)(r))`
pressure drop=`DeltaP=(32mu underset(avg)(v)(L-Le))/(underset(h)(D)^2)`
here,
`underset(h)(D)`=Diameter
Entry Length=`underset(e)(L)=0.06*Re*underset(h)(D)`
Length of pipe=L=`0.05+underset(e)(L)`
Kinematic pressure=`(DeltaP)/(rho)`
sheer stress=`tau=-mu(du)/(dy)`
Governing equations
Pipe Geometry
For pipe flow applications, it’s often important to be able to determine whether a given flow condition is laminar flow or turbulent flow. Different equations or methods of analysis often apply to laminar flow and turbulent flow conditions.The typical criterion for whether pipe flow is laminar or turbulent is the value of the Reynolds Number.
The Reynolds number for pipe flow is defined as
`Re=DVρ/μ = uL/ν` where
D is the pipe diameter,
V is the average flow velocity in the pipe,
ρ is the density of the flowing fluid and
μ is the dynamic viscosity of the flowing fluid.
`nu` is the kinematic viscosity =`ρ/μ`=`1/nu`.
Re is a dimensionless number.
Any consistent set of units can be used for D, V, ρ, and, will result in Reynolds number being dimensionless. The generally accepted criteria for laminar flow and turbulent flow in terms of Re are as follows:
For Re < 2100, it is laminar flow.
For Re > 4000, it is turbulent flow.
For 2100 < Re < 4000, the flow is in the transition region.
For flow in the transition region, the flow may be either laminar or turbulent, depending upon the nature of the entrance to the pipe and the pipe wall roughness.
A pressure drop between the beginning and the end of a pipe is the drive for a flow. However, the pressure gradient that forms over the length of the pipe has two tasks, so to speak. The pressure drop must not only compensate for the frictional force during the flow of the fluid but must accelerate the fluid to achieve the characteristic flow profile in the first place.
Therefore, the total pressure drop can be divided into a part that is used to accelerate the fluid and a part that is used to compensate for friction. In the above equations, the pressure gradient dp/dx only refers to that part of the total pressure gradient which is necessary to overcome the friction.
To determine this pressure gradient, one must not simply take the pressure difference between the beginning and end of the pipe and then divide it by the length of the pipe. This is why the above equations always referred to a pipe section where it was implicitly assumed that the parabolic flow profile had already been fully established.
The equation of Hagen-Poiseuille does not apply in the inlet of a pipe (inlet flow), within which the pressure gradient must accelerate the fluid to the typical parabolic profile in addition to overcoming friction.
Under certain circumstances, it may nevertheless be justified to describe flow by simply determining the pressure gradient in the above equations from the quotient of pressure drop and pipe length. Namely, whenever the work of acceleration is negligible compared to the friction work. This will always be the case if the pipe is relatively long compared to its diameter, or if the flow velocities are relatively low. In such cases, the frictional work usually clearly outweighs the relatively low work of acceleration due to the long pipe length.
Solver: To solve the incompressible form we will solve it by icoFoam solver of OpenFOAM.
Boundary Types
patch
The basic patch type for a patch condition contains no geometric or topological information about the mesh (with the exception of the wall), e.g. an inlet or an outlet.
wall
There are instances where a patch that coincides with a wall needs to be identifiable as such, particularly where specialist modelling is applied at wall boundaries. A good example is wall turbulence modelling where a wall must be specified with a wall patch type so that the distance from the wall of the cell centres next to the wall are stored as part of the patch.
symmetry plane
For a symmetry plane.
empty
While OpenFOAM always generates geometries in 3 dimensions, it can be instructed to solve in 2 (or 1) dimensions by specifying a special empty condition on each patch whose plane is normal to the 3rd (and 2nd) dimension for which no solution is required.
wedge
For 2 dimensional axisymmetric cases, e.g. a cylinder, the geometry is specified as a wedge of small angle (e.g. 
) and 1 cell thick running along the plane of symmetry, straddling one of the coordinate planes. The axisymmetric wedge planes must be specified as separate patches of wedge type. The details of generating wedge-shaped geometries using blockmesh.
cyclic
Enables two patches to be treated as if they are physically connected; used for repeated geometries, e.g. heat exchanger tube bundles. One cyclic patch is linked to another through a neighbourPatch keyword in the boundary file. Each pair of connecting faces must have a similar area within a tolerance given by the matchTolerance keyword in the boundary file. Faces do not need to be of the same orientation.
processor
If a code is being run in parallel, on a number of processors, then the mesh must be divided up so that each processor computes on roughly the same number of cells. The boundaries between the different parts of the mesh are called processor boundaries.
Creating blocks with fewer than 8 vertices:
It is possible to collapse one or more pair(s) of vertices onto each other in order to create a block with fewer than 8 vertices. The most common example of collapsing vertices is when creating a 6-sided wedge-shaped block for 2-dimensional axisymmetric cases that use the wedge patch type. The process is best illustrated by using a simplified version of our example block shown in Figure below. Let us say we wished to create a wedge-shaped block by collapsing vertex 7 onto 4 and 6 onto 5. This is simply done by exchanging the vertex number 7 by 4 and 6 by 5 respectively so that the block numbering would become: hex (0 1 2 3 4 5 5 4).
Creating a wedge-shaped block with 6 vertices
The same applies to the patches with the main consideration that the block face containing the collapsed vertices, previously (4 5 6 7) now becomes (4 5 5 4). This is a block face of zero area which creates a patch with no faces in the polyMesh, as the user can see in a boundary file for such a case. The patch should be specified as empty in the blockMeshDict and the boundary condition for any fields should consequently be empty also.
The edges
Each edge joining 2 vertex points is assumed to be straight by default. However, any edge may be specified to be curved by entries in a list named edges. The list is optional; if the geometry contains no curved edges, it may be omitted.
edges
(
arc 1 5 (1.1 0.0 0.5)
);
Main code
% ********************************************** %
%************************************************%
clear all
close all
clc
%% Problem Initialization
%Material: Water @25 deg C inside
%Outside: Ambient conndition
Re = 2100; %Reynolds Number
rho = 1000; % Density@ 25 deg C : Unit: kg/m3
mu = 0.89e-3; %0.00089; %nu*rho; %Dynamic Viscosity@25 deg C: Unit: Pa.s
nu = mu/rho; %kinematic viscosity@25 deg C: Unit: m2/s
D = 0.02; % Dia of the pipe in meters
A_c = (pi/4)*D^2; % Cross sectional area
p = pi*D; % wetted perimeter
D_h = 4*A_c/p; %-> Hydraulic Dia of pipe
r_h = D_h/2; %Hydralic radius of pipe
v_avg = (Re*mu)/(rho*D_h);% Avg velocity
v_max = 2*v_avg; %Fully developed V max
tau_max = 2*mu*v_max/r_h;
if Re < 2300
L_e = 0.06*Re*D_h;
elseif Re == 2300 %Ref: Cengel Fluid Dynamics
L_e = 115*D_h;
elseif Re > 4000
L_e = 1.359*D_h*Re^(1/4);
end
L =L_e+0.18; %Meters --> length of Pipe
%% Pressure Loss %Daracy Friction factor --> for Laminar Circular pipe flow
f = 65/Re;
%Haigen Poiseuille Eqn
P = (32*mu*(L-L_e)*v_avg)/(D^2)
g = 9.81;
h_loss = P/(rho*g); %head loss
P_kinematicPress = P/rho;
theta = 3; %3 degree wedge angle
r=r_h;
%% Print out
haagen(r,L,theta)Function code for blockMesh file
function haagen(r,L,theta)
%% creating blockmesh dict file%% Opening Part
h1='/*--------------------------------*- C++ -*----------------------------------*';
h2=' ========= |';
h3=' / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4=' / O peration | Website: https://openfoam.org';
h5=' / A nd | Version: 9';
h6=' / M anipulation |';
h7='*---------------------------------------------------------------------------*/';
h8='FoamFile';
h9='{';
h10=' format ascii;';
h11=' class dictionary;';
h12=' object blockMeshDict;';
h13='}';
h14='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
f1=fopen('blockMeshDict.txt','wt');
fprintf(f1,'%sn',h1);
fprintf(f1,'%sn',h2);
fprintf(f1,'%sn',h3);
fprintf(f1,'%sn',h4);
fprintf(f1,'%sn',h5);
fprintf(f1,'%sn',h6);
fprintf(f1,'%sn',h7);
fprintf(f1,'%sn',h8);
fprintf(f1,'%sn',h9);
fprintf(f1,'%sn',h10);
fprintf(f1,'%sn',h11);
fprintf(f1,'%sn',h12);
fprintf(f1,'%sn',h13);
fprintf(f1,'%sn',h14);
%%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 0 0];
v4 = [L r*cosd(theta/2) r*sind(theta/2)];
v5 = [L r*cosd(theta/2) -r*sind(theta/2)];
mid1 = [0 r 0];
mid2 = [L r 0];
%%Write vertices
h15='convertToMeters 1;';
fprintf(f1,'%snn',h15);
h16='vertices';
fprintf(f1,'%sn',h16);
h17='(';fprintf(f1,'%sn',h17);
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);h18=');';
fprintf(f1,'%sn',h18);
%% write Blocks
h19='blocks';
h20='(';
h21=' hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 0.2 1)';
h22=');';
fprintf(f1,'n%sn',h19);
fprintf(f1,'%sn',h20);
fprintf(f1,'n%sn',h21);
fprintf(f1,'%sn',h22);
%% write edges
h23='edges';
h24='(';
h25=' arc 1 2';
h26=' arc 4 5';
h27=');';
fprintf(f1,'n%sn',h23);
fprintf(f1,'%sn',h24);
fprintf(f1,'%s (%d %d %d)n',h25,0, r, 0);
fprintf(f1,'%s (%d %d %d)n',h26,L, r, 0);
fprintf(f1,'%sn',h27);
%% write boundary
inlet = [0 1 2 0];
outlet=[3 5 4 3];
surface=[2 1 4 5];
axis=[0 3 3 0];
front=[0 3 4 1];
back=[0 2 5 3];
fprintf(f1,'n%sn','boundary');
fprintf(f1,'%sn','(');
%inlet patch
fprintf(f1,'t%sn','inlet');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type patch;');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',inlet);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%outlet patch
fprintf(f1,'nt%sn','outlet');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type patch;');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',outlet);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%surface
fprintf(f1,'nt%sn','surface');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wall;');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',surface);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%axis
fprintf(f1,'nt%sn','axis');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type empty;');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',axis);
fprintf(f1,'tt%sn',');');
fprintf(f1,'nt%sn','}');
%front
fprintf(f1,'nt%sn','front');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wedge; //symmetry');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',front);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%back
fprintf(f1,'nt%sn','back');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wedge; //symmetry');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',back);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
fprintf(f1,'%sn',');');
fprintf(f1,'n%sn','mergePatchPairs');
fprintf(f1,'%sn','(');
fprintf(f1,'%sn',');');
fprintf(f1,'n%sn','// ************************************************************************* //');
endBlockmeshDict file
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 9
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
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)
(2.700000e+00 0 0)
(2.700000e+00 9.996573e-03 2.617695e-04)
(2.700000e+00 9.996573e-03 -2.617695e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.700000e+00 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
surface
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
front
{
type wedge; //symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type wedge; //symmetry;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //U (initial)
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 9
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
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.0935 0 0);
}
outlet
{
type zeroGradient;
}
surface
{
type noSlip;
//type fixedValue;
//value uniform(0,0,0)
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //P(initial)
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 9
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
surface
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //controlDict File
It is important to provide sufficient time for the flow to develop and pass through the entire cross-section of pipe.
`"Time"="Distance"/"Velocity"`
`"Time"=2.7/0.0935=28.87s`
Hence we are taking the sufficient time of 50s for the flow to pass through pipe and `DeltaT=0.01s`
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 9
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //transportProperties
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 9
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-07;
// ************************************************************************* //1. In the first step we have defined the geometry parameters and then passed the radius, length of pipe and wedge angle to the haagen() function which basically printed the blockMeshDict file.
2. We have set the file structure of the openFOAM for the icoFoam solver.
3. Afterwards we have typed the following keys on the ubuntu terminal to simulate using icoFoam
4. Post-processing was done using paraView.
Undeveloped velocity at 0.05m.
| RESULTS |
Max velocity m/s |
Max Pressure loss Pa |
Max Sheer Stress Pa |
| Numerical | 0.1854 | 0.00121065 | 0.03279 |
| Theoretical | 0.1869 | 0.0012 | 0.0333 |