Week 11 - Simulation of Flow through a pipe in OpenFoam

Simulation of Flow through a Pipe in OpenFoam

Aim:

Simulate an axi-symmetric flow by applying the wedge boundary condition in OpenFoam and write a code in Matlab to automate the generation of blockMeshDict file.

• Validate Hydro-dynamic length with numerical result.
• Validate the fully developed flow velocity profile with its analytical profile.
• Validate maximum velocity and pressured drop for fully developed flow.
• Post-process Shear stress and validate wall shear stress for fully developed flow.

Introduction:

OpenFOAM (Open Field Operation and Manipulation) is an open-source software for computational fluid dynamics (CFD) modeling. This report provides a detailed analysis of the simulation of flow through a pipe in OpenFOAM. The simulation will be performed to study the flow characteristics and its behavior in a pipe.

Problem Description:

The problem considered in this simulation is the flow of fluid through a pipe of length L and diameter D. Here we will simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold’s number 2100.

The fluid is assumed to be incompressible and laminar with a constant viscosity. So we will perform steady state simulation and validate the results with Hagen-poiseuill’s equation.

The velocity profile at the inlet of the pipe is considered to be zerogradient, and the pressure at the outlet is considered to be atmospheric. The simulation is performed to study the following result.

• Validating Hydro-dynamic length with the numerical result.
• Validating the fully developed flow velocity profile with its analytical profile.
• Validating maximum velocity and pressured drop for fully developed flow.
• Post-process Shear stress and validate wall shear stress for fully developed flow.

Hagen-Poiseuil’s equation:

Hagen-Poiseuille equation is a mathematical expression that describes the laminar flow of an incompressible fluid in a cylindrical pipe. The equation is derived from the Navier-Stokes equations, which describe the motion of fluid. The Hagen-Poiseuille equation predicts the pressure drop and flow rate through a pipe as a function of fluid viscosity, pipe diameter, length, and the pressure difference between the ends of the pipe.

The equation is widely used in fluid dynamics and is particularly useful in the design of fluid transport systems, such as pipelines and medical devices. It is only valid for laminar flow conditions, where the fluid moves in smooth and parallel layers, and not for turbulent flow conditions, where the fluid has chaotic and irregular movement.

The Hagen-Poiseuille equation is expressed as:

Q = (π * r^4 * ΔP) / (8 * μ * L)

where:

• Q is the flow rate through the pipe (m^3/s)
• r is the radius of the pipe (m)
• ΔP is the pressure drop across the pipe (Pa)
• μ is the dynamic viscosity of the fluid (Pa * s)
• L is the length of the pipe (m)

Note that this equation assumes that the fluid is incompressible, laminar, and Newtonian, and that the velocity of the fluid is uniform along the radial direction and is a maximum at the center of the pipe. It also assumes that the flow is fully developed, meaning that the velocity profile reaches a steady state before the fluid reaches the outlet of the pipe.

Hydrodynamic Length:

Hydrodynamic length, also known as the entrance length, is the length of the pipe required for the flow of fluid to reach a fully developed velocity profile. In other words, it is the distance required for the fluid to adjust to the pipe geometry and reach a steady state.

The hydrodynamic length can be calculated using the following equation:

Lh = 0.05 * Re * D

where:

• Lh is the hydrodynamic length (m)
• Re is the Reynolds number, a dimensionless number that describes the flow regime (laminar, transitional or turbulent)
• D is the diameter of the pipe (m)

Shear Stress:

The equation to calculate the shear stress at the wall of a pipe is:

τ = μ * dv/dy

where:

• τ is the shear stress at the wall (Pa or N/m^2)
• μ is the dynamic viscosity of the fluid (Pa * s)
• dv/dy is the velocity gradient perpendicular to the wall, which represents the rate of change of velocity in the direction perpendicular to the wall.

To validate the wall shear stress for fully developed flow, one can compare the calculated shear stress with the theoretical value for fully developed flow in a circular pipe.

Methodology:

