Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective – To simulate laminar flow through a part of a pipe and obtain the fully developed velocity profile and shear stress profile by creating a MATLAB script that creates the BlockMeshDict file for a given wedge angle.
Laminar flow through pipe
The nature of flow of a fluid is characterized by Reynolds number, a non-dimensional quantity that measures the ratio of the inertial and viscous forces.
The categories of the Reynolds number are as follows-
Re<2000 – Laminar
2000<Re<4000 -Transitional
Re>4000- Turbulent
Therefore, in this problem a Reynolds number of 1800 is chosen to ensure laminar flow.
For a Reynolds number lesser than 2000, the viscous forces will dampen out any turbulence that enters the flow and the fluid flows like sheets rubbing over each other.
For a flow through a cylindrical pipe the velocity of the fluid in contact with the walls of the pipe is assumed zero due to shear stresses, a condition known as NO-SLIP.
The velocities vary in a parabolic sense from the wall up to the axis of pipe. This is called a fully developed flow and the distance form the inlet (uniform profile) at which such a profile forms is called entrance length.
Hagen Poiseuille Law
In fluid dynamics the Hagen Poiseuille equation is used for obtaining the pressure drop for laminar flow in cylindrical pipe of constant cross section.
The pressure drop is given as..
`DeltaP=(32*mu*V*L)/(D*D)`
`DeltaP`- Pressure drop
`mu`-Absolute viscousity
V- Velocity
L-length from inlet
D-diameter of pipe
Assumptions-
1.The liquid must be a streamline flow witht the streamlines parallel to the axis of the cylinder.
2.The pressure must be constant over the circular cross section.
3.Velocity of the fluid must be zero at the wall, and the flow should strictly be laminar.
4.The fluid obeys Newton's law of viscousity - `tau=mu*delU/dely`
Formulas used-
1.Entrance length , Le = 0.05.Re.D
2.V(r)=`underset(max)V*(1 - r^2/R^2)`
This shows that in the fully developed region the velocity varies parabolically from the center(max) to the circumference (zero).
3.Velocity(V)(inlet)=`Re.mu/(rho.D)`
4.Kinematic pressure = `P/rho`
Data:
Liquid - water
Density - 997 kg/m3
Absolute Viscousity-0.00089 Pas
Diameter of pipe 0.004m
Length = 1m
Reynolds Number -1800
Pressure drop (kinematic) =0.7172 Pa
Inlet Velocity = 0.4017 m/s uniformly distributed
Max velocity- at center of fully developed = 0.8034 m/s
Max value of shear stress = 0.715 Pa
CFD software used - OpenFOAM.
Solver -icoFOAM
Mesh structure - nx:500 , ny:50 ,nz:1.
Geometry
The geometry used in this challenge is a wedge shape , a part of a cylinder . The main determinant of this geometry is the wedge angle . The angle is taken as 3 degrees.
The boundry surfaces are inlet , outlet , top , two sides and the axis.
Boundary Conditions-
Inlet:
velocity - 0.4017m/s
Pressure- zeroGradient.
Outlet:
Pressure -0.7172 Pa
Velocity Zerogradient
Mesh obtained-
Simulation Results-
Velocity profiles in time steps
Time elapsed -0s
Time elapsed -0.1s
Time elapsed-0.2s
Time elapsed-0.3s
Time elapsed-0.4s
Time elapsed -0.5s
Time elapsed -10s
Time elasped -25s
Shear stress profile in start and end time steps
Start-
End-
Plots Obtained-
Velocity distribution along different points -
Shear stress distributions along different points-
Velocity Distribution cumulative-
Shear Stress distribution cumulative-
Conclusion-
1. The velocity maximum value
Simulation - 0.8029
Analytical - 0.8034
2.The shear stress maximum value
Simulation - 0.6860
Analytical - 0.715
3. At about the middle section of the pipe at from 0.5 seconds the laminar flow velocity profile is obtained in which the velocity distribution is parabolic and the shear stress variation is a linear straight line.
MATLAB SCRIPT
clear all
close all
clc
reyno=1800; %Reynolds Number%
rho=997;
mu=0.00089;
d=0.004;
L = 1;
L_act = 0.05*reyno*d;
V=(reyno*mu/rho/d);
r=d/2;
pressure_drop = (32*mu*V*L/d/d/rho);
tetha=3;
v_max=2*V
tau_max=2*mu*2*V/r;
f1 = fopen('blockMeshDict.txt','w');
line1='/*--------------------------------*- C++ -*----------------------------------*\';
line2=' ========= |';
line3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4=' \\ / O peration | Website: https://openfoam.org';
line5=' \\ / A nd | Version: 7'
line6=' \\/ M anipulation '
line7='\*---------------------------------------------------------------------------*/'
fprintf(f1,'%s\n',line1);
fprintf(f1,'%s\n',line2);
fprintf(f1,'%s\n',line3);
fprintf(f1,'%s\n',line4);
fprintf(f1,'%s\n',line5);
fprintf(f1,'%s\n',line6);
fprintf(f1,'%s\n',line7);
fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n',' version 2.0;');
fprintf(f1,'%s\n',' format ascii;');
fprintf(f1,'%s\n',' class dictionary;');
fprintf(f1,'%s\n',' object blockMeshDict;');
fprintf(f1,'%s\n','}');
fprintf(f1,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\n');
fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t(0 0 0)\n');
fprintf(f1,'\t(0 %d %d)\n',r*cosd(tetha/2),r*sind(tetha/2));
fprintf(f1,'\t(0 %d %d)\n',r*cosd(tetha/2),(-1)*r*sind(tetha/2));
fprintf(f1,'\t(%d 0 0)\n',L);
fprintf(f1,'\t(%d %d %d)\n',L,r*cosd(tetha/2),r*sind(tetha/2));
fprintf(f1,'\t(%d %d %d)\n',L,r*cosd(tetha/2),(-1)*r*sind(tetha/2));
fprintf(f1,');');
fprintf(f1,'\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(');
fprintf(f1,'\t hex (%d %d %d %d %d %d %d %d) (%d %d %d) simpleGrading (%d %d %d)',0,3,5,2,0,3,4,1,500,50,1,1,1,1);
fprintf(f1,'\n\n);');
fprintf(f1,'\nedges');
fprintf(f1,'\n(');
fprintf(f1,'\t arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1,'\t arc %d %d (%d %d %d)\n',4,5,L,r,0);
fprintf(f1,'\n);');
fprintf(f1,'\nboundary\n(');
fprintf(f1,'\n\t inlet \n\t { \n\t\t type patch; \n\t\t faces \n\t\t ( \n\t\t\t (0 1 2 0) \n\t\t ); \n\t }');
fprintf(f1,'\n\t outlet \n\t { \n\t\t type patch; \n\t\t faces \n\t\t ( \n\t\t\t (3 5 4 3) \n\t\t ); \n\t }');
fprintf(f1,'\n\t top \n\t { \n\t\t type wall; \n\t\t faces \n\t\t ( \n\t\t\t (1 4 5 2) \n\t\t ); \n\t }');
fprintf(f1,'\n\t frontside \n\t { \n\t\t type wedge; \n\t\t faces \n\t\t ( \n\t\t\t (0 3 4 1) \n\t\t ); \n\t }');
fprintf(f1,'\n\t backside \n\t { \n\t\t type wedge; \n\t\t faces \n\t\t ( \n\t\t\t (0 2 5 3) \n\t\t ); \n\t }');
fprintf(f1,'\n\t centerline \n\t { \n\t\t type empty; \n\t\t faces \n\t\t ( \n\t\t\t (0 3 3 0) \n\t\t ); \n\t }');
fprintf(f1,'\n);');
fprintf(f1,'\n\nmergerPatchPairs\n(\n);');
fprintf(f1,'\n//************************************************************************* //');
fclose(f1);
MATLAB graph plotting
figure(1)
plot(UA.arc_length,UA.U_Magnitude)
hold on
plot(UB.arc_length,UB.U_Magnitude)
hold on
plot(UC.arc_length,UC.U_Magnitude)
hold on
plot(UD.arc_length,UD.U_Magnitude)
hold on
plot(UE.arc_length,UE.U_Magnitude)
set(gcf,'position',[10,10,2000,400]);
legend('At distance 0.1 m','At distance 0.25 m','At distance 0.5 m','At distance 0.75 m','At distance 0.95 m');
title('Velocity Distribution from the center axis to the surface');
xlabel('Radial distance form central axis');
ylabel('Velocity');
figure(2)
plot(UD.arc_length,UD.U_Magnitude)
xlabel('Radial distance form central axis');
ylabel('Velocity');
title('Velocity Distribution at 0.75 m distance');
set(gcf,'position',[10,10,2000,400]);
figure(3)
plot(UE.arc_length,UE.U_Magnitude)
xlabel('Radial distance form central axis');
ylabel('Velocity');
title('Velocity Distribution at 0.95 m distance');
set(gcf,'position',[10,10,2000,400]);
figure(4)
plot(TA.arc_length,TA.tau_Magnitude)
hold on
plot(TB.arc_length,TB.tau_Magnitude)
hold on
plot(TC.arc_length,TC.tau_Magnitude)
hold on
plot(TD.arc_length,TD.tau_Magnitude)
hold on
plot(TE.arc_length,TE.tau_Magnitude)
set(gcf,'position',[10,10,2000,400]);
legend('At distance 0.1 m','At distance 0.25 m','At distance 0.5 m','At distance 0.75 m','At distance 0.95 m');
title('Variation of Shear Stress along radial direction');
xlabel('Radial distance form central axis');
ylabel('Shear stress');
figure(5)
plot(TD.arc_length,TD.tau_Magnitude)
set(gcf,'position',[10,10,2000,400]);
title('Variation of Shear Stress along radial direction at 0.75 m form inlet');
xlabel('Radial distance form central axis');
ylabel('Shear stress');
figure(6)
plot(TE.arc_length,TE.tau_Magnitude)
set(gcf,'position',[10,10,2000,400]);
title('Variation of Shear Stress along radial direction at 0.95 m form inlet');
xlabel('Radial distance form central axis');
ylabel('Shear stress');
figure(7)
plot(UA.arc_length,UA.U_Magnitude)
xlabel('Radial distance form central axis');
ylabel('Velocity');
title('Velocity Distribution at 0.1 m distance');
set(gcf,'position',[10,10,2000,400]);
figure(8)
plot(UC.arc_length,UC.U_Magnitude)
xlabel('Radial distance form central axis');
ylabel('Velocity');
title('Velocity Distribution at 0.5 m distance');
set(gcf,'position',[10,10,2000,400]);
figure(9)
plot(TA.arc_length,TA.tau_Magnitude)
set(gcf,'position',[10,10,2000,400]);
title('Variation of Shear Stress along radial direction at 0.1 m form inlet');
xlabel('Radial distance form central axis');
ylabel('Shear stress');
figure(10)
plot(TC.arc_length,TC.tau_Magnitude)
set(gcf,'position',[10,10,2000,400]);
title('Variation of Shear Stress along radial direction at 0.5 m form inlet');
xlabel('Radial distance form central axis');
ylabel('Shear stress');
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 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 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)
);
}
frontside
{
type wedge;
faces
(
(0 3 4 1)
);
}
backside
{
type wedge;
faces
(
(0 2 5 3)
);
}
centerline
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergerPatchPairs
(
);
// ************************************************************************* //
controlDict
/*--------------------------------*- 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 25;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Pressure boundry conditions
/*--------------------------------*- 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.7172;
}
frontside
{
type wedge;
}
backside
{
type wedge;
}
top
{
type zeroGradient;
}
centerline
{
type empty;
}
}
// ************************************************************************* //
Velocity boundry 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.4017 0 0);
}
outlet
{
type zeroGradient;
}
frontside
{
type wedge;
}
backside
{
type wedge;
}
top
{
type noSlip;
}
centerline
{
type empty;
}
}
// ************************************************************************* //