Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective
The objective of this project is to numerically simulate the Hagen Poiseuille flow through a pipe in OpenFOAM by accomplishing the following objectives.
- Writing a program in MATLAB to automate the generation of the computational mesh for any wedge angle and grading schemes. Wedge angle should be less than 5 degrees.
- For a baseline mesh, show that the velocity profile matches with the Hagen Poiseuille Equation.
- Show that the velocity profile is fully developed.
- Post process shear stress developed using Paraview
Project Methodology
This project aims to simulate the Hagen-Poiseuille Flow through a pipe of circular cross section, by considering all of the assumptions. Since the flow is axisymmetric in nature (symmetric along the axis of the cylindrical pipe), a reduced model, i.e. only a section of the pipe (wedge) is modelled through which fluid flow is then simulated. This is done to reduce complexity of the problem, and the computational resources required extensively. The wedge is created by considering a specific wedge angle. A MATLAB Script is created to automatically generate the computational mesh by accepting the wedge angle as the input. The blockMeshDict file is created from the output of the MATLAB Script. The icoFoam Solver which is used to solve the incompressible, laminar, Navier Stokes Equation is used to simulate the flow through the pipe. The results are then post-processed using ParaView, to observe the velocity profile for the fully developed flow, and at the Hydrodynamic Entrance Length of the pipe. The variation of Shear Stress developed along the length of the pipe (axial direction) is also plotted. The data from the CFD simulation is validated by comparing against the analytical solutions obtained from the Hagen-Poiseuille Equations to ensure the fidelity of results.
Hagen-Poiseuille Flow
The Hagen-Poiseuille Flow is a type of fluid flow based on the Hagen-Poiseuille Law, which 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. Hagen-Poiseuille Flow or simply Poiseuille Flow is the guiding principle behind a variety of applications like, air flow through the lung alveoli, fluid flow through a drinking straw or through a hypodermic needle. In essence, Poiseuille's Equation is used to describe the pressure drop due to the viscosity of the fluid.
The assumptions of the HP Equation is that the fluid is incompressible and Newtonian, and the flow is laminar through a pipe of constant cross section, with a length that is sufficiently larger than the diameter of the pipe, and that there is no acceleration of the fluid through the pipe.
Consider a Cylindrical pipe of circular cross section as shown above. Let the diameter of the pipe be D and the Length of the Pipe be L. Assume a fluid flowing through the pipe with an inlet velocity of v1 and let pressure at the outlet of the pipe be p2. For this project, the diameter of the pipe is assumed to be 0.02m or 20mm, and the length of the pipe is taken as 3m. The length of the pipe has been decided based on the Hydrodynamic Entrance Length (Le) obtained, which will be explained in upcoming sections. Let water be the working fluid flowing through the pipe. Thus, the characteristic properties of water such as density (`rho`) and viscosity (`mu`) are taken into account to simulate the fluid flow. Since Hagen-Poiseuille Flow is a laminar flow, a Reynold's Number pertaining to the Laminar Flow regime (Re = 2100) is considered to simulate the fluid flow and calculate the fluid velocity. In order to minimse the computational resources required for simulating the problem, a section of the pipe or a wedge with a specific wedge angle (less than 5 degrees) is considered, and the symmetry boundary conditions are applied. Considering the aforementioned input values, a MATLAB Script is created to automate the generation of the computational mesh using blockMeshDict to be used in OpenFOAM.
Input Parameters
- Diameter of the pipe (D) = `0.02m`
- Wedge Angle (`theta`) = 4 (must be less than 5)
- Length of the Pipe (L) = `3 m`
- Density of Fluid (Water) (`rho`)=`997 kg//m^3`
- Viscosity of Water at 298K (`mu`)= 8.9x`10^-4` `Pa.s`
- Reynold's Number for Laminar Flow (`Re`) = 2100
Description of Flow Physics
As the working fluid (water in this case) enters the pipe with a constant inlet velocity, the fluid comes in contact with the pipe walls. Because of the no-slip boundary condition, the fluid particles in the layer in contact with the wall come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually due to friciton between the layers (fluid viscosity). To make up for this velocity reduction, the velocity of the fluid at the mid-section of the pipe has to increase so that the mass flow rate through the pipe remains constant. As a result of this, a velocity gradient develops along the axial direction or length of the pipe.
The region of flow in which the effects of viscous shearing forces caused by the fluid viscosity are significantly felt is called the velocity boundary layer or simply the boundary layer. The boundary surface divides the flow in a pipe into two regions: the boundary layer region, in which the viscous effects and velocity changes are significant, and the irrotational (core) flow region (or) the inviscid flow region, in which the friction effects are negligible and the velocity essentially remains constant in the radial direction.
The thickness of the boundary layer gradually increases in the flow direction, until it reaches the mid-section of the pipe, at which point it fills the entire pipe, and the velocity profile becomes fully developed. The region of the pipe from the pipe inlet, to the point at which the flow becomes fully developed is called as the hydrodynamic entrance region, and the length of this region is called the hydrodynamic entrance length (`L_e`). The flow in this region is called hydrodynamically developing flow since this is the region where the velocity profile develops and keeps changing. The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged across the length of the pipe, is called the hydrodynamically fully developed region. The flow in this region is said to be a fully developed flow and remains constant across the pipe. The velocity profile in the fully developed region is parabolic for laminar flow as shown in the figure below.
Hydrodynamically fully developed flow: `(delU_(r,x))/(delx) = 0` `->` U = U(r)
Hydrodynamic Entrance Length: For Laminar Flow, the expression for calculating the Hydrodynamic Entrance Length is given as follows,
`L_e = 0.05*Re*D`
We know that for Newtonian Fluids, the shear stress is directly proportional to the velocity gradient. Thus, the shear stress at the pipe wall (`tau_w`) is also related to the slope of the velocity profile. Since the velocity profile remains unchanged in the fully developed region, the wall shear stress also remains constant in this region. Also, the wall shear stress is highest at the inlet where the thickness of the boundary layer is the smallest, and decreases gradually to the fully developed value as shown in the figure below.
Other commonly used forms of the equation are described below.
`Deltap = (8muLV_(avg))/(R^2)` and `Deltap = (32muLV_(avg))/(D^2)`
where,
`Deltap` = `P_1-P_2` Difference in pressures between the inlet and outlet of the pipe, and all the other quantities are as described earlier.
From the above equation, it is clear that the pressure drop across the pipe is directly proportional to the viscosity of the fluid flowing through the pipe. If there were no viscosity (i.e. no friction), the pressure loss would be zero. Hence, the pressure drop in the above case is entirely due to the fluid viscosity.
Deriving Analytical Solutions for Flow parameters from HP Equation
In order to validate the CFD Results obtained, and to ensure their fidelity, they should be compared with analytically obtained values. These analytical values are obtained through different analytical equations, describing different flow parameters.
- Reynold's Number `=>``Re = (rhoV_(avg)D)/(mu)` `=>``2100 = (997xxV_(avg)xx0.02)/(8.92xx10^-4` `therefore``V_(avg) = 0.0937 m//s`
- Maximum Velocity at Centre line `=> V_max = 2.V_(avg)``=> 2xx0.0937 = 0.1875m//s`
- Hydrodynamic Entrance Length `=> L_h = 0.06xxDxxRe = 0.06xx0.02xx2100 = 2.52m`
- Pressure drop across the pipe from HP Equation `=> Delta p = (32muLV_(avg))/(D^2) = (32xx8.9xx10^-4xx3xx0.0937)/(0.02^2) = 20.012 Pa`
- Kinematic Pressure Drop `==> Deltap_k = (Deltap)/(rho) = 20.021/997 = 0.0201Pa`
- Velocity Profile at any location in the fully developed region `=> U(r) = 2*V_(avg)*(1-(r^2)/(R^2))`, where 'r' represents the radial distance from the centre of the pipe
The analytical values and equations described above will be used for comparing against the CFD results to ensure their validity. A MATLAB Script is written to evaluate these values for the flow parameters, and to generate the computational mesh automatically through the blockMeshDict file for usage in OpenFOAM.
The Geometry under consideration for fluid flow simulation is a reduced model that has been designed to ensure minimal use of computational resources. This geometry is shown in the figure below. Since the geometry of the pipe is axisymmetric along the x axis, a section of the pipe is modelled and the boundary conditions are provided accordingly. Also, since it is essential to capture the boundary layer developed at the topwall of the wedge, a grading factor is given to the mesh to obtain a locally refined mesh. This locally refined mesh is used to sufficiently model and capture the boundary layer developed during the fluid flow, and to determine the velocity profile developed during the flow.
MATLAB Code
The MATLAB Script file is written by assuming initial values for flow parameters as described earlier. Using these values, the analytical solutions for flow parameters is obtained by solving the Hagen-Poiseuille Equation. Next, the code for automatically generating the computational mesh is written by considering a specific wedge angle for the pipe. The MATLAB Code for Mesh Generation and Evaluation of Analytical Solutions is given below.
clear all
close all
clc
%BlockMesh file for simulating flow through a pipe
Re = 2100; %Reynolds Number for Laminar Flow
rho = 997; %Density of Water
mu = 0.89*1e-3; %Dynamic viscosity of water
D = 0.02; %Diameter of the pipe
R = D/2; %radius of the pipe
L_e = 0.06*D*Re; %Hydrodynamic Entrance Length of Pipe
L = 3; %Total Length of the pipe
v_avg = (Re*mu)/(rho*D); %Average Velocity from Reynolds Number
v_max = 2*v_avg; %Maximum Velocity at the centre of pipe
delta_p = (32*mu*L*v_avg)/(D^2); %pressure drop from Hagen-Poiseulle Equation
kin_p = (delta_p)/(rho);
theta = input('Enter the wedge angle = ');
%Mesh parameters
%No. of Mesh points
mesh_x = 200;
mesh_y = 20;
mesh_z = 1;
%Mesh grading factors
grade_x = 1;
grade_y = 0.2;
grade_z = 1;
%Determining the vertices of wedge
v0 = [0,0,0];
v1 = [0,(R*cosd(theta/2)),(R*sind(theta/2))];
v2 = [0, (R*cosd(theta/2)), (-R*sind(theta/2))];
v3 = [L,0,0];
v4 = [L,(R*cosd(theta/2)),(R*sind(theta/2))];
v5 = [L, (R*cosd(theta/2)), (-R*sind(theta/2))];
%Parsing BlockMeshDict file
f1 = fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n========= |');
fprintf(f1,'\n\\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1, '\n\\ / O peration | Website: https://openfoam.org');
fprintf(f1, '\n\\ / A nd | Version: 7');
fprintf(f1, '\n\\/ M anipulation |');
fprintf(f1, '%s','\n\*---------------------------------------------------------------------------*/');
fprintf(f1, '\nFoamFile');
fprintf(f1, '\n{');
fprintf(f1,'\nversion 2.0');
fprintf(f1, '\nformat ascii');
fprintf(f1, '\nclass dictionary');
fprintf(f1,' object blockMeshDict');
fprintf(f1,'}');
fprintf(f1,'\n%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\nconvertToMeters 1;');
fprintf(f1,'\n\nvertices');
fprintf(f1,'\n(');
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point0',v0(1),v0(2),v0(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point1',v1(1),v1(2),v1(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point2',v2(1),v2(2),v2(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point3',v3(1),v3(2),v3(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point4',v4(1),v4(2),v4(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point5',v5(1),v5(2),v5(3));
fprintf(f1,'\n);');
fprintf(f1,'\n\nblocks');
fprintf(f1,'\n(');
fprintf(f1,'\n\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (%f %f %f)',mesh_x,mesh_y,mesh_z,grade_x,grade_y,grade_z);
fprintf(f1,');');
fprintf(f1,'\n\nedges\n');
fprintf(f1,'(');
fprintf(f1,'\n\t arc 1 2 (%f %f %f)',0,R,0);
fprintf(f1,'\n\t arc 4 5 (%f %f %f)',L,R,0);
fprintf(f1,'\n);');
fprintf(f1,'\n\nboundary');
fprintf(f1,'\n(');
fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\ttopwall\n\t{');
fprintf(f1,'\n\t\ttype wall;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(1 4 5 2)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tfront\n\t{');
fprintf(f1,'\n\t\ttype wedge;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 3 4 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tback\n\t{');
fprintf(f1,'\n\t\ttype wedge;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 2 5 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\taxis\n\t{');
fprintf(f1,'\n\t\ttype empty;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n);');
fprintf(f1,'\n\nmergePatchPairs');
fprintf(f1,'\n(');
fprintf(f1,'\n);');
fprintf(f1,'\n\n%s','// ************************************************************************* //');The code used above generates a blockMeshDict text file that is used by OpenFOAM to create the geometry of the pipe (pipe section as a wedge) with the computational mesh. The blockMeshDict file generated by the above code is attached below.
BlockMeshDict file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000) //point0
(0.000000 0.009994 0.000349) //point1
(0.000000 0.009994 -0.000349) //point2
(3.000000 0.000000 0.000000) //point3
(3.000000 0.009994 0.000349) //point4
(3.000000 0.009994 -0.000349) //point5
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1.000000 0.200000 1.000000));
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (3.000000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
topwall
{
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 command blockMesh is used to execute the blockMeshDict file and generate the computational mesh. Next, the initial and boundary conditions for the flow parameters (Velocity and Pressure), the time-step value, and the transport properties (i.e. fluid viscosity) are to be defined. The generated computational mesh using blockMeshDict is shown below.
The figures shown above depict the geometry and the mesh generated from the blockMeshDict file, from different view angles.
Specifying Initial Conditions
The initial conditions for pressure and velocity (for the timestep 0) are to be described in the 0 directory in OpenFOAM. The dictionary files for Pressure and Velocity Initial conditions are given below.
Initial Velocity
The average velocity calculated from the MATLAB Program is given as the inlet velocity,
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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.0937 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
topwall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //Initial Pressure
The outlet pressure is given as 0, and the inlet pressure is given as zeroGradient, to calculate the pressure drop across the pipe
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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;
}
topwall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //ControlDict file - Specifying Time Step Value
After the blockMeshDict file and initial conditions have been defined, the time step value should be defined to perform the simulation. The time-step value for the simulation is determined by the CFL Criterion, given below.
`CFL = (U_max*dt)/(dx)`, where CFL Represents the CFL Number or Courant Number, 'dt' represents the time step size, and 'dx' represents the mesh sizing.
To ensure stability of the numerical scheme used, Courant Number should be less than 1 (`CFL <1`)
Using the above expression and considering the maximum velocity, and mesh sizing for the problem, a time step size of 0.01s is given to simulate the flow. Also, in order to ensure that the fluid flow completely covers the pipe and attains a steady-state value, a sufficient amount of time must be provided. Keeping this in mind, an end time value of 30s is given to simulate the flow. These values are given as inputs to the controlDict dictionary file in OpenFOAM.
controlDict file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Simulation of Fluid flow using Solver
Once the physical problem has been set up with the appropriate initial and boundary conditions, the problem setup should be solved using the solvers provided in OpenFOAM. The solver used for simulating the flow is selected depending upon the type of problem to be solved. In this case, since the Hagen-Poiseuille Flow is the flow to be simulated, the characteristics of the Hagen-Poiseuille Flow should align with the solver characteristics. Since Hagen-Poiseuille Flow is an incompressible, Newtonian, and Laminar Flow, the incompressible solver icoFoam is used to simulate the flow phenomena. icoFoam solves the incompressible laminar Navier-Stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions. The icoFOAM code can take mesh non-orthogonality into account with successive non-orthogonality iterations. The number of PISO corrections and non-orthogonality corrections are controlled through user input.
The solver is invoked using the keyword icoFoam in the command terminal. The next step after the generation of results by the solver, is to post-process the obtained results and analyze them to obtain useful data.
Post-Processing of Results
The post-processing of results is done using ParaView, which is an open-source application for data analysis and visualization, provided as part of OpenFOAM. ParaView is invoked by using the command paraFoam in the command terminal. The main flow parameters to be analyzed are the velocity profile of the flow along the length of the pipe especially in the fully developed region, the pressure drop across the length of the pipe, the wall shear stress distribution along the axial and radial distances of the pipe.
Velocity Distribution along the Pipe
Velocity Distribution along the pipe (Front View)
Velocity Distribution at Pipe Inlet - Start of Boundary Layer Formation
Velocity Distribution at Pipe Outlet -Fully Developed Flow Region and Boundary Layer
Velocity Distribution at Pipe Outlet - Isometric View
Velocity Distribution at Hydrodynamic Entrance Region
Pressure Distribution along the Pipe
Pressure Distribution at Pipe Inlet
Pressure Distribution at Pipe Outlet
Plot of Velocity Distribution along the Radial Distance at different points along the Pipe Length
Velocity Profile at x=0m (Pipe Inlet)
Velocity Profile at x=0.2m
Velocity Profile at x=0.4m
Velocity Profile at x=0.8m
Velocity Profile at x=1m
Velocity Profile at x=2m
Velocity Profile at x=2.5m
Velocity Profile at x=2.8m
Velocity Profile at x=3m (Pipe Outlet)
- From the above plots, it is clear that as the fluid enters the pipe with a constant velocity (`V_(avg) = 0.0937m//s`) (x=0m) and progresses through the pipe, the velocity profile along the radial distance keeps changing.
- As the fluid approaches the hydrodynamically fully developed region, the velocity profile becomes progressively more parabolic, and the maximum velocity of the fluid also keeps increasing. (x=0.2,0.4,0.8,1,2)
- The hydrodynamic entrance length of the flow is `L_e = 2.52m`. After the fluid crosses this length, the flow becomes fully developed and remains constant throughout the remaining length of the pipe.
- Thus, at the point x=2.8m, and x=3m, the flow is fully developed, and the velocity profile is completely parabolic, as can be seen in the last two plots shown above. The maximum velocity obtained from the plots in the fully developed region is given as, `V_(max) = 0.1855`
Plot of Pressure Distribution along the Pipe Length
- The above plot shows the pressure drop across the length of the pipe. The pressure at outlet was given as 0 Pa, and the pressure at inlet from the plot is shown to be 0.0256 Pa.
- Thus, the pressure drop across the pipe from CFD Data, can be calculated as, `Delta P = P1 - P2 = 0.0256-0=0.0256Pa`
Plot of Wall Shear Stress Distribution along the Radial Distance at x=3m
- The above plot shows the variation of shear stress along the radial distance (y axis) of the pipe at x=3m (fully developed flow region)
- In the hydrodynamically fully developed flow regime, the wall shear stress varies linearly from the centre of the pipe (i.e. r=0) to the pipe walls (r=0.01m)
- The wall shear stress value reaches a maximum approximately near the pipe walls.
- This distribution is in agreement with the analytical results, which states that the wall shear stress produced during the flow is directly proportional to the velocity gradient along the y axis. (`tau = mu(delu)/(dely)`)
- The maximum shear stress obtained from the plot is given as, `tau_(max) = 0.032Pa`.
Comparison of Analytical Vs. CFD Results - Velocity Profile along Radial Distance
- The above plot shows the comparison of CFD Results obtained from the simulation, with the analytical results calculated from the MATLAB Code.
- From the plot, it is clear that the CFD Results obtained from the flow simulation is in good agreement with the results obtained from the analytical equations, thus ensuring the fidelity of the results and the simulation.
- The minor errors that occur in the plots can be attributed to discretization and approximation errors that occur during the solution process. Thus, the CFD Results have been validated with the analytical solutions.
| Flow Parameter | CFD Value | Analytical Value | % Error |
| Maximum Velocity along Radial Distance | 0.1855 m/s | 0.1875 m/s | 1.06% |
From the table above, it is clear that the obtained CFD Solution is in close agreement with the analytical solution. Hence, the validity of the reduced model (i.e. wedge) has been established.
Inference and Conclusion
- The Hagen-Poiseuille flow was computationally simulated successfully and the results are found to be in good agreement with the analytical values derived from the Hagen-Poiseuille Equations.
- The Hagen-Poiseuille Flow was simulated by using a reduced model in the form of a wedge with a specific angle, instead of a fully cylindrical pipe with circular cross-section. This was done to reduce the amount of computational time and resources required to simulate the flow.
- The CFD results obtained for the reduced model were found to be in good agreement with the analytical results, hence confirming the validity of the reduced model as a replacement for the full pipe geometry.
- The wall shear stress distribution along the radial distance of the pipe was found and the maximum wall shear stress obtained from the CFD Results was similar to the theoretical maximum, calculated from analytical equations.
- The pressure drop across the pipe due to the viscosity of the fluid, obtained from CFD Results, was also found to be in good agreement with the result obtained from the Hagen-Poiseuille Equation.
- Hence, it can be concluded that the CFD simulation with correct assumptions and approximations yields results similar to and in close agreement with results obtained from close-form analytical equations.