CFD simulation of flow through a pipe in OpenFoam using symmetry and wedge boundary conditions for icoFoam solver
Abstract
The given problem is described to demonstrate a flow simulation over a pipe to observe the boundary layer formation at different angles of theta. Symmetry boundary conditions are introduced and they are compared with the wedge boundary conditions. The velocity profile is captured at three different zones which are at start of entry length, at the end of entry length and after the entry length for the two different boundary conditions. Similarly, pressure drop is also observed at the fully developed region. The computational numerical data are validated against the analytical results to support the simulation algorithm, which efficiently predicted the physics of fluid flow along the entire geometrical domain.
Theory:
When a real fluid flows past a stationary solid boundary, a layer of fluid which comes in contact with the boundary surfaces adheres to it (on account of viscosity) and condition of no slip occurs. Thus, the layer of fluid which does not slips away undergoes retardation, this retarded layer further causes retardation in the vicinity of boundary layers in which the velocity of the fluid increases rapidly from zero at the boundary surface to maximum velocity at the main stream. The layer adjacent to the boundary is known as boundary layer. The formation is given below:
Geometry
Initial conditions
Length of the pipe= 1.5m
Density=1000 kg/m3
Diameter=0.01 m
Dynamic viscosity= 8.9*10-4
Mesh specification
No of cells in x direction = 20
No of cells in y direction = 1
No of cells in z direction = 300
Mesh grading in x direction=0.2
Mesh grading in y direction=1
Mesh grading in z direction=1
Vertices of geometry
Block vertices
(0 2 1 0 3 5 4 3)
Workflow diagram
Physics behind simulation
- The angle is varied at 10, 25, 45 degree to capture physics of flow from the entry region to the exit region. This allowed to visualize the boundary flow formation.
- The flow dynamics is captured at different zone to predict the numerical maximum velocity.
- Utilized both symmetry and wedge boundary conditions to predict the difference in boundary flow regions.
- The pressure gradient flow and shear stress at the fully developed region is numerically predicted to validate against the analytical solutions.
MATLAB CODE
BLOCKMESHDICT FILE GENERATION
clc
clear all
close all
% %------------simulating pipe flow-------------------------------------%
% % specify the condition for the flow as well as the dimension of domain
Re=2100;
rho=1000;
mu=8.9e-4;
d=0.01;
r=d/2;
l=5;
th=10;
area=pi/4*d^2;
l_entry=0.05*d*Re;
v_avg=Re*mu/(rho*d);
v_max=2*v_avg;
Q=v_avg*area;
pressure=32*mu*v_avg*l/d^2;
kinematic_pressure=pressure/rho;
% specifying the number of mesh cells in x,y,z direction
nx=20;
ny=1;
nz=300;
% specifying the grading
gx=0.2;
gy=1;
gz=1;
%specifying geometry of the wedges
w_0=[0 0 0];
w_1=[0 r*cosd(th/2) r*sind(th/2)];
w_2=[0 r*cosd(th/2) -r*sind(th/2)];
w_3=[5 0 0];
w_4=[5 r*cosd(th/2) r*sind(th/2)];
w_5=[5 r*cosd(th/2) -r*sind(th/2)];
arc_1=[0.0 r 0.0];
arc_2=[5 r 0.0];
%Using file parsing technique to create blockmeshdict file
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 | CodeDeveloper: SOUDIP HAZRA ';
line7='*-----------------------------------------------------------------*/';
line8='FoamFile ';
line9=' { ';
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%sn',line1);
fprintf(f1,'%sn',line2);
fprintf(f1,'%sn',line3);
fprintf(f1,'%sn',line4);
fprintf(f1,'%sn',line5);
fprintf(f1,'%sn',line6);
fprintf(f1,'%sn',line7);
fprintf(f1,'n');
fprintf(f1,'%sn',line8);
fprintf(f1,'%snn',line9);
fprintf(f1,'%16s %10s n','version','2.0; ');
fprintf(f1,'%15s %12s n','format' ,'ascii; ');
fprintf(f1,'%14s %11s n','class' ,'dictionary; ');
fprintf(f1,'%15s %11s n','object' ,'blockMeshDict;');
fprintf(f1,'n } nn');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //n');
fprintf(f1,' convertToMeters %d;n n',1);
% defining geometry%
fprintf(f1,'tvertices n');
fprintf(f1,'tt( n');
fprintf(f1,'ttt(%f %f %f)n',w_0);
fprintf(f1,'ttt(%f %f %f)n',w_1);
fprintf(f1,'ttt(%f %f %f)n',w_2);
fprintf(f1,'ttt(%f %f %f)n',w_3);
fprintf(f1,'ttt(%f %f %f)n',w_4);
fprintf(f1,'ttt(%f %f %f)n',w_5);
fprintf(f1,'tt);n n');
fprintf(f1,'ttblocks n');
fprintf(f1,'tt( n');
fprintf(f1,'ttt hex (0 2 1 0 3 4 5 3) (%d %d %d) simplegrading (%d %d %d)ntt); n n',nx,ny,nz,gx,gy,gz);
% defining sides of geometry%
% edge of pipe%
fprintf(f1,'ttedges n');
fprintf(f1,'tt( n');
fprintf(f1,'tttarc 1 2 (%f %f %f) n',arc_1);
fprintf(f1,'tttarc 4 5 (%f %f %f) n',arc_2);
fprintf(f1,'tt);n n');
% defining boundary conditions%
fprintf(f1,'ttboundaryntt(n');
% inlet section of pipe
fprintf(f1,'%16sn','inletWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type patch');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(0 1 2 0)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
% outlet section of pipe
fprintf(f1,'%16sn','outletWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type patch');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(3 5 4 3)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
%front side of pipe
fprintf(f1,'%16sn','frontWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type symmetry');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(0 3 4 1)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
%back side of pipe
fprintf(f1,'%16sn','backWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type symmetry');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(2 5 3 0)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
%top side of wall
fprintf(f1,'%16sn','topWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type wall');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(1 4 5 2)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
%bottom of pipe
fprintf(f1,'%16sn','bottomWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type empty');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(0 3 3 0)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
fprintf(f1,');n n');
%merge patch pair
fprintf(f1,'%20sn','mergePatchPairs');
fprintf(f1,'%6sn','(');
fprintf(f1,'%5sn',');');
fprintf(f1,'%sn','//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//');
fclose(f1)
control dict 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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 100;
deltaT 1e-2;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
pressure boundary 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
{
inletWall
{
type zeroGradient;
}
outletWall
{
type fixedValue;
value uniform 0.2661;
}
topWall
{
type zeroGradient;
}
bottomWall
{
type empty;
}
frontWall
{
type symmetry;
}
backWall
{
type symmetry;
}
}
// ************************************************************************* //
velocity boundary 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.1869 0 0);
boundaryField
{
inletWall
{
type fixedValue;
value uniform (0.1869 0 0);
}
outletWall
{
type zeroGradient;
}
topWall
{
type fixedValue;
value uniform (0 0 0)
}
bottomWall
{
type empty;
}
frontWall
{
type symmetry;
}
backWall
{
type symmetry;
}
}
// ************************************************************************* //
MATLAB code for validation between numerical and analytical solution
% numerical vs analytical solution at fully developed region
clc
clear all
T10=readtable('velocity10w.csv');
T25=readtable('velocity25w.csv');
T45=readtable('velocity45w.csv');
TS_10=readtable('velocity10symmetry.csv');
TS_25=readtable('velocity25symm.csv');
TS_45=readtable('velocity45symm.csv');
% r_10=T10{:,1}';% radial distance
r_10=linspace(0,0.005,101);
r_25=T25{:,1}';
r_45=T45{:,1}';
rs_10=TS_10{:,1}';% radial distance
rs_25=TS_25{:,1}';
rs_45=TS_45{:,1}';
numerical_velocity_10w=T10{:,2}';% numerical velocity data
numerical_velocity_25w=T25{:,2}';
numerical_velocity_45w=T45{:,2}';
numerical_velocity_10s=TS_10{:,2}';% numerical velocity data
numerical_velocity_25s=TS_25{:,2}';
numerical_velocity_45s=TS_45{:,2}';
%---------extracting the maximum numerical velocity----------------
max_velo_10w=max(numerical_velocity_10w);
max_velo_25w=max(numerical_velocity_25w);
max_velo_45w=max(numerical_velocity_45w);
max_velo_10s=max(numerical_velocity_10s);
max_velo_25s=max(numerical_velocity_25s);
max_velo_45s=max(numerical_velocity_45s);
%-------------------------------------------------------------------
R=0.005;% radius of pipe
v_max=0.3738;% anlytical max velocity
M10=readtable('pressure10w.csv');
M25=readtable('pressure25w.csv');
M45=readtable('pressure45w.csv');
MS_10=readtable('pressure10symmetry.csv');
MS_25=readtable('pressure25symm.csv');
MS_45=readtable('pressure45symm.csv');
pressure_10=M10{:,2}';
pressure_25=M25{:,2}';
pressure_45=M45{:,2}';
pressure_10s=MS_10{:,2}';
pressure_25s=MS_25{:,2}';
pressure_45s=MS_45{:,2}';
max_pressure_10w=max(pressure_10);
max_pressure_25w=max(pressure_25);
max_pressure_45w=max(pressure_45);
max_pressure_10s=max(pressure_10s);
max_pressure_25s=max(pressure_25s);
max_pressure_45s=max(pressure_45s);
min_pressure_10w=min(pressure_10);
min_pressure_25w=min(pressure_25);
min_pressure_45w=min(pressure_45);
min_pressure_10s=min(pressure_10s);
min_pressure_25s=min(pressure_25s);
min_pressure_45s=min(pressure_45s);
% %-----------wedge numerical pressure difference-----------
numerical_pressure_difference_10w=max_pressure_10w-min_pressure_10w;%numerical hagen poiseulli pressure
numerical_pressure_difference_25w=max_pressure_25w-min_pressure_25w;
numerical_pressure_difference_45w=max_pressure_45w-min_pressure_45w;
% %-----------symmetry numerical pressure difference-----------
numerical_pressure_difference_10s=max_pressure_10s-min_pressure_10s;%numerical hagen poiseulli pressure
numerical_pressure_difference_25s=max_pressure_25s-min_pressure_25s;
numerical_pressure_difference_45s=max_pressure_45s-min_pressure_45s;
for i=1:length(r_10)
analytical_velocity_10w=v_max*(1-(r_10/R).^2);%analytical velocity data
analytical_velocity_25w=v_max.*(1-(r_10./R).^2);
analytical_velocity_45w=v_max.*(1-(r_10./R).^2);
analytical_velocity_10s=v_max.*(1-(r_10./R).^2);%analytical velocity data
analytical_velocity_25s=v_max.*(1-(r_10./R).^2);
analytical_velocity_45s=v_max.*(1-(r_10./R).^2);
wedge_error_10=abs(analytical_velocity_10w-numerical_velocity_10w);
wedge_error_25=abs(analytical_velocity_25w-numerical_velocity_25w);
wedge_error_45=abs(analytical_velocity_45w-numerical_velocity_45w);
symmetry_error_10=abs(analytical_velocity_10s-numerical_velocity_10s);
symmetry_error_25=abs(analytical_velocity_25s-numerical_velocity_25s);
symmetry_error_45=abs(analytical_velocity_45s-numerical_velocity_45s);
end
figure(1)
suptitle('Analytical velocity vs Numerical velocity at fully developed region for wedge boundary condition');
subplot(3,1,1)
plot(r_10,analytical_velocity_10w,'-D','linewidth',1.5)
hold on
plot(r_10,numerical_velocity_10w,'-*','linewidth',1.5)
xlabel('Radial distance');
ylabel('velocity at 10 degree');
legend('analytical velocity','CFD numerical velocity')
subplot(3,1,2)
plot(r_25,analytical_velocity_25w,'-D','linewidth',1.5)
hold on
plot(r_25,numerical_velocity_25w,'-*','linewidth',1.5)
xlabel('Radial distance');
ylabel('velocity at 25 degree');
legend('analytical velocity','CFD numerical velocity')
subplot(3,1,3)
plot(r_45,analytical_velocity_45w,'-D','linewidth',1.5)
hold on
plot(r_45,numerical_velocity_45w,'-*','linewidth',1.5)
xlabel('Radial distance');
ylabel('velocity at 45 degree');
legend('analytical velocity','CFD numerical velocity')
figure (2)
suptitle('Analytical velocity vs Numerical velocity at fully developed region for symmetry boundary condition');
subplot(3,1,1)
plot(rs_10,analytical_velocity_10s,'-D','linewidth',1.5)
hold on
plot(rs_10,numerical_velocity_10s,'-*','linewidth',1.5)
xlabel('Radial distance');
ylabel('velocity at 10 degree');
legend('analytical velocity','CFD numerical velocity')
subplot(3,1,2)
plot(rs_25,analytical_velocity_25s,'-D','linewidth',1.5)
hold on
plot(rs_25,numerical_velocity_25s,'-*','linewidth',1.5)
xlabel('Radial distance');
ylabel('velocity at 25 degree');
legend('analytical velocity','CFD numerical velocity')
subplot(3,1,3)
plot(rs_45,analytical_velocity_45s,'-D','linewidth',1.5)
hold on
plot(rs_45,numerical_velocity_45s,'-*','linewidth',1.5)
xlabel('Radial distance');
ylabel('velocity at 45 degree');
legend('analytical velocity','CFD numerical velocity')
figure(3)
suptitle('Error plot for wedge boundary conditions')
subplot(3,1,1)
bar(wedge_error_10,'facecolor',[0.12 0.4 0.8]);
xlabel('velocity')
ylabel('error at 10 degree')
subplot(3,1,2)
bar(wedge_error_25,'facecolor',[0.12 0.4 0.8]);
xlabel('velocity')
ylabel('error at 25 degree')
subplot(3,1,3)
bar(wedge_error_45,'facecolor',[0.12 0.4 0.8]);
xlabel('velocity')
ylabel('error at 45 degree')
figure(4)
suptitle('Error plot for symmetry boundary conditions')
subplot(3,1,1)
bar(symmetry_error_10,'facecolor',[0.12 0.4 1]);
xlabel('velocity')
ylabel('error at 10 degree')
subplot(3,1,2)
bar(symmetry_error_25,'facecolor',[0.12 0.4 1]);
xlabel('velocity')
ylabel('error at 25 degree')
subplot(3,1,3)
bar(symmetry_error_45,'facecolor',[0.12 0.4 1]);
xlabel('velocity')
ylabel('error at 45 degree')
%--------------boundary layer formation for symmetry boundary condition---------
A=readtable('0.01s.csv');
B=readtable('1s.csv');
C=readtable('2.5s.csv');
D=readtable('4s.csv');
v_0=A{:,2}';
v_1=B{:,2}';
v_2=C{:,2}';
v_3=D{:,2}';
radial_distance=A{:,1}';
figure(5)
plot(radial_distance,v_0,'-D')
hold on
plot(radial_distance,v_1,'-D')
hold on
plot(radial_distance,v_2,'-D')
hold on
plot(radial_distance,v_3,'-D')
xlabel('radial distance')
ylabel('velocity profile')
legend('at 0.01m','at 1m','at 2.5m','at 4m')
title('bondary region formation at different velocity profile')
%-------------computing and validating shear stress-------------------
S=readtable('tau.csv');
numerical_shear_stress=S{:,2}';
for j=2:101
analytical_shear_stress(j)=0.00089.*(analytical_velocity_10s(j-1)-analytical_velocity_10s(j))/(radial_distance(j)-radial_distance(j-1));
end
figure(6)
plot(radial_distance,numerical_shear_stress,'linewidth',1.5)
hold on
plot(radial_distance,analytical_shear_stress,'linewidth',1.5)
xlabel('radial distance')
ylabel('tau')
legend('numerical shear stress','analytical shear stress')
Results and Analysis
1.for wedge boundary conditions
at `theta`=10 degree
geometrical mesh
velocity simulation and plots
at start of entry length
at end of entry length
at fully developed region
pressure simulation and plots
at `theta`=25 degree
geometrical mesh
velocity simulation and plots
at start of entry length
at end of entry length
at fully developed region
pressure simulation and plots
at `theta`=45 degree
geometrical mesh
velocity simulation and plots
at start of entry length
at end of entry length
at fully developed region
pressure simulation and plots
2.for symmetry boundary conditions
at `theta`=10 degree
geometrical mesh
velocity simulation and plots
at start of entry length
at end of entry length
at fully developed region
pressure simulation and plots
at `theta`=25 degree
geometrical mesh
velocity simulation and plots
at start of entry length
at end of entry length
at fully developed region
pressure simulation and plots
at `theta`=45 degree
geometrical mesh
velocity simulation and plots
at start of entry length
at end of entry length
at fully developed region
pressure simulation and plots
from the given simulation plot at different wedge angle for both synmmetry and wedge boundary condition the results for maximum velocity correlates with the analytical velocity
comparison between numerical velocity and analytical velocity
wege condition
symmetry condition
error plot
wege condition
symmetry condition
The error margin for both the case is seen to get reduce down which predicts the accuracy of the CFD model.
Boundary layer formation for symmetry conditions
shear stress plot
numerical shears stress vs analytical shear stress
Comparison between wedge and symmetry boundary conditions
Conclusions:
- From courant number point of view, it is seen that symmetry boundary conditions gives more stable and lower values than that of wedge boundary conditions.
- For wedge angle of 10, 25 degree the error in pressure difference and velocity is lower for both symmetry and wedge boundary conditions.
- From computational time point of view the wedge boundary condition took less time than symmetry boundary condition.
- The numerical computation is in good agreement with the analytical solutions.
- From all the above points, it is concluded that symmetry boundary conditions is the best for simulating and capturing the physics of flow through a pipe.