Simulation of Flow through a pipe in OpenFoam part-2
Objectives :-
- Write a Matlab program that takes an angle as input and generates a blockMesh file for the given angle.
- Angles to test 10,25,45
- Show that the velocity profile is fully developed.
- Compare the above results and discuss.
Given Inputs :-
- Reynolds number(Re) = 2100
- Working fluid = Water
- Solver = icoFoam
- Density = 997 kg/m^2
- Kinematic viscosity = 1.004*10^-6 m^2/s
- Dynamic viscosity = 1.002*10^-3 N.s/m^2
Note :- We are already calculated all the flow variable in previous project at an angle 2 degree \"OpenFOAM pipe flow challenge - part 1\"
Calculated result from above input :-
- Diameter of pipe (D) = 0.05 m
- Radius (R) = D/2
- Length of pipe (L) = 0.05*Re*D = 5.25 m ~ 5.5 m
- Average Velocity (Vavg) = Re*mu/rho*D = 0.042168 m/s
- Maximum velocity (Vmax) = 2*Vavg = 0.084336 m/s
- Pressure difference(del_p) = 32*mu*V*L/D^2 = 2.8393 N/m^2
- Kinematic pressure = del_p/rho = 0.00284784 m^/s^2
- Maximum Shear Stress = 2*mu*Vmax/R = 0.676037 N/m^2
Assumption :-
- Fluid is incompressible.
- The flow is laminar through the pipe.
- Follow the newtonion law of fluid.
- Constant circular cross section throughout.
- No acceleration of fluid in the pipe.
Matlab code that can generate the computational mesh for any wedge angle and grading factor for symmetry boundary condition :-
% Code to create blockMeshDict file
clear all
close all
clc
L = 5.5; % Length of the pipe
theta = input(\'Wedge angle = \'); % Angle
D = 0.05; % Diameter of the pipe
r = D/2; % Radius of the pipe
nx = 200; % grid point along x-axis
ny = 30; % grid point along y-axis
nz = 1; % grid point along z-axis
% Creating blockMesh file
l1 = \'/*--------------------------------*- C++ -*----------------------------------*\\\';
l2 = \'| ========= | |\';
l3 = \'| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\';
l4 = \'| \\\\ / O peration | Version: 4.1 |\';
l5 = \'| \\\\ / A nd | Web: www.OpenFOAM.org |\';
l6 = \'| \\\\/ M anipulation | |\';
l7 = \'\\*---------------------------------------------------------------------------*/\';
Block = fopen(\'blockMeshDict\',\'w\');
fprintf(Block,\'%s\\n\',l1);
fprintf(Block,\'%s\\n\',l2);
fprintf(Block,\'%s\\n\',l3);
fprintf(Block,\'%s\\n\',l4);
fprintf(Block,\'%s\\n\',l5);
fprintf(Block,\'%s\\n\',l6);
fprintf(Block,\'%s\\n\',l7);
fprintf(Block,\'FoamFile\\n{\\n\');
fprintf(Block,\'%12s \\t 2.0;\\n\',\'version\');
fprintf(Block,\'%11s \\t ascii;\\n\',\'format\');
fprintf(Block,\'%10s \\t dictionary;\\n\',\'class\');
fprintf(Block,\'%11s \\t blockMeshDict;\\n\',\'object\');
fprintf(Block,\'}\\n\');
l8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
fprintf(Block,\'%s\\n\\n\',l8);
fprintf(Block,\'convertToMeters 1;\\n\\n\');
fprintf(Block,\'vertices\\n\');
fprintf(Block,\'(\\n\');
fprintf(Block,\'\\t (%f %f %f)\\n\',0,0,0); % point 0
fprintf(Block,\'\\t (%f %f %f)\\n\',0,r*cosd(theta/2),r*sind(theta/2)); % point 1
fprintf(Block,\'\\t (%f %f %f)\\n\',0,r*cosd(theta/2),-r*sind(theta/2)); % point 2
fprintf(Block,\'\\t (%f %f %f)\\n\',L,0,0); % Pont 3
fprintf(Block,\'\\t (%f %f %f)\\n\',L,r*cosd(theta/2),r*sind(theta/2)); % point 4
fprintf(Block,\'\\t (%f %f %f)\\n\',L,(r*cosd(theta/2)),-r*sind(theta/2)); % point 5
fprintf(Block,\');\\n\\n\');
fprintf(Block,\'blocks\\n\');
fprintf(Block,\'(\\n\');
fprintf(Block,\'\\t hex ( 0 3 5 2 0 3 4 1) (%d %d %d)\',nx,ny,nz);
fprintf(Block, \' SimpleGrading (1 0.1 1)\\n\\n\');
fprintf(Block,\');\\n\\n\');
fprintf(Block,\'edges\\n(\\n\\t arc %d %d (%f %f %f)\\n\',1,2,0,r,0);
fprintf(Block,\'\\t arc %d %d (%f %f %f)\\n\',4,5,L,r,0);
fprintf(Block,\');\\n\\n\');
fprintf(Block,\'boundary\\n(\\n\');
fprintf(Block,\'\\t inlet\\n\');
fprintf(Block,\'\\t{\\n\');
fprintf(Block,\'\\t\\t type patch;\\n\');
fprintf(Block,\'\\t\\t faces\\n\');
fprintf(Block,\'\\t\\t (\\n\');
fprintf(Block,\'\\t\\t\\t (0 1 2 0)\\n\');
fprintf(Block,\'\\t\\t );\\n\');
fprintf(Block,\'\\t}\\n\');
fprintf(Block,\'\\t outlet\\n\');
fprintf(Block,\'\\t{\\n\');
fprintf(Block,\'\\t\\t type patch;\\n\');
fprintf(Block,\'\\t\\t faces\\n\');
fprintf(Block,\'\\t\\t (\\n\');
fprintf(Block,\'\\t\\t\\t (3 5 4 3)\\n\');
fprintf(Block,\'\\t\\t );\\n\');
fprintf(Block,\'\\t}\\n\');
fprintf(Block,\'\\t top \\n\');
fprintf(Block,\'\\t{\\n\');
fprintf(Block,\'\\t\\t type wall;\\n\');
fprintf(Block,\'\\t\\t faces\\n\');
fprintf(Block,\'\\t\\t (\\n\');
fprintf(Block,\'\\t\\t\\t (1 4 5 2)\\n\');
fprintf(Block,\'\\t\\t );\\n\');
fprintf(Block,\'\\t}\\n\');
fprintf(Block,\'\\t back\\n\');
fprintf(Block,\'\\t{\\n\');
fprintf(Block,\'\\t\\t type symmetry;\\n\');
fprintf(Block,\'\\t\\t faces\\n\');
fprintf(Block,\'\\t\\t (\\n\');
fprintf(Block,\'\\t\\t\\t (0 2 5 3)\\n\');
fprintf(Block,\'\\t\\t );\\n\');
fprintf(Block,\'\\t}\\n\');
fprintf(Block,\'\\t front\\n\');
fprintf(Block,\'\\t{\\n\');
fprintf(Block,\'\\t\\t type symmetry;\\n\');
fprintf(Block,\'\\t\\t faces\\n\');
fprintf(Block,\'\\t\\t (\\n\');
fprintf(Block,\'\\t\\t\\t (0 3 4 1)\\n\');
fprintf(Block,\'\\t\\t );\\n\');
fprintf(Block,\'\\t}\\n\');
fprintf(Block,\'\\t axis\\n\');
fprintf(Block,\'\\t{\\n\');
fprintf(Block,\'\\t\\t type empty;\\n\');
fprintf(Block,\'\\t\\t faces\\n\');
fprintf(Block,\'\\t\\t (\\n\');
fprintf(Block,\'\\t\\t\\t (0 3 3 0)\\n\');
fprintf(Block,\'\\t\\t );\\n\');
fprintf(Block,\'\\t}\\n\');
fprintf(Block,\');\\n\\n\');
fprintf(Block,\'mergePatchPairs\\n(\\n\');
fprintf(Block,\');\\n\');
line9=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
fprintf(Block,\'%s\\n\',line9);
fclose(Block);blockMeshDict file generated by Matlab for an angle 2 degree with wedge BC:-
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.024996 0.000436)
(0.000000 0.024996 -0.000436)
(5.250000 0.000000 0.000000)
(5.250000 0.024996 0.000436)
(5.250000 0.024996 -0.000436)
);
blocks
(
hex ( 0 3 5 2 0 3 4 1) (200 30 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0.000000 0.025000 0.000000)
arc 4 5 (5.250000 0.025000 0.000000)
);
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)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //blockMeshDict file generated by Matlab for an angle 10 degree with symmetry BC:-
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.024905 0.002179)
(0.000000 0.024905 -0.002179)
(5.500000 0.000000 0.000000)
(5.500000 0.024905 0.002179)
(5.500000 0.024905 -0.002179)
);
blocks
(
hex ( 0 3 5 2 0 3 4 1) (200 30 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0.000000 0.025000 0.000000)
arc 4 5 (5.500000 0.025000 0.000000)
);
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)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //blockMeshDict file generated by Matlab for an angle 25 degree with symmetry BC:-
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.024407 0.005411)
(0.000000 0.024407 -0.005411)
(5.500000 0.000000 0.000000)
(5.500000 0.024407 0.005411)
(5.500000 0.024407 -0.005411)
);
blocks
(
hex ( 0 3 5 2 0 3 4 1) (200 30 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0.000000 0.025000 0.000000)
arc 4 5 (5.500000 0.025000 0.000000)
);
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)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //blockMeshDict file generated by Matlab for an angle 45 degree with symmetry BC:-
/*--------------------------------*- C++ -*----------------------------------*\\
| ========= | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 4.1 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\/ M anipulation | |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.023097 0.009567)
(0.000000 0.023097 -0.009567)
(5.500000 0.000000 0.000000)
(5.500000 0.023097 0.009567)
(5.500000 0.023097 -0.009567)
);
blocks
(
hex ( 0 3 5 2 0 3 4 1) (200 30 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0.000000 0.025000 0.000000)
arc 4 5 (5.500000 0.025000 0.000000)
);
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)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Velocity Boundary condition :-
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ M anipulation |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
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.042168 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //Pressure Boundary condition :-
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ 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.0028393;
}
top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //For the simulation purpose we have to edit controlDict file accordingly that courant number less than 1. So the solution is converge.
Total time for simulation is 50 sec.
controlDict file :-
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ M anipulation |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location \"system\";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Transport properties - Kinematic viscosity :-
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 6
\\\\/ M anipulation |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location \"constant\";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1.004e-6;
// ************************************************************************* //Note :- All the fluid properties are same as the previous challenges \"OpenFOAM pipe flow challenge - part 1\". In this challenge we are comparing our result on different wedge angle.
Result and plot of velocity profile at different wedge angle :-
Fig1 : Mesh generation on the profile with grading 0.1.
Fig2 : Velocity distribution profile at starting of the pipe.
Fig3 : Velocity distribution profile at end of the pipe.
Fig4 : Pipe flow model with an angle 2 degree.
Fig5 : Pipe flow model with an angle 10 degree.
Fig6 : Pipe flow model with an angle 25 degree.
Fig7 : Pipe flow model with an angle 45 degree.
Fig8 : Fully developed velocity profile for an angle 2 degree.
Fig9 : Fully developed velocity profile for an angle 10 degree.
Fig10 : Fully developed velocity profile for an angle 25 degree.
Fig11 : Fully developed velocity profile for an angle 45 degree.
Table : Comparision the result between the different wedge angle.
Note :- In this table wedge angle 2 degree is with wedge BC and rest of all 10,25 and 45 degree with symmetry BC.
From the above table it shown that as the wedge angle increases the error and the difference between theoritical and simulation result will decreases but the problem in greater angle is that it decreases the radial distance of the profile.
Conclusion :-
- Matlab code is for creating the model of the pipe and generating mesh in both x and y direction (x direction means along the pipe and y direction means radially).
- There will be velocity and pressure boundary condition for the input and output of the fluid in the pipe.
- Editing the controlDict file accordingly that courant number less than 1. So the solution is converge.
- In this we have to edit the transport properties file or the value of kinematic viscosity according the fluid properties.
- Fig 2 and 3 shows velocity distribution profile at start and end of the pipe.
- Fig 4,5,6 and 7 shows the pipe flow model at different wedge angle.
- Fig 8,9,10 and 11 shows fully developed velocity profile for different wedge angle.
- So, as the wedge angle increases error will decreases.
- As the wedge angle increases, velocity profile is getting srink in radially (see Fig11) 0.025 to in between 0.025-0.02. So it not give correct result for the simulation.
- So by the above comparision it will find that, at wedge angle 2 degree with wedge BC velocity profile should be fully developed and it take lower simulation time.