BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

Simulation of Laminar Viscous Flow through a regular Pipe using Wedge Boundary conditions in OpenFOAM and post-processed in ParaView

SShivanshu Kumar21 Sept 2020
CFDMATLAB

Problem Specification:

  • Write a program in MATLAB that can generate the computational mesh automatically for any wedge angle (wedge angle should be less than 5 degrees) and grading schemes.
  • For a baseline mesh, show that the velocity profile matches with the Hagen Poiseuille's equation.
  • Show that the velocity profile is fully developed.
  • Post-process velocity, pressure and shear stress.

Introduction:

In this project, we’re going to simulate the incompressible laminar flow through a cylindrical pipe using the icoFoam solver in OpenFoam. Here, we’re also going to validate numerical results with the analytical solution obtained using the Hagen-Poiseuille equation.

Simulation of flow through an entire pipe would be computationally demanding, and since the geometry of the pipe is axisymmetric we can use a wedge section of pipe for the simulation.

NOTE: Post-processing is done in ParaView.

Theory:

We’re assuming an incompressible, laminar, steady-flow through a smooth cylindrical pipe. As shown in the above diagram, when the fluid enters the pipe, we can see viscous effects because of the flow resistance characteristic of the pipe wall. The viscous effects are important at the boundary layer. The boundary layer is a dynamic phenomenon, and it is seen that it grows from the wall to such an extent that they all merge on the centerline of the pipe. That is when the flow is said to be a fully developed flow. The length required for the flow to reach a fully developed flow condition is called the hydrodynamic entrance length `(Le)`.

For laminar flow, `(Le)/D=0.06 Re`

where `Le=` Hydrodynamic entrance length

`D=` Diameter of pipe

`Re=` Reynolds number

The fully developed pipe flow is generally known as Hagen-Poiseuille flow. The velocity profile of a fully developed laminar horizontal circular pipe flow is given by:

`V=1/(4μ) (∂p)/(∂x)(y^2-R^2 )`

And, the dynamic pressure difference in the pipe is given by:

`∆p=(32μLV)/D^2`

where `∆p=` dynamic pressure (Pa)

`μ=` dynamic viscosity (Pa.s)

`y=` vertical distance from the pipe axis (m)

`R=` radius of the pipe (m)

`D=` diameter of the pipe (m)

`V=` average velocity

Pipe Geometry:

We considered a smooth cylindrical pipe of diameter 8mm and length 1500mm for this case study.

`L=1500 mm or 1.5m`

`D=8mm or 0.008m`

Preliminary Data:

Fluid: Water

Temperature: `298K`

Dynamic Viscosity @ 298K and 1 atm: `μ=0.001003 Pa.s`

Density of water: `ρ=998.2 (kg)/m^3`

Kinematic Viscosity: `ν=1.0048e-06 m^2/s`

Reynold’s number: `Re=2100`

Analytical Solution: (Hagen-Poiseuille)

Average Velocity, `V=(μ*Re)/(ρ*D)=(0.001003*2100)/(998.2*0.008)=0.2637 m/s`

Maximum Velocity, `V_max=2*V_(avg)=2*0.2637=0.5274 m/s`

Hydrodynamic entrance length, `Le=0.06*Re*D=0.06*2100*0.008=1.008 m`

Total length of the pipe, `L=1.5 m`

Pressure Difference for 1.1m to 1.5m, `∆p=(32*μ*L*V)/D^2 =(32*0.001003*0.4*0.2637)/(0.008*0.008)=52.89822Pa`

NOTE: For constant density problems we can simplify the system of equations. The incompressible solvers transform the pressure from static pressure, `p_s` [Pa] to kinematic pressure `p_k`, i.e. `p_k=p_s/ρ` `[m^2/s^2]`. Note that the kinematic pressure has the name `p` in OpenFOAM solvers, the same as the name for the static pressure `p`. To interpret the results from incompressible solvers, the kinematic pressure should be multiplied by the reference density (constant) value.

Therefore, Kinematic Pressure= Static Pressure/Density=52.89822/998.2=0.0529 or 0.053 `m^2/s^2` (approx.)

Pre-Processing:

Geometry and Mesh generation: Flow through a cylinder can be very challenging because of high resource requirements. To minimize the requirement of high computing resources and considering the axisymmetric property of pipe, only a wedge section of the domain was simulated.

A MATLAB code was written which will auto-generate the computational mesh and geometry for any length, radius, and wedge angle.

