Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim:Simulation of Flow through a pipe in OpenFoam.
Objectives:
- write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.
- Simulation of Flow through a pipe in OpenFoam using icoFoam solver.
- Post process the velocity profile and shear stress and compare the simulation results with analytical results obtained using the Hagen Poiseuille's equation.
Assumptions:
- The steady laminar and incompressible flow of a fluid is taken.
- Reynolds number for the laminar flow is 2100.
- The wedge shape pipe is considered for simulation because the flow is axis- symmetric.
- The equation governing the fluid flow is described by the Hagen Poiseuille's equation.
Hagen Poiseuille's equation :
The assumptions of the equations are the fluid is incompressible and Newtonian, the flow is laminar through a pipe of a constant cross section and its length is longer than the diameter of pipe.
Observations based on Hagen Poiseuille's equation :
- Diameter of pipe D=0.02 m
- Hydrodynamic entrance length =`0.06*R_e*D`=2.52 m
- Wedge angle `theta`=`4 ^0`
- Dynamic viscosity of water `mu=0.89e^-3 pascal-sec`
- Density of water `rho`=997 `(kg)/m^3`
- Kinematic Viscosity `nu`=`mu/rho=``(0.89e^-3)/997`=`8.9268e^-07 m^2/sec`
- Average velocity `v_(avg)`=`(mu*R_e)/(rho*D)=0.0937 m/sec`
- Velocity profile `u(r)=``2*v_(avg)*(1-r^2/R^2)`
- The velocity profile is maximum at the center of the pipe and becomes zero at the pipe wall due to no slip condition at the pipe wall.
- Maximum velocity `u_max=2*v_(avg)=0.1875 m/sec`
- Pressure drop for the fluid in laminar flow is `DeltaP=(32*mu*v_(avg)*L_c)/(D^2)`=0.0188 kPa
- where `L_c` is the total length of pipe we assume `L_c=L_e+0.3`=2.82m
- shear stress `tau =``(2*mu*u_max*r)/(R^2)`
- where r= Radial distance 0 to R , R = Radius of pipe ,
Wedge profile of the pipe used in blockMeshDict file :
Matlab code for generation of blockMeshDict script file:
% simulation of a flow in a pipe
clear all
close all
clc
%Defining the inputs
%For laminar flow
Re=2100;
% Diameter of pipe
D=0.02;
% Radius of pipe
R=D/2;
% Wedge angle
theta=4;
% Hydrodynamic entry length L
L=0.06*Re*D;
% Assumed characteristic length
L_c=L+0.3;
% Dynamic viscosity of water at 25*c
mu=0.89e-3;
% Density of water at 25*c
rho=997;
% Kinematic viscosity
nu=mu/rho;
% Calculating the average velocity
V_avg=(mu*Re)/(rho*D);
% Calculating the maximum velocity at the centerline of pipe
U_max=2*V_avg;
% Darcy friction factor
f=64/Re;
% Pressure drop for the laminar flow
delta_p=(32*mu*V_avg*L_c)/D^2;
% Finding Vertices
v0=[0 0 0];
v1=[0 R*cosd(0.5*theta) -R*sind(0.5*theta)];
v2=[0 R*cosd(0.5*theta) R*sind(0.5*theta)];
v3=[L_c 0 0];
v4=[L_c R*cosd(0.5*theta) -R*sind(0.5*theta)];
v5=[L_c R*cosd(0.5*theta) R*sind(0.5*theta)];
% Writing the script for blockMeshDict file
% Create a permission to write the file
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*----------------------------*-C++-*------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'======== | \n');
fprintf(f1,'%s','\\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / A nd | Version: 9');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / M anipulation |');
fprintf(f1,'\n');
fprintf(f1,'%s','\*------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'FoamFile\n');
fprintf(f1,'{ \n');
fprintf(f1,' format ascii;\n');
fprintf(f1,' class dictionary;\n');
fprintf(f1,' object blockMeshDict;\n');
fprintf(f1,'} \n');
fprintf(f1,'%s','//* * * * * * * * * * * * * * * * * * * * * * * * * * * * *//');
fprintf(f1,'\n\nconvertToMeters 1;\n\n');
% To form vertices
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v0(1),v0(2),v0(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v1(1),v1(2),v1(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v2(1),v2(2),v2(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v3(1),v3(2),v3(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v4(1),v4(2),v4(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v5(1),v5(2),v5(3));
fprintf(f1,')\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,'%s %s',' ','hex (0 1 2 0 3 4 5 3)');
fprintf(f1,blanks(1));
fprintf(f1,'%s','(50 1 300)');
fprintf(f1,'%s %s','','simpleGrading (0.4 1 1');
fprintf(f1,')\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','edges');
fprintf(f1,'\n');
fprintf(f1,'( \n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 1 2(');
fprintf(f1,'%f %f %f',0,R,0);
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 4 5(');
fprintf(f1,'%f %f %f',L_c,R,0);
fprintf(f1,')\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','boundary');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,'%s %s',' ','inlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','outlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(3 4 5 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','Front');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','Back');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 1 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','Top');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','Axis');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','mergePatchPairs');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','//* * * * * * * * * * * * * * * * * * * * * * * * * * * * *//');
fclose(f1);
Creating the geometry of the wedge:
BlockMesh Dict file:
Control Dict file :
Transport properties file:
Pressure File:
Velocity file:
Result: Generated geometry and mesh:
Velocity distribution result by applying the initial and boundary condition:
Velocity profile [Graphical representation]:
Pressure profile [Graphical representation]:
Shear stress
Shear stress can be calculated as under `tau=mu*(du)/(dr)`
where
`tau=`Shear stress in pascal
`mu`=Dynamic viscosity of water in Pascal-sec
`(du)/(dr)`=Velocity gradient `
In Newtonian fluid dynamic viscosity is constant thus shear stress is `prop` velocity gradient `(du)/(dr)`
Calculating Shear stress numerically:
Matlab code:
% Calculating shear stress
% Dynamic viscosity of water at 25*c
mu=0.89e-3;
%Diameter of pipe
D=0.02;
%Radius of pipe
R=D/2;
%Average velocity
V_avg=0.0937;
%Maximum velocity
U_max=2*V_avg;
% Radial distance
r=linspace(-0.01,0.01,100);
for i= 1:length(r)
Tau(i)=(2*mu*U_max*abs (r(i)))/(R^2);
end
plot(r,Tau,'linewidth',3)
grid on
xlabel('Radial distance in m')
ylabel('Shear stress in Pascal')
title('Shear stress for fully developed flow at x=2.52 m ')
legend('Shear stress of the profile')
Shear stress graph from Matlab :
Shear stress graph from OpenFoam:
The OpenFoam results are tabulated as under
Results Hagen Poiseuille equation Wedge angle `4^0`
CFL Number mean - 0.1020
CFL Number max - 0.3761
Execution Time - 387.99 sec
Clock Time - 390 sec
Max. Velocity 0.1875 m/sec 0.19 m/sec
Conclusion:
- The velocity profile is plotted at various distances from the inlet of the pipe .
- The velocity at a distance x=0.01m from the inlet.from the graph we can conclude that the velocity is gradually developing and suddenly become zero at the wall because of the no slip condition at wall.
- The velocity profile at a distance of x=0.4 , x=1m, x=2.52m, x=2.82m from the observation of graph we can conclude that the nature of graph is parabolic because the velocity profile moves to the hydrodynamically fully developed region.
- The shear stress plot is at distance of x= 0.01m, 2.2 m,2.52m. is plotted.
- From the above table shows that the maximum velocity of simulation is approximately equal to the HP equation velocity.