Simulation of Flow through a pipe and Comparison between Wedge and Symmetrical Boundary Condition in terms of Wedge angle, Computation time and CFL number 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.
Objective:
- Write a Matlab program that takes an angle as input and generates a blockMesh file for the given angle. Angles to test 10,25,45
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 = input('Enter the value of wedge angle='); %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,'frontSymmetry');
fprintf(f1,'%s %sn',s3,'{');
fprintf(f1,'%s %sn',s3,'type symmetry;');
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,'backSymmetry');
fprintf(f1,'%s %sn',s3,'{');
fprintf(f1,'%s %sn',s3,'type symmetry;');
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:
For Theta = 10
/*--------------------------------*- 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.00498 0.000435)
(0 0.00498 -0.000435)
(1.5 0 0)
(1.5 0.00498 0.000435)
(1.5 0.00498 -0.000435)
);
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)
);
}
frontSymmetry
{
type symmetry;
faces
(
(0 1 4 3)
);
}
backSymmetry
{
type symmetry;
faces
(
(3 5 2 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************ ************* //
For Theta = 25 degree
/*--------------------------------*- 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.00488 0.00108)
(0 0.00488 -0.00108)
(1.5 0 0)
(1.5 0.00488 0.00108)
(1.5 0.00488 -0.00108)
);
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)
);
}
frontSymmetry
{
type symmetry;
faces
(
(0 1 4 3)
);
}
backSymmetry
{
type symmetry;
faces
(
(3 5 2 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************ ************* //
For Theta = 45 degree.
/*--------------------------------*- 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.00461 0.00191)
(0 0.00461 -0.00191)
(1.5 0 0)
(1.5 0.00461 0.00191)
(1.5 0.00461 -0.00191)
);
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)
);
}
frontSymmetry
{
type symmetry;
faces
(
(0 1 4 3)
);
}
backSymmetry
{
type symmetry;
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;
}
frontSymmetry
{
type symmetry;
}
backSymmetry
{
type symmetry;
}
}
// ************************************************************ ************* //
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;
}
frontSymmetry
{
type symmetry;
}
backSymmetry
{
type symmetry;
}
}
// ************************************************************ ************* //
6) Postprocessing Results:
Velocity Profile
For Theta = 10 degree
At X = 0.1 At X=0.9
For Theta = 25 degree
At X = 0.1 At X=0.9
For Theta = 45 degree
At X = 0.1 At X=0.9
At the outlet (Fully developed profile):
For Theta = 10 degree
At X = 1.1 At X=1.5
For Theta = 25 degree
At X = 1.1 At X=1.5
For Theta = 45 degree
At X = 1.3 At X=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
For Entire domain
For Theta =10,25,45 degrees
For Fully developed region
For Theta = 10 degree.
For Theta = 25 degree.
For Theta = 45 degree.
Inference from the Pressure drop profile:
1) As the flow enters in fully develop region, the change in pressure over the length of 0.3m is 0.01623 Pa.
2) As the pressure is decreasing aceoss 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
Plot of Analytical shear Stress:
Inference from the Plot: Shear Stress is maximum at the surface of the pipe due to Viscosity, whereas it is minimum at the center of the Pipe.
Comparison between analytical and numerical Results for Symmetry Boundary condition at different Wedge angles:
|
Wedge angle (degree) |
Maximum Velocity |
Pressure drop |
Shear Stress |
|
Theta = 10 |
0.360 |
0.01616 |
0.1271 |
|
Theta= 25 |
0.360 |
0.01688 |
0.1280 |
|
Theta = 45 |
0.3641 |
0.01910 |
0.129 |
Comparison between Wedge and Symmetrical Boundary Condition in terms of Computation time and CFL number:
|
Parameter |
Symmetry |
Wedge |
|
CFL Number |
0.375 |
0.258 |
|
Computation Time |
104.3 seconds |
300.4 seconds |
Conclusion:
1) Overall Wedge Boundary condition gives more accurate results than Symmetry Boundary Condition.
2) However, the Symmetry Boundary condition takes less computation time than the Wedge Boundary condition.