Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective: To simulate the flow through the pipe using OpenFoam and to generate an automated blockMeshdict file using Matlab.
Given:
- Axisymmetric flow
- Laminar flow
- Incompressible flow
- Constant cross section area
Theory:
Consider water flowing through a pipe of diameter,D and length,L such that Reynolds number is below 2100.
Due to viscosity ,no slip condition prevails at the solid boundary. As a result, a zone near the solid boundary develops where velocity gradient exists.This narrow region is called the boundary layer. After a certain distance the boundary layer merges , this distance is called the entrance length. Inside the developing flow region the velocity changes and hence it is an accelerating flow. The velocity profile gradually takes the shape of parabola as flow moves towards developed region.
The shear stress in the pipe is developed at the boundary due to velocity gradient.In general the shear force that gets developed on a pipe due to viscous flow depends upon following factors
- Relative velocity of flow
-Mass density of fluid
-Viscousity of fluid
-Diameter/Length of Pipe
Now to simulate the pipe flow in openFoam we need to define the geomerty. To ease of our , we will tale the sector of the pipe or wedge and write the matlab code for the same. The sector is shown below.
Matlab code:
clear all
close all
clc
%% GIVEN PARAMETERS
rho=1000;
v=0.1246; %Velocity inlet (m/s)
mu=8.90e-4; %DYNAMIC VISCOUSITY (Pa-s)
Re=1000; %REYNOLDS NUMBER
D=0.015;
r=D/2;
theta= 4;
L_e=0.05*Re*D;
L=L_e+0.5;
dP=(32*mu*v)/(D^2);
%% vertices
v_0=[0 0 0];
v_1=[r*sind(theta/2) r*cosd(theta/2) 0];
v_2=[-r*sind(theta/2) r*cosd(theta/2) 0];
v_3=[0 0 L];
v_4=[r*sind(theta/2) r*cosd(theta/2) L];
v_5=[-r*sind(theta/2) r*cosd(theta/2) L];
%% BLOCKS
H_0='hex (0 3 4 1 0 3 5 2) (100 50 1) simpleGrading (1 1 1)';
%% edges
edge1=' arc 1 2 (0 0.0075 0)';
edge2=' arc 4 5 (0 0.0075 1.25)';
%% Faces
inlet=' (0 2 1 0)'
outlet=' (3 4 5 3)'
wedge1=' (0 1 4 3)'
wedge2=' (0 3 5 2)'
noSlipwall=' (1 2 5 4)'
axis=' (3 0 0 3)'
f1=fopen('blockMeshDict.txt','w')
text1='/*------------------------------*-C++-*----------------------------------- *\'
text2='| ========== | |'
text3='| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |'
text4='| \\ / O peration | Version: v2112 |'
text5='| \\ / A nd | Website: www.openfoam.com |'
text6='| \\/ M anipulation | |'
text15='\*------------------------------------------------------------------------*/'
text8= 'FoamFile'
text9='{'
text10= 'version 2.0;'
text11= 'format ascii; '
text12= 'class dictionary;'
text13= 'object blockMeshDict; '
text14='}'
text16='scale 1;'
text17='vertices'
text18='('
text19=');'
text20='blocks'
text21='('
text22=');'
text23='edges'
text24='('
text25=');'
text26='boundary'
text27='('
text28=' inletWall'
text29=' {'
text30=' type patch;'
text31=' faces'
text32=' ('
text33=' );'
text34=' }'
text35=' outletWall'
text36=' {'
text37=' type patch;'
text38=' faces'
text39=' ('
text40=' );'
text41=' }'
text42=' front'
text43=' {'
text44=' type wedge;'
text45=' faces'
text46=' ('
text47=' );'
text48=' }'
text49=' Back'
text50=' {'
text51=' type wedge;'
text52=' faces'
text53=' ('
text54=' );'
text55=' }'
text56=' noSlipwall'
text57=' {'
text58=' type wall;'
text59=' faces'
text60=' ('
text61=' );'
text62=' }'
text63=' axis'
text64=' {'
text65=' type empty;'
text66=' faces'
text67=' ('
text68=' );'
text69=' }'
text70=');'
text71='//****************************************************//'
fprintf(f1,'%s\n' ,text1);
fprintf(f1,'%s\n' ,text2);
fprintf(f1,'%s\n' ,text3);
fprintf(f1,'%s\n' ,text4);
fprintf(f1,'%s\n' ,text5);
fprintf(f1,'%s\n' ,text6);
fprintf(f1,'%s\n' ,text15);
fprintf(f1,'%s\n' ,text8);
fprintf(f1,'%s\n' ,text9);
fprintf(f1,'%s\n' ,text10);
fprintf(f1,'%s\n' ,text11);
fprintf(f1,'%s\n' ,text12);
fprintf(f1,'%s\n' ,text13);
fprintf(f1,'%s\n\n' ,text14);
fprintf(f1,'%s\n\n' ,text16);
fprintf(f1,'%s\n' ,text17);
fprintf(f1,'%s\n\n' ,text18);
fprintf(f1,'(%f %f %f) \n' ,v_0);
fprintf(f1,'(%f %f %f) \n' ,v_1);
fprintf(f1,'(%f %f %f) \n' ,v_2);
fprintf(f1,'(%f %f %f) \n' ,v_3);
fprintf(f1,'(%f %f %f) \n' ,v_4);
fprintf(f1,'(%f %f %f) \n\n' ,v_5);
fprintf(f1,'%s\n\n',text19);
fprintf(f1,'%s\n',text20);
fprintf(f1,'%s\n',text21);
fprintf(f1,'%s\n',H_0);
fprintf(f1,'%s\n\n',text22);
fprintf(f1,'%s\n' ,text23);
fprintf(f1,'%s\n' ,text24);
fprintf(f1, '%s\n %s\n' ,edge1,edge2 );
fprintf(f1,'%s\n\n' ,text25);
fprintf(f1,'%s\n' ,text26);
fprintf(f1,'%s\n' ,text27);
fprintf(f1,'%s\n' ,text28);
fprintf(f1,'%s\n' ,text29);
fprintf(f1,'%s\n' ,text30);
fprintf(f1,'%s\n' ,text31);
fprintf(f1,'%s\n' ,text32);
fprintf(f1,'%s\n',inlet);
fprintf(f1,'%s\n' ,text33);
fprintf(f1,'%s\n\n' ,text34);
fprintf(f1,'%s\n' ,text35);
fprintf(f1,'%s\n' ,text36);
fprintf(f1,'%s\n' ,text37);
fprintf(f1,'%s\n' ,text38);
fprintf(f1,'%s\n' ,text39);
fprintf(f1,'%s\n' ,outlet);
fprintf(f1,'%s\n' ,text40);
fprintf(f1,'%s\n\n' ,text41);
fprintf(f1,'%s\n',text42);
fprintf(f1,'%s\n',text43);
fprintf(f1,'%s\n',text44);
fprintf(f1,'%s\n',text45);
fprintf(f1,'%s\n',text46);
fprintf(f1,'%s\n',wedge1);
fprintf(f1,'%s\n',text47);
fprintf(f1,'%s\n\n',text48);
fprintf(f1,'%s\n',text49);
fprintf(f1,'%s\n',text50);
fprintf(f1,'%s\n',text51);
fprintf(f1,'%s\n',text52);
fprintf(f1,'%s\n',text53);
fprintf(f1,'%s\n',wedge2);
fprintf(f1,'%s\n',text54);
fprintf(f1,'%s\n\n',text55);
fprintf(f1,'%s\n',text56);
fprintf(f1,'%s\n',text57);
fprintf(f1,'%s\n',text58);
fprintf(f1,'%s\n',text59);
fprintf(f1,'%s\n',text60);
fprintf(f1,'%s\n',noSlipwall);
fprintf(f1,'%s\n',text61);
fprintf(f1,'%s\n\n',text62);
fprintf(f1,'%s\n',text63);
fprintf(f1,'%s\n',text64);
fprintf(f1,'%s\n',text65);
fprintf(f1,'%s\n',text66);
fprintf(f1,'%s\n',text67);
fprintf(f1,'%s\n',axis);
fprintf(f1,'%s\n',text68);
fprintf(f1,'%s\n\n',text69);
fprintf(f1,'%s\n',text70);
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n%s\n','(',');');
fprintf(f1,'%s\n',text71);
fclose(f1);
%% Calculation of Shear Stress
v_max=2*v;
Tau_max=(2*mu*v_max)/r;
%% Haigen-Poiseuille equation
dP=(32*mu*v)/(D^2);
blockMeshdict file
/*------------------------------*-C++-*----------------------------------- *\
| ========== | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2112 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
scale 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000262 0.007495 0.000000)
(-0.000262 0.007495 0.000000)
(0.000000 0.000000 1.250000)
(0.000262 0.007495 1.250000)
(-0.000262 0.007495 1.250000)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (100 50 1) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.0075 0)
arc 4 5 (0 0.0075 1.25)
);
boundary
(
inletWall
{
type patch;
faces
(
(0 2 1 0)
);
}
outletWall
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
Back
{
type wedge;
faces
(
(0 3 5 2)
);
}
noSlipwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type empty;
faces
(
(3 0 0 3)
);
}
);
mergePatchPairs
(
);
//****************************************************//
Transport Properties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-7;
// ************************************************************************* //Controldict file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 17;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Results
Plot of velocity at x= 0.2m
Plot of velocity at x=0.8m
Plot of Shear stress at x=0.8m
In above graph we can see that shear stress goes on increasing towards the wall.The shear stress is maximum at pipe wall and is 0 or minimum at center . The shear stress is maximum at boundary because velocity gradient is maximum at wall. As we have assumed newtonian fluid , the dynamic viscousity is constant and hence shear stress variation is linear, as expected.
Also as per theorotical result the velocity profile and shear stress variation obtained is parabolic in fully developed region and linear respectively.