SIMULATING LAMINAR FLOW THROUGH A PIPE USING ICOFOAM SOLVER AND WEDGE BOUNDARY CONDITION IN OPENFOAM (PART 1)
Aim :
- To calculate the flow variables
- To generate the blockMeshDict file using wedge boundary condition.
- To create a blockMeshDict file using Matlab such that the mesh is automatically generated for any wedge angle
- To simulate the pipe flow in openfoam using icoFoam solver.
- postprocess the velocity profile and shear stress and compare this simulation results with analytical results obtained using the Hagen-poiseuille equation.
Entrance length :
The entrance length is the distance a flow travels after entering a pipe before the flow becomes fully developed. Entrance length refers to the length of the 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 change with increased distance along the pipe.
Hydrodynamic entrance length :
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. This region is characterized by a non-uniform flow. The fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile. For the conservation of mass to hold true, the velocity of layers of the fluid in the center of the pipe increases to compensate for the reduced velocities of the layers of fluid near the pipe surface. This develops a velocity gradient across the cross section of the pipe.
The length of the hydrodynamic entry region along the pipe is called the hydrodynamic entry length. It is a function of Reynolds number of the flow. In case of laminar flow, this length is given by:
Where, Re is the Reynold's number and D is the diameter of the pipe.
In this experimental analysis a simple Hagen-poisuelli flow is simulated for Re=2100 in circular pipe using water. To reduce the calculation a wedge of pipe cross section was taken instead of using whole circular cross section
Data :
- Working fluid = water
- Density = 997 Kg/m^3
- Dynamic viscosity = 8.9*10^-4 Kg/ms
- Renolyds number = 2100
- Wedge angle = 5 degree
- Pipe diameter = 0.01 m
- pipe length from calculations = 1.56 m
- Number of cells along X axis = 100
- Number of cells along radially outward direction =20
Hagen-poiseuille's law :
In non ideal fluid dynamics this law gives pressure drop in an incompressible and newtonian fluid ina laminar flow through a long pipe of uniform cross section.
Average velocity
Pressure difference
Kinematic pressure difference
shear stress
Sample geometry :
Matlab code for generation of blockMesh dict file for any wedge angle and computing the analytical data :
close all
clear all
clc
%The Renoylds number
Re=2100;
%Diameter of the pipe in meters
D=0.01;
%The entry length of the pipe
Le=0.06*Re*D;
%Total length of the pipe
L=Le+0.3;
%Radius of the pipe
R=D/2;
%Angle of the wedge in degrees
Theta=5;
%Dynamic viscosity
mu=0.89e-3;
%Density in kg/m^3
rho=997;
% Haigen poisuelle flow equations
%Avg velocity
Vavg=(Re*mu)/(rho*D);
%Maximum velocity
Vmax=2*Vavg;
%Piezo metric pressure drop
Delp=32*mu*Vavg*L/D^2;
%Kinematic pressure
Kinematic_p=Delp/rho;
%creating set of Diameters for plotting the shear stress
Rvar=linspace(-R,R,1000);
%Calculating shear stress
for i=1:length(Rvar)
Tau(i)=2*mu*(Rvar(i))*Vmax/R^2;
v(i)=Vmax*(R^2-Rvar(i)^2)/R^2;
end
%Maximum shear stress value
Tau_max=max(Tau)
%plotting
figure(1)
plot(Tau,Rvar)
hold on
axis([-0.5 0.5 -R R])
xL = xlim;
yL = ylim;
line(xL,[0 0],'linestyle','--','color','r','linewidth',1); %x-axis
ylabel('Radius')
xlabel('Shear stress')
title('Analytical variation of shear stress along the radial distance')
grid on
legend('Shear stress','Axis of the pipe')
figure(2)
plot(v,Rvar)
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1); %x-axis
grid on
xlabel('Velocity')
ylabel('Radius')
legend('Velocity','Axis of the pipe')
title('Analytical variation of velocity along the radial distance')
%creating the blockmesh file
header1='/*---------------------------*- C++ -*--------------------------------*/'
header2=' ========== | '
header3=' \ / F ield | OpenFOAM: The open source CFD Toolbox '
header4=' \ / O perration | website: https://openfoam.org '
header5=' \ / A nd | version : 8 '
header6=' \ / M anipulation | '
header7='*------------------------------------------------------------------*/'
f1=fopen('blockMeshDict.c','w')
fprintf(f1,'%sn',header1)
fprintf(f1,'%sn',header2)
fprintf(f1,'%sn',header3)
fprintf(f1,'%sn',header4)
fprintf(f1,'%sn',header5)
fprintf(f1,'%sn',header6)
fprintf(f1,'%sn',header7)
fprintf(f1,'FoamFilen{n')
fprintf(f1,'%12s t 2.0;n','version')
fprintf(f1,'%12s t ascii;n','format')
fprintf(f1,'%12s t dictionary;n','class')
fprintf(f1,'%12s t blockMeshDict;n','object')
fprintf(f1,'}n')
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1,'%snn',header8)
fprintf(f1,'convertToMeters 1;nn')
fprintf(f1,'verticesn')
fprintf(f1,'(n')
fprintf(f1,'t (%d %d %d)n',0,0,0)
fprintf(f1,'t (%d %d %d)n',0,R*cosd(Theta/2),R*sind(Theta/2))
fprintf(f1,'t (%d %d %d)n',0,R*cosd(Theta/2),-R*sind(Theta/2))
fprintf(f1,'t (%d %d %d)n',L,0,0)
fprintf(f1,'t (%d %d %d)n',L,R*cosd(Theta/2),R*sind(Theta/2))
fprintf(f1,'t (%d %d %d)n',L,R*cosd(Theta/2),-R*sind(Theta/2))
fprintf(f1,');nn')
fprintf(f1,'blocksn')
fprintf(f1,'(n')
fprintf(f1,'t hex (0 2 1 0 3 5 4 3) (%d %d %d)',20,1,100)
fprintf(f1, 'simpleGrading (0.2 1 1)nn')
fprintf(f1,');nn')
fprintf(f1,'edgesn')
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,');nn')
fprintf(f1,'boundaryn')
fprintf(f1,'(n')
fprintf(f1,'t inlet n')
fprintf(f1,'t {n')
fprintf(f1,'tt type patch;n')
fprintf(f1,'tt facesn')
fprintf(f1,'tt (n')
fprintf(f1,'ttt (0 1 2 0)n')
fprintf(f1,'tt );n')
fprintf(f1,'t }nn')
fprintf(f1,'t outlet n')
fprintf(f1,'t {n')
fprintf(f1,'tt type patch;n')
fprintf(f1,'tt facesn')
fprintf(f1,'tt (n')
fprintf(f1,'ttt (3 5 4 3)n')
fprintf(f1,'tt );n')
fprintf(f1,'t }nn')
fprintf(f1,'t top n')
fprintf(f1,'t {n')
fprintf(f1,'tt type wall;n')
fprintf(f1,'tt facesn')
fprintf(f1,'tt (n')
fprintf(f1,'ttt (1 4 5 2)n')
fprintf(f1,'tt );n')
fprintf(f1,'t }nn')
fprintf(f1,'t front n')
fprintf(f1,'t {n')
fprintf(f1,'tt type wedge;n')
fprintf(f1,'tt facesn')
fprintf(f1,'tt (n')
fprintf(f1,'ttt (0 3 4 1)n')
fprintf(f1,'tt );n')
fprintf(f1,'t }nn')
fprintf(f1,'t back n')
fprintf(f1,'t {n')
fprintf(f1,'tt type wedge;n')
fprintf(f1,'tt facesn')
fprintf(f1,'tt (n')
fprintf(f1,'ttt (0 2 5 3)n')
fprintf(f1,'tt );n')
fprintf(f1,'t }nn')
fprintf(f1,'t axisn')
fprintf(f1,'t {n')
fprintf(f1,'tt type empty;n')
fprintf(f1,'tt facesn')
fprintf(f1,'tt (n')
fprintf(f1,'ttt (0 3 3 0)n')
fprintf(f1,'tt );n')
fprintf(f1,'t }nn')
fprintf(f1,');n')
fprintf(f1,'mergePatchPairsn')
fprintf(f1,'(n')
fprintf(f1,');n')
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1,'%sn',header9)
fclose(f1)BlockMeshDict file :
/*---------------------------*- C++ -*--------------------------------*/
========== |
\ / F ield | OpenFOAM: The open source CFD Toolbox
\ / O perration | website: https://openfoam.org
\ / A nd | version : 8
\ / M anipulation |
*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 4.995241e-03 2.180969e-04)
(0 4.995241e-03 -2.180969e-04)
(1.560000e+00 0 0)
(1.560000e+00 4.995241e-03 2.180969e-04)
(1.560000e+00 4.995241e-03 -2.180969e-04)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (20 1 100)simpleGrading (0.2 1 1)
);
edges
(
arc 1 2 (0 5.000000e-03 0)
arc 4 5 (1.560000e+00 5.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)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
controlDict file :
/*--------------------------------*- 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 10;
deltaT 0.05;
writeControl timeStep;
writeInterval 10;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Velocity boundary condition :
/*--------------------------------*- 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.1875 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Pressure boundary condition :
/*--------------------------------*- 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;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //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.9268e-7;
// ************************************************************************* //
Results :
Mesh
Inlet velocity :
Outlet velocity :
Pressure variation
Inlet
Outlet
Various velocity profiles for 5 degree wedge angle
at x= 0.01
at x= 0.05
at x=0.1
at x=0.5
at x=1
at x=1.56
Pressure variation along the length :
post processing
- After opening paraview click on apply to see the geometry
- press on play button to post process to the end time
- click on filters and select plot over line and plot velocity magnitude and pressure wrt to radial distance
- Export the velocity data in .csv format
- For calculating the shear stress import the .csv data in matlab
- select the suitable column to plot shear stress
- plot the result for shear stress
Matlab code for calculating shear stress by importing . csv files from openfoam :
clear all
close all
clc
%Dynamic viscosity
mu=0.89e-3;
%The Renoylds number
Re=2100;
%Diameter of the pipe in meters
D=0.01;
%The entry length of the pipe
Le=0.06*Re*D;
%Total length of the pipe
L=Le+0.3;
%Radius of the pipe
R=D/2;
%For 5 degree wedge angle
%At location 1.5 meter
data1=csvread('velocityatx=1.5.csv',1,0);
u1=data1(:,2);
y1=data1(:,1);
for i=1:length(y1)
y1_n(i)=-y1(length(y1)+1-i);
y1_neg=y1_n';
u1_n(i)=u1(length(u1)+1-i);
u1_neg=u1_n';
end
y1_complete=[y1_neg;y1(2:end)];
u1_complete=[u1_neg;u1(2:end)];
plot(u1_complete,y1_complete,'r','linewidth',2)
u1_maximum=max(u1_complete)
%At location 1 meter
data2=csvread('velocityatx=1.csv',1,0);
u2=data2(:,2);
y2=data2(:,1);
for i=1:length(y2)
y2_n(i)=-y2(length(y2)+1-i);
y2_neg=y2_n';
u2_n(i)=u2(length(u2)+1-i);
u2_neg=u2_n';
end
y2_complete=[y2_neg;y2(2:end)];
u2_complete=[u2_neg;u2(2:end)];
u2_maximum=max(u2_complete)
hold on
plot(u2_complete,y2_complete,'g','linewidth',2)
%At location 0.5 meter
data3=csvread('velocityatx=0.5.csv',1,0);
u3=data3(:,2);
y3=data3(:,1);
for i=1:length(y3)
y3_n(i)=-y3(length(y3)+1-i);
y3_neg=y3_n';
u3_n(i)=u3(length(u3)+1-i);
u3_neg=u3_n';
end
y3_complete=[y3_neg;y3(2:end)];
u3_complete=[u3_neg;u3(2:end)];
u3_maximum=max(u3_complete)
hold on
plot(u3_complete,y3_complete,'b','linewidth',2)
%At location 0.1 meter
data4=csvread('velocityatx=0.1.csv',1,0);
u4=data4(:,2);
y4=data4(:,1);
for i=1:length(y4)
y4_n(i)=-y4(length(y4)+1-i);
y4_neg=y4_n';
u4_n(i)=u4(length(u4)+1-i);
u4_neg=u4_n';
end
y4_complete=[y4_neg;y4(2:end)];
u4_complete=[u4_neg;u4(2:end)];
u4_maximum=max(u4_complete)
hold on
plot(u4_complete,y4_complete,'y','linewidth',2)
%At location 0.05 meter
data5=csvread('velocityatx=0.05.csv',1,0);
u5=data5(:,2);
y5=data5(:,1);
for i=1:length(y5)
y5_n(i)=-y5(length(y5)+1-i);
y5_neg=y5_n';
u5_n(i)=u5(length(u5)+1-i);
u5_neg=u5_n';
end
y5_complete=[y5_neg;y5(2:end)];
u5_complete=[u5_neg;u5(2:end)];
u5_maximum=max(u5_complete)
hold on
plot(u5_complete,y5_complete,'k','linewidth',2)
%At location 0.01 meter
data6=csvread('velocityatx=0.01.csv',1,0);
u6=data6(:,2);
y6=data6(:,1);
for i=1:length(y6)
y6_n(i)=-y6(length(y6)+1-i);
y6_neg=y6_n';
u6_n(i)=u6(length(u6)+1-i);
u6_neg=u6_n';
end
y6_complete=[y6_neg;y6(2:end)];
u6_complete=[u6_neg;u6(2:end)];
u6_maximum=max(u6_complete)
hold on
plot(u6_complete,y6_complete,'c','linewidth',2)
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1); %x-axis
xlabel('velocity')
ylabel('Radial distance')
legend('At length=1.5','At length=1','At length=0.5','At length=0.1','At length=0.05','At length=0.01','axis')
title('velocity variation along the radial distance for 5 degree wedge angle')
%For 2 degree wedge angle and at the end of the pipe
data7=csvread('2degreex=1.55.csv',1,0);
u7=data7(:,2);
y7=data7(:,1);
for i=1:length(y7)
y7_n(i)=-y7(length(y7)+1-i);
y7_neg=y7_n';
u7_n(i)=u7(length(u7)+1-i);
u7_neg=u7_n';
end
y7_complete=[y7_neg;y7(2:end)];
u7_complete=[u7_neg;u7(2:end)];
u7_maximum=max(u7_complete)
figure(2)
plot(u7_complete,y7_complete,'c','linewidth',2.5)
%For 3 degree wedge angle and at the end of the pipe
data8=csvread('3degreex=1.55.csv',1,0);
u8=data8(:,2);
y8=data8(:,1);
for i=1:length(y8)
y8_n(i)=-y8(length(y8)+1-i);
y8_neg=y8_n';
u8_n(i)=u8(length(u8)+1-i);
u8_neg=u8_n';
end
y8_complete=[y8_neg;y8(2:end)];
u8_complete=[u8_neg;u8(2:end)];
u8_maximum=max(u8_complete)
hold on
plot(u8_complete,y8_complete,'g','linewidth',2)
%For 4 degree wedge angle and at the end of the pipe
data9=csvread('4degreex=1.55.csv',1,0);
u9=data9(:,2);
y9=data9(:,1);
for i=1:length(y9)
y9_n(i)=-y9(length(y9)+1-i);
y9_neg=y9_n';
u9_n(i)=u9(length(u9)+1-i);
u9_neg=u9_n';
end
y9_complete=[y9_neg;y9(2:end)];
u9_complete=[u9_neg;u9(2:end)];
u9_maximum=max(u9_complete)
u9_avg=u9_maximum/2;
hold on
plot(u9_complete,y9_complete,'k','linewidth',1)
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1); %x-axis
xlabel('velocity')
ylabel('Radial distance')
legend('2 degree wedge angle','3 degree wedge angle','4 degree wedge angle','axis')
title('velocity variation at length =1.55 meter along the radial distance')
%calculation of shear stress for 4 degree wedge angle
for i=1:length(y9_complete)
Tau(i)=2*mu*(y9_complete(i))*u9_maximum/R^2;
end
Tau_max=max(Tau)
figure(3)
plot(Tau,y9_complete,'k','linewidth',2)
xlabel('shear stress')
ylabel('Radial distance')
axis([-0.15 0.15 -R R])
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1); %x-axis
legend('shear stress','axis')
title('shear stress variation along the radius for 4 degree wedge angle')
Comparison :
| Variables | Analytical | 2 degree wedge | 3 degree wedge | 4 degree wedge |
| Maximum velocity in m/sec | 0.3749 | 0.3714 | 0.3714 | 0.3714 |
| Wall shear stress in N/m^2 | 0.1335 | 0.1322 | 0.1322 | 0.1321 |
| pressure drop in pa | 83.2873 |
0.105*997 Kinematic pressure *density =104.68 |
104.68 | 104.68 |
Conclusion :
- The wedge angle should not be greater than 5 degree as per openfoam
- The simulated results does not vary much under 5 degree