MATLAB code for generating blockMeshDict file:

This MATLAB code creates the blockMeshDict file for wedge angle 3deg, radius 0.004m, and total length 1.5m. The blockMeshDict file is used to generate the geometry and computational mesh for the simulation using blockMesh command in OpenFoam.

clear all
close all
clc

%Inputs

%Total length (m)
L=1.5; 
%Radius of the pipe (m)
r=0.004; 
%Wedge angle (deg)
theta=3;
%No. of cells along x-axis
nx=500;
%no. of cells along y-axis
ny=40;
%no. of cells along z-axis
nz=1;

line_1='/*--------------------------------*- C++ -*----------------------------------*';
line_2='| =========                 |                                                 |';                
line_3='| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |';
line_4='|  \    /   O peration     | Website:  https://openfoam.org                  |';
line_5='|   \  /    A nd           | Version:  8                                     |';
line_6='|    \/     M anipulation  |                                                 |';
line_7='*---------------------------------------------------------------------------*/';


%blockMeshDict
f1=fopen('blockMeshDict','w')
fprintf(f1,'%sn',line_1);
fprintf(f1,'%sn',line_2);
fprintf(f1,'%sn',line_3);
fprintf(f1,'%sn',line_4);
fprintf(f1,'%sn',line_5);
fprintf(f1,'%sn',line_6);
fprintf(f1,'%sn',line_7);

%FoamFile Generation
fprintf(f1,'FoamFilen{n');
fprintf(f1,'%12st 2.0;n','version');
fprintf(f1,'%11st ascii;n','format');
fprintf(f1,'%10st dictionary;n','class');
fprintf(f1,'%11st blockMeshDict;n','object');
fprintf(f1,'}n');
line_8='// ************************************************************************* //';
fprintf(f1,'%snn',line_8);

fprintf(f1,'convertToMeters 1;nn');

%Vertices
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'t (%d %d %d)n',0,0,0); %Point 0
fprintf(f1,'t (%d %d %d)n',0,r*cosd(theta/2),r*sind(theta/2)); %Point 1
fprintf(f1,'t (%d %d %d)n',0,r*cosd(theta/2),-r*sind(theta/2)); %Point 2
fprintf(f1,'t (%d %d %d)n',L,0,0); %Point 3
fprintf(f1,'t (%d %d %d)n',L,r*cosd(theta/2),r*sind(theta/2)); %Point 4
fprintf(f1,'t (%d %d %d)n',L,r*cosd(theta/2),-r*sind(theta/2)); %Point 5
fprintf(f1,');nn');

%Blocks
fprintf(f1,'blocksn');
fprintf(f1,'(n');
fprintf(f1,'     hex (0 3 5 2 0 3 4 1) (%d %d %d) ',nx,ny,nz);
fprintf(f1,'simpleGrading (1 0.1 1)nn');
fprintf(f1,');nn');

%Edges
fprintf(f1,'edgesn(nt arc %d %d (%d %d %d)n',1,2,0,r,0);
fprintf(f1,'t arc %d %d (%d %d %d)n',4,5,L,r,0);
fprintf(f1,');nn');

%Boundary
fprintf(f1,'boundaryn(n');

fprintf(f1,'t inletn'); %Inlet
fprintf(f1,'t{n');
fprintf(f1,'tt type patch;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt(n');
fprintf(f1,'ttt (0 1 2 0)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t outletn'); %Outlet
fprintf(f1,'t{n');
fprintf(f1,'tt type patch;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (3 5 4 3)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t Noslipwalln'); %Top-Noslipwall
fprintf(f1,'t{n');
fprintf(f1,'tt type wall;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (1 4 5 2)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t backn'); %Back
fprintf(f1,'t{n');
fprintf(f1,'tt type wedge;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 2 5 3)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t frontn'); %Front
fprintf(f1,'t{n');
fprintf(f1,'tt type wedge;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 3 4 1)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');

fprintf(f1,'t axisn'); %Axis
fprintf(f1,'t{n');
fprintf(f1,'tt type empty;n');
fprintf(f1,'tt facesn');
fprintf(f1,'tt (n');
fprintf(f1,'ttt (0 3 3 0)n');
fprintf(f1,'tt );n');
fprintf(f1,'t}n');
fprintf(f1,');nn');

%mergePatchPairs
fprintf(f1,'mergePatchPairsn(n');
fprintf(f1,');n');

