Week 11 - Simulation of Flow through a pipe in OpenFoam
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.
Assumptions:-
Fluid obey's the neuton's law of viscosity(τ=μ⋅(dudy)τ=μ⋅(dudy))
No slip of the fluid particles at the boundary.
Flow variables:-
The working fluid is water.
Density =1000 kg⋅m−3kg⋅m-3
Dynamic viscosity = 8.90⋅10−48.90⋅10-4pa⋅spa⋅s
Reynolds number= 2100
diameter of pipe=0.004 m
Formulae:-
V(r)=Vmax⋅(1−r2R2)V(r)=Vmax⋅(1-r2R2)
τ(r)=2⋅μ⋅Vmax⋅rR2τ(r)=2⋅μ⋅Vmax⋅rR2
Length of pipe=0.05⋅Re⋅DLength of pipe=0.05⋅Re⋅D
Velocity(v)=Re⋅μρ⋅DVelocity(v)=Re⋅μρ⋅D
Vmax=2⋅VVmax=2⋅V
Pressure(P)=32⋅μ⋅V⋅LD2Pressure(P)=32⋅μ⋅V⋅LD2
Kinematic pressure=PρKinematic pressure=Pρ
Maximum shear stress(τmax)=2⋅μ⋅VmaxRMaximum shear stress(τmax)=2⋅μ⋅VmaxR
Calculated results:-
L=0.42 m
V=0.4686 m/s
Vmax=0.9373msVmax=0.9373ms
Kinematic pressure = 0.3493161
τmax=0.831705Nm2τmax=0.831705Nm2
since the minimum length required for fully developed flow is 0.42m, hence we take the pipe length as 1m.
Hydrodynamic Entrance Length :-
It is the distance a flow travels in the pipe before it becomes fully developed. Entrance length refers to the length of entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer. When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer changes with increased distance along the pipe.
For laminar flow, the hydrodynamic entrance length is given as-
L=0.05⋅Re⋅DorL=1.05mL=0.05⋅Re⋅DorL=1.05m
Hence the length of the pipe for simulation is taken as 1.2m.
Analytical solution by Hagen-Poiseuille equation:-
Average velocity=Vavg=Re⋅μρ⋅D=0.46725msVavg=Re⋅μρ⋅D=0.46725ms
maximum velocity=Vmax=2⋅Vavg=0.9345msecVmax=2⋅Vavg=0.9345msec
pressure drop =△P=32⋅μ⋅L⋅VavgD2=349.3161△P=32⋅μ⋅L⋅VavgD2=349.3161
Kinematic pressure drop=△Pk=△Pρ=0.3493161△Pk=△Pρ=0.3493161
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: 4.1 |';
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∘3∘
/*----------------------------------------*- 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 dictonary;
object blockMeshDict;
}
//**************************************************************************************//
convertToMeters 1;
vertices
(
(0 0 0)
(0 1.999315e-03 5.235390e-05)
(0 1.999315e-03 -5.235390e-05)
(1 0 0)
(1 1.999315e-03 5.235390e-05)
(1 1.999315e-03 -5.235390e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (500 50 1) simpleGrading(1 0.1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (1 2.000000e-03 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: 7
\\/ 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.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: 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.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∘3∘
Results for the wedge angle=3∘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
τ=μ⋅dudyPaτ=μ⋅dudyPa
Pressure Profile:-