Simulation of Flow through a pipe in Matlab OpenFoam
velocity profile varies according to the spatial distance in the pipe
shear stress variation is the spatial distance
Entry length of laminar flow is given by
`L=0.05*Re*D_h`
Hydraulic diameter is given by
`D_h=4A/P`
`A=π/4 D^2`
`P=πD`
By substituting A,P in `D_h`
`D_h=D`
`Re=vD/ϑ`
Assuming Diameter=1/8”
Where Inner Diameter=8.563mm
D=8.563/1000=0.008563
Where
Given condition as water
At T=30c (density)rho=995.65 kg/m3
Re=2100
the wedge angle is assumed as 4degree
by calculated from Hagen poiseuille\'s equation
Hydraulic pressure drop
`dP=(32*mu*L*v_(avg))/d^2;`
average velocity
`v_(avg)=(Re*nu)/d;`
maximum velocity
`v_max=2*v_(avg);`
kinematic pressure
`P_k=(dP)/rho;`
Creating a Matlab Program for OpenFoam
clear all;
close all;
clc;
% % Diameter of pipe in inch to meter
% % P=3.14*d;
% % A=(3.14/4)*(d^2);
% % D=4*A/P=d;
d=8.563/1000; %% which is simply ID of 1/8\" pipe
r=d/2;
theta=4; % wedge angle 10,25,45
% % Arc length
% % property of water at 30c
al=2*3.14*r*(theta/360);
rho=995.65;
mu=0.798*10^-3;
nu=mu/rho;
Re=2100; %% laminar flow
L_ent=0.05*Re*d;
L=L_ent+0.1009; %% make L_ent<L to capture velocity profile at upto fully develop ehere there is no changes after L_ent
% % average velocity
v_avg=Re*nu/d;
% % maximum velocity
v_max=2*v_avg;
% % Hydrolic pressur drop
dP=32*mu*L*v_avg/d^2;
% % kinamatic pressure
dP_k=dP/rho;
dx=L/(15*2000);
% % calculating vertex
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)];
% % creating blockmesh
header1=\"/*--------------------------------*- C++ -*----------------------------------*\\\";
header2=\" ========= | \";
header3=\" \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \";
header4=\" \\\\ / O peration | Website: https://openfoam.org \";
header5=\" \\\\ / A nd | Version: 7 \";
header6=\" \\\\/ M anipulation | \";
header7=\"\\*---------------------------------------------------------------------------*/\";
D1=\" version 2.0;\";
D2=\" format ascii;\";
D3=\" class dictionary;\";
D4=\" object blockMeshDict;\";
D5=\"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\";
c0=\"convertToMeters 1;\";
blocks1=\" hex (0 1 2 0 3 4 5 3) (15 1 2000) simpleGrading (1 1 1)\";
footer=\"// ************************************************************************* //\";
b4=blanks(4);
b8=blanks(8);
b12=blanks(12);
% % printing file as a file format
fp=fopen(\'blockMeshDict.txt\',\'wt\');
fprintf(fp,\'%s\\n\',header1);
fprintf(fp,\'%s\\n\',header2);
fprintf(fp,\'%s\\n\',header3);
fprintf(fp,\'%s\\n\',header4);
fprintf(fp,\'%s\\n\',header5);
fprintf(fp,\'%s\\n\',header6);
fprintf(fp,\'%s\\n\',header7);
fprintf(fp,\'FoamFile\\n{\\n\');
fprintf(fp,\'%s\\n\',D1);
fprintf(fp,\'%s\\n\',D2);
fprintf(fp,\'%s\\n\',D3);
fprintf(fp,\'%s\\n}\\n\',D4);
fprintf(fp,\'%s\\n\',D5);
fprintf(fp,\'\\n%s\\n\',c0);
% % printing vertices
fprintf(fp,\'\\nvertices\\n(\\n\');
fprintf(fp,\'\\t(%d %d %d)\\n\',v0);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v1);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v2);
fprintf(fp,\'\\t(%g %g %g)\\n\',v3);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v4);
fprintf(fp,\'\\t(%g %0.4f %0.4f)\\n\',v5);
% % printing blocks
fprintf(fp,\');\\n\\nblocks\\n(\\n\')
fprintf(fp,\'%s\\n);\\n\\n\',blocks1)
% % printing edges
fprintf(fp,\'edges\\n(\\n\');
fprintf(fp,\'%sarc 1 2 (0 %0.3f 0)\\n\',b4,r);
fprintf(fp,\'%sarc 4 5 (%g %0.3f 0)\\n\',b4,L,r);
% % printing boundary
fprintf(fp,\');\\n\\nboundary\\n(\\n\');
fprintf(fp,\'%saxis\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype empty;\\n%sfaces\\n%s(\\n%s(0 3 3 0)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(fp,\'%swall\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype wall;\\n%sfaces\\n%s(\\n% (2 5 4 1)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(fp,\'%sinlet\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype patch;\\n%sfaces\\n%s(\\n%s(0 0 2 1)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(fp,\'%soutlet\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype patch;\\n%sfaces\\n%s(\\n%s(3 4 5 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(fp,\'%sFront\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 3 5 2)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(fp,\'%sBack\\n%s{\\n\',b4,b4);
fprintf(fp,\'%stype wedge;\\n%sfaces\\n%s(\\n%s(0 1 4 3)\\n%s);\\n%s}\\n\',b8,b8,b8,b12,b8,b4);
fprintf(fp,\');\\n\\nmergePatchPairs\\n(\\n);\\n\\n\');
fprintf(fp,\'%s\',footer);
OpenFoam Output File
System folder
BlockMeshDict
/*--------------------------------*- 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 1;
vertices
(
(0 0 0)
(0 0.0043 -0.0001)
(0 0.0043 0.0001)
(1.00002 0 0)
(1.00002 0.0043 -0.0001)
(1.00002 0.0043 0.0001)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (15 1 2000) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.004 0)
arc 4 5 (1.00002 0.004 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
Front
{
type wedge;
faces
(
(0 3 5 2)
);
}
Back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //controlDict
adjust CFL number maximum as 0.5 to run the simulation for all simulation and run up to dt=500
folder 0
P
/*--------------------------------*- 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.0685;
}
wall
{
type zeroGradient;
}
axis
{
type empty;
}
Front
{
type wedge;
}
Back
{
type wedge;
}
}
// ************************************************************************* //
U
/*--------------------------------*- 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.1966 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1966 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type noSlip;
}
axis
{
type empty;
}
Front
{
type wedge;
}
Back
{
type wedge;
}
}
// ************************************************************************* //
Folder Constant
transportProperties
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu [0 2 -1 0 0 0 0] 8.014e-07;
// ************************************************************************* //
turbulenceProperties
/*--------------------------------*- 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 turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType laminar;
RAS
{
// Tested with kEpsilon, realizableKE, kOmega, kOmegaSST, v2f,
// ShihQuadraticKE, LienCubicKE.
RASModel kEpsilon;
turbulence off;
printCoeffs off;
}
// ************************************************************************* //Output Result
convergency graph
Velocity Magnitude contour
Front View
side View
Pressure contour
Front View
Side View
Plots
Velocity Magnitude
centre of pipe entre velocity profile variation
variable distance Radial Velocity profile Variation
Pressure
Pressure variation at the centre of the pipe to entire span length
shear stress
from Newtonian flow
`τ=μ*(∂U)/(∂r)`
the value of U is calculated from OpenFOAM result at fully develope flow region
the value is discretized in 1000 nodes
clear all;
close all;
clc;
mu=0.798*10^-3;
resw=xlsread(\'s0.95.csv\');
vw=resw(:,1);
raw=resw(:,2);
ress=xlsread(\'w0.95.csv\');
vs=(ress(:,1));
ras=(ress(:,2));
thuw=zeros(length(vw),1);
thus=zeros(length(vs),1);
for i=2:length(vw)
thuw(i)=mu.*((vw(i-1)-vw(i))/(raw(i-1)-(raw(i))));
end
plot(abs(thuw),abs(raw),\'color\',\'b\');
title(\'shear stress Vs spatial Radial at Fully Develope region\');
xlabel(\'shear stress\');
ylabel(\'spatial Radial\');Conclusion
From Hagen poiseuille\'s equation, kinamic pressure is calculated and given input for transient simulation and CFL is maintained below 0.5 & iterated till converged in OpenFOAM, the Aspect ratio of the mesh is maintained below 5 for better accuracy of the mesh quality where the physical property of profile variation is obtained in graphical format & contour. After fully develop region there is no velocity profile variation.