Simulation of Flow through a pipe in OpenFoam and MATLAB
Objective:
In this challenge, you will simulate an axisymmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. You can do a transient or steady-state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.
You will also write code in Matlab to automate the generation of the blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc, and output a fully functional blockMeshDict file.
Expected Results
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressure drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
Hint
- You can create a wedge from a hex block by collapsing its bottom face (by merging both sets of adjacent vertices)
- The wedge angle should be less than 5 degrees
- You can assume the pipe diameter
Solutions:
Laminar flow, type of fluid (gas or liquid) flow in which the fluid travels smoothly or in regular paths, in contrast to turbulent flow, in which the fluid undergoes irregular fluctuations and mixing. In laminar flow, sometimes called streamline flow, the velocity, pressure, and other flow properties at each point in the fluid remain constant. Laminar flow over a horizontal surface may be thought of as consisting of thin layers, or laminae, all parallel to each other. The fluid in contact with the horizontal surface is stationary, but all the other layers slide over each other. A deck of new cards, as a rough analogy, may be made to “flow” laminarly.
Laminar flow in a straight pipe may be considered as the relative motion of a set of concentric cylinders of fluid, the outside one fixed at the pipe wall and the others moving at increasing speeds as the center of the pipe is approached. Smoke rising in a straight path from a cigarette is undergoing laminar flow. After rising a small distance, the smoke usually changes to turbulent flow, as it eddies and swirls from its regular path.
Laminar flow is common only in cases in which the flow channel is relatively small, the fluid is moving slowly, and its viscosity is relatively high. Oil flow through a thin tube or blood flow through capillaries is laminar. Most other kinds of fluid flow are turbulent except near solid boundaries, where the flow is often laminar, especially in a thin layer just adjacent to the surface. See fluid mechanics.
Hydrodynamic entrance length:
For internal flow regime an entrance region is typical. In this region a nearly inviscid upstream flow converges and enters the tube. To characterize this region the hydrodynamic entrance length is introduced and is approximately equal to:
The maximum hydrodynamic entrance length, at ReD,crit = 2300 (laminar flow), is Le = 138d, where D is the diameter of the pipe. This is the longest development length possible. In turbulent flow, the boundary layers grow faster, and Le is relatively shorter. For any given problem, Le / D has to be checked to see if Le is negligible when compared to the pipe length. At a finite distance from the entrance, the entrance effects may be neglected, because the boundary layers merge and the inviscid core disappears. The tube flow is then fully developed.
In principle, the speed of the flow within the pipe is not constant over the cross-section. There is friction between the fluid and the pipe, which is why the flow velocity near the wall is lower than in the middle of the pipe. At the wall, the fluid even adheres to the wall due to the adhesive forces. Thus, there is no relative speed between fluid and wall. This is also called the no-slip condition. But the individual fluid layers also have friction between each other due to the viscosity of the fluid. This leads to the formation of a certain velocity profile, the course of which is yet to be determined. It should be noted, however, that the flow velocity is maximum in the middle of the pipe and then drops to zero as far as the wall.
The flow resistance to be overcome, which is necessary to shear the individual fluid layers against each other, depends on the viscosityof the fluid. The flow resistance is expressed by the shear stress τ, i.e. as the force per unit area required to move a fluid layer. This shear stress depends on how much the fluid layers are displaced against each other. This in turn is expressed by the velocity gradient perpendicular to the flow direction, i.e. by the slope of the velocity profile in the radial direction dv(r)/dr (also called shear rate). The mathematical relationship between the two quantities is described by the viscosity η (Newton’s law of fluid motion):
Mean flow velocity
With its parabolic velocity profile, a flow inside a pipe has no characteristic velocity in principle. Therefore, a mean flow velocity c is often defined to characterize such flows. The mean flow velocity is defined as a constant velocity over the entire cross-section which provides the same flow rate as the real velocity profile. With A=π⋅R2 as the inner pipe cross-section, the mean flow velocity c can be determined by the volume flow rate V*
In the above derivation, it was used that the quotient of traveled distance Δx and time Δt just corresponds to the mean flow velocity. Equating the equations finally provides the mean flow velocity
The following relationship applies between the maximum flow velocity Vmax and the mean flow velocity
Pressure loss
As already explained, the pressure gradient dp/dx in the upper equations represents the drive for the flow of the fluid. According to the equation, this directly influences the flow rate. Conversely, this means that a pressure gradient must exist for the desired flow rate. This requires a pump that ultimately generates this pressure gradient (at this point the negative sign is omitted since only the amounts are decisive anyway)
For a friction-free flow, on the other hand, no such pressure drop would have to be applied in the steady-state to keep the fluid flowing. Once set in motion, the fluid would flow through the pipe without loss of speed. The application of the pressure is therefore due to the friction loss, which is equivalent to a pressure drop. In the steady case, the pressure created by a pump must fully compensate for the pressure loss in order to keep the fluid flowing. The above formula is therefore equivalent to the pressure loss Δpl, which occurs over the pipe length ΔL:
Now the following calculation are done for analytical solutions of the Hagen Poiseuille law.
- Reynolds Number = 2100
- Kinematuc viscosity = 1e-6
- Dynamic Viscosity = 0.00089 N-s/m2
- Density of water = 1000Kg/m3
- Diameter of pipe = 20mm
- Hydrodynamic length Le = 00.06ReD = 2.53m
- Average velocity Vavg = Re*v/D
= 0.09345ms-1
- Maximum velocity Vmax = 2Vavg
0.01869ms-1
- Pressure
p = 8VavgmuL/R2 = 18.763Nm-2
In OpenFOAM, we have to provide the pressure boundary condition in terms of kinematic pressure and for that to begin we are going to write code to develop a Block mesh file for OpenFOAM to execute.
Pk = DeltaP/p
= 0.018m2s-2
Matlab Code to create blockmesh dict file
%Laminar incomprssible flow over channel
%MATLAB code to creat a blockMeshDict file for Openfoam
%simulation of incomprssible flow and then analysis of Prssure and shear
%stress of Hagen-Posuilles Equation
clear all
close all
clc
%simulation parameters of laminar flow
Re = 2100;
rho = 1000;
D = 0.02;
R = D/2;
L = 0.06*Re*D;
mu = 0.00089;
nu = mu/rho;
theta = 2;
%Defining the Boundary condition
v_avg = Re*nu/D;
v_max = 2*v_avg;
l = L + 0.3;
deltap = (8*v_avg*mu*l)/ R^2;
%vertices and co-ordinates in the local coordinate system
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)];
grading = 0.2;
b1 = 50;
b2 = 1;
b3 = 200;
%creating BlockMesh file code for Openfoam in its official format
%the first section of OpenFOam decription
h1 =' /*--------------------------------*- C++ -*----------------------------------*';
h2 ='| ========= | |';
h3 ='| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
h4 ='| \ / O peration | Website: https://openfoam.org |';
h5 ='| \ / A nd | Version: 8 |';
h6 ='| \/ M anipulation | |';
h7 =' *---------------------------------------------------------------------------*/';
h8 ='|// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%writing and creating the Blockmesh file which is writable.
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 n',h7);
%second section of the file where the parameters are initilaized and
%simulated with time steps
fprintf(f1, 'Foamfile n {n');
fprintf(f1,'%16s %10sn','version','2.0;');
fprintf(f1,'%16s %12sn','format','ascii');
fprintf(f1,'%16s %17sn','class','dictionary;');
fprintf(f1,'%16s %20sn','object','blockMeshDict');
fprintf(f1,'%2sn','}');
%the third section of the paramters
fprintf(f1,'%2snn','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%2snn','converToMeters 1;');
fprintf(f1, 'vertices n(n');
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);
fprintf(f1, '%2snn',');');
%defining the flow parametrs in the file
fprintf(f1, 'blocks n (n');
fprintf(f1, 't %s (%d %d %d) %s (%d 1 1) n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3, 'simpleGrading',grading);
fprintf(f1, '%2snn', ');');
fprintf(f1, 'edges n (n');
fprintf(f1, 't %s (0 %d 0)n','arc 1 2',R);
fprintf(f1, 't %s (%d %d 0)n','arc 4 5',L,R);
fprintf(f1,' %2snn',');')
fprintf(f1,'boundary n (n');
fprintf(f1,'t %5s n','axis','{');
fprintf(f1,'t %13s %sn','type','empty;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 3 3 0)',');','}');
fprintf(f1,'t %5s n','wall','{');
fprintf(f1,'t %13s %sn','type','wall;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(2 5 4 1)',');','}');
fprintf(f1,'t %5s n','inlet','{');
fprintf(f1,'t %13s %sn','type','patch;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 0 2 1)',');','}');
fprintf(f1,'t %5s n','outlet','{');
fprintf(f1,'t %13s %sn','type','patch;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(3 4 5 3)',');','}');
fprintf(f1,'t %5s n','wedgeFront','{');
fprintf(f1,'t %13s %sn','type','wedge;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 3 5 2)',');','}');
fprintf(f1,'t %5s n','wedgeBack','{');
fprintf(f1,'t %13s %sn','type','wedge;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn %snn','(0 1 4 3)',');','}',');');
fprintf(f1,'mergePatchPairs n (n );n');
fprintf(f1,'%2snn', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
%finally closing the file
fclose(f1);
Below is how you find the directory in your terminal command window
skill_lync@skillync:~$ cd $FOAM_TUTORIALS
skill_lync@skillync:/opt/openfoam8/tutorials$ ll
total 84
drwxr-xr-x 17 root root 4096 Apr 18 12:45 ./
drwxr-xr-x 11 root root 4096 Apr 18 12:45 ../
-rwxr-xr-x 1 root root 504 Mar 16 18:22 Allclean*
-rwxr-xr-x 1 root root 3217 Mar 16 18:22 Allrun*
-rwxr-xr-x 1 root root 6048 Mar 16 18:22 Alltest*
drwxr-xr-x 5 root root 4096 Mar 8 18:19 basic/
drwxr-xr-x 9 root root 4096 Mar 8 18:19 combustion/
drwxr-xr-x 6 root root 4096 Mar 8 18:19 compressible/
drwxr-xr-x 4 root root 4096 Mar 8 18:19 discreteMethods/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 DNS/
drwxr-xr-x 4 root root 4096 Mar 8 18:19 electromagnetics/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 financial/
drwxr-xr-x 5 root root 4096 Mar 8 18:19 heatTransfer/
drwxr-xr-x 13 root root 4096 Mar 8 18:19 incompressible/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 IO/
drwxr-xr-x 9 root root 4096 Mar 8 18:19 lagrangian/
drwxr-xr-x 8 root root 4096 Mar 8 18:19 mesh/
drwxr-xr-x 13 root root 4096 Mar 8 18:19 multiphase/
drwxr-xr-x 4 root root 4096 Mar 8 18:19 resources/
drwxr-xr-x 4 root root 4096 Mar 8 18:19 stressAnalysis/
skill_lync@skillync:/opt/openfoam8/tutorials$ cd incompressible/
skill_lync@skillync:/opt/openfoam8/tutorials/incompressible$ ll
total 52
drwxr-xr-x 13 root root 4096 Mar 8 18:19 ./
drwxr-xr-x 17 root root 4096 Apr 18 12:45 ../
drwxr-xr-x 3 root root 4096 Mar 8 18:19 adjointShapeOptimizationFoam/
drwxr-xr-x 6 root root 4096 Mar 8 18:19 boundaryFoam/
drwxr-xr-x 4 root root 4096 Mar 8 18:19 icoFoam/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 nonNewtonianIcoFoam/
drwxr-xr-x 5 root root 4096 Mar 8 18:19 pimpleFoam/
drwxr-xr-x 5 root root 4096 Mar 8 18:19 pisoFoam/
drwxr-xr-x 5 root root 4096 Mar 8 18:19 porousSimpleFoam/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 shallowWaterFoam/
drwxr-xr-x 13 root root 4096 Mar 8 18:19 simpleFoam/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 SRFPimpleFoam/
drwxr-xr-x 3 root root 4096 Mar 8 18:19 SRFSimpleFoam/
skill_lync@skillync:/opt/openfoam8/tutorials/incompressible$ cd ll
bash: cd: ll: No such file or directory
skill_lync@skillync:/opt/openfoam8/tutorials/incompressible$ cd simpleFoam/
skill_lync@skillync:/opt/openfoam8/tutorials/incompressible/simpleFoam$ ll
total 52
drwxr-xr-x 13 root root 4096 Mar 8 18:19 ./
drwxr-xr-x 13 root root 4096 Mar 8 18:19 ../
drwxr-xr-x 5 root root 4096 Apr 18 12:45 airFoil2D/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 mixerVessel2D/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 motorBike/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 pipeCyclic/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 pitzDaily/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 pitzDailyExptInlet/
drwxr-xr-x 6 root root 4096 Apr 18 12:45 roomResidenceTime/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 rotorDisk/
drwxr-xr-x 6 root root 4096 Apr 18 12:45 T3A/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 turbineSiting/
drwxr-xr-x 5 root root 4096 Apr 18 12:45 windAroundBuildings/
skill_lync@skillync:/opt/openfoam8/tutorials/incompressible/simpleFoam$ cp -r pitzDaily/ $FOAM_RUN/laminartestThe following are the Laminar flow properties in a file which is in the folder of FOAM_RUN
which are BlockMeshFile, U properties, and P properties with control diet and Transport properties.
BlockMesh
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Note: this file is a Copy of $FOAM_TUTORIALS/resources/blockMesh/pitzDaily
converToMeters 1;
vertices
(
(0 0 0)
(0 9.823953e-03 -1.745241e-04)
(0 9.998477e-03 1.745241e-04)
(2.520000e+00 0 0)
(2.520000e+00 9.823953e-03 -1.745241e-04)
(2.520000e+00 9.998477e-03 1.745241e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.520000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchpairs
(
);
// ************************************************************************* //
ControlDict
/*--------------------------------*- 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 simpleFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
cacheTemporaryObjects
(
kEpsilon:G
);
functions
{
#includeFunc streamlines
#includeFunc writeObjects(kEpsilon:G)
}
// ************************************************************************* //
U properties
/*--------------------------------*- 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.09345 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type noSlip;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
P properties
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.018;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
Transport Properties
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu [0 2 -1 0 0 0 0] 1e-05;
// ************************************************************************* //
Output:
The parafoam view of the profile with surface mesh.
velocity profile.
Pressure:
Shear stress:
Conclusion:
1. The error gradually increases with an increase in wedge length whereas it decreases as well with finer mesh o the profile.
2. Haigen-Poiseuille equation satisfies the numerical analysis of the equation.
References:
Derivation of Hagen-Poiseuille equation for pipe flows with friction - tec-science
https://www.thermal-engineering.org/what-is-hydrodynamic-entrance-length-definition/