Simulation of Flow through a pipe in OpenFoam
Simulation of Flow through a pipe in OpenFoam
Aim:-
- Write a program in Matlab that can generate the computational mesh automatically for any wedge angle.
- Provide proper grading scheme for block meshing of pipe.
- Post process velocity profile and shear stress.
- Show that the velocity profile is fully developed.
- Compare velocity profile calculated by numerical method with the velocity profile calculated with Hagen_poiseuille equation.
Observation:-
- First, understand Hagen-Poiseulle equation for laminar incompressible fluid.
- Then, start with defining problem statement for the fluid flowing through pipe. Mention all the inlet parameters of the working fluid.
- Open openFoam tutorials and find the the solver to solve incompressible-laminar-viscous flow problem.
- Next step is to copy required file in solver and paste it in OpenFOAM run folder so that the official file of OpenFoam remain untouched.
- Program block mesh script file in Matlab to create domain space and provide mesh specification.
- Find ‘controlDict’ file to edit start and end time for simulation.
- In initial folder ‘0’ edit initial parameter for velocity and pressure according to faces of the given domain.
- Use blockMesh command to generate mesh and checkMesh to verify it.
- Open simulation file in paraView and plot the velocity profile at fully developed point from inlet.
- Next, calculate pressure difference over a length of file.
- Also, calculate shear stress and plot its graph in fully developed region.
- Extract all the data from paraFoam to read it in Matlab program.
- To compare data calculate velocity and pressure using Hagen Poiseuill’s equation.
- Plot velocity graph of fully developed flow for both numerical and analytical method.
- We can also compare pressure difference plot of complete length of pipe by both numerical and analytical method.
- Write a detailed report explaining all the command and codes used in project.
Hagen-Poiseuille equation
In non fluid dynamics, the Hagen-Poiseuille equation is a physical law that gives the pressure drop in an incopressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
The assumption of the equation are
- The fluid is incompressible and Newtonian, which means, the viscosity property of fluid follows Newton’s law of viscosity, that is, `tau=mu(du)/(dy)`
- There is no relative motion between fluid particles and solid boundaries, that is no slip of fluid at the solid boundary.
- The flow is laminar through a pipe of constant circular cross section that is substantially longer than its diameter
- There is no acceleration of fluid in the pipe.
Poiseuille’s Equation describes the pressure drop due to the viscosity of the fluid.
`DeltaP = (8muLQ)/(piR^4) = (8pimuLQ)/(A^2)`
Where
`DeltaP` is the pressure difference between the two ends
L is the length of pipe
`mu` is the dynamic viscosity
Q is the volumetric flow rate
R is the pipe radius
A is the cross section of pipe
There are some limitation for Hagen-Poiseuille’s equations
- The equation dose not hold close to the pipe entrance.
- The equation fails in the limit of low viscosity, wide and/or short pipe. Low viscosity or a wide pipe may result in turbulent flow, making it necessary to use more complex models, such as the Darcy-Weisbach equation.
As fluid moves downstream from the entrance to a channel, its speed distribution or velocity profile changes in shape until it becomes fully developed. The length of an arrow indicates the magnitude of the flow velocity at the location of the arrow.
The Hagen-Poiseuille equation is valid only for the condition of fully developed flow, where the parabolic velocity profile is no longer changing with distance along the cylinder. As shown in above figure. For some distance downstream of the entrance to the cylinder, the velocity profile will be changing with distance. The entrance length usually is defined as the distance from the entrance to the location where the velocity at the center of the fully developed flow downstream. The entrance length, L_E may be estimated easily as
`L_E = 0.058 R_e d,`
Where,
`R_e` is the dimensionless Reynolds number,
% Entrance Length
Le = 0.058*Re*D;`R_e = (rho DV)/mu`
Where,
D is the internal diameter of the cylinder(m)
`rho` is density of the fluid (kg/`m^3`)
V is the average velocity of the fluid inside the cylinder(m/s)
Therefore, the entrance length will increase with flow rate.
Problem statement
For flow in pipe of diameter D, experimental observations shows that for “fully developed” flow, laminar flow occurs when `Re_D` < 2300. Therefore, let us, consider Re for flow is 2000. Assume diameter of pipe (D)=0.02m and length (L) = 2.5 m. also consider working fluid as water at temperature 25oC so the values of density (`rho`)=997kg/`m^3` and Dynamic Viscosity (`mu`)=8.90X10-4Pa.sec.
Pre-processing
In order to start the simulation process, First start terminal in Linux Operating system using ctrl+alt+T buttons on key board or right click -> Open Terminal. After starting terminal yo will see following window.
Now, we have to simulate our problem in openFOAM and we are going to solve our problem using tutorials present in openFOAM, therefore, just type ‘tut’ to open OpenFOAM tutorial directory. And type ‘ls’ to identify what kind of tutorials present inside this directory.
We have to solve incopressible flow simulation for this problem so select incompressible directory by typing 'cd incompressible' and then again type 'ls' to explore folders in incopressible folders.
In above picture, we can see the plethora of solver folders. For this particular solver we need solver for laminar-viscous flow. So here we choose 'icoFoam'
icoFoam is the solver in openFoam which is use as transient solver for incopressible, laminar flow of Newtonian fluids.
(Newtonian fluid is a fluid in which the viscous stress arising from its flow, at every point, are linearly correlated to the local strain rate - the rate of change of its deformation over time. Newtonian fluids are the simplest mathematical models of fluids that account for viscosity.)
Therefore, we will select icoFoam folder by typing 'cd icoFoam' and by typing 'ls' we explore this folder in terminal.
Here, we have, again to option, cavity and elbow. They both are part of the icoFoam solver and cavity section use to create required domain and black mesh so select 'cavity' folder.
In cavity, we are, again, going to see cavity folder, where we have to define our fluid domain and generate block mesh.
But we are not going to make any changes here in tutorial section. First of all, OpenFoam does no let anyone make any changes in their official tutorial and secondly, any changes with their official tutorials my corrupt openFoam. An in a bid to correct it, we may need to reinstall open foam all over again.
So, we are going to copy this cavity folder using 'cp -r cavity' command. Next, rename it as ‘Pipe’ and paste it in openFoam run folder by typing 'cp -r cavity $FOAM_RUN/Pipe'.
And then open run folder by typing 'cd $FOAM_RUN'. where we can see 'Pipe' file.
Now, we can make any changes to the Pipe file, which will never affect the openFoam official files. Open Pipe folder. Here, we can provide specification for our block mesh in system folder.
As we can see in above picture, there are 4 script file system folder and by opening blockMeshDict we can create domain and provide required specification for block mesh.
FILE WRITING (in Matlab)
But, this time, we are going to create blockMeshDict file using Matlab coding by file parsing techniques. Prior, lest define all the given data first in our main program
Condition of working fluid according to problem statement.
Reynolds number (`R_e`) = 2000
Diameter of the pipe (D) = 0.02m
Length of the pipe (L) = 2.5m
Density of working fluid (`rho`) = 997 kg/`m^3`
Dynamic viscosity of working fluid (`mu`) = 8.9X10^-4 Pa.sec
%% Input parameters
% Reynold's Number
Re = 2000;
% Total Length of pipe(m)
L = 2.5;
% Diameter of pipe in meters
D = 0.02;
% Wedge Angle in degrees
th = input('<0> Wedge Angle in degrees - ');
% Density (kg/m^3)
rho = 997;
% Dynamic Viscosity (Pa.sec)
mu = 8.90e-4;
% Radius of pipe(m)
R = D/2;Now, let us create separate function file to write blockMeshDict file in Matlab program. Create blockMeshDict file and start writing it using following ‘fopen’ command.
function[f1] = blockMesh(R,th,L,nx,ny,nz)
%% Writing BlockMeshDict file for openFoam
f1 = fopen('blockMeshDict','w');Before defining any data we are going to state with printing header information in the form of comment which is not executable in OpenFoam
fprintf(f1,"/*--------------------------------*- C++ -*----------------------------------*\\ \n");
fprintf(f1,' ========= | \n');
fprintf(f1,' \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \n');
fprintf(f1,' \\\\ / O peration | Website: https://openfoam.org \n');
fprintf(f1,' \\\\ / A nd | Version: 8 \n');
fprintf(f1,' \\\\/ M anipulation |\n');
fprintf(f1,"\\*---------------------------------------------------------------------------*/ \n");Next, we set units in meter, by converting all type of length in to meter.
fprintf(f1,'convertToMeters 1; \n\n');Vertices:-
Afterward, in codes, we provide list of vertices by providing x,y,and z coordinates locations; its always start from vertex (0,0,0). so we have to create following geometry.
Mesh generation
The mesh of the Pipe case includes a new feature: a block with fewer than 8 vertices will be used. There is an axi-symmetry in the case so it is not necessary to mesh the whole pipe but only a wedge of it. Since `del/(deltheta) = 0` for the whole domain, the results of the simulation will be the same while the number of elements of the mesh will be lower. Following figure represents a schematic view of the geometry that it is going to be created:
First of all, it is important to comprehend the collapsing of the vertices. It can be seen that vertices 0 and 3 are repeated in an adequate position within the parentheses to indicate that the two vertices that would be occupying this position have been collapsed.
Secondly, in the patches definition, the same vertices 0 and 3 are joined creating an empty patch forming the axis of the pipe. The patch type wedge is used in two faces to indicate that an axi-symmetry exists and there are no physical walls.
fprintf(f1,'vertices\n(\n');
fprintf(f1,' (0 0 0)\n');
fprintf(f1,(' (%f %f 0)\n'),R*cosd(th/2),R*sind(th/2));
fprintf(f1,(' (%f %f %f)\n'),R*cosd(th/2),R*sind(th/2),L);
fprintf(f1,(' (0 0 %f)\n'),L);
fprintf(f1,(' (%f %f 0)\n'),R*cosd(th/2),-R*sind(th/2));
fprintf(f1,(' (%f %f %f)\n);\n\n'),R*cosd(th/2),-R*sind(th/2),L);Next, we are going to define block for meshing. If we look at he mesh format first letter, in block parenthesis, is ‘hex’ which define type of mesh which is hexahedral mesh.
Then, again in parenthesis, we provide vertices number according to order which follow right hand thumb rule, meaning, thumb is the z axis and order of vertices follow rotational order of our other finger.
Then comes, number of mesh grid that we want to divide along x, y, and z length. We define it as (nx,ny,nz) format in main program.
Given mesh specification:
- Number of cells along the zdirection(longer dimension) = 500
- Number of cells along the longer xdirection = 20
- Number of cell for y direction =1
% Number of cells in x,y,z Directions
nx = 20;
ny = 1;
nz = 500;Lastly, we give simple grading which is cell expansion ratios
Grading factor or expansion ratios:-
The expansion ratios enables the mesh to be graded, or refined, in specified directions. The ratio is that of the width of the end cell `delta_e` along one edge of a block to the width of the start cell `delta_s` along the edge as shown in below figure.
In simple format,
Expansion ratio = (size of end cell) / (size of starting cell)
The simple description specifies uniform expansion in the local `x_1`,`x_2`and `x_3` directions respectively with only 3 expansion ratio. So here, we have 1 wedge shape block so we have to define hex mesh for 1 block.
There is a grading towards the wall to compute a more exact value of wall shear stress and towards the inlet where the flow is not fully developed and therefore the velocity profile changes.
fprintf(f1,('hex (0 4 1 0 3 5 2 3) (%d %d %d) simpleGrading (0.2 1 1)\n);\n\n'),nx,ny,nz);In this case, for edge 4-1 and 5-2 we have radius curve so we will define this in edges section following way.
fprintf(f1,'edges\n(\n');
fprintf(f1,' arc 4 1 (0 %f 0)\n',R);
fprintf(f1,' arc 5 2 (%f %f 0)\n);\n\n',L,R);Next, define boundary of the problem, in this we have to define faces of the our 3D block this should be in order according to right hand thumb rule. So set faces for this problem following way.
fprintf(f1,'boundary\n(\n');
fprintf(f1,' inlet\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 4 1 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' outlet\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (3 5 2 3)\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,' (0 1 2 3)\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,' (0 3 5 4)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' pipewall\n');
fprintf(f1,' {\n');
fprintf(f1,' type wall;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (1 4 5 2)\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,'// ************************************************************************* // \n');
fclose(f1);After running above codes, we will get our BlockMeshDict file as output.
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0.009990 0.000436 0)
(0.009990 0.000436 2.500000)
(0 0 2.500000)
(0.009990 -0.000436 0)
(0.009990 -0.000436 2.500000)
);
blocks
(
hex (0 4 1 0 3 5 2 3) (20 1 500) simpleGrading (0.2 1 1)
);
edges
(
arc 4 1 (0 0.010000 0)
arc 5 2 (2.500000 0.010000 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 4 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 2 3)
);
}
front
{
type wedge;
faces
(
(0 1 2 3)
);
}
back
{
type wedge;
faces
(
(0 3 5 4)
);
}
pipewall
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* // Copy above blockMeshDict file and replace it in the system directory of the pipe folder and by running this block mesh file we will get following mesh:
Furthermore, we are going to edit ‘controlDict’ file to provides time limits - start time and end time. Also to represent deltat. While selection of deltat CFL value shoud be less than 1. otherwise, the solution could stopped forcibly or it may diverge at the time of converging. So make following changes in ‘controlDict’ file.
/*--------------------------------*- 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 startTime;
startTime 0;
stopAt endTime;
endTime 2;
deltaT 0.0001;
writeControl timeStep;
writeInterval 10;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Boundary and initial conditions
Then, we have to define initial value for pressure and velocity. In order to do that, we have to go back in cavity_copy directory using ‘cd ..’ command and select folder named ‘0’, which stand for initial values. In this case initial value for pressure and velocity.
When solving circulating internal flow, the most common boundary conditions are to fix an inlet velocity and an outlet pressure.
Initial setting for pressure
/*--------------------------------*- 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 0;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
pipewall
{
type zeroGradient;
}
}
// ************************************************************************* //In the U dictionary one can find out a new type of boundary condition: inlet and outlet. It switches U and p between fixedValue and zeroGradient: the first when the flow is in-going and the second when it is outgoing.
Example: in circular pipe case, the outlet velocity will always be adapted to the conditions according to zeroGradient except if there is an inflow at the outlet patch. If this happens, to this inflow velocity it is assigned the value specified under the inletOutlet instruction and so the flow will always be outgoing.
Initial setting for Velocity
/*--------------------------------*- 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 0.0893);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
pipewall
{
type noSlip;
}
}
// ************************************************************************* //Next, we have to mention dynamic viscosity in transposeProperty in constant directory of the Pipe.
/*--------------------------------*- 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.90e-4;
// ************************************************************************* //After deploying mesh, now we are ready to run the simulation. Simulation that we selected was icoFoam. Hence, just type ‘icoFoam’ and simulation will start.
Post-processing
Now, let us view the results in paraView by using ‘paraFoam’ command
RESULTS OF THE SIMULATION
The velocity at the inlet of the pipe (non-developed flow)
In above figure we can notice that the shear stress along wall surface is starting to increase as fluid moving forward.
The velocity field at the outlet of the pipe (developed flow)
And at the end of the pipe, we can see, developed wall shear stress and velocity remain stable.
PRESSURE FIELD ALONG THE PIPE
Pressure at inlet of pipe
For pressure, we can see, the pressure at the beginning is high and as fluid is moving forward through pipe its decreasing.
Pressure at outlet of the pipe
This is the lowest pressure at the end of pipe as depicted in above picture.
Velocity at various distance from the inlet of pipe and comparison of it with fully developed velocity flow at distance 2.32 from inlet of the pipe
As we can see above graph, its portrays the velocity profile along radius at various distance from the inlet. So we can observe that as water moving forward through pipe its velocity profile along the radius of pipe is stabilizing and there is a point within a length of point where it does not change that point is call fully developed flow region. And it comes after entrance length which we calculated in above section of this project.
Pressure change along the length of pipe
If we check the pressure plot along the length of the pipe its gradually decreasing as the velocity increasing.
Shear stress along the radius at developed region of the pipe
Next, we will calculate shear stress
Plotting shear stress in paraView is not a straight forward task. The reason for this is that there is no direct option to view the shear-stress tensor field inside the domain unlike Pressure and Velocity fields which can be accessed directly at a single click. So by using following techniques we compute the shear stress during the post-processing of simulation in ParaView.
1] Computing the magnitude of the vector field U using the Calculator tool.
By clicking on the calculator tool from extreme left of the common toolbar create the magnitude of velocity “U”.
2] Computing the gradients for U-magnitude using
Now, Press ‘Ctrl+space’ or click Filters -> search in paraview. A search widow will pop up for the mentioned filter. Type gradient of Unstructured Data set and select it.
The property window changes as follows. Select scalar array as mag_u i.e. magnitude of velocity and give required result array name and click apply which will give us a gradient of the magnitude of velocity.
3] Compute the shear-stress quantity
Shear stress formula is given as
`Tau = mu (du)/(dy)`
Where,
`mu` is Dynamic viscosity
`tau` is Shear stress = `F/A`
`(du)/(dy)` is Rate of shear deformation
Therefore, we can obtain shear stress by multiplying the velocity gradient with the dynamic viscosity of the fluid under study.
So, in this project, we have considered the water as working fluid at 25oC room temperature which dynamic viscosity is 8.90e-4 N/`m^2`.
Then, again click on the Calculator tool and same way create tau array using product of velocity gradient and dynamic viscosity.
After plotting shear stress we will get following plot as our result.
If we focus on shear stress along radius of the pipe, in fully developed flow region, we can observe there is very less shear stress at the center as there is no wall is present to restrict the flow of fluid. But from center to wall the, we can see, the shear stress increasing gradually and its value near the wall is maximum.
Next, plotting above results we can export above data from paraView to read it in matlab
Click on the spread sheet to open spread sheet of all data related to simulation.
After opening spread sheet it will look like below picture. And by clicking upper right corner export data button we can export our data in .csv format.
FILE READING
For file reading we are going to make separate function program.
function [L_n,V_n,L_nn,P_n,Tau_n,L_nnn] = velocity_pressure_file_reading()First to read file we are going to keep that particular .csv file, that we have to read in matlab, into same directory where our program is located.
Start reading file using following fopen command
%% Reading Velocity File Extracted from paraView
f2 = fopen('fully_developed_flow_data.csv','r');We can use fgetl to read lines or skip lines sequence wise so by typing fgetl we will skip first line and then create loop to read file line by line and store velocity and length data in V_n and L_n.
Lastly, we use str2double to convert string data to number data.
A=fgetl(f2);
for i = 1:1000
C = fgetl(f2);
D = strsplit(C,',');
V_n(i) = D(3);
L_n(i) = D(6);
end
V_n=str2double(V_n');
L_n=str2double(L_n');Similar way, we will read and store pressure data in Matlab and end the function.
%% Reading Pressure File Extracted from paraView
f3 = fopen('Pressure.csv','r');
A=fgetl(f3);
for i = 1:1000
C = fgetl(f3);
D = strsplit(C,',');
P_n(i) = D(4);
L_nn(i) = D(6);
end
P_n=str2double(P_n');
L_nn=str2double(L_nn');
endTo compare above stored numerical data with Analytical data we first need to calculate our analytical results.
ANALYTICAL CALCULATIONS
Average Velocity:-
`V_(avg) = (R_e mu)/(rho D)`
% Average Velcity
V_avg = (Re*mu)/(rho*D);Velocity at various radius(r) near fully developed flow region
`V_r =2V_(avg)(1-(r^2)/(R^2))`
% Velocity at every radius in a fully developed region
r = linspace(0,R,1000);
Vr=2*V_avg*(1-(r.^2/R^2));
% Velocity at center
V_max = 2*V_avg;Pressure drop across the length of pipe
`DeltaP = (32 mu L V_(avg))/D^2`
%% Pressure Drop
% Pressure Drop in across the length of pipe (Pa)
L_an = linspace(0,L,1000);
P_drop = (32*mu.*L_an.*V_avg)/D^2;
P_drop = fliplr(P_drop);Wall Shear stress
To calculate wall shear stress first calculate Darcy friction factor (f)
`f = (64)/(R_e)`
`Tau_w = f/4.(rho V_(avg)^2)/2`
%% Shear Stress
% Darcy Friction Factor
f = 64/Re;
% wall Shear stress
tau_w = (f/4)*((rho*V_avg^2)/2);Now, let us call the function which store our Numerical results and plot the data
%% Velocity and Pressure file reading
[L_n,V_n,L_nn,P_n] = velocity_pressure_file_reading();Velocity comparison of Numerical and analytical data
%% [1]-Plotting - Velocity comparison from analytical method and numerical method
figure(1)
set(gcf, 'Position', [100, 100, 800, 800])
plot(r,Vr,'r','linewidth',1)
hold on
plot(L_n,V_n,'--','color','k','linewidth',2.5)
xlabel('<---- Radius of pipe in m ---->','fontweight','bold','fontsize',14)
ylabel(' Velocity in m/s ------>','fontweight','bold','fontsize',14)
legend(' Analytical method',' Numerical method')
grid onResult Plot:-
To validate our results, we can compare our numerical results from finite volume method with analytical method, which depicted in above figure for velocity along the radius. We can notice, both velocity profile, for Numerical method and for Analytical method, are almost similar which validate our results.
After calculating maximum error in both velocity results we get to know that max error is 0.1774, which is very less; In other words negligible, as shown in below picture of command window.
Pressure comparison of Numerical and analytical data
%% [2]-Plotting - Pressure comparison from analytical method and numerical method
figure(2)
set(gcf, 'Position', [100, 100, 800, 800])
plot(L_an,P_drop,'r','linewidth',1)
hold on
plot(L_nn,P_n,'--','color','k','linewidth',2.5)
xlabel('<---- Length of pipe in m ---->','fontweight','bold','fontsize',14)
ylabel(' Pressure in Pa ------>','fontweight','bold','fontsize',14)
title('Pressure variation over the length of pipe','fontsize',18)
legend(' Analytical method',' Numerical method')Result Plot:-
Above line graph depicts results for the pressure variation over the length of pipe for both method Analytical method and Numerical method. and we can notice, we get almost similar result for pressure with both methods.
Conclusion:-
- By using Matlab, successfully written codes for OpenFoam blockMeshDict file to create fluid domain and meshing.
- Understood the Hagen-Poiseuille equation to calculate our analytical or exact results of laminar viscus flow through pipe.
- Using icoFoam we also calculate same results by Numerical method using finite volume method.
- Finally, results are validated successfully by comparing both analytical and Numerical results.
References:-
[1] FLUID MECHANICS by Prof. Dr. A. BULU - https://web.itu.edu.tr/~bulu/fluid_mechanics_files/lecture_notes_07.pdf
[2] Hagen-Poiseuille Equation - https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation
[3] Application of the Hagen-Poiseuille Equation to Fluid Feeding through Short Tubes -
Complete Function Program for blockMeshDict File Writing
function[f1] = blockMesh(R,th,L,nx,ny,nz)
%% Writing BlockMeshDict file for openFoam
f1 = fopen('blockMeshDict','w');
fprintf(f1,"/*--------------------------------*- C++ -*----------------------------------*\\ \n");
fprintf(f1,' ========= | \n');
fprintf(f1,' \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \n');
fprintf(f1,' \\\\ / O peration | Website: https://openfoam.org \n');
fprintf(f1,' \\\\ / A nd | Version: 8 \n');
fprintf(f1,' \\\\/ M anipulation |\n');
fprintf(f1,"\\*---------------------------------------------------------------------------*/ \n");
fprintf(f1,'FoamFile \n');
fprintf(f1,'{ \n');
fprintf(f1,' version 2.0;\n');
fprintf(f1,' format ascii;\n');
fprintf(f1,' class dictionary;\n');
fprintf(f1,' object blockMeshDict;\n');
fprintf(f1,'} \n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n');
fprintf(f1,'convertToMeters 1; \n\n');
fprintf(f1,'vertices\n(\n');
fprintf(f1,' (0 0 0)\n');
fprintf(f1,(' (%f %f 0)\n'),R*cosd(th/2),R*sind(th/2));
fprintf(f1,(' (%f %f %f)\n'),R*cosd(th/2),R*sind(th/2),L);
fprintf(f1,(' (0 0 %f)\n'),L);
fprintf(f1,(' (%f %f 0)\n'),R*cosd(th/2),-R*sind(th/2));
fprintf(f1,(' (%f %f %f)\n);\n\n'),R*cosd(th/2),-R*sind(th/2),L);
fprintf(f1,'blocks\n(\n');
fprintf(f1,('hex (0 4 1 0 3 5 2 3) (%d %d %d) simpleGrading (0.2 1 1)\n);\n\n'),nx,ny,nz);
fprintf(f1,'edges\n(\n');
fprintf(f1,' arc 4 1 (0 %f 0)\n',R);
fprintf(f1,' arc 5 2 (%f %f 0)\n);\n\n',L,R);
fprintf(f1,'boundary\n(\n');
fprintf(f1,' inlet\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 4 1 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' outlet\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (3 5 2 3)\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,' (0 1 2 3)\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,' (0 3 5 4)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n');
fprintf(f1,' pipewall\n');
fprintf(f1,' {\n');
fprintf(f1,' type wall;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (1 4 5 2)\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,'// ************************************************************************* // \n');
fclose(f1);
endComplete Function Program for velocity and Pressure File Reading
function [L_n,V_n,L_nn,P_n,Tau_n,L_nnn] = velocity_pressure_file_reading()
%% Reading Velocity File Extracted from paraView
f2 = fopen('fully_developed_flow_data.csv','r');
A=fgetl(f2);
for i = 1:1000
C = fgetl(f2);
D = strsplit(C,',');
V_n(i) = D(3);
L_n(i) = D(6);
end
V_n=str2double(V_n');
L_n=str2double(L_n');
%% Reading Pressure File Extracted from paraView
f3 = fopen('Pressure.csv','r');
A=fgetl(f3);
for i = 1:1000
C = fgetl(f3);
D = strsplit(C,',');
P_n(i) = D(4);
L_nn(i) = D(6);
end
P_n=str2double(P_n');
L_nn=str2double(L_nn');
endComplete main program
clear all
close all
clc
%% Input parameters
% Reynold's Number
Re = 2000;
% Total Length of pipe(m)
L = 2.5;
% Diameter of pipe in meters
D = 0.02;
% Wedge Angle in degrees
th = input('<0> Wedge Angle in degrees - ');
% Density (kg/m^3)
rho = 997;
% Dynamic Viscosity (Pa.sec)
mu = 8.90e-4;
% Radius of pipe(m)
R = D/2;
% Entrance Length
Le = 0.058*Re*D;
% Number of cells in x,y,z Directions
nx = 20;
ny = 1;
nz = 500;
%% BlockMeshDict file writting
[f1] = blockMesh(R,th,L,nx,ny,nz);
%% Calculations for Analytical methods
%% Velocity
% Analytical results
% Average Velcity
V_avg = (Re*mu)/(rho*D);
% Velocity at every radius in a fully developed region
r = linspace(0,R,1000);
Vr=2*V_avg*(1-(r.^2/R^2));
% Velocity at center
V_max = 2*V_avg;
%% Pressure Drop
% Pressure Drop in across the length of pipe
L_an = linspace(0,L,1000);
P_drop = (32*mu.*L_an.*V_avg)/D^2;
P_drop = fliplr(P_drop);
%% Wall Shear Stress
% Darcy Friction Factor
f = 64/Re;
% wall Shear stress
tau_w = (f/4)*((rho*V_avg^2)/2);
%% Velocity and Pressure file reading
[L_n,V_n,L_nn,P_n] = velocity_pressure_file_reading();
%% [1]-Plotting - Velocity comparison from analytical method and numerical method
figure(1)
set(gcf, 'Position', [100, 100, 800, 800])
plot(r,Vr,'r','linewidth',1)
hold on
plot(L_n,V_n,'--','color','k','linewidth',2.5)
xlabel('<---- Radius of pipe in m ---->','fontweight','bold','fontsize',14)
ylabel(' Velocity in m/s ------>','fontweight','bold','fontsize',14)
legend(' Analytical method',' Numerical method')
grid on
%% [2]-Plotting - Pressure comparison from analytical method and numerical method
figure(2)
set(gcf, 'Position', [100, 100, 800, 800])
plot(L_an,P_drop,'r','linewidth',1)
hold on
plot(L_nn,P_n,'--','color','k','linewidth',2.5)
xlabel('<---- Length of pipe in m ---->','fontweight','bold','fontsize',14)
ylabel(' Pressure in Pa ------>','fontweight','bold','fontsize',14)
title('Pressure variation over the length of pipe','fontsize',18)
legend(' Analytical method',' Numerical method')
grid on
%Error
E_max = max(max(abs(Vr-V_n)))