1. Geometry and Mesh Generation: The pipe geometry is created using OpenFOAM's built-in CAD tools. The mesh file is generated using Matlab. The mesh quality is checked using checkMesh.
2. Case Setup: The simulation case is set up in OpenFOAM by defining the governing equations, boundary conditions, and initial conditions. The flow is modeled using the Navier-Stokes equations.
3. Solver Selection: A suitable solver is selected from OpenFOAM's library of solvers to solve the governing equations. The icoFoam solver is used in this simulation.
4. Running the Simulation: The simulation is run using the OpenFOAM run script. The solution is monitored and checked for convergence.
5. Post-Processing: The results of the simulation are post-processed using ParaView, an open-source visualization tool. The velocity and pressure distribution along the pipe length are analyzed.

Input Values:

Diameter of the Pipe (d)= 0.04

Wedge angle (theta) = 4 (in degree)

Renold’s number(Re) = 2100

Density of the fluid(rho) = 997 kg/m^3

Dynamic Viscosity of water (mu) = 0.00089 Ns/m^2 or (Pa.s)

Kinematic Viscosity (nu) =(mu)/rho

Formulas:

Hydrodynamic Length(Lh) = 0.05 * Re *d = 4.2 m

Total Length of the pipe (L) = Lh + 0.7 = 4.9 m

Average Velocity(V_avg) = Re*(nu/d) = 0.0468 m/s

Maximum Velocity(V_max) = 2*V_avg = 0.0937 m/s

Pressure Drop (delatap) = (8*V_avg*mu*L)/r^2 = 4.081896 Pa

Radius of the Pipe (r) = 0.02 m

%creating blockMeshDict file for OpenFoam

d=0.04 %in meters, assumed

r=d/2;

theta=4

Re=2100

rho=997; % in kg/m^3

L_etr=0.06*Re*d;

L=L_etr+0.7;

mu=0.00089; % viscosity of water

nu=mu/rho; % kinetic viscosity

% boundary conditions

v_avg= Re*nu/d

v_max=2*v_avg;

deltap= (8*v_avg*mu*L)/r^2;

% calculating shear stress

R = linspace(0,(d/2),1000);

for i = 1:length(R)

tau(i) = 2*mu*v_max*(abs(R(i)))/r^2

end

plot(R,tau)

xlabel('Radius')

ylabel('Shear stress')

grid on

legend ('Shear stress');

h1='/*----------------------------------*- C++ -*------------------------------------*\\';

h2=' ========= | ';

h3=' \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox ';

h4=' \\\\ / O peration | Website: https://openfoam.org ';

h5=' \\\\ / A nd | Version: 7 ';

h6=' \\\\/ M anipulation | ';

h7='\\*----------------------------------------------------------------------------------*/';

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)]

f1=fopen('blockMeshDict.txt','w')

fprintf(f1,'%s\n',h1);

fprintf(f1,'%s\n',h2);

fprintf(f1,'%s\n',h3);

fprintf(f1,'%s\n',h4);

fprintf(f1,'%s\n',h5);

fprintf(f1,'%s\n',h6);

fprintf(f1,'%s\n',h7);

fprintf(f1, 'FoamFile\n');

fprintf(f1, '{\n');

fprintf(f1, ' format ascii;\n');

fprintf(f1, ' class dictionary;\n');

fprintf(f1, ' object blockMeshDict;\n');

fprintf(f1, '}\n');

fprintf(f1, '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n');

fprintf(f1, 'convertToMeters 1;\n\n')

fprintf(f1, 'vertices\n(\n');

fprintf(f1, '\t (%d %d %d)\n', v0);

fprintf(f1, '\t (%d %d %d)\n', v1);

fprintf(f1, '\t (%d %d %d)\n', v2);

fprintf(f1, '\t (%d %d %d)\n', v3);

fprintf(f1, '\t (%d %d %d)\n', v4);

fprintf(f1, '\t (%d %d %d)\n', v5);

fprintf(f1,');\n\n');

fprintf(f1,'blocks\n(\n');

fprintf(f1,'\t hex(0 1 2 0 3 4 5 3) (300 300 1) simpleGrading (1 0.1 0)\n);\n\n');

fprintf(f1,'edges\n(\n');

