Simulation of Flow through a pipe in OpenFoam
Your Answers
problem specification:
let us assume a cylindrical pipe of 0.01m diameter and length 2m
Fluid flowing with a velocity of 15m/sec with a kinematic viscosity of 1e-4 sqm/sec .
d = 0.01m.
nu = 1e-4 sqm/sec.
v = 15m/sec.
l = 2m.
rho = 1000 kg/cubicm.
mu = 10 Ns/sqm.
outlet is open to atmospheric pressure.
The Renolds number will be 1500
Hagen poiseuille's flow
In nonideal fluid dynamics the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross-section. It can be successfully applied to airflow in lung alveoli, or the flow through a drinking straw or through a hypodermic needle.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of the constant circular cross-section that is substantially longer than its diameter, and there is no acceleration of fluid in the pipe. For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent, leading to larger pressure drops than calculated by the Hagen–Poiseuille equation.
The specific calculation of the Reynolds number, and the values where laminar flow occurs, will depend on the geometry of the flow system and flow pattern. The common example flows through a pipe, where the Reynolds number is defined as
velocity distribution in the radial direction for fully developed laminar flow in a pipe of constant diameter:
shear stress distribution in the radial direction for fully developed laminar flow in a pipe of constant diameter:
which is linear, zero at a centre and maximum at walls as shown in the figure above.
The fully developed flow starts after the entrance length which can be calculated using the equation.
To get correct velocity and shear stress profiles we should measure the profiles after this entrance length.
1. Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes:
clear all
close all
clc
%% MATLAB program for auto-generation of blockmesh file in OpenFoam %%
d = 0.01; % diameter
L_t = 2; %total length
r = d/2; % radius
theta_actual = 4;
theta = 2; %degree(input according to the need)
%% Calculating the vertices of the wedge
v0 = [ 0 0 0 ];
v1 = [ L_t 0 0 ];
v2 = [ 0 -r*sind(theta) r*cosd(theta)];
v3 = [ L_t -r*sind(theta) r*cosd(theta)];
v4 = [ L_t r*sind(theta) r*cosd(theta)];
v5 = [ 0 r*sind(theta) r*cosd(theta)];
%% String format specifying for Blockmesh file %%
head1 = '/*--------------------------------*- C++ -*----------------------------------*';
head2 = '| ========= | |';
head3 = '| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
head4 = '| \ / O peration | Version: 7 |';
head5 = '| \ / A nd | Website: https://openfoam.org |';
head6 = '| \/ M anipulation | |';
head7 = '*---------------------------------------------------------------------------*/';
b4 = blanks(4);
b8 = blanks(8);
b12 = blanks(12);
b16 = blanks(16);
dict1 = ' version 2.0;';
dict2 = ' format ascii;';
dict3 = ' class dictionary;';
dict4 = ' object blockMeshDict;';
dict5 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
conv = 'convertToMeters 1;';
vertex1 = 'vertices';
blocks1 = 'hex (0 1 1 0 2 3 4 5) (100 1 5) simpleGrading (1 1 0.2)';
footer = '// ************************************************************************* //';
%% Creting blockMesh dict for OepnFOAM %%
f1 = fopen('blockMeshdict.txt','w')
% Printing the headder %%
fprintf(f1,'%sn' ,head1);
fprintf(f1,'%sn' ,head2);
fprintf(f1,'%sn' ,head3);
fprintf(f1,'%sn' ,head4);
fprintf(f1,'%sn' ,head5);
fprintf(f1,'%sn' ,head6);
fprintf(f1,'%sn' ,head7);
fprintf(f1,'FoamFilen{n' );
fprintf(f1,'%s%sn', b4,dict1);
fprintf(f1,'%s%sn', b4,dict2);
fprintf(f1,'%s%sn', b4,dict3);
fprintf(f1,'%s%sn', b4,dict4);
fprintf(f1,'}n%s%snnn', b4,dict5);
fprintf(f1,'%sn', conv);
%% Printing the vertices %%
fprintf(f1,'verticesn(n');
fprintf(f1,'%s(%f %d %d)n', b4,v0(1),v0(2),v0(3));
fprintf(f1,'%s(%f %d %d)n', b4,v1(1),v1(2),v1(3));
fprintf(f1,'%s(%f %f %f)n', b4,v2(1),v2(2),v2(3));
fprintf(f1,'%s(%f %f %f)n', b4,v3(1),v3(2),v3(3));
fprintf(f1,'%s(%f %f %f)n', b4,v4(1),v4(2),v4(3));
fprintf(f1,'%s(%f %f %f)n', b4,v5(1),v5(2),v5(3));
fprintf(f1,');n');
%% Printing the blocks %%
fprintf(f1,'nblocksn(n');
fprintf(f1,'%s%sn', b4,blocks1);
fprintf(f1,');n');
%% Printing edges %%
fprintf(f1,'nedgesn(n');
fprintf(f1,'%sarc 2 5 (0 0 %0.3f)n', b4,r);
fprintf(f1,'%sarc 3 4 (%0.2f 0 %0.3f)n', b4,L_t,r);
fprintf(f1,');n');
%% Printing boundary %%
fprintf(f1,'nboundaryn(n');
fprintf(f1,'%sinletn%s{n', b4,b4);
fprintf(f1,'%stype patch;n', b8);
fprintf(f1,'%sfacesn%s(n', b8,b8);
fprintf(f1,'%s(0 2 5 0)n%s);n%s}n', b12,b8,b4);
fprintf(f1,'%saxisn%s{n', b4,b4);
fprintf(f1,'%stype empty;n', b8);
fprintf(f1,'%sfacesn%s(n', b8,b8);
fprintf(f1,'%s(0 1 1 0)n%s);n%s}n', b12,b8,b4);
fprintf(f1,'%swalln%s{n', b4,b4);
fprintf(f1,'%stype wall;n', b8);
fprintf(f1,'%sfacesn%s(n', b8,b8);
fprintf(f1,'%s(2 3 4 5)n%s);n%s}n', b12,b8,b4);
fprintf(f1,'%soutletn%s{n', b4,b4);
fprintf(f1,'%stype patch;n', b8);
fprintf(f1,'%sfacesn%s(n', b8,b8);
fprintf(f1,'%s(1 3 4 1)n%s);n%s}n', b12,b8,b4);
fprintf(f1,'%saxi_symm-fn%s{n', b4,b4);
fprintf(f1,'%stype wedge;n', b8);
fprintf(f1,'%sfacesn%s(n', b8,b8);
fprintf(f1,'%s(0 1 4 5)n%s);n%s}n', b12,b8,b4);
fprintf(f1,'%saxi_symm-fn%s{n', b4,b4);
fprintf(f1,'%stype wedge;n', b8);
fprintf(f1,'%sfacesn%s(n', b8,b8);
fprintf(f1,'%s(0 1 3 2)n%s);n%s}n);n', b12,b8,b4);
fprintf(f1,'nmergePatchPairsn(n);n');
fprintf(f1,'%s', footer);
fclose(f1)The blockmesh file generated is:
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: 7 |
| \ / A nd | Website: https://openfoam.org |
| \/ M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0 0)
(2.000000 0 0)
(0.000000 -0.000174 0.004997)
(2.000000 -0.000174 0.004997)
(2.000000 0.000174 0.004997)
(0.000000 0.000174 0.004997)
);
blocks
(
hex (0 1 1 0 2 3 4 5) (200 1 20) simpleGrading (1 1 0.2)
);
edges
(
arc 2 5 (0 0 0.005)
arc 3 4 (2.00 0 0.005)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 5 0)
);
}
wall
{
type wall;
faces
(
(2 3 4 5)
);
}
outlet
{
type patch;
faces
(
(1 3 4 1)
);
}
axi_symm-f
{
type wedge;
faces
(
(0 1 4 5)
);
}
axi_symm-f
{
type wedge;
faces
(
(0 1 3 2)
);
}
axis
{
type empty;
faces
(
(0 1 1 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //Inlet conditions:
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2006 |
| \ / A nd | Website: www.openfoam.com |
| \/ 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;
}
wall
{
type zeroGradient;
}
axi_symm-f
{
type wedge;
}
axi_symm-r
{
type wedge;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2006 |
| \ / A nd | Website: www.openfoam.com |
| \/ M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (15 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (15 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
wall
{
type noSlip;
}
axi_symm-f
{
type wedge;
}
axi_symm-r
{
type wedge;
}
}
// ************************************************************************* //
transportProperties:
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: 7 |
| \ / A nd | Website: www.openfoam.org |
| \/ M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu [0 2 -1 0 0 0 0] 1e-04;
// ************************************************************************* //
Post-process velocity and shear stress:
The graphs are plotted between the location (1.9 0 0) and (1.9 0 0.004997) in order to avoid entrance length or for a fully developed flow.
The plot is plotted in paraview using the instruction given in fresh desk.
Validation of the above two graphs using analytic solutions :
As in graph 1, the wall shear stress is nearly 120000N/sqm.
Also, we observe that maximum shear stress or wall shear stress is at the walls(location at 0.005m in x-axis) and zero at the centre (location at 0m in x-axis) due to no-slip condition.
The profile follows a linear path as in the graph.
As in graph 4, the velocity is nearly 30 m/sec.
Also, we observe that maximum velocity is at the centre (location at 0m in x-axis) and zero at the walls (location at 0.005m in x-axis) due to no-slip condition.
The profile follows a parabolic path as in the graph.