Week 11 - Simulation of Flow through a pipe in OpenFoam
In nonideal flueid dynamics,the higen_poiseuille equation, also known as Hagen-poiseuille law, iis a pysical law that gives the pressuree drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section
EQUATION:
In standard fluid-kinetics notation.

where:
Δp is the pressure dsifference between the two ends
L is the length of pipe
μ is the dynamic viscosity
Q is the volumetric flow rate
R is the pipe radius
A is the cross section of pipe
The equation does not hold close to p[ipe entrance.
Given parameters.
- Wedge angle less than 5 degrees.
- Flow is axi symetrical.
- Fluid is incompressible and the flow is laminar.
Expected Velocity Profile.

Assuming the following:
- The distance the fluid has to travel to become fully developed is given by the Entry Length formula.
- Given Reynolds number is 2100.This implies flow is laminar.
- Hence Entry length formula for Laminar is given by:L/d=0.06.Re
- where L is effective length; d is dia of pipe and Re is reynolds number.
- Assuming the dia of pipe as 10mm;substituting the values,we get L=1.2m;
Steps to follow
1. First, the geometry is created in the blockmesh Dict file, this is done by creating the file through the matlab code.
2. The mesh is generated according to requirement.
3. The velocity and preassure bc are applied
4. The controlDict file is edited for Simulation run time and timestep;
5. After the geometry is input, the simulation can be started by calling 'icoFoam' command in the terminal.
6. 'paraFoam' command helps to open paraView
7. The shear stress diagram is obtained by post-processing the velocity results obtained from Openfoam.
Matlab programm for automatically generate blockmeshDict file for a wedge angle provided by user.
1. Main Code
clear all
close all
clc
% problem initialization
Re=2100;
rho=1000;
mu=0.89e-3;
nu=mu/rho;
D=0.02;
A_c=(pi/4)*D^2;
p=pi*D;
D_h=4*A_c/p;
r_h=D_h/2;
v_avg=(Re*mu)/(rho*D_h);
v_max=2*v_avg;
tau_max=2*mu*v_max/r_h;
if Re<2300
L_e=0.06*Re*D_h;
elseif Re==2300
L_e=115*D_h;
elseif Re>4000
L_e=1.359*D_h*Re^(1/4);
end
L=L_e+0.18; % length of pipe
%% Pressure loss
f=65/Re;
% Haigen Poiseuille Eqn
P=(32*mu*(L-L_e)*v_avg)/(D^2);
g=9.81;
h_loss=P/(rho*g);
P_kinematicPress=P/rho;
theta=3;
r=r_h;
%% Print out
system(r,L,theta)
2. system()
function system(r,L,theta)
%% creating blockmesh dict file
%%opening part
h1='/=------------------------*-C++ -*------------------------*';
h2=' ======== |';
h3=' \ / F ield | opennFoam: The open source CFD toolbox';
h4=' \ / O peration | website: http://openfoam.org';
h5=' \ / A nd | Version: 9';
h6=' \/ M anipulation |';
h7='--------------------------------------*/';
h8='FoamFile';
h9='{';
h10=' format ascii;';
h11=' class dictionary;';
h12=' object blockMeshDict;';
h13='}';
h14='// *********************************//';
f1=fopen('blockMeshDict.txt','wt');
fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'%s\n',h12);
fprintf(f1,'%s\n',h13);
fprintf(f1,'%s\n',h14);
%%vertices
v0=[0 0 0];
v1=[0 r*cosd(theta/2) r*sind(theta/2)];
v2=[0 r*cosd(theta/2) -r*sind(theta/2)];
v3=[L 0 0];
v4=[L r*cosd(theta/2) r*sind(theta/2)];
v5=[L r*cosd(theta/2) -r*sind(theta/2)];
mid1=[0 r 0];
mid2=[L r 0];
%%write vertices
h15='convertToMeters 1;';
fprintf(f1,'%s\n',h15);
h16='vertices';
fprintf(f1,'%s\n',h16);
h17='(';
fprintf(f1,'%s\n',h17);
fprintf(f1,'(%d %d %d)\n',v0);
fprintf(f1,'(%d %d %d)\n',v2);
fprintf(f1,'(%d %d %d)\n',v3);
fprintf(f1,'(%d %d %d)\n',v4);
fprintf(f1,'(%d %d %d)\n',v5);
h18=');';
fprintf(f1,'%s\n',h18);
%% write blocks
h19='blocks';
h20='(';
h21=' hex (0 3 4 1 0 3 5 2) (200 50 1) simpleGrading (1 0.2 1)';
h22=');';
fprintf(f1,'%s\n',h19);
fprintf(f1,'%s\n',h20);
fprintf(f1,'%s\n',h21);
fprintf(f1,'%s\n',h22);
%% write edges
h23='edges';
h24='(';
h25=' arc 2 1';
h26=' arc 5 4';
h27=');';
fprintf(f1,'%s\n',h23);
fprintf(f1,'%s\n',h24);
fprintf(f1,'%s (%d %d %d)\n',h25,0, r, 0);
fprintf(f1,'%s (%d %d %d)\n',h26,L, r, 0);
fprintf(f1,'%s\n',h27);
%% write boundary
inlet= [0 1 2 0];
outlet=[3 5 4 3];
surface=[2 1 4 5];
axis=[0 3 3 0];
front=[1 0 3 4];
back=[0 2 5 3];
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
%inlet patch
fprintf(f1,'%s\n','inlet');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n','type patch');
fprintf(f1,'%s\n','faces');
fprintf(f1,'%s\n','(');
fprintf(f1,'(%d %d %d)\n',inlet);
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','}');
%outlet patch
fprintf(f1,'%s\n','outlet');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n','type patch');
fprintf(f1,'%s\n','faces');
fprintf(f1,'%s\n','(');
fprintf(f1,'(%d %d %d)\n',outlet);
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','}');
%surface
fprintf(f1,'%s\n','surface');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n','type wall');
fprintf(f1,'%s\n','faces');
fprintf(f1,'%s\n','(');
fprintf(f1,'(%d %d %d)\n',surface);
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','}');
%axis
fprintf(f1,'%s\n','axis');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n','type empty');
fprintf(f1,'%s\n','faces');
fprintf(f1,'%s\n','(');
fprintf(f1,'(%d %d %d)\n',axis);
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','}');
%front
fprintf(f1,'%s\n','front');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n','type wedge');
fprintf(f1,'%s\n','faces');
fprintf(f1,'%s\n','(');
fprintf(f1,'(%d %d %d)\n',front);
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','}');
%back
fprintf(f1,'%s\n','back');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n','type wedge');
fprintf(f1,'%s\n','faces');
fprintf(f1,'%s\n','(');
fprintf(f1,'(%d %d %d)\n',back);
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','}');
fprintf(f1,'%s\n',')');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',')');
fprintf(f1,'%s\n','// ******************************************* //');
end
OpenFoam Files:
1. System Folder:
a. blockMeshDict:
/*--------------------------------*- 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.996573e-03 -2.617695e-04)
(0 9.996573e-03 2.617695e-04)
(2.7 0 0)
(2.7 9.996573e-03 -2.617695e-04)
(2.7 9.996573e-03 2.617695e-04)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 1e-02 0)
arc 4 5 (2.7 1e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
surface
{
type wall;
faces
(
(2 5 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(0 3 5 2)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
controlDict:
/*--------------------------------*- 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 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 4;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
velocity (U):
/*--------------------------------*- 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0935 0 0);
}
outlet
{
type zeroGradient;
}
surface
{
type noSlip;
//type fixedValue;
//value uniform (0 0 0);
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
pressure (P):
/*--------------------------------*- 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
surface
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
Results and graph:
For wedge angle 3degree
1.Velocity line plot
velocity line plot across the radial distance X=2.5m

velocity line plot across the radial distance X=0.05m

maximum shear stress plot:
