Week 11 - Simulation of Flow through a pipe in OpenFoam

Aim:-

To simulate flow through a pipe using OpenFoam and validate the results.

 

Objective:-

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

 

Methodology:-

In this project, we are considering laminar flow of Incompressible fluid of constant properties through a straight circular pipe with fully developed flow. Generally, simulating an entire pipe in Openfoam takes a lot resources and time. Since, we are considering axi-symmetric flow, it will be easier to capture the physics. To achieve this, we apply a wedge boundary condition to the problem. For the pre-processing part, we are writing a code in Matlab to obtain the BlockMeshDict file which contains the Mesh parameters and properties. This code is written is such a way where just giving the input parameters of the problem yields us with a functional BlockMeshDict file. After this step, other input files are initialized like U,p, controlDict which contain the Boundary conditions for Velocity, pressure and for setting up the simulation. We will go in detail when we come across this step. After the simulation in done, we will be post-processing the results using Paraview and Matlab validating our results with the analytical solution.

 

Physics of the problem:-

 If you have ever observed smoke plume from a candle flame, you would have seen a rising smooth plume of smoke for some centimeters then starts to flutter in all directions as it rises above. Similarly, when fluid is observed in a pipe, we see a similar trend where the fluid is streamlined at low velocities and becomes chaotic as the velocity hits and increases from a critical value. The first flow regime is known as Laminar, as it characterises a smooth streamlines and highly ordered flow and the second flow regime is known as Turbulrent, as it is charecterised by a very chaotic and uneven streamlines and highly disordered flow. But, the transition from laminar flow to turbulent flow dosen't occur immedeately but rather slowly where the laminar flow region after some distance has some flow fluctuations and slowly becomes turbulent. This transition depends on various factors like geometry, roughness of surface, etc. The flow regime mainly depends on the ratio of interial to viscous forces. This ratio is known as Reynolds Number (Re) which is a dimensionless number. It is given as follows

$Re=\frac{\text{Inertial Forces}}{\text{Viscous forces}}$

$Re=\frac{V_{avg}\cdot D}{\nu }$

$Re=\frac{\rho \cdot V_{avg}\cdot D}{\mu }$

where $V_{avg}$ is the average flow velocity (m/s)

$D$ is the characteristic length of the geometry (Diameter in this case) and

$\nu =\frac{\mu }{\rho }$ is the kinematic viscosity (m^2/s) of the fluid with $\mu$ being Dynamic viscosity.

 

The Reynolds number at which the flow becomes turbulent is known as Critical Reynolds number $\left(R{e}_{cr}\right)$. As we are considering flow in a circular pipe, the generally accepted value of Critical Reynolds number is 2300.

Entrance Length:

Consider a fluid entering a pipe with a uniform velocity. The fluid particles at the wall come to a complete stop due to the high shear stress present, this is known as no-slip condition. This condition affects the adjacent fluid layers to slow down moderately as a result of friction. As this process of velocity reduction moves onto the mid-section or at the axis of the pipe, the velocity increases to keep the mass flow rate constant throughout the pipe.

The region at which the effects of the viscous shearing forces are felt due to the viscosity of the fluid is known as velocity boundary layer or just boundary layer. This hypothetical surface line divides the flow into two regions namely boundary region where the viscous effects and the velocity changes ae significant and Irrotational (core) flow region where the velocity remains constant in the radial direction due to the frictional effects being negligible. The thickness of the boundary layer keeps increasing as it moves towards the center of the pipe and thus covers the entire pipe. The region from the pipe inlet to the fully developed flow is called the Hydrodynamic entrance region and the length of this region is known as the Hydrodynamic entry length. The flow in this region is still in the devveloping phase as it ends, beyond this region, we have fully developed velocity profile and it remains constant. This region is known as Hydrodynamically fully developed region. This development of the velocity boundary layer process is illustrated below

          

                                                    Development of Velocity Boundary Layer

 The shear stress at the pipe wall is related to the velovity profile at the surface. Keeping that in mind, the velocity profile remains unchanged in the hydrodynamic fully developed region as we can see in the above illustration. The wall shear stress also remians constant in that region. We can see opposite effets happening at the hydrodynamic entrance region.

The Hydrodynamic entry length is generally taken from the pipe entrance to where the shear stress if the wall reaches within 2 percent of the fully developed value. The hydrodynammic entry length for a laminar flow is taken as $L_h\cong 0.05\cdot Re$.

 

Hagen-Poiseuille Equation:

The Hagen-Poiseuille equation is used to determine the pressure loss in a steady, hydrodynamically fully-developed, no-swirl flow in a circular pipe. It is given as follows

$\Delta p=\frac{8\cdot \mu \cdot V_{avg}\cdot L}{R^2}$

where $\mu$is the Dynamic viscosity,

$V_{avg}$ is the average velocity,

$L$ is the total length of the pipe and

$R$ is the radius of the pipe

 

Input and Given parameters:

Pipe Data:

Diameter of the pipe $D=$0.0025 meters with Radius $R=\frac{D}{2}$

Reynold's Number $Re=$2100 (given)

Hydrodynaimc entry length $L_e=0.05\cdot Re\cdot D$

Total length of the pipe $L=0.3+L_e$

Angle of the wedge $\theta =3^o$

 

Fluid Data:

Type: Water at ${20}^o\text{Celsius}$

Density $\rho =997\left(\frac{kg}{m^3}\right)$

Dynamic Viscosity $\mu =0.000890\left(\frac{N\cdot s}{m^2}\right)$

Kinematic Viscosity $\nu =\frac{\mu }{\rho }$

 

Other calculated parameters:

Average velocity $V_{avg}=\frac{Re\cdot \mu }{\rho \cdot D}$

Maximum velocity $u_{max}=2\cdot V_{avg}$

Kinematic pressure $k_p=\frac{\Delta p}{\rho }$

 

Pre-Processing:-

The pre-processing part contains setting up the block mesh, boundary and initial conditions, physical properties and the simulation time. For the current simulation, we set up the mesh using BlockMeshDict file, the boundary and initial conditions using U, p files namely for velocity and pressure , physical properties using transportProperties file and the control (simulation parameters) using the ControlDict file.

The BlockMeshDict file is the most important and the primary step as it contains the block (mesh) parameters. To simplify editing the BlockMeshDict file, we have written a Matlab code which takes minimal parameters and yields the entire block mesh file. The Matlab code and the generated BlockMeshDict file is given below:

Matlab Code:

clear all
close all
clc

% Input Data (meters)
D = 0.0025; % Diameter
R = D/2; % Radius
Re = 2100; % Reynolds Number
L_e = 0.05*Re*D; % Hydrodynamic Entry length
L = 0.3 + L_e; % Total Length
theta = 3; % Wedge Angle

scale = 1; % All values referenced to 1 meter scale

% Vertices
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)];

% Block Parameters
cellx = 50; % Cells along x-axis
celly = 1; % Cells along y-axis
cellz = 400; % Cells along z-axis
gradx = 1; % Grading along x-axis
grady = 1; % Grading along y-axis
gradz = 1; % Grading along z-axis

% Arc inputs
arc1 = [1,2];
arc2 = [4,5];

