Week 11 - Simulation of Flow through a pipe in OpenFoam
%simulation for flow through a pipe
clear all
close all
clc
%parameters
D = 0.025;
Re = 1200;
L = (0.05*Re*D)+0.75;
theta = 2;
z = (D/2)*sind(theta/2);
y = (D/2)*cosd(theta/2);
g_x = 1;
g_y = 0.1;
g_z = 1;
v_avg = (Re*0.89e-3)/(1000*D)
v_max = 2*(v_avg)
r = linspace(0,(D/2),1000)
for i = 1:length(r)
tau(i) = 2*(0.89e-3)*v_max*(abs(r(i)))/(D/2)^2
end
plot(r,tau)
xlabel('Radius')
ylabel('Shear stress')
grid on
legend ('Shear stress')
f1 = fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'%s',' ========= |');
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: 6');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\/ M anipulation |');
fprintf(f1,'\n');
fprintf(f1,'%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'%s','FoamFile');
fprintf(f1,'\n');
fprintf(f1,'%s','{');
fprintf(f1,'\n');
fprintf(f1,'%s',' version 2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s',' format ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s',' class dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s',' object blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'%s','}');
fprintf(f1,'\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','vertices');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,' (%d %d %d)\n',0,0,0);
fprintf(f1,' (%d %d %d)\n',0,y,z);
fprintf(f1,' (%d %d %d)\n',0,y,-z);
fprintf(f1,' (%d %d %d)\n',L,0,0);
fprintf(f1,' (%d %d %d)\n',L,y,z);
fprintf(f1,' (%d %d %d)\n',L,y,-z);
fprintf(f1,');\n\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,' hex (0 3 5 2 0 3 4 1) (400 60 1) simpleGrading (%d %d %d)\n',g_x,g_y,g_z);
fprintf(f1,');\n\n');
fprintf(f1,'edges\n');
fprintf(f1,'(\n');
fprintf(f1,' arc 1 2 (0 %d 0)\n',D/2);
fprintf(f1,' arc 4 5 (%d %d 0)\n',L,D/2);
fprintf(f1,');\n\n');
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
fprintf(f1,' inletwall\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (2 0 0 1)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' outletwall\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (4 3 3 5)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' Noslipwall\n');
fprintf(f1,' {\n');
fprintf(f1,' type wall;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (5 2 1 4)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' front\n');
fprintf(f1,' {\n');
fprintf(f1,' type wedge;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (1 0 3 4)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' Back\n');
fprintf(f1,' {\n');
fprintf(f1,' type wedge;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (2 5 3 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' axis\n');
fprintf(f1,' {\n');
fprintf(f1,' type empty;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 3 3 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n\n')
fprintf(f1,'// ************************************************************************* //');
fclose(f1);/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 1.249810e-02 2.181551e-04)
(0 1.249810e-02 -2.181551e-04)
(2.500000e+00 0 0)
(2.500000e+00 1.249810e-02 2.181551e-04)
(2.500000e+00 1.249810e-02 -2.181551e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (400 60 1) simpleGrading (1 1.000000e-01 1)
);
edges
(
arc 1 2 (0 1.250000e-02 0)
arc 4 5 (2.500000e+00 1.250000e-02 0)
);
boundary
(
inletwall
{
type patch;
faces
(
(2 0 0 1)
);
}
outletwall
{
type patch;
faces
(
(4 3 3 5)
);
}
noslipwall
{
type wall;
faces
(
(5 2 1 4)
);
}
front
{
type wedge;
faces
(
(1 0 3 4)
);
}
back
{
type wedge;
faces
(
(2 5 3 0)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* ///*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.007;
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: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.893e-6;
// ************************************************************************* ///*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inletwall
{
type zeroGradient;
}
outletwall
{
type fixedValue;
value uniform 1;
}
noslipwall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* ///*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inletwall
{
type fixedValue;
value uniform (0.05 0 0);
}
outletwall
{
type zeroGradient;
}
noslipwall
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //In order to simulate the flow through a pipe, we have to write a blockMeshDict file using matlab code inside which we will be specifying all the parameters necessary to run the simulation. We first specify the given parameters, then open a .txt file, and then write all the steps inside the file and save it. Then, we open ubuntu, tutorials and inside tutorials incompressible inside which the solve icoFoam is chosen.
Inside icoFoam, we get cavity folder which we copy to our $FOAM_RUN folder using cp -r cavity $FOAM_RUN/cavity2 so that the original files remain untouched. then after copying, we replace all the data inside the blockMeshDIct file with the file we earlier created using matlab. Then, we specify all the necessary parameters such as the simulation run time, the time steps, the starting time and purge. Then we open the pressure file inide which we define the faces and their types, and similarly modify the velocity file as well. After entering all the required data, we give the value for nu inside the physicalproperties folder and svae it.
AFter saving all the files, we use blockMesh command to run the mesh and then use checkMesh command to check for the appropriateness of mesh. Then we can use icoFoam command to run the simulation for the stipulated time and use paraview using paraFoam command to view the results.
We check all the parameters we need to see and then we view the velocity and pressure profiles correspondingly and use plot by line command to plot the line graph for velocity and pressure at last time step. Then, we can use the plot command in matlab to get plot for the shear stress as well.
It is evident that the velocity profile at the entrance is not fully developed. This is due to the no slip condition specified. At the outlet the velocity profile is fully developed. Also, the comparison between the anlytical and the openfoam values for the velocity is near enough and the velocity is maximun at the center and decreases as a form or parabola towards the wall. SHear stress increases linearly towards the wall and has a value of 0 at the centre.
Velocity profile at inlet:
Velocity profile at the outlet:
Velocity at final time step:
Pressure at final time step:
Shear stress:
