Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim: Simulation-of-Flow-through-a-pipe-in-OpenFoam.
Objective
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
Preprocessing
BlockMeshDict preparation
%clear all
close all
clc
Re=2100 %Reynolds number
rho=1000; %density water
mu=0.89e-3; %dynamic viscosity water at 25degC (mPa)
nu=mu/rho; %kinematic viscosity =mu/rho
R=0.01; %radius of pipe (m)
D=2*R;
L_e=0.05*Re*D %entrance length for laminar flow
theta=2.5; %axis symmetric max 5deg
Vavg=Re*nu/D % approximated from Reynolds equation for RE=2100
vmax=2*Vavg %max velocity
dp=8*mu*L_e*Vavg/(R^2); %pressure drop from Hagen's
dp_k=dp/rho %kinematic pressure difference for comparison against OpenFoam
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot preparation entites
r=linspace(-R,R,1001);
r1=linspace(0,R,1001);
dx=r1(2)-r1(1)
v=2*Vavg.*(1-(r.^2/R.^2)); %vel=f(r)
v_s=2*Vavg.*(1-(r1.^2/R.^2)); %vel=f(r)
%v_ss=flip(v_s)
vmax1=linspace(0,vmax,1001);
%plot analytical v=f(pipe radius)
figure(1)
plot(v,r);
ylabel('distance along pipe diameter (m)')
xlabel('Analytical Max Velocity (m/s)')
title('Analytical velocity: 0m radius=pipe centre');
%import max velocity profile from OpenFoam
figure(3)
plot(UProf(:,2),UProf(:,1))
ylabel('pipe radius (m)')
xlabel('Openfoam Max Velocity (m/s)')
title('OpenFoam velocity 0m radius=wall');
%Analytical Shear stress
ud=diff(v_s)
dr=diff(r1)
dtau=mu.*ud./dx;
figure(2)
plot(r1(1,2:end),dtau)
xlabel('pipe radius (m)')
ylabel('Analytical Shear Stress (Pa)')
title('Analytical Shear stress: 0m radius=pipe centre');
%Openfoam Shear stress
% ud_o=diff(flip(UProf(:,2)))
% dr_o=diff(r1)
% dtau=mu.*ud_o./dx;
% figure(4)
% plot(UProf(2:end,1),dtau)
% xlabel('pipe radius (m)')
% ylabel('Openfoam Shear Stress (Pa)')
% title('Openfoam Shear stress: 0m radius=wall'); %defining the verties for blockMesh input file generation
v1=[0 0 0];
v2=[-R*sind(theta) R*cosd(theta) 0];
v3=[+R*sind(theta) R*cosd(theta) 0];
v4=[0 0 L_e];
v5=[-R*sind(theta) R*cosd(theta) L_e];
v6=[+R*sind(theta) R*cosd(theta) L_e];
% mesh discretisation
nx=20
ny=1; % 1 cell thick according to userguide pg U-39
nz=500;
%blockMeshdict file generation
f=fopen('blockMeshDict','w')
fprintf(f,'%s\n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f,'%s\n','========= |');
fprintf(f,'%s\n','\\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f,'%s\n','\\ / O peration | Website: https://openfoam.org');
fprintf(f,'%s\n','\\ / A nd | Version: 8');
fprintf(f,'%s\n','\\/ M anipulation |');
fprintf(f,'%s\n','\*---------------------------------------------------------------------------*/');
fprintf(f,'%s\n','FoamFile');
fprintf(f,'%s\n','{');
fprintf(f,'%s\n','version 2.0;');
fprintf(f,'%s\n','format ascii;');
fprintf(f,'%s\n','class dictionary;');
fprintf(f,'%s\n','object blockMeshDict;');
fprintf(f,'%s\n','}');
fprintf(f,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f,'\n');
fprintf(f,'%s\n','convertToMeters 1;');
fprintf(f,'\n');
fprintf(f,'%s\n','vertices');
fprintf(f,'%s\n','(');
fprintf(f,'\t (%f %f %f)\n',v1(1), v1(2), v1(3));
fprintf(f,'\t (%f %f %f)\n',v2(1), v2(2), v2(3));
fprintf(f,'\t (%f %f %f)\n',v3(1), v3(2), v3(3));
fprintf(f,'\t (%f %f %f)\n',v4(1), v4(2), v4(3));
fprintf(f,'\t (%f %f %f)\n',v5(1), v5(2), v5(3));
fprintf(f,'\t (%f %f %f)\n',v6(1), v6(2), v6(3));
fprintf(f,'%s\n',');');
fprintf(f,'\n');
fprintf(f,'%s\n','blocks');
fprintf(f,'%s\n','(');
fprintf(f,'\t hex (0 2 1 0 3 5 4 3)\n');
fprintf(f,'\t (%d %d %d)\n',nx,ny,nz);
fprintf(f,'\t %s\n','simpleGrading (1 1 1)');
fprintf(f,'%s\n',');');
fprintf(f,'\n');
fprintf(f,'%s\n','edges');
fprintf(f,'%s\n','(');
fprintf(f,'%s\n',');');
fprintf(f,'\n');
fprintf(f,'%s\n','boundary');
fprintf(f,'%s\n','(');
fprintf(f,'\t inlet \n');
fprintf(f,'\t { \n');
fprintf(f,'\t \t type patch; \n');
fprintf(f,'\t \t faces \n');
fprintf(f,'\t \t ( \n');
fprintf(f,'\t \t (0 1 2 0) \n');
fprintf(f,'\t \t ); \n');
fprintf(f,'\t } \n');
fprintf(f,'\t outlet \n');
fprintf(f,'\t { \n');
fprintf(f,'\t \t type patch; \n');
fprintf(f,'\t \t faces \n');
fprintf(f,'\t \t ( \n');
fprintf(f,'\t \t (3 4 5 3) \n');
fprintf(f,'\t \t ); \n');
fprintf(f,'\t } \n');
fprintf(f,'\t upperwall \n');
fprintf(f,'\t { \n');
fprintf(f,'\t \t type wall; \n');
fprintf(f,'\t \t faces \n');
fprintf(f,'\t \t ( \n');
fprintf(f,'\t \t (1 4 5 2) \n');
fprintf(f,'\t \t ); \n');
fprintf(f,'\t } \n');
fprintf(f,'\t frontwall \n');
fprintf(f,'\t { \n');
fprintf(f,'\t \t type wedge; \n');
fprintf(f,'\t \t faces \n');
fprintf(f,'\t \t ( \n');
fprintf(f,'\t \t (0 3 4 1) \n');
fprintf(f,'\t \t ); \n');
fprintf(f,'\t } \n');
fprintf(f,'\t backwall \n');
fprintf(f,'\t { \n');
fprintf(f,'\t \t type wedge; \n');
fprintf(f,'\t \t faces \n');
fprintf(f,'\t \t ( \n');
fprintf(f,'\t \t (0 2 5 3) \n');
fprintf(f,'\t \t ); \n');
fprintf(f,'\t } \n');
fprintf(f,'%s\n',');');BlockMeshDict
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(-0.000436 0.009990 0.000000)
(0.000436 0.009990 0.000000)
(0.000000 0.000000 2.100000)
(-0.000436 0.009990 2.100000)
(0.000436 0.009990 2.100000)
);
blocks
(
hex (0 2 1 0 3 5 4 3)
(20 1 500)
simpleGrading (1 1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
upperwall
{
type wall;
faces
(
(1 4 5 2)
);
}
frontwall
{
type wedge;
faces
(
(0 3 4 1)
);
}
backwall
{
type wedge;
faces
(
(0 2 5 3)
);
}
);
Velocity Vector Field
/*--------------------------------*- 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 0 0.0934);
}
outlet
{
type zeroGradient;
}
upperwall
{
type noSlip;
}
backwall
{
type wedge;
}
frontwall
{
type wedge;
}
}
// ************************************************************************* //
Pressure Field
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0140;
}
upperwall
{
type zeroGradient;
}
backwall
{
type wedge;
}
frontwall
{
type wedge;
}
}
// ************************************************************************* //
Control Settings
/*--------------------------------*- 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 icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20;
deltaT 0.01;
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: 8
\\/ 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; //kinematic viscosity
// ************************************************************************* //
Analysis
We use Hagen-Poiseuille’s equations to characterise the flow in a circular pipe for Reynolds number Re=2100. It predicts the pressure drop in an incompressible Newtonian fluid in a long cylinder the length of which is determine by the entrance number Le=0.05*Re*D where D is the pipe diameter. This is correlated to a fully developed velocity profile as shown below.
The Hagen-Poiseuille equation is given as
The following parameters are given/calculated for the purpose of comparing the analytical against openfoam values.
Re=2100 %Reynolds number
rho=1000; %density water
mu=0.89e-3; %dynamic viscosity water at 25degC (mPa)
nu=mu/rho; %kinematic viscosity =mu/rho
R=0.01; %radius of pipe (m)
D=2*R;
L_e=0.05*Re*D=2.1 %entrance length for laminar flow
theta=2.5; %axis symmetric max 5deg
Vavg=Re*nu/D=0.0934 % approximated from Reynolds equation for RE=2100
vmax=2*Vavg=0,1869%max velocity
dp=8*mu*L_e*Vavg/(R^2)=13.9726; %pressure drop from Hagen's
dp_k=dp/rho:0.0140 %kinematic pressure difference for comparison against OpenFoam
Results/Conclusion
The following 2 plots show a comparison of the analytical against that predicted by Openfoam.
This plot shows the variation of the velocity
- At the wall in which velocity is 0 due to no slip condition -> max shear stress
- At the pipe centre, the veloctiy is max -> min shear stress
- Convergence of the velocity profile occurs approximate > 17secs. Total simulation time =20secs
- The calculated max velocity correlates with the analytical value.
This plots shows the analytical pressure drop which was applied as a boundary condition on the outlet patch confirming the model is set up correctly.
The fluid shear stress is defined as the rate of change of the velocity times the kinematic viscosity divided by the distance along the pipe radius. This is illustrated below.
Conditions of shear stress at the wall/pipe centre.
- The flow at the wall is minimum -> the gradient du/dx will be max -> shear stress will be max
- The flow at the centre is max -> gradient du/dx will be min -> shear stress will be min.
The following plot shows the analytical shear stress.
At this point in time the correct openfoam postprocessing procedure is missing to compare against this analytical curve.