fprintf(f1,'\t arc 1 2 (0 %.6e 0)\n', r);

fprintf(f1,'\t arc 4 5 (%d %.6e 0)\n', L, r);

fprintf(f1,');\n\n');

fprintf(f1,['boundary\n(\n\t axis\n\t\t{\n\t\t type empty;\n\t\t faces\n\t\t (\n']);

fprintf(f1,'\t\t\t (0 3 3 0)\n\t\t );\n\t\t}\n');

fprintf(f1,['\t wall\n\t\t{\n\t\t type wall;\n\t\t faces\n\t\t (\n']);

fprintf(f1,'\t\t\t (2 5 4 1)\n\t\t );\n\t\t}\n');

fprintf(f1,['\t inlet\n\t\t{\n\t\t type patch;\n\t\t faces\n\t\t (\n']);

fprintf(f1,'\t\t\t (0 0 2 1)\n\t\t );\n\t\t}\n');

fprintf(f1,['\t outlet\n\t\t{\n\t\t type patch;\n\t\t faces\n\t\t (\n']);

fprintf(f1,'\t\t\t (3 4 5 3)\n\t\t );\n\t\t}\n');

fprintf(f1,['\t wedgeFront\n\t\t{\n\t\t type wedge;\n\t\t faces\n\t\t (\n']);

fprintf(f1,'\t\t\t (0 3 5 2)\n\t\t );\n\t\t}\n');

fprintf(f1,['\t wedgeBack\n\t\t{\n\t\t type wedge;\n\t\t faces\n\t\t (\n']);

fprintf(f1,'\t\t\t (0 1 4 3)\n\t\t );\n\t\t}\n);\n');

fprintf(f1, 'mergePatchPairs')

fprintf(f1,'\n(\n);\n')

fprintf(f1,'// ************************************************************************* //')

% % fprintf(fileID, '\t arc 1 2 (0 2.000000e-03 0)\n');

% fprintf(fileID, '\t arc 1 2 (0 %.6e 0)\n', r);

fclose(f1);

Shear Stress Plot along the Radius of the Pipe

blockMeshDict file

/*--------------------------------*- C++ -*----------------------------------*\

========= |

\\ / F ield | OpenFOAM: The Open Source CFD Toolbox

\\ / O peration | Website: https://openfoam.org

\\ / A nd | Version: 9

\\/ M anipulation |

\*---------------------------------------------------------------------------*/

FoamFile

