Week 11 - Simulation of Flow through a pipe in OpenFoam
AIM:-
1. Calculation of flow variables.
2. Generation of the blockMeshDict using the wedge boundary condition. create the blockMeshDict file in MATLAB such that the mesh is autoatically created for various wedge angles.
3. Simulating the pipe flow in openFoam using icoFoa solver.
4. Post process the velocity profile and shear strass and comparfe the simulation results with analytical results obtained using the hagen poiseuille's equation.Post processing is carried out in the paraVIEW.
Hagen Poiseuuiles Equation:-
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille 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.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter; and there is no acceleration of fluid in the pipe. For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent, leading to larger pressure drops than calculated by the Hagen–Poiseuille equation.
In standard fluid-kinetics notation
`Deltap=(8*mu*L*Q)/A^2`
where,
`Deltap` is the pressure difference between the two ends,
`L`is the length of pipe,
`mu`is the dynamic viscosity,
`Q`is the vloumetric flow rate.
`R` is the pipe radius.
`A`is the cross section of pipe.
The equation does not hold close to the pipe entrance.
This is the diagramatic representation of flow through pipe.This is what the fluid profile looks like.
The layer in which the shearing viscous forces are significant, is called the boundry layer. This boundary layer is a hypothetical concept. It divides the flow in pipe into two regions:
- Boundary layer region: The region in which viscous effects and the velocity changes are significant.
- The irrotational (core) flow region: The region in which viscous effects and velocity changes are negligible, also known as the inviscid core.
When the fluid just enters the pipe, the thickness of the boundary layer gradually increases from zero moving in the direction of fluid flow and eventually reaches the pipe center and fills the entire pipe. This region from the entrance of pipe to the point where the boundary layer covers the entire pipe is termed as the hydrodynamic entrance region and length of the pipe in this region is termed as the hydrodynamic entry length. In this region the velocity profile develops and thus the flow is called the hydrodynamically developing flow. After this region, the velocity profile is fully developed and continues unchanged. This region is called the hydrodynamically fully developed region.
Assumptions:-
Fluid obey's the neuton's law of viscosity(`tau=mu*((du)/dy)`)
No slip of the fluid particles at the boundary.
Given:-
- Reynlods number :2100
- Working fluid :water .
- Fluid temperature:25 `overset(circ)C`
- Diameter of the pipe =0.02m
- Dynamic viscosity(`nu`)=`0.89times10^-3`
- Kinematic Viscosity =`0.89times10^-6``m^2s^-1`
- Also `DeltaP= (32*mu*L*V_("avg"))/(D^2)`
- `theta=4^o`
- L=2.52m
Formulae:-
`V(r)=V_max *(1-r^2/R^2)`
`tau(r)=(2*mu*V_max*r)/R^2`
`"Length of pipe"=0.05*Re*D`
`"Velocity(v)"=(Re*mu)/(rho*D)`
`V_max =2*V`
`"Pressure(P)"=(32*mu*V*L)/D^2`
`"Kinematic pressure"=P/rho`
`"Maximum shear stress"(tau_max)=(2*mu*V_max)/R`
Analytical Solution:-
- Average Velocity : `V_("avg") = (R_e*mu)/(rho*D)`= 0.09373`ms^-1`
- Maximum velocity :`V_max = 2*V_("avg")`= 0.18746`ms^-1`
- Pressure Drop = 16.817pa
- Kinematic Pressure Drop : `(DeltaP)/(rho` =0.016817m.
- We use the incompressible folder and we navigate to then use icoFoam for solving the problem.we copy the folder cavity and then bring it torun folder.and make the following changes for different files.
- set of commands:
- cp -r cavity $FOAM_RUN/simulation
- cd $FOAM_RUN/simulation/cavity
- gedit -filename.
- finally ,blockMesh and icoFoam for creating mesh and running simulation respectively.
- Then open paraView for plotting and other purpose of running the simulation.
Matlab code for automatic generation of the blockMeshDict file:-
clear all
close all
clc
L=1;
D=0.004;
r=D/2;
theta=3;
nx=500;
ny=50;
nz=1;
line1='/*----------------------------------------*- C++ -*-------------------------------------*\';
line2='| ========= | |';
line3='| \\ / F ield | openFoam: The open source CFD Toolbox |';
line4='| \\ / O peration | Version: 9 |';
line5='| \\ / A nd | Web: www.openFoam.org |';
line6='| \\/ M anipulation | |';
line7='\*--------------------------------------------------------------------------------------*/';
bmd=fopen('blockMeshDict','w')
fprintf(bmd,'%s\n',line1);
fprintf(bmd,'%s\n',line2);
fprintf(bmd,'%s\n',line3);
fprintf(bmd,'%s\n',line4);
fprintf(bmd,'%s\n',line5);
fprintf(bmd,'%s\n',line6);
fprintf(bmd,'%s\n',line7);
fprintf(bmd,'Foamfile\n{\n');
fprintf(bmd,'%12s\t 2.0;\n','version');
fprintf(bmd,'%11s\t ascii;\n','format');
fprintf(bmd,'%10s\t dictionary;\n','class');
fprintf(bmd,'%11s\t blockMeshDict;\n','object');
fprintf(bmd,'}\n');
line8='//**************************************************************************************//';
fprintf(bmd,'%s\n\n',line8);
fprintf(bmd,'convertToMeters 1;\n\n');
fprintf(bmd,'vertices\n')
fprintf(bmd,'(\n')
fprintf(bmd,'\t (%d %d %d)\n',0,0,0); %point 0
fprintf(bmd,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2)); %point 1
fprintf(bmd,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2)); %point 2
fprintf(bmd,'\t (%d %d %d)\n',L,0,0); %point 3
fprintf(bmd,'\t (%d %d %d)\n',L,r*cosd(theta/2),r*sind(theta/2)); %point 4
fprintf(bmd,'\t (%d %d %d)\n',L,r*cosd(theta/2),-r*sind(theta/2)); %point 5
fprintf(bmd,');\n\n');
fprintf(bmd,'blocks\n');
fprintf(bmd,'(\n');
fprintf(bmd,' hex (0 3 5 2 0 3 4 1) (%d %d %d) ',nx,ny,nz);
fprintf(bmd,'simpleGrading(1 0.1 1)\n\n');
fprintf(bmd,');\n\n');
fprintf(bmd,'edges\n(\n\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(bmd,'\tarc %d %d (%d %d %d)\n',4,5,L,r,0);
fprintf(bmd,');\n\n');
fprintf(bmd,'boundary\n(\n');
fprintf(bmd,'\t inlet\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type patch;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t(\n');
fprintf(bmd,'\t\t\t (0 1 2 0)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');
fprintf(bmd,'\t outlet\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type patch;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t (\n');
fprintf(bmd,'\t\t\t (3 5 4 3)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');
fprintf(bmd,'\t top\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type wall;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t (\n');
fprintf(bmd,'\t\t\t (1 4 5 2)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');
fprintf(bmd,'\t back\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type wedge;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t (\n');
fprintf(bmd,'\t\t\t (0 2 5 3)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');
fprintf(bmd,'\t front\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type wedge;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t (\n');
fprintf(bmd,'\t\t\t (0 3 4 1)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');
fprintf(bmd,'\t axis\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type empty;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t (\n');
fprintf(bmd,'\t\t\t (0 3 3 0)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');
fprintf(bmd,');\n\n')
fprintf(bmd,'mergePatchPairs\n(\n');
fprintf(bmd,');\n')
fprintf(bmd,'%s\n',line8);
fclose(bmd)blockMeshDict file generated for wedge angle =`3^@`
/*----------------------------------------*- C++ -*-------------------------------------*\
| ========= | |
| \\ / F ield | openFoam: The open source CFD Toolbox |
| \\ / O peration | Version: 9 |
| \\ / A nd | Web: www.openFoam.org |
| \\/ 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)
(1 0 0)
(1 9.996573e-03 2.617695e-04)
(1 9.996573e-03 -2.617695e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (500 50 1) simpleGrading(1 0.1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (1 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
//**************************************************************************************//
Velocity boundary condition:-
/*--------------------------------*- 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.46725 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
back
{
type wedge;
}
front
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Pressure boundary condition:-
/*--------------------------------*- 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.3493161;
}
top
{
type zeroGradient;
}
back
{
type wedge;
}
front
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //OpenFoam case setUp:-
A wedge boundary condition is used in this case as it is cylinder and flow is axi-symmetric along the axis.
After updating the blockMeshDict and pressure and velocity boundary conditions file simulation is ran using the icoFoam solver which is used to simulate the turbulent, incompressible flows.
Base mesh for wedge angle=`3^@`
Results for the wedge angle=`3^@`
Velocity contours:-
Pressure contours:-
Velocity profile At 0.01m:-
the velocity profile near the entrance of the pipe is flat. It means the flow is still developing.
Velocity profile At 0.1m:-
velocity profile At 0.43m:-
velocity profile At 0.75m:-
A parabolic velocity profile is obtained after a certain distance in the pipe. This profiles remains unchanged thereafter. it concludes that the length of the pipe as calculated theoritically is suffiecient to get a fully developed flow at the outlet of the pipe.
The shear stress is calculated from the given formula
`tau =mu*(du)/(dy) Pa`