Simulation using OpenFOAM in Circular Pipe
In this project, we will be observing the flow of a fluid in a pipe using OpenFOAM.
Brief description of the physical scenario
Boundary Layer - When a real fluid flows past a solid body or solid wall, the fluid particles stick to the boundary and the condition is known as the no-slip condition. This means that the velocity of the layer which is adjacent to the stationary body will be the same as that of the boundary. If the boundary is stationary the velocity of the fluid layer will also be zero.
Further away from the boundary the fluid will be experiencing less frictional force as compared to the layer near the boundary. A velocity gradient `frac{du}{dy}`is formed as we move away from the boundary. This variation of velocity from zero to free stream velocity(`U`) in the direction normal to the surface of the boundary takes place in a small region of the boundary. This region is separated by a layer which is termed as the boundary layer. And the region is termed as boundary layer region.
In the case of flow over flat plates, the formation of the boundary layer starts at the edge of the plate. However, in the case of a pipe, the flow of fluid develops over a certain length of pipe before attaining fully developed flow state. This length is termed as entrance length, `L_e`.
Development of fully developed flow in a circular pipe is shown in fig.
Mathematically, entrance length is a function of Reynold's number and for a laminar flow in a circular pipe it is given as
`L_e = 0.05ReD`
Shear Stress of fully developed velocity profile - The fluid adjacent to the boundary will be experiencing the maximum shear stress while the fluid layer flowing at the centre of the pipe will be facing minimalistic resistance. Hence, at the centre, the velocity of the flow will be maximum and the resistance or the shear stress will be minimum.
Also, the wall shear stress is maximum at the inlet of the pipe and decreases along the length of the pipe. This is the reason for the highest pressure drop at the inlet of the pipe.
Hagen - Poiseuille Equation: This 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.
Let's assume `D = 0.01m = 10mm`
We know, `Reynold's Number, Re = frac{Dv_(avg)}{nu}`
For water at `27^circ C`
Kinematic Viscosity `nu = 8.532 times 10^-7 m^2/s`
Dynamic Viscosity `mu = 8.502 times 10^-4 Pa.s`
Density `rho = 997 kg /m^3`
For maintaining laminar flow `Re = 2100`
`rArr v_(avg) = frac{nuRe}{D}`
`rArr v_(avg) = frac{8.532 times 10^-7 times 2100}{0.01}`
`rArr v_(avg) = 0.179 m/s`
`v_(max) = 2v_(avg) = 0.358m/s`
Entrance Length , `L_e = 0.05ReD = 0.05 times 2100 times 0.01 = 1.05m`
Considering the length of pipe, `L = 1.5m`
Using Hagen - Poiseuille formula , `Delta P = frac{32v_(avg)muL}{D^2} = 73.04 Pa.s`
Kinematic Pressure difference `= frac{Delta P}{rho} = 0.0732`
Time taken `t = frac{1.5}{0.179} = 8.37 s`
Case Setup in OpenFOAM
Generating a mesh for a circular pipe and solving it computational can be resource-heavy. In order to study the flow, we will be simulating a section/wedge of the circular pipe and will be observing the values there. The same wedge is shown in the figure below:
wedge boundary condition:
The wedge boundary condition defines an axis-symmetric boundary condition. The model and flow must be axis-symmetric along a central line such that all physical variables of the flow have same value and distribution at a radius and angle.
The axis-symmetric ‘wedge’ boundary is specified by two planes that must be selected on separate sides of the domain on surfaces running along the central line.
Mesh:
Wedge Angle `alpha = 4^circ`
Simulation Time = `11s`
Time Step to maintain Courant No less than 1 = `0.00586s`
Shear Stress is calculated using Newton's law of viscosity `tau = mu frac{du}{dy}`
In this case `tau = mu frac{du}{dr}`
Program to generate a blockMeshDict file for a given wedge angle
close all;
clear all;
clc;
d = 0.01;
L = 1.5;
alpha = input("Enter wedge angle: ");
y_dist = (d/2) * sin(0.5*alpha*pi/180);
x_dist = (d/2) * cos(0.5*alpha*pi/180);
g_x = input("Enter Grading factor for x-axis: ");
g_y = input("Enter Grading factor for y-axis: ");
g_z = input("Enter Grading factor for z-axis: ");
foamFileString = "FoamFile {version 2.0;format ascii;class dictionary;object blockMeshDict;}";
scaleString = "convertToMeters 1;";
verticesString = sprintf("vertices (\
(0 0 0) \
(%d %d 0) \
(%d %d %d) \
(0 0 %d) \
(%d -%d 0) \
(%d -%d %d));",x_dist, y_dist,x_dist, y_dist, L, L,x_dist, y_dist, x_dist, y_dist, L);
blockString = sprintf("blocks (hex (0 4 1 0 3 5 2 3) (20 1 100) simpleGrading (%d %d %d) );", g_x,g_y,g_z);
edgesString = "edges ();";
if alpha < 5
conditionType = "wedge";
else
conditionType = "symmetry";
end
boundaryString = sprintf("boundary (front {type %s; faces ((0 1 2 3));} back {type %s; faces ((0 3 5 4));} tankWall {type wall;faces ((1 4 5 2));} inlet {type patch; faces ((0 4 1 0));} outlet {type patch;faces ((3 2 5 3));} axis {type empty;faces ((0 3 3 0));} );", conditionType, conditionType);
mergeString = "mergePatchPairs();";
blockMeshDictString = strcat(foamFileString,"\n", scaleString ,"\n" , verticesString , "\n" , blockString , "\n" , edgesString , "\n" , boundaryString , "\n" , mergeString);
filename = 'blockMeshDict';
fid = fopen(filename, 'w');
fputs(fid, blockMeshDictString);
fclose(fid);
disp('Success!! blockMeshDict file generated');blockMeshDict file
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0.0049970 0.00017450 0)
(0.0049970 0.00017450 1.5)
(0 0 1.5)
(0.0049970 -0.00017450 0)
(0.0049970 -0.00017450 1.5)
);
blocks
(
hex (0 4 1 0 3 5 2 3) (20 1 100) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
front //movingWall
{
type wedge;//wall;
faces
(
(0 1 2 3)
);
}
back //fixedWalls
{
type wedge;//wall;
faces
(
(0 3 5 4)
);
}
tankWall //frontAndBack
{
type wall;
faces
(
(1 4 5 2)
);
}
inlet
{
type patch;
faces
(
(0 4 1 0)
);
}
outlet
{
type patch;
faces
(
(3 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
controlDict file
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 11;
deltaT 0.00586;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //U
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ 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
{
front
{
type wedge;
}
back
{
type wedge;
}
tankWall
{
type noSlip;
}
inlet
{
type fixedValue;
value uniform (0 0 0.179); //0.179m/s
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //p
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
front
{
type wedge;
}
back
{
type wedge;
}
tankWall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0732;
}
axis
{
type empty;
}
}
// ************************************************************************* //transportProperties File
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.532e-7;
// ************************************************************************* //Program to generate Velocity Profiles and Shear Stress Profile
close all;
clear all;
clc;
#Reading values from csv
profile_0_2 = dlmread('length_0.2.csv', ',' , 1,0);
profile_1_05 = dlmread('length_1.05.csv', ',' , 1,0);
profile_1_2 = dlmread('length_1.2.csv', ',' , 1,0);
profile_1_4 = dlmread('length_1.4.csv', ',' , 1,0);
# Assigning values to variables
radial_distance = profile_0_2(:, 1);
vel_at_0_2 = profile_0_2(:, 2);
vel_at_1_05 = profile_1_05(:, 2);
vel_at_1_2 = profile_1_2(:, 2);
vel_at_1_4 = profile_1_4(:, 2);
f1 = figure(1);
# Plotting values for +ve x-axis
plot(radial_distance, vel_at_0_2,'k');hold on;
plot(radial_distance, vel_at_1_05, 'r');hold on;
plot(radial_distance, vel_at_1_2,'g');hold on;
plot(radial_distance, vel_at_1_4,'b');hold on;
legend('Profile at distance 0.2 from inlet', 'Profile at distance 1.05 from inlet','Profile at distance 1.2 from inlet', 'Profile at distance 1.4 from inlet');hold on;
# Plotting values for -ve x-axis
plot(-radial_distance, vel_at_0_2,'k');hold on;
plot(-radial_distance, vel_at_1_05,'r');hold on;
plot(-radial_distance, vel_at_1_2,'g');hold on;
plot(-radial_distance, vel_at_1_4,'b');hold on;
xlabel('Radial Distance');
ylabel('Velocity');
title('Velocity Profile at various points from inlet');
#Calculating shear stress at various points of dy
mu = 8.502e-4;
dr = diff(radial_distance);
du = abs(diff(vel_at_1_4));
tau = mu * du./dr;
f2 = figure(2);
#Plotting shear stress
plot(radial_distance(1:end-1),tau, 'k');hold on;
plot(-radial_distance(1:end-1),tau, 'k');hold on;
title('Shear Stress profile at a distance 1.4m from inlet');
xlabel('Radial Distance');
ylabel('Shear Stress');
#Saving the plots
saveas(f1, 'VelocityProfiles','jpg');
saveas(f2, 'ShearStressProfile','jpg');Velocity Profiles
Shear Stress Profile
Flow development from Inlet to Outlet of the pipe
Conclusion:
- From velocity profiles, it can be observed till the entrance length the profile is developing and after the entrance length(i.e. `L_e = 1.05`), we get a fully developed profile. With `L gt 1.05`, the variation in velocity is minimum. This confirms our analytical assertion.
- The velocity of flow calculated using H-P equation and the velocity of flow observed using simulation matches for fully developed flow as seen in velocity profiles.
- The shear stress profile also matches with our analytical deduction i.e. the value of shear stress is maximum near the boundary/wall and decrease as we move towards the centre.
- For the given wedge angle (`alpha = 4^circ`) and initial conditions, the solution remains stable with time step = 0.00586s.
- Analytical `v_(max)` = 0.358 m/s , Simulation of `v_(max)` = 0.355563 m/s obtained near outlet. Percentage Error = 0.68073%