{

format ascii;

class dictionary;

object blockMeshDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices

(

(0 0 0)

(0 1.998782e-02 -6.979899e-04)

(0 1.998782e-02 6.979899e-04)

(5.740000e+00 0 0)

(5.740000e+00 1.998782e-02 -6.979899e-04)

(5.740000e+00 1.998782e-02 6.979899e-04)

);

blocks

(

hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (1 1 1)

);

edges

(

arc 1 2 (0 2.000000e-02 0)

arc 4 5 (5.740000e+00 2.000000e-02 0)

);

boundary

(

axis

{

type empty;

faces

(

(0 3 3 0)

);

}

fixedWalls

{

type wall;

faces

(

(2 1 4 5)

);

}

front

{

type wedge;

faces

(

(0 3 5 2)

);

}

back

{

type wedge;

faces

(

(0 1 4 3)

);

}

inlet

{

type patch;

faces

(

(0 1 2 0)

);

}

outlet

{

type patch;

faces

(

(3 5 4 3)

);

}

);

mergePatchPairs

(

);

// ************************************************************************* //

Pressure File

/*--------------------------------*- 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

{

inlet

{

type zeroGradient;

}

outlet

{

type fixedValue;

value uniform 1;

}

fixedWalls

{

type zeroGradient;

}

front

{

type wedge;

}

back

{

type wedge;

}

axis

{

type empty;

}

}

// ************************************************************************* //

Velocity File

/*--------------------------------*- 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

{

inlet

{

type fixedValue;

value uniform (0.0937312 0 0);

}

outlet

{

type zeroGradient;

}

fixedWalls

{

type noSlip;

}

front

{

type wedge;

}

back

{

type wedge;

}

axis

{

type empty;

}

}

// ************************************************************************* //

ControlDict file

/*--------------------------------*- 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.01;

writeControl timeStep;

writeInterval 200;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;

// ************************************************************************* //

Physical Properties File

/*--------------------------------*- 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] 1e-06;

// ************************************************************************* //

Matlb Code to find the Velocity Profile

clear all

close all

clc

D = 0.04; % pipe diameter [m]

U = 0.0937312; % fluid velocity [m/s]

mu = 0.00089; % fluid dynamic viscosity

R = D/2; % pipe radius [m]

rho = 997 % Density of water

nu =mu/rho % Kinematic Viscosity

Re = 2100; % Reynolds number

d = linspace(0,D,50); % radial coordinate

%

% % % calculate velocity

u = zeros(size(d));

for i = 1:length(d)

u(i) = 4.081896*(22/7)*(1/16)*d(i)^4/(128*mu*4.9)

end

plot(u, d)

title('Velocity Profile')

xlabel('Velocity')

ylabel('Diametre of the Pipe')

Velocity Profile taken at the length of 4.2 m from inlet

X axis- Radius of the Pipe

Y_axis- Velocity of the Pipe

So the maximum velocity reached at the centre of the pipe is approximately 0.017622 m/s

Mirror Results in ParaView:

Step 1: Choose the Filter

Go to Filters, Alphabetical and Search. If the filter is disabled/not visible it can be accessed using Filters and search and then typing name Reflect in the search window.

Step 2: Select the Symmetry of the Plane

We want to reflect about the central axis which is the bottom edge of the wedge. It is important to make note of the coordinate axes before selecting the symmetry plane. In this case, the bottom edge of the wedge is aligned to the central axis of the flow and since that is axis of symmetry, we want to choose the Y-plane that passes through the symmetry axis.

And as the Y- coordinate is minimum for the bottom edge of the wedge, we choose the symmetry plane as Y Min.

Step 3: Click on the Apply button

By clicking on the apply button, we can see the complete 2D solution for our problem. The solution has been mirrored about the central axis. Now, we can use the plot over line tool to visualize the velocity profile across the complete diameter. The velocity profile appears as shown below.

Velocity Profile taken along the Diameter of the pipe

X-axis- Diameter of the Pipe

Y-axis – Velocity of the Pipe

Creating Shear stress plot inside Paraview

Step 1 : Computing the Magnitude of vector field U using Calculator Tool

1. Make sure that the common toolbar under view tab is enabled.

2. Click on the Calculator tool that can be found in the extreme left of the commom toolbar.

3. A calculator will appear in the properties box on the left pane of the Paraview window.
4. In the calculator, we want to compute the magnitude of the U vector, so type mag(U) inside the calculation box . Vectors button can also be used to select U. The result of the computation is stored inside that name so that it can be accessed later. In the above picture, mag_ U is the result array name.

5. On clicking on the apply button, the magnitude of the velocity vector U is computed and applied in the Field Variable selection box.

Step 2: Computing the Gradients for U_magnitude using the ‘Gradients ‘ tool.

1. Now we need to compute the gradients of U_magnitude inside the entire domain, since velocity gradient is what determines the shear stress developed inside the fluid.
2. This can be done using the “ Gradients “ tool inside Filter.

3. The properties windows shows as follows.

4. In the Scalar Array box, select the array name where the U_magnitude values were previously stored. Since, we stored the results inside mag_U, that quantity is selected as the Scalar Array in the above picture.
5. Type the result array name as desired. That is where the gradients of U_magnitude values will be stored, so that we can access them in the future. In the above picture grad_mag_U is the array name chosen for that purpose.
6. Clicking on the Apply button will create a new quantity of the specified name inside the Field Variable Selection box.

Step 3: Computing the Shear-stress quantity

1. In the previous step, we computed the velocity gradients data. Now we need to use the formula for shear stress to calculate it.

The equation to calculate the shear stress at the wall of a pipe is:

τ = μ * dv/dy

where:

• τ is the shear stress at the wall (Pa or N/m^2)
• μ is the dynamic viscosity of the fluid (Pa * s)
• dv/dy is the velocity gradient perpendicular to the wall, which represents the rate of change of velocity in the direction perpendicular to the wall.

2. Therefore, we can obtain shear stress by multiplying the velocity gradient with the dynamic viscosity of the fluid under study.
3. Assuming water is our working fluid, and that is at the room temperature 25 C, the dynamic viscosity is 0.00089 Ns/m^2.
4. Now click on the calculator, and select our velocity gradient vector array by clicking on the vectors button on the bottom-right and then multiply that by dynamic viscosity as shown in the picture below.

We have used the variable tau to store shear stress date.

5. Clicking on the apply button will create a new quantity inside the Field Variable selection box as shown below.

Step-4: Plotting Shear-stress

We can create the Plot Over Line tool located in the common toolbar.

Results

1. Validating Hydrodynamic Length with the numerical result.

Hydrodynamic length = 0.05 * Re * d = 0.05 *2100 * 0.04 => 4.2 m

2. And the fully developed velocity profile obtained at the distance of 4.2 m from the inlet in ParaView. Maximum Velocity obtained at the centre of the Pipe is approximately 17622 m/s

Maximum Velocity(V_max) = 2*V_avg = 0.187 m/s

X- axis – Diameter of the Pipe

Y – axis – Velocity along the diameter

Velocity Profile- Matlab

3. Pressure Drop for Fully Developed Flow

Diameter of the Pipe (d)= 0.04

Renold’s number(Re) = 2100

Density of the fluid(rho) = 997 kg/m^3

Dynamic Viscosity of water (mu) = 0.00089 Ns/m^2 or (Pa.s)

Kinematic Viscosity (nu) =(mu)/rho

Formulas:

Hydrodynamic Length(Lh) = 0.05 * Re *d = 4.2 m

Total Length of the pipe (L) = Lh + 0.7 = 4.9 m

Average Velocity(V_avg) = (Re*nu)/(density*diametre) = 0.093m/s

Maximum Velocity(V_max) = 2*V_avg = 0.187 m/s

Pressure Drop (delatap) = (8*V_avg*mu*L)/r^2 = 4.081896 Pa

Radius of the Pipe (r) = 0.02 m

Pressure is constant along the cross section for a flow through a circular pipe due to the symmetry of the pipe. The pressure exerted by the fluid on the walls of the pipe is normal to the wall at every point, which implies that there is no variation in the pressure in the radial direction of the pipe. Therefore, the pressure remains constant along any cross section of the pipe.

This result is a consequence of the continuity equation, which states that the mass flow rate of a fluid through a pipe is conserved. In order for the mass flow rate to be conserved, the velocity of the fluid must vary inversely with the cross-sectional area of the pipe. This inverse variation of velocity and area is what leads to the constancy of pressure along the cross section of the pipe.

This phenomenon is known as Poiseuille's law, and it is widely used to model fluid flow through circular pipes.

Shear stress

X-axis- Radius of the Pipe

Y-axis – Shear stress

Shear stress increases along the Radial direction and reaches maximum at the wall.

The wall shear stress can be calculated using the following formula:

τw = μ(dv/dy),

where τw is the wall shear stress, μ is the dynamic viscosity of the fluid, and dv/dy is the velocity gradient at the wall of the surface. The velocity gradient can be calculated using the following formula:

dv/dy = (u(x,y+Δy) - u(x,y)) / Δy,

where u(x,y) is the velocity of the fluid at a point (x,y) on the surface, and Δy is the distance between the point and the wall.

To calculate the wall shear stress in practice, we would typically need to obtain velocity data for the fluid near the wall using experimental techniques such as particle image velocimetry (PIV) or hot wire anemometry. Once the velocity data is obtained, we can calculate the velocity gradient at the wall and then use the formula above to calculate the wall shear stress. In cases where the flow is laminar and the velocity profile is known, the wall shear stress can also be calculated using the Hagen-Poiseuille equation. As we are simulating for laminar flow the wall shear stress would be similar to the shear stress calculated in the above graph.