fprintf(f1,'%sn',line_8);

fclose(f1);

blockMeshDict:

/*--------------------------------*- 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 3.998629e-03 1.047078e-04)
	 (0 3.998629e-03 -1.047078e-04)
	 (1.500000e+00 0 0)
	 (1.500000e+00 3.998629e-03 1.047078e-04)
	 (1.500000e+00 3.998629e-03 -1.047078e-04)
);

blocks
(
     hex (0 3 5 2 0 3 4 1) (500 40 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 4.000000e-03 0)
	 arc 4 5 (1.500000e+00 4.000000e-03 0)
);

boundary
(
	 inlet
	{
		 type patch;
		 faces
		(
			 (0 1 2 0)
		 );
	}
	 outlet
	{
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	}
	 Noslipwall
	{
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	}
	 back
	{
		 type wedge;
		 faces
		 (
			 (0 2 5 3)
		 );
	}
	 front
	{
		 type wedge;
		 faces
		 (
			 (0 3 4 1)
		 );
	}
	 axis
	{
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

Initial and Boundary Conditions:

Velocity, U

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

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.2637 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    Noslipwall
    {
        type            noSlip;
    }

    back
    {
        type            wedge;
    }
    
    front
    {
        type            wedge;
    }
    
    axis
    {
        type            empty;
    }
}

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

Pressure, p

/*--------------------------------*- 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;
    }
    
    Noslipwall
    {
        type            zeroGradient;
    }

    back
    {
        type            wedge;
    }
    
    front
    {
        type            wedge;
    }
    
    axis
    {
        type            empty;
    }
}

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

controlDict

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

fvSchemes

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

fvSolution

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

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


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

transportProperties: (Kinematic Viscosity=1.0048e-06)

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


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

Post-Processing:

Velocity Contour: Hydrodynamic entrance region (near Inlet)

Velocity Contour: Fully developed flow region

Velocity Plots: Hydrodynamic entrance region

at x=0.15m

at x=0.35m

at x=0.55m

at x=0.75m

Velocity Plots: Fully developed flow region

at x=1.2m

at x=1.3m

at x=1.4m

at x=1.5m

Pressure Contour:

Pressure-drop Plots:

NOTE: In OpenFoam, the kinematic pressure is calculated, and in the below plots y-axis is denoting the kinematic pressure.

x=1.1m to x=1.5m (Fully developed flow region)

x=0 to x=1.5 (over entire length)

Shear stress (tau) Plot:

MATLAB code to generate shear stress plot (analytical solution):

Shear stress, `tau=mu(delu)/(dely)` 

clear all
close all
clc

%Inputs
r=0.004; %radial length
y=linspace(0,r,100);%radius array
Y=linspace(0,r,99); %raidus array
u_max=0.5274; %maxm. velocity
mu=0.001003; %dynamic viscosity

for i=1:length(y)
    u(i)=u_max*(1-((y(i)/r)^2));
end

for i=1:length(y)-1
    du=u(i+1)-u(i);
    dy=y(i+1)-y(i);
    
    tau (i)=mu*(du/dy);
end

figure (1)
plot (Y,tau)
xlabel ('Radial Distance (m)')
ylabel ('Shear-stress (analytical)')
set (gca,'YDIR','reverse')
grid on

Comparison between Analytical and Computational results:

Parameter Analytical Computational Error %
Maximum Velocity (m/s) 0.5274 0.53 0.49
Kinematic Pressure (`m^2/s^2`) 0.0529 0.0532 0.56
Shear Stress (Pa) 0.260 0.261 0.38

Conclusion:

1. From the velocity plots, we can observe that: for the hydrodynamic entrance region the velocity profile changes, but as the flow reaches to fully developed region, the velocity profile doesn't change and remains constant; which satisfies the Haigen-Poiseuille statement.

2. The maximum velocity is at the center of the pipe, whereas the minimum velocity is at the pipe surface because of viscosity.

3. A parabolic velocity profile is obtained after the hydrodynamic entrance length region, and this profile remains unchanged thereafter. By this, we can conclude that the length of the pipe as calculated analytically suffices to get a fully developed flow.

4. The pressure difference is maximum at the inlet of the pipe, whereas it gets 0 (boundary condition for pressure) at the outlet of the pipe.

5. Shear stress is maximum at the surface of the pipe, whereas its value plummets as it reaches the center of the pipe.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.