f1 = fopen('BlockMeshDict.txt','w');
fprintf(f1,'%sn','/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'%sn','  =========                 |');
fprintf(f1,'%sn','  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'%sn','   \    /   O peration     | Website:  https://openfoam.org');
fprintf(f1,'%sn','    \  /    A nd           | Version:  6');
fprintf(f1,'%sn','     \/     M anipulation  |');
fprintf(f1,'%sn','*---------------------------------------------------------------------------*/');
fprintf(f1,'%sn','FoamFile');
fprintf(f1,'%sn','{');
fprintf(f1,'%sn','    version     2.0;');
fprintf(f1,'%sn','    format      ascii;');
fprintf(f1,'%sn','    class       dictionary;');
fprintf(f1,'%sn','    object      blockMeshDict;');
fprintf(f1,'%sn','}');
fprintf(f1,'%snn','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%s %d;nn','convertToMeters',scale);

% Defining Vertices
fprintf(f1,'%sn','vertices');
fprintf(f1,'%sn','(');
fprintf(f1,'t(%d %d %d)n',v0(1),v0(2),v0(3));
fprintf(f1,'t(%d %d %d)n',v1(1),v1(2),v1(3));
fprintf(f1,'t(%d %d %d)n',v2(1),v2(2),v2(3));
fprintf(f1,'t(%d %d %d)n',v3(1),v3(2),v3(3));
fprintf(f1,'t(%d %d %d)n',v4(1),v4(2),v4(3));
fprintf(f1,'t(%d %d %d)n',v5(1),v5(2),v5(3));
fprintf(f1,'%snn',');');

% Defining Blocks
fprintf(f1,'%sn','blocks');
fprintf(f1,'%sn','(');
fprintf(f1,'thex (0 1 2 0 3 4 5 3)(%d %d %d) %s (%d %d %d)n',cellx,celly,cellz,'simpleGrading',gradx,grady,gradz);
fprintf(f1,'%snn',');');

% Defining Edges
fprintf(f1,'%sn','edges');
fprintf(f1,'%sn','(');
fprintf(f1,'tarc %d %d (%d %d %d)n',arc1(1),arc1(2),0,R,0);
fprintf(f1,'tarc %d %d (%d %d %d)n',arc2(1),arc2(2),L,R,0);
fprintf(f1,'%sn',');');

% Defining Boundaries
fprintf(f1,'%sn','boundary');
fprintf(f1,'%sn','(');
fprintf(f1,'tinletnt{n');
fprintf(f1,'tttype patch;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 1 2 0)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'toutletnt{n');
fprintf(f1,'tttype patch;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(3 4 5 3)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'tfrontnt{n');
fprintf(f1,'tttype wedge;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 3 4 1)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'tbacknt{n');
fprintf(f1,'tttype wedge;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 3 5 2)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'ttopnt{n');
fprintf(f1,'tttype wall;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(1 2 5 4)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'taxisnt{n');
fprintf(f1,'tttype empty;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 3 3 0)n');
fprintf(f1,'tt);nt}n');
fprintf(f1,');nn');
fprintf(f1,'mergePatchPairsn(n);nn');
fprintf(f1,'// ************************************************************************* //');
fclose(f1);

 

BlockMeshDict:

/*--------------------------------*- 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;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
	(0 0 0)
	(0 1.249572e-03 -3.272119e-05)
	(0 1.249572e-03 3.272119e-05)
	(5.625000e-01 0 0)
	(5.625000e-01 1.249572e-03 -3.272119e-05)
	(5.625000e-01 1.249572e-03 3.272119e-05)
);

blocks
(
	hex (0 1 2 0 3 4 5 3)(50 1 400) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.250000e-03 0)
	arc 4 5 (5.625000e-01 1.250000e-03 0)
);
boundary
(
	inlet
	{
		type patch;
		 faces 
		(
			(0 1 2 0)
		);
	}
	outlet
	{
		type patch;
		 faces 
		(
			(3 4 5 3)
		);
	}
	front
	{
		type wedge;
		 faces 
		(
			(0 3 4 1)
		);
	}
	back
	{
		type wedge;
		 faces 
		(
			(0 3 5 2)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

We use the 'blockMesh' command to generate the mesh and 'checkMesh' command to lookout for any discrepancies in our desired mesh parameters or any errors.

Below are the snapshots of the generated mesh

According to the given parameters we have given, we have obtained a mesh with simple grading and 50 cells in x-direction, 1 cell in y-direction and finally 400 cells in z-direction which can corroborated in the above mesh snapshots.

 

U file:

The U file is located in the '0' folder and is used for specifying the inlet and boundary conditions for velocity. Here, we have given inlet as uniform and a value of 0.7498 m/sec which is fixed in the x-direction. The outlet is specified as zeroGradient which is a Neumann Boundary condition as we have a pressure outlet and we will calculate the velocity at the outlet as part of the time marching simulation. Moving on, the front and back are specified as wedge boundary condition. The top (wall) is specified with noSlip conditon. The pipe wall imparts a large shear stress onto the fluid particles that are in contact with it causing them to have zero velocity, this is called the noSlip condition. Finally, the axis specified as empty to regulate the calculation of fluid characteristics at the axis domain.

/*--------------------------------*- 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.7498 0 0);

boundaryField
{

    inlet
    {
        type            fixedValue;
		value			uniform (0.7498 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

    top
    {
        type            noSlip;
    }

    axis
    {
        type            empty;
    }
}

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

 

p file:

This file is also located in the '0' folder and is used for specifying the inlet and boundary conditions for pressure. The inlet is given as zeroGradient which is a Neumann Boundary Condition and will be calculated during the time marching part of the simulation. The outlet is a fixed value and specified as uniform 0 indicating the guage pressure which is related to the atmospheric pressure in this case 1 atm. Similar to the U file, front and back are specified as wedge boundary conditions and axis as empty to regulate the calculation of fluid characteristics at the axis domain. top is specified as zeroGradient as it will be calculated during the time marching part of the simulation. 

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

    top
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
}

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

 

Note:

While setting up the BlockMeshDict, U and p files, it is very important to maintain same boundary parameters (inlet, outlet, front, etc..), if not OpenFoam throws an error called "cannot find patchfield entry" due to mismatched parameters.

 

controlDict file:

The OpenFOAM solvers begin all runs by setting up a database that controls Input/Output and, since output of data is usually requested at intervals of time during the simualtion run, as time is an inextricable part of the database. The controlDict dictionary stands for time and data Input/Output control and sets the I/O parameters for the creation of the database. This dictionary contains the start, end time for the simulation including the time step size. Here, we have given end time of 2 with time step of 0.0001, so we have about 20000 iterations. It also contains a parameter called writeInterval at which data is saved into the solution database which means it will take data for every 5 iterations. These parameters are mandatory and have default values stored in sample cases which can be tweaked according to the user's specification and the nature of the simulation.

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

writeControl    timeStep;

writeInterval   5;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

 

transportProperties:

In OpenFoam, solvers which do not include energy/heat parameters, include libraries of models for viscosity $\left(\nu \right)$. These models generally relate viscosity to strain rate $\left(\stackrel{.}{\gamma }\right)$ and are specified by the user in the transportProperties dictionary.

/*--------------------------------*- 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.92e-07;


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

 

At this point, we are done with the pre-processing part of this project. To run the simulation, we use the 'icoFoam' command to invoke the icofoam solver which is an incompressible laminar Navier-Stokes equation solver using the PISO algorithm. The simulation will start and the residuals are displayed accordingly and will run as specified in the controlDict file until the residuals converge. After the simulation is done, we step into the post-processing phase of our project.

 

Post-Processing:

Once we are done with the simulation run, we move onto the post-processing phase which is data visualization. Here, we will be using the Paraview platform to visualize our flow properties. We will first visualize the velocity distribution at the entry and the exit of the pipe (shown in wedge form).

 

Velocity Distribution at the entrance of the pipe:

In the above snapshot, we can see the velocity distribution at the entrance of the pipe and along the its length. There is a boundary layer which is gradually forming at the hydrodynamic entry region.

 

Velocity Distribution at the exit of the pipe:

 

 Here we see the velocity distribution at the end of the pipe and the uniform velocity gradient at the hydrodynamically fully developed flow region.

 

We will now be plotting the velocity profile at different locations of the x domain to understand how the flow is gradually formed.

Velocity profile at x = 0.001 m:

Here, we can see the inception of the boundary layer formation. The velocity is plotted from the axis to the wall of the pipe. So, the profile starts from the axis with velocity $\approx$0.791 m/s (maximum velocity) to 0 m/s at the wall. The velocity at the wall is 0 due to high shear stress. This is affected onto the flow to the adjacent fluid layers which are almost staitonary due to the shear stress. As a result, the velocity gradient is non-uniform. Since this is near the hydrodynamic entry region, the boundary flow is slowly being initiated.

 

Velocity profile at x = 0.1 m:

 

The above velocity profile at 0.1 meters belongs to the fully developed hydrodynamic flow region. Here, we can see the boundary layer forming and we have a unform velocity profile. We can easily observe that the profile is parabolic in nature. And the maximum velocity attained here is 1.334 m/s. Similar to the 0.001 m snapshot, here we have negligible velocity at the wall due to the shear stress. We also observe that there is a lot of viscosity changes between the fluid layers.

 

Velocity profile at x = 0.28 m:

Now, we have the velocity profile at x = 0.28 meters. We observe similar trend when compared to the 0.1 m plot with the exception of more uniformity in the gradient. Here, the maximum velocity achieved is 1.481 m/s. The location taken here is from the analytical value of the hydrodynamic entry length. We will see the fully developed flow near the exit of the pipe in the next snapshot.

 

Velocity profile at x = 0.5 m:

Above, we have a plot of the velocity profile at x-domain location of 0.5 meters. Here, the maximum velocity achieved is 1.497 m/s which is at the axis of the pipe. This plot was taken to observe and compare the velocity of fully developed flow near the exit of the pipe to the above plot of 0.2 m. The maximum velocity in 0.2 meters plot is 1.481 m/s and in the above plot at 0.5 meters is 1.497 m/s which is very close.

 

Pressure profile along pipe length:

Now we observe the pressure profile captured along the pipe. And it is clearly evident that the pressure is decreasing as we move along the domain. This tells us that the pressure is high at the pipe entrance due to the non-uniformity in the velocity gradient due to high shear stress as we have observed in the above velocity profiles. The pressure naturally decreases as the velocity moves towards more uniformity where the boundary layers are uniform in nature. We can see the non-uniform pressure gradient till the hydrodynamic entry length. As the flow moves on, we see constant and uniform decrease in the pressure gradient compensating for the viscous forces in the flow. So, we can concur that the pressure drop is due to the viscous forces of the fluid. Here, the maximum pressure captured is about 2.282 pascals.

 

Velocity along Radial distance at different x-domain locations:

In the above plot, we have all the velocity profiles taken at different locations from the inlet. These profiles show how the velocity distribution is starting and slowly developing as it moves forward in the domain.

 

Validation Studies:

Now, we will extract the numerical solution data from Paraview and process them in Matlab and compare them with the analytical data. This is done so that we can understand how accurate our numerical solution is when compared with the analytical solution. The Matlab code is given below

 

Post-process Matlab code:

clear all
close all
clc

% Input Data (meters)
D = 0.0025; % Diameter
R = D/2; % Radius
Re = 2100; % Reynolds Number
L_e = 0.05*Re*D; % Hydrodynamic Entry length
L = 0.3 + L_e; % Total Length
theta = 3; % Wedge Angle
mu = 0.000890; % (Pa.s) Dynamic Viscosity
rho = 997; % (kg/m^3) Density of Water
nu = mu/rho; % (m^2/s) Kinematic viscosity

% Reading and Extracting Velocity Data from Paraview
f1 = fopen('velocity.csv','r');
A = fgetl(f1);
for i = 1:1000
    B = fgetl(f1);
    C = strsplit(B,',');
    Vel_n(i) = C(2);
    Len_n(i) = C(1);
end
Vel_n = str2double(Vel_n');
Len_n = str2double(Len_n');

%Reading and Extracting Pressure Data from Paraview
f2 = fopen('Pressure.csv','r');
A = fgetl(f2);
for i = 1:1000
    B = fgetl(f2);
    C = strsplit(B,',');
    Press_n(i) = C(2);
    Len_n2(i) = C(1);
end
Press_n = str2double(Press_n');
Len_n2 = str2double(Len_n2');

%Reading and Extracting Shear Stress Data from Paraview
f3 = fopen('tau.csv','r');
A = fgetl(f3);
for i = 1:1000
    B = fgetl(f3);
    C = strsplit(B,',');
    tau_n(i) = C(2);
    Len_n3(i) = C(1);
end
tau_n = str2double(tau_n');
Len_n3 = str2double(Len_n3');

r = linspace(0,R,1000); % Radial Distance
% Hagen poiseuille's Solution
V_avg = (Re*mu)/(rho*D); % Average velocity
u_max = 2*V_avg; % Maximum Velocity
V_r = u_max*(1-(r.^2./R^2));


% Hagen-Poiseuille Equations for Pressure Drop
P_len = linspace(0,L,1000); % Pipe Length
for i = 1:length(P_len)
    delta_p(i) = (8*mu*P_len(i)*V_avg)/R^2; % Pressure Drop
    k_press(i) = delta_p(i)/rho; % Kinematic Pressure
end

k_press = fliplr(k_press);

% Analytical Shear Stress calculation
for i=1:length(r)
   tau_an(i)=(2*mu*u_max*r(i))/R^2;
end


figure(1)
plot(r,V_r,'r','linewidth',1) % Anaytical Data
hold on
plot(Len_n,Vel_n,'--','color','k','linewidth',1) % Numerical Data
title('Velocity: Numerical vs Analytical Data')
xlabel('Radial Length (m)')
ylabel('Velocity in m/s')
legend('Analytical Method','Numerical Method')
grid on

figure(2)
plot(P_len,k_press,'r','linewidth',1) % Analytical Data
hold on
plot(Len_n2,Press_n,'--','color','k','linewidth',1) % Numerical Data
title('Pressure: Numerical vs Analytical Data')
xlabel('Radial Length (m)')
ylabel('Pressure (Pa)')
legend('Analytical Method','Numerical Method')
grid on

figure(3)
plot(r,tau_an,'r','linewidth',1) % Analytical Data
hold on
plot(Len_n3,tau_n,'--','color','k','linewidth',1) % Numerical Data
title('Shear Stress: Numerical vs Analytical Data')
xlabel('Radial Length (m)')
ylabel('Shear Stress (N/m^2)')
legend('Analytical Method','Numerical Method')
grid on

In the above code, we are reading the extracted data files from paraview and parse them to plot with the analytical solution.

The formulae used in the analytical solution are given below

Average Velocity:

$V_{avg}=\frac{Re\cdot \mu }{\rho \cdot D}$

Radius at various radial lengths at the fully developed region:

$V_r=2\cdot V_{avg}\cdot (1-\frac{r^2}{R^2})$

Pressure drop across the pipe length:

$\Delta p=\frac{8\cdot \mu \cdot V_{avg}\cdot L}{R^2}$

Kinematic Pressure:

$K_{press}=\frac{\Delta P}{\rho }$

Shear Stress:

$\tau =\frac{2\cdot \mu \cdot u_{max}\cdot r}{R}$

 

All of the following plots show the comparison of Numerical to Analytical studies for Velocity, Pressure and Shear stress respectively

 

In the above plot, we can see that both data are almost similar.

 

In the above plot, we can see that the data is almost similar and is following the pressure drop trend.

 

In the above plot, we can see that the data is almost similar.

 

Conclusion:-

We have automated the process of creating the blockMeshdict file for our simulations.

We have seen that the data obtained from the numerical simulation is very close to the analytical solutions.

 

References:-

Fluid Mechanics Fundamentals and Applications by Yunus Cengel, John Cimbala

https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation

https://cfd.direct/openfoam/user-guide