Week 11 - Simulation of Flow through a pipe in OpenFoam
The entrance length is the distance that the fluid need to travel through before the full development of the flow. see the following image
the entry length for laminar flow is
L = 0.06*Re*D
the study variable values that I chose
Reyonlds number, re = 2100
Diameter, d = 0.004
radius = d/2 = 0.002
theta = 4
Density, rho = 997
Dynamic viscosity, mu = 0.00089
According to Hagen pioseuille's law we can calculate the following variable
Avarage velocity V_avg= re*mu/(rho*d)
Maximum velocity = 2*V_avg
Pressure differemce deltap = 8*mu*L*V_avg/r^2 where r is the distance from the center
Kinematic pressure difference deltap_k = deltap/rho
Shear stress tau = 2*mu*V_max*radius/r^2
to plot the shear stress the velocity field should be defined
for i= 1:length(r)
U(i) = u_max*abs(1-r(i)^2/radius^2);
end
then the related shear stress values could be calculated
for i= 1:length(r)-1
tau1(i) = mu*(U(i+1)-U(i))/dr;
end
the result will be as following graph
The Blockmesh programm:
First defining all requiewd variables then open a file to write
f1 = fopen('blockMeshDict', 'w');
after that using fprintf command to write in that file
fprintf(f1, '%s\n', header1);
where '\n' means new line. '\t' means tab, %s means string variable and %d means double variable
vertices file part is writen using repeated set of command according to the number of vertices
r1 = r*cosd(theta/2);
r2 = r*sind(theta/2);
fprintf(f1, 'vertices\n');
fprintf(f1, '(\n');
fprintf(f1, '\t (%d %d %d)\n',0,0,0);
fprintf(f1, '\t (%d %d %d)\n',0,r1,r2);
fprintf(f1, '\t (%d %d %d)\n',0,r1,-r2);
fprintf(f1, '\t (%d %d %d)\n',l_pipe,0,0);
fprintf(f1, '\t (%d %d %d)\n',l_pipe,r1,r2);
fprintf(f1, '\t (%d %d %d)\n',l_pipe,r1,-r2);
fprintf(f1, ');\n\n');
then define the blockes with the mesh number of elements and the scaling factor
fprintf(f1, 'blocks\n');
fprintf(f1, '(\n');
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1);
fprintf(f1, ' SimpleGrading (1 0.1 1)\n\n');
fprintf(f1, ');\n\n');
after that the edges of the structure
fprintf(f1, 'edges\n');
fprintf(f1, '(\n');
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,l_pipe,r,0);
fprintf(f1, ');\n\n');
at the end the boundary conditions for inlet, outlet, top, front, back and axis
fprintf(f1, 'boundary\n');
fprintf(f1, '(\n');
fprintf(f1, '\t inlet \n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (0 1 2 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
then close the file
fclose(f1)
The blockMesh file is generated as following
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeter 1;
vertices
(
(0 0 0)
(0 1.998782e-03 6.979899e-05)
(0 1.998782e-03 -6.979899e-05)
(1.204000e+00 0 0)
(1.204000e+00 1.998782e-03 6.979899e-05)
(1.204000e+00 1.998782e-03 -6.979899e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (1.204000e+00 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
the other controll files
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.4686 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.4686 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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 1.0073;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9256e-7;
// ************************************************************************* //
the velocity result near the wall is almost zero and increases when moving towards the center. check the next photo
the preesure and velocity profile over the pipe X axis
To show the development of the velocity profile over Y axis, I plot the velocity on line parallel to Y axis on different X values
We can estimate the entry length from previous figures as 0.4m which comparable to theoretical value 0.5m.
For post shear stress processing:
shear stress in fluid developes as a result of particles relative movement to each other.
I ploted the velocity on Y axis then I used grad function to calculate shear stress in "paraFoam" after that I exported the file and ploted the result with theoretical shear stress. There is error between the two figures because of the discretization error. this error reduces when having finer mesh.
clear all
close all
clc
% the Reynold number
re = 2100;
% the diamnter of the pipe
d = 0.004;
% the entry length of pipe
l = 0.06*re*d;
% the total length of pipe
l_pipe = l+0.7;
% pipe radius
r = d/2;
% angle of the wedge
theta = 4;
% viloscity
mu = 0.00089;
% density
rho = 997;
% Hagen poisuille flow equation
v_mean = (re*mu)/(rho*d);
v_max = 2*(v_mean);
p_del = (8*mu*l_pipe*v_mean)/(r^2);
p_kine = p_del/rho;
header1 = '/*--------------------------------*- C++ -*----------------------------------*\';
header2 = ' ========= |';
header3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
header4 = ' \\ / O peration | Website: https://openfoam.org';
header5 = ' \\ / A nd | Version: 8';
header6 = ' \\/ M anipulation |';
header7 = '\*---------------------------------------------------------------------------*/';
f1 = fopen('blockMeshDict', 'w');
fprintf(f1, '%s\n', header1);
fprintf(f1, '%s\n', header2);
fprintf(f1, '%s\n', header3);
fprintf(f1, '%s\n', header4);
fprintf(f1, '%s\n', header5);
fprintf(f1, '%s\n', header6);
fprintf(f1, '%s\n', header7);
fprintf(f1, 'FoamFile\n{\n');
fprintf(f1, '%12s \t 2.0;\n', 'version');
fprintf(f1, '%12s \t ascii;\n', 'format');
fprintf(f1, '%12s \t dictionary;\n', 'class');
fprintf(f1, '%12s \t blockMeshDict;\n', 'object');
fprintf(f1, '}\n');
header8 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1, '%s\n\n', header8);
fprintf(f1, 'convertToMeter 1;\n\n');
r1 = r*cosd(theta/2);
r2 = r*sind(theta/2);
fprintf(f1, 'vertices\n');
fprintf(f1, '(\n');
fprintf(f1, '\t (%d %d %d)\n',0,0,0);
fprintf(f1, '\t (%d %d %d)\n',0,r1,r2);
fprintf(f1, '\t (%d %d %d)\n',0,r1,-r2);
fprintf(f1, '\t (%d %d %d)\n',l_pipe,0,0);
fprintf(f1, '\t (%d %d %d)\n',l_pipe,r1,r2);
fprintf(f1, '\t (%d %d %d)\n',l_pipe,r1,-r2);
fprintf(f1, ');\n\n');
fprintf(f1, 'blocks\n');
fprintf(f1, '(\n');
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1);
fprintf(f1, ' SimpleGrading (1 0.1 1)\n\n');
fprintf(f1, ');\n\n');
fprintf(f1, 'edges\n');
fprintf(f1, '(\n');
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,l_pipe,r,0);
fprintf(f1, ');\n\n');
fprintf(f1, 'boundary\n');
fprintf(f1, '(\n');
fprintf(f1, '\t inlet \n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (0 1 2 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t outlet\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (3 5 4 3)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t top\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wall;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (1 4 5 2)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t front\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wedge;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (0 3 4 1)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t back\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wedge;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (0 2 5 3)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t axis\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type empty;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t(\n');
fprintf(f1, '\t\t\t (0 3 3 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, ');\n');
fprintf(f1, 'mergePatchPairs\n');
fprintf(f1, '(\n');
fprintf(f1, ');\n');
header9 = '// ************************************************************************* //';
fprintf(f1, '%s\n', header9);
fclose(f1);
clear all
close all
clc
d = 0.02;
radius = d/2;
re = 2100;
rho = 997.13;
mu = 8.9e-4;
nu = mu/rho;
v = re*nu/d;
u_max = 2/v;
u_all= csvread('tau.csv',1);
tau2 = u_all((1:1000),2)';
r = u_all((1:1000),1)';
U = u_max*abs(1-(r.^2./radius^2));
for i= 1:length(r)-1
tau1(i) = -mu*(U(i+1)-U(i))/(r(i+1)-r(i));
end
figure(1)
plot(r(1:999),tau1);
xlabel('Distance in m');
ylabel('Shear stress in N/m^2');
title('shear stress on radial distance');
grid minor
hold on
plot(r,tau2);
xlabel('Distance in m');
ylabel('Shear stress in N/m^2');
title('shear stress on radial distance OpenFOAM');
grid minorShear stress result file
|