Simulation of Flow through a pipe in OpenFoam
Objective- To use the icoFOAM solver and simulate the flow through a pipe by creating wedge.
To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes
Dimensions of Pipe:
Diameter(D)=10mm
To calculate length of pipe entry length formula is used
`Le=0.06⋅Re⋅D`
Re=2100
=0.06*2100*10=1260mm=12.6m
So the length is taken as
L=1500mm=1.5m
Calculation for inlet velocity:
`Re=(VD)/ν`
`ν=8.9xx10^(−7`
`V = (2100xx8.9xx10^-7) /.01`
Wedge angle=2
Mesh specification
- Number of cells along the x direction (longer dimension) = 400
- Number of cells along the y direction = 50
blockMeshDict 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;
vertices
(
(0 0 0)
(0 4.99924 -0.087262)
(0 4.99924 0.087262)
(1500 0 0)
(1500 4.99924 -0.087262)
(1500 4.99924 0.087262)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 400) simpleGrading (0.1 1 10)
);
edges
(
arc 1 2 (0 5 0)
arc 4 5 (1500 5 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
upperWall
{
type wall;
faces
(
(1 2 5 4)
);
}
lowerWall
{
type empty;
faces
(
(0 3 3 0)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
Back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //Automated blockmesh generation:
close all;
clear all;
clc;
theta=2;
L=1500;
dia=10;
r=dia/2;
y=r*cosd(theta/2);
z=r*sind(theta/2);
f1=fopen('blockmeshdict.txt','wt');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
fprintf (f1,'========= |n')
fprintf (f1,' / F ield | OpenFOAM: The Open Source CFD Toolboxn');
fprintf (f1,' / O peration | Website: https://openfoam.orgn');
fprintf (f1,' / A nd | Version: 7n');
fprintf (f1,' / M anipulation |n');
fprintf(f1, '')
fprintf(f1,'*---------------------------------------------------------------------------*/n');
fprintf(f1,'FoamFilen{n')
fprintf(f1,' version 2.0;n');
fprintf(f1,' format ascii;n');
fprintf(f1,' class dictionary;n');
fprintf(f1,' object blockMeshDict;n}n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //n');
fprintf(f1,'convertToMeters 0.001;nn');
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,' (0 0 0)n')
fprintf(f1,' (0 %g -%g)n',y,z)
fprintf(f1,' (0 %g %g)n',y,z)
fprintf(f1,' (%g 0 0)n',L)
fprintf(f1,' (%g %g -%g)n',L,y,z)
fprintf(f1,' (%g %g %g)n',L,y,z)
fprintf(f1,');n');
fprintf(f1,'blocksn');
fprintf(f1,'(n');
fprintf(f1,' hex (0 1 2 0 3 4 5 3) (50 1 400) simpleGrading (0.1 1 10)n')
fprintf(f1,');n');
fprintf(f1,'edgesn');
fprintf(f1,'(n');
fprintf(f1,' arc 1 2 (0 %g 0)n',r)
fprintf(f1,' arc 4 5 (%g %g 0)n',L,r)
fprintf(f1,');n');
fprintf(f1,'boundaryn');
fprintf(f1,'(n');
fprintf(f1,'inletn {n');
fprintf(f1,' type patch;n');
fprintf(f1,' facesn(n');
fprintf(f1,' (0 2 1 0)n);n }n');
fprintf(f1,'outletn{n');
fprintf(f1,' type patch;n');
fprintf(f1,' facesn (n');
fprintf(f1,' (3 4 5 3)n);n}n');
fprintf(f1,'upperWalln{n');
fprintf(f1,' type wall;n');
fprintf(f1,' facesn (n');
fprintf(f1,' (1 2 5 4)n);n}n');
fprintf(f1,'lowerWalln{n');
fprintf(f1,' type empty;n');
fprintf(f1,' facesn(n');
fprintf(f1,' (0 3 3 0)n);n}nn');
fprintf(f1,'frontn{n');
fprintf(f1,' type wedge;n');
fprintf(f1,' facesn(n');
fprintf(f1,' (0 3 5 2)n);n}nn');
fprintf(f1,' Backn{n');
fprintf(f1,' type wedge;n');
fprintf(f1,' facesn(n');
fprintf(f1,' (0 1 4 3)n);n}n);n');
fprintf(f1,'mergePatchPairsn(n);n');
fprintf(f1,'// ************************************************************************* //');
controlDict file
In this file time step and time of simulation is provided.The time-step was taken as 1e-3 which gives a courant number below 1.
The simulation was run for 9 seconds.
/*--------------------------------*- 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 9;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Transport Properties:
/*--------------------------------*- 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.9e-7;
// ************************************************************************* //fvShemes
The discretization scheme is provided here.
/*--------------------------------*- 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 fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //fvSolution
The iterative solver is included here.
/*--------------------------------*- 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 fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //Initial and Boundary conditions:
Initial and boundary conditions are given in p and U file located in the 0 folder.
The initial condition is indicated by the name Internal Field and is set to 0,which means pressure is 0 at the start.
Boundary condition is given in the boundary field
Intlet velocity is given as 0.1869 m/s and pressure is extrapolated from inner nodes
Oultet velocity is extrapolated from inner nodes and pressure is set to 0
No-slip boundary condition at top wall and empty at bottom walls.
Front and back wedge is set to wedge boundary conditions
p 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;
}
upperWall
{
type zeroGradient;
}
lowerWall
{
type empty;
}
front
{
type wedge;
}
Back
{
type wedge;
}
}
// ************************************************************************* //U 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.1869 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1869 0 0);
}
outlet
{
type zeroGradient;
}
upperWall
{
type noSlip;
}
lowerWall
{
type empty;
}
front
{
type wedge;
}
Back
{
type wedge;
}
}
// ************************************************************************* //Output:
Velocity contour
Velocity magnitude at 0.085m from inlet
Velocity profile at different locations:
Velocity profile during under developed flow:
At the inlet of pipe flow is completely under developed
Here also the flow is under developed
Slowly Velocity profile is changing
At 1m the flow is turning fully developed
As per entry length calulation flow will be fully developed after 1.26m and it can be seen at x=1.5m
Shear Stress distribution:
Pressure Drop:
Comparision with Haigen poissule equation:
Analytical solution,
According to Haigen poissule equation
`Δp=(32μLV)/(D^2)`
=`(32xx8.9xx10^(−4)xx0.1869xx1.50)/(.01^2)`
`=.0798m^2/(s^2)`
And
Vmax=2*Vavg=2*0.1869=0.3738 m/s
Simulation results,
`Deltap=0.1 m^2/s^2`
And
Vmax=0.375 m/s
So the simulation results matches closely to analytical haigen poissule eqn