Simulation of Flow through a pipe in OpenFoam
Simulation of Flow through a pipe in OpenFoam
This project objective is to be able to accurately simulate the laminar incompressible flow through a pipe. This project will be broken into two different parts. This part of the project focused on the simulation set-up. This problem can be computationally expensive, and that is why the symmetry condition will be considered. Instead of simulating all the pipe, it will also have to be determined what is the minimum wedge angle that should be analyzed to realistically simulate the flow.
Objective:
Based on the desired Reynolds number and the pipe diameter, the inlet velocity and maximum velocity will be calculated. The working fluid use is water.
Then the length of the pipe using the entry length formula for laminar flow through a pipe will be calculated.
Then the entry length is sufficient to produce a fully developed flow and the velocity profiles at different positions along the length of the pipe will be shown.
Solver use is icoFoam.
A program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
The wedge angle will be less than 5 degrees.
The number of cells along the axis will be symmetric.
For a baseline mesh, the velocity profile matches with the Hagen poiseuille's equation.
Show that the velocity profile is fully developed.
Post process velocity and shear stress
The above figure shows that, the flow in a pipe takes some distance to adjust the sudden flow change from the pipe entrance. This is called the entrance or hydrodynamically entrance length. It is a function of Reynolds number.
1) Preliminary Data
The first step of the project is to establish the properties of flow and pipe simulated. The properties assumed for this projects are:
1. Reynolds number (Re) = 2100
2. Diameter of the pipe(D) = 0.01m
3. Dynamic viscosity(mu) = 8.9e-4 Pa.S, for water at 25C.
4. Density = 997 kg/m^3
5. Kinematic viscosity(nu) = 8.92e-7 m^2/s
2) Analytical Solution – Hagen Pouiseuille
Assuming the fully developed flow, the Navier Stoke Integral equations lead to the following results (Hagen Poiseuille equations based on these assumptions.)
1. Velocity = (Re*mu)/(rho*D) = 0.187 m/s
2. Maximum velocity = 2*Velocity = 0.375 m/s
3. Pressure difference = (32*mu*L*Velocity)/D^2
4. Entry Length(for laminar flow) = 0.05*Re*D = 1.05m
5. Total length of the pipe = 1.5m
3) Geometry and Meshing
4) Boundary and Initial conditions:
For internal flow through a pipe, we need to fix the inlet velocity condition and the outlet pressure. The inlet velocity has been calculated from the Reynolds number and the fluid properties. The outlet pressure must be ambient pressure (0) and the pressure drop in the fully developed region will be calculated numerically and will be compared with the analytical pressure drop.
5) Codes:
Coding by Matlab to generate BlockMeshDict file:
clear all
close all
clc
%Inputs
Re = 2100; %Reynolds number
d = 0.01; %diameter of pipe
r = d/2; %radius of pipe
Le = 1.05; %Entrance length
Lt = 1.5; %Total length
rho = 997; %Density of water
v_avg = 0.187; %average velocity
v_max = 0.375; %max velocity
p_k = 0.056; %kinematic pressure
mu = 0.00089; % dynamic viscocity
theta = 4; %wedge angle
h1 = '/*--------------------------------*- C++ -*----------------------------------*';
h2 = '| ========= | |';
h3 = '| \ / F ield | foam-extend: Open Source CFD |';
h4 = '| \ / O peration | Version: 3.1 |';
h5 = '| \ / A nd | Web: http://www.extend-project.de |';
h6 = '| \/ M anipulation | |';
h7 = '*---------------------------------------------------------------------------*/';
h8 = 'FoamFile'
h9 = '{';
h10 = ' version 2.0';
h11 = ' format ascii';
h12 = ' class dictionary';
h13 = ' object blockMeshDict';
h14 = '}';
h15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h16 = 'mergePatchPairs';
h17 = '(';
h18 = ')';
h19 = '// ************************************************************************* //';
conv = 'convertToMeters 1;';
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 = [Lt 0 0];
v4 = [Lt r*cosd(theta/2) r*sind(theta/2)];
v5 = [Lt r*cosd(theta/2) -r*sind(theta/2)];
%creating blank spaces for vertices
s1 = blanks(4);
s2 = blanks(11);
s3 = blanks(7);
s4 = blanks(3);
%Creating block
block = 'hex (0 2 1 0 3 5 4 3) (10 1 300) simpleGrading(1 1 1)';
edge1 = 'arc 1 2 (0 0.005 0)';
edge2 = 'arc 4 5 (1.5 0.005 0)';
%Writing blockMeshDict file
f1 = fopen('blockMeshDict.txt','w');
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);
fprintf(f1,'%sn',h15);
fprintf(f1,'%sn',conv);
fprintf(f1,'n');
fprintf(f1,'%sn','vertices');
fprintf(f1,'%sn','(');
fprintf(f1,'%s(%d %d %d)n',s1,v0(1),v0(2),v0(3));
fprintf(f1,'%s(%d %d %d)n',s1,v1(1),v0(2),v0(3));
fprintf(f1,'%s(%d %d %d)n',s1,v2(1),v0(2),v0(3));
fprintf(f1,'%s(%d %d %d)n',s1,v3(1),v0(2),v0(3));
fprintf(f1,'%s(%d %d %d)n',s1,v4(1),v0(2),v0(3));
fprintf(f1,'%s(%d %d %d)n',s1,v5(1),v0(2),v0(3));
fprintf(f1,'%sn',')');
fprintf(f1,'n');
fprintf(f1,'%sn','blocks');
fprintf(f1,'%sn','(');
fprintf(f1,'%sn',s1,block);
fprintf(f1,'%sn',');');
fprintf(f1,'n');
fprintf(f1,'%sn','egdes');
fprintf(f1,'%sn','(');
fprintf(f1,'%s',s1);
fprintf(f1,'%sn',edge1);
fprintf(f1,'%s',s1);
fprintf(f1,'%sn',edge2);
fprintf(f1,'%sn',');');
fprintf(f1,'n');
fprintf(f1,'%sn','boundary');
fprintf(f1,'%sn','(');
fprintf(f1,'%sn',h14);
fprintf(f1,'%s',s1);
fprintf(f1,'%sn','inlet');
fprintf(f1,'%s',s1);
fprintf(f1,'%sn','{')
fprintf(f1,'%s %sn',s3,'type patch;');
fprintf(f1,'%s %sn',s3,'faces');
fprintf(f1,'%s %sn',s3,'(');
fprintf(f1,'%s %sn',s4,'(0 2 1 0)');
fprintf(f1,'%s %sn',s3,');');
fprintf(f1,'%s %sn',s1,'}');
fprintf(f1,'n');
fprintf(f1,'%s %sn',s1,'outlet');
fprintf(f1,'%s %sn',s1,'{');
fprintf(f1,'%s %sn',s3,'type patch;');
fprintf(f1,'n');
fprintf(f1,'%s %sn',s3,'faces');
fprintf(f1,'%s %sn',s3,'(');
fprintf(f1,'%s %sn',s4,'(3 4 5 3)');
fprintf(f1,'%s %sn',s3,');');
fprintf(f1,'%s %sn',s1,'}');
fprintf(f1,'n');
fprintf(f1,'%s %sn',s1,'topWall');
fprintf(f1,'%s %sn',s3,'{');
fprintf(f1,'%s %sn',s3,'type wall;');
fprintf(f1,'%s %sn',s3,'faces');
fprintf(f1,'%s %sn',s3,'(');
fprintf(f1,'%s %sn',s4,'(1 2 5 4)');
fprintf(f1,'%s %sn',s3,');');
fprintf(f1,'%s %sn',s1,'}');
fprintf(f1,'n');
fprintf(f1,'%s %sn',s1,'bottomWall');
fprintf(f1,'%s %sn',s3,'{');
fprintf(f1,'%s %sn',s3,'type empty;');
fprintf(f1,'%s %sn',s3,'faces');
fprintf(f1,'%s %sn',s3,'(');
fprintf(f1,'%s %sn',s4,'(0 3 3 0)');
fprintf(f1,'%s %sn',s3,');');
fprintf(f1,'%s %sn',s1,'}');
fprintf(f1,'n');
fprintf(f1,'%s %sn',s1,'frontWedge');
fprintf(f1,'%s %sn',s3,'{');
fprintf(f1,'%s %sn',s3,'type wedge;');
fprintf(f1,'%s %sn',s3,'faces');
fprintf(f1,'%s %sn',s3,'(');
fprintf(f1,'%s %sn',s4,'(0 1 4 3)');
fprintf(f1,'%s %sn',s3,');');
fprintf(f1,'%s %sn',s1,'}');
fprintf(f1,'n');
fprintf(f1,'%s %sn',s1,'backWegde');
fprintf(f1,'%s %sn',s3,'{');
fprintf(f1,'%s %sn',s3,'type wedge;');
fprintf(f1,'%s %sn',s3,'faces');
fprintf(f1,'%s %sn',s3,'(');
fprintf(f1,'%s %sn',s4,'(3 5 2 0)');
fprintf(f1,'%s %sn',s3,');');
fprintf(f1,'%s %sn',s1,'}');
fprintf(f1,'%sn',');');
fprintf(f1,'n');
fprintf(f1,'%sn',h16);
fprintf(f1,'%sn',h17);
fprintf(f1,'%sn',h18);
fprintf(f1,'%sn',h19);
fclose(f1)
blockMeshDict file:
/*--------------------------------*- 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.00499 0.000174)
(0 0.00499 -0.000174)
(1.5 0 0)
(1.5 0.00499 0.000174)
(1.5 0.00499 -0.000174)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (10 1 300) simpleGrading (0.1 1 1)
);
edges
(
arc 1 2 (0 0.005 0)
arc 4 5 (1.5 0.005 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topWall
{
type wall;
faces
(
(1 2 5 4)
);
}
bottomWall
{
type empty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
backWedge
{
type wedge;
faces
(
(3 5 2 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************ ************* //
controlDict File:
/*--------------------------------*- 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.01;
writeControl timeStep;
writeInterval 20;
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: 7
\/ M anipulation |
*---------------------------------------------------------- -----------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.92e-7;
// ************************************************************ ************* //
Pressure file:
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
topWall
{
type zeroGradient;
}
bottomWall
{
type empty;
}
frontWedge
{
type wedge;
}
backWedge
{
type wedge;
}
}
// ************************************************************ ************* //
Velocity file:
/*--------------------------------*- 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.187 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.187 0 0);
}
outlet
{
type zeroGradient;
}
topWall
{
type noSlip;
}
bottomWall
{
type empty;
}
frontWedge
{
type wedge;
}
backWedge
{
type wedge;
}
}
// ************************************************************ ************* //
6) Postprocessing Results:
Velocity Profile
At the inlet
At X = 0.1, 0.4, 0.9
At the outlet (Fully developed profile):
At X = 1,4, 1.5
Analytical Velocity Profile for Fully developed region:
Inference from the Graph: From the above velocity profiles, we can see that, for entrance condition the velocity profile changes, but fr fully develop region the velocity profile remain constant, which satisfy Haigen Poisseuile statement.
Pressure Drop
at Inlet
at outlet
At X = 1.2-1.5(For fully developed region) & x=0-1.5 (Over entire length
Inference from the Pressure drop profile:
1) As the flow enters in fully develop region, the pressure over the length of 0.3m is 0.01623 Pa.
2) As the pressure is decreasing across the pipe , Hence at the outlet pressure condition our pressure difference is equal 0, equal to ambient pressure.
Shear Stress
By Computational
By Analytical
![]()
Coding by Matlab to generate Velocity and shear stress profile analytically:
%Code to calculate velocity profile
clear all
close all
clc
%Inputs
D = 0.01; %diameter of pipe
R = 0.005; %radius
r = linspace(0,R,100) ; %radius array
k = linspace(0,R,99) ; %radius array
u_max = 0.375; %maximum velocity
mu = 8.9e-4; %Dynamic viscosity
for i=1:length(r)
U(i) = u_max*(1-((r(i)/R)^2))
end
for i = 1:length(r)-1
dU = U(i+1)-U(i);
dR = r(i+1)-r(i);
tau(i) = mu*(dU/dR);
end
%Plotting
figure(1)
plot(r,U)
xlabel('Radial distance')
ylabel('Velocity Profile')
legend('Velocity profile wrt radial distance')
title('analytical velocity profile')
grid
figure(2)
plot(k,tau)
xlabel('Radial distance')
ylabel('Shear stress')
legend('Shear stress wrt radial distance')
title('analytical Shear stress')
set(gca,'YDIR','reverse')
grid
Analytical Plot:
Inference: Shear stress is maximum at the surface of the pipe, while it is minimum at the center of the pipe, due to viscosity.
Comparison between Analytical and Computational results
|
Parameter |
Analytical |
Computational |
Error(%) |
|
Maximum Velocity |
0.375 |
0.371 |
0.4% |
|
Pressure Drop |
0.01602 |
0.01623 |
0.02% |
|
Shear stress |
0.137 |
0.134 |
0.3% |
Conclusion:
1) The maximum velocity is at the center of the pipe, whereas minimum velocity is at the pipe surface due to viscosity.
2) There is a change in velocity in the Hydrodynamic entrance region, but it gets constant in a fully developed region.
3) The pressure difference is maximum at the inlet of the pipe, whereas it gets 0 at the exit of the pipe.
4) Shear stress is maximum at the surface of the pipe, while it is minimum at the center of the pipe.