Pipe flow simulation in OpenFOAM (Part 1/2: Wedge boundary conditions)

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM:

[PART 1/2]: WEDGE BOUNDARY CONDITION:

OBJECTIVES:
1. Calculation of Pre-eliminary quantinites related to the flow.
2. Creating the Mesh script (blockMeshDict file) for specifying the geometry and boundary condition.
3. Simulation Results and post processing the velocity profile at various locations of the pipe for fully developed flow.

[PART 2/2]: SYMMETRY BOUNDARY CONDITION:

OBJECTIVES:
1. Creating the Mesh-Script (blockMeshDict file) with symmetry boundary condition.
2. Simulation Results and post processing the velocity profile at different points of the pipe for fully developed flow.
3. Comparison of Wedge Boundary Condition and Symmetry Boundary Condition against Hagen Poiseuille Equation.

____________________________________________________________________________

[PART 1/2]:

[I] INTRODUCTION:

In this project we aim to simulate a section of pipe (Wedge shaped geometry) to understand various aspects of flow through a pipe. The fluid chosen is water and the flow is laminar and incompressible. Also we choose to neglect any changes in cross section. The specified Reynolds Number `Re` is specified as 2100.

Internal flow occurs when the fluid is filled entirely in a conduit and is driven primarily by a pressure difference.

Fluid velocity in a pipe changes from zero at the surface because of the no-slip boundary condition to a maximum at the pipe center. The average velocity `V_(avg)` remains constant in an incompressible flow when the cross section of the pipe is constant.

                                      pipe

Reynold's Number `Re`: The tranistion from Laminar flow to Turbulent flow occurs from Laminar (`Re <= 2100`) to Transition (`2100 <= Re <= 4000`) to Turbulent (`Re >= 4000`) for a smooth pipe with no viscous dissipation. it depends on the geometry, face roughness, flow velocity, surface temperature etc. Reynold's number is an important non-dimensional number which is the ratio of Inertial Forces to Viscous Forces and is expressed for internal flow in a circular pipe as: `Re = (V_(avg)D)/nu = (rho*V_(avg)*D)/mu`

Hydrodynamic entry length: The region from the pipe inlet to the point at which the velocity boundary layer merges at the centreline is called the hydrodynamic entrance region. The region beyond the entrance region in which the velocity profile is fully developed is called the fully developed region.

pipe 2

In laminar flow the entry length is given approximately as `L_(laminar) ~= 0.06*Re*D`

Velocity profile for laminar flow: The velocity profile for a fully developed laminar flow in a pipe is parabolic with a maximum at the centreline and minimum at R = 0.

`u(r) = R^2/(4mu)*((dp)/dx)(1-r^2/R^2)` It is obtained for a fluid element at a radial distance `r` from the center. The average velocity can be written as: `V_(avg) = R^2/(8mu)*((dp)/dx)` hence the final expression for the velocity profile can be written as:

`u(r) = 2V_(avg)(1-r^2/R^2)`

Hagen Poiseuille Flow: It states that ''for a specified flowrate, the pressure dop and thus the required pumping power is inversely proportional to the length of the pipe and the viscosity of the fluid, but it is inversely proportional to the fourth power of the radius (or diameter) of the pipe''.

Or in other words the flowrate: `Q = V_(avg)A_c = (DeltaPpiD^4)/(128muL)`
Thus the average velocity is `V_(avg) = (DeltaPD^2)/(32muL)`

[II] SPECIFICATIONS OF THE PROBLEM:

1. Given Reynolds Number = 2100

2. The working fluid is water at `25^oC`

3. The entry length is `L_(entry) = 0.06ReD` or `L_(entry) = 2.52m`

4. The length of the pipe taken for simulation is 2.8m

5. The Density of water `rho = 997.0 kgm^-3`

6. The dynamic viscosity of water is `mu = 0.89*10^-3 Pa.S`

7. The kinematic viscosity of water is `nu = 8.927*10^-7 m^2s^-1`

CALCULATED VALUES: (Analytical Solution)

1. Average Velocity `V_(avg) = (Re*mu)/(rho*D) = 0.09373 ms^-1`

2. Maximum Velocity `V_(max) = 2*V_(avg) = 0.1875 ms^-1`

3. Pressure drop `DeltaP = (32muLV_(avg))/D^2 = 16.817 Pa`

4. Kinematic Pressure drop `DeltaP_(kinematic) = (DeltaP)/rho = 0.01687 m`

[III] MATLAB PROGRAM:

In this project we have used MATLAB to write a program that will generate the mesh file (blockMeshDict file) for OpenFOAM. The geometry used is wedge made out of a single block and for different wedge angles our MATLAB program creates the required blockMeshDict file.

%% MATLAB program to Auto-generate the blockMeshDict file for a Wedge geometry:

clear all;
close all;
clc;

%% Input Parameters
Re = 2100; % Given Reynold's Number
d = 0.02; % Diameter of the pipe
r = d/2; % Radius of the pipe
L_entrance = 0.06*Re*d; % Length of the pipe
fprintf('Enter 2, 3, 4, 5 for wedge boundary condition nEnter 10, 25, 45 for symmetry boundary conditionn');
theta = input('Enter the angle in angles:n '); % Angle of wedge
mu = 0.89e-3; % Viscosity of water at 25 deg C
rho = 997.0; % Density of water at 25 deg C
nu = mu/rho; % Kinematic Viscosity of water at 25 deg C

%% Calculation of pressure and velocity:
L = L_entrance + 0.3; % Length of the Wedge
v_avg = Re*mu/(d*rho); % Average Velocity
v_max = 2*v_avg; % Max Velocity
del_P = (32*mu*v_avg*L)/(d^2); % Pressure Drop 
kin_P = del_P/rho; % KinematicPressure Drop

%% Calculating the vertices of the 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)];

%% String format specifying for BlockMesh File:

head1 = "/*--------------------------------*- C++ -*----------------------------------*";
head2 = "  =========                 |                                                  ";
head3 = "  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox            ";
head4 = "   \    /   O peration     | Website:  https://openfoam.org                   ";
head5 = "    \  /    A nd           | Version:  6                                      ";
head6 = "     \/     M anipulation  |                                                  ";
head7 = "*---------------------------------------------------------------------------*/";

b4 = blanks(4);
b8 = blanks(8);
b12 = blanks(12);
b16 = blanks(16);

dict1 = "    version     2.0;";
dict2 = "    format      ascii;";
dict3 = "    class       dictionary;";
dict4 = "    object      blockMeshDict;";
dict5 = "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //";

conv = "convertToMeters 1;";

vertex1 = "vertices";

blocks1 = b4 + "hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)";

footer = "// ************************************************************************* //";

%% Creating BlockMeshDict file for OpenFOAM:

f1 = fopen('blockMeshDict_45.txt','wt');
% Printing the Header File
fprintf(f1,'%sn',head1);
fprintf(f1,'%sn',head2);
fprintf(f1,'%sn',head3);
fprintf(f1,'%sn',head4);
fprintf(f1,'%sn',head5);
fprintf(f1,'%sn',head6);
fprintf(f1,'%sn',head7);

fprintf(f1,'FoamFilen{n');
fprintf(f1,'%sn',dict1);
fprintf(f1,'%sn',dict2);
fprintf(f1,'%sn',dict3);
fprintf(f1,'%sn}n',dict4);
fprintf(f1,'%sn',dict5);

fprintf(f1,'n%sn',conv);

% Printing the vertices
fprintf(f1,'nverticesn(n');
fprintf(f1,'%s(%d %d %d)n',b4,v0(1),v0(2),v0(3));
fprintf(f1,'%s(%d %f %f)n',b4,v1(1),v1(2),v1(3));
fprintf(f1,'%s(%d %f %f)n',b4,v2(1),v2(2),v2(3));
fprintf(f1,'%s(%f %d %d)n',b4,v3(1),v3(2),v3(3));
fprintf(f1,'%s(%f %f %f)n',b4,v4(1),v4(2),v4(3));
fprintf(f1,'%s(%f %f %f)n',b4,v5(1),v5(2),v5(3));
% Printing the blocks
fprintf(f1,');nnblocksn(n');
fprintf(f1,'%sn);nn',blocks1);
fprintf(f1,'edgesn(n');
% Printing the edges
fprintf(f1,'%sarc 1 2 (0 %0.3f 0)n',b4,r);
fprintf(f1,'%sarc 4 5 (%0.4f %0.3f 0)n',b4,L,r);
% Printing the boundaries
fprintf(f1,');nnboundaryn(n');

fprintf(f1,'%saxisn%s{n',b4,b4);
fprintf(f1,'%stype empty;n%sfacesn%s(n%s(0 3 3 0)n%s);n%s}n',b8,b8,b8,b12,b8,b4);

fprintf(f1,'%swalln%s{n',b4,b4);
fprintf(f1,'%stype wall;n%sfacesn%s(n%s(2 5 4 1)n%s);n%s}n',b8,b8,b8,b12,b8,b4);

fprintf(f1,'%sinletn%s{n',b4,b4);
fprintf(f1,'%stype patch;n%sfacesn%s(n%s(0 0 2 1)n%s);n%s}n',b8,b8,b8,b12,b8,b4);

fprintf(f1,'%soutletn%s{n',b4,b4);
fprintf(f1,'%stype patch;n%sfacesn%s(n%s(3 4 5 3)n%s);n%s}n',b8,b8,b8,b12,b8,b4);

fprintf(f1,'%swedgeFrontn%s{n',b4,b4);
fprintf(f1,'%stype wedge;n%sfacesn%s(n%s(0 3 5 2)n%s);n%s}n',b8,b8,b8,b12,b8,b4);

fprintf(f1,'%swedgeBackn%s{n',b4,b4);
fprintf(f1,'%stype wedge;n%sfacesn%s(n%s(0 1 4 3)n%s);n%s}n',b8,b8,b8,b12,b8,b4);

fprintf(f1,');nnmergePatchPairsn(n);nn');
fprintf(f1,'%s',footer);



 This code would generate the blockMeshDict file for Wedge angles (2, 3 and 4 degrees). The blockMeshDict file for Wedge angle = `3^o` is shown below:

/*--------------------------------*- 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 0.009997 -0.000262)
    (0 0.009997 0.000262)
    (2.820000 0 0)
    (2.820000 0.009997 -0.000262)
    (2.820000 0.009997 0.000262)
);

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

edges
(
    arc 1 2 (0 0.010 0)
    arc 4 5 (2.8200 0.010 0)
);

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

mergePatchPairs
(
);

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

PRESSURE BOUNDARY CONDITION:

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

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    axis
    {
        type            empty;
    }

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
	value		uniform 0.0168682;
    }

    wall
    {
        type            zeroGradient;
    }

    wedgeFront
    {
        type            wedge;
    }

    wedgeBack
    {
        type            wedge;
    }
}

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

VELOCITY BOUNDARY CONDITION:

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

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.0937312 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

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

    outlet
    {
        type            zeroGradient;
    }
    
    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    wedgeFront
    {
 	type		wedge;
    }
  
    wedgeBack
    {
 	type		wedge;
    }
}

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

[IV] MESH PROFILE:

 3 deg mesh 1

 wedge 3 deg 2

Wedge Boundary Condition: The wedge boundary condition defines an axis-symmetric boundary condition. The model and flow must be axis-symmetric along a central line such that all physical variables of the flow have same value and distribution at a radius and angle.The axis-symmetric ‘wedge’ boundary is specified by two planes that must be selected on separate sides of the domain on surfaces running along the central line.

[V] SIMULATION RESULTS:

1. WEDGE ANGLE: 3

Simulation time: 50s

A. Velocity Profile Near the Entrance: (Hydrodynamic Boundary Layer Formation)

wedge ent

B. Velocity Profile near the exit (Fully Developed flow)

wedge fully developed

C. Velocity Profile Curves:

i. Velocity Profile at x = 0.02m from the entrance:

wedge 0.02m 3deg

The velocity profile near the entrance shows three distinct regions:

A. Viscous Sublayer: The effect of velocity boundary layer is seen in formation in this region of the pipe.
B. Transitional Buffer layer: This layer gradually decreases with the viscous sublayer merges with the turbulent core with increasing distance from the entrance.
C. Turublent Core: The velocity profile is almost flat in this region and is comprised of the turbulent component of the flow.

ii. Velocity Profile at x = 0.1m

wedge 3 deg 0.1m

This plot shows us that the velocity boundary layer increases in size at the cost of the shrinking turbulent core and the transitional buffer layer.

iii. Velocity profile at x = 0.25m

wedge 3 0.25

iv. Velocity Profile at x = 0.5m

wedge 3 0.5m

v. Velocity Profile at x = 1.0m

wedge 3 1.0m

Near the entrance length, the velocity profile almost resembles the fully developed laminar flow with a very small turbulent component at the center line.

vi. Velocity Profile for fully developed laminar flow (at `x = L_(entrance) = 2.52m`)

wedge 3 2.52m

Fig: Fully developed flow at entrance length showing one uni-directional flow along the x axis and no radial components.

Shear Stress Distribution for fully developed flow:

shear stress

The shear stress profile matches with the literature except for slight non-linearity near the edges for unknown reasons. However it is zero at the the centreline and maximum at the walls. The linear least square fit of the data is given below.

shear stress fit

(The following curves were obtained using OriginPro 2018 graphing software)

2. VELOCITY PROFILES FOR OTHER WEDGE ANGLES:

WEDGE ANGLE: 2

A. Mesh Profile:

wedge 2 mesh

B. Velocity Profile near entrance:

2deg ent

C. Velocity Profile near the exit (fully formed flow) region:

2deg exit

WEDGE ANGLE: 4

A. Mesh Profile:

 4 deg mesh

B. Velocity Profile near entrance:

WEDGE 4 ENT

C. Velocity Profile near the exit (fully formed flow) region:

wedge 4 fully formed

[VI] CONCLUSIONS AND REMARKS:

Comparison of maximum velocities: `V_(max)`

We compare the maximum velocities obtained for various wedge angles and compare them against the Hagen Poiseuille Equation:

vel compare

Comparison of execution times and courant numbers for the simulation is shown below:

c compare


Projects by Priyotosh Bairagya

LAMINAR INCOMPRESSIBLE FLOW SIMULATION THROUGH A PIPE IN OPENFOAM: [PART 2/2]: SYMMETRY BOUNDARY CONDITION: OBJECTIVES:1. Creating the Mesh-Script (blockMeshDict file) with symmetry boundary condition.2. Simulation Results and post processing the velocity profile at Read more

I. Interpolation Schemes in Finite Volume Method: The approximation of surface and volume integrals may require values of the variable at locations other than the computational nodes of the CV. Values at these locations are obtained using interpolation formulae. Some o Read more

ANALYSIS OF NUMERICAL STABILITY OF VARIOUS ITERATIVE SOLVERS FOR TRANSIENT 2D HEAT CONDUCTION: [Part: 3/3] INTRODUCTION: The criterion of stability of a numerical scheme is determined by the way the errors propagate while the solution moves from one time-step to th Read more

MESH GENERATION AND ANALYSIS USING BLOCKMESH FOR FLOW OVER A BACKWARD FACING STEP: The purpose of the following project is to generate the geometry for a variation of the incompressible cavity flow problem in OpenFOAM. For this purpose we have modified the lid-driven c Read more

ANALYSIS OF VARIOUS ITERATIVE SCHEMES FOR THE SOLUTION OF A 2D HEAT CONDUCTION PROBLEM: [PART: 2/3] In the previous part we had explained the problem statement and the MATLAB Program in detail. In this Part we are going to explain the outputs from the 2D Heat Conduct Read more

NUMERICAL SOLUTION OF 1D SUPERSONIC NOZZLE FLOW SIMULATION BY MAC-CORMACK METHOD PROJECT OBJECTIVES: i. Numerical solution of the governing equations in both conservative and non-conservative forms. ii. Creating user defined functions for calculating the flow quan Read more

Simulation of a 2D Heat Conduction problem in steady and unsteady/transient forms using iterative methods. Project Objectives: 1. Solving the 2 Dimensional Heat conduction equation in the generalized form using various iterative techniques: i. Explicit Solver (for Read more

UNDERSTANDING LINEAR SYSTEMS(ANALYSIS OF VARIOUS ITERATIVE SCHEMES TO SOLVE A SYSTEM OF LINEAR EQUATIONS TO FIND THE EIGEN VALUES AND SPECTRAL RADIUS) (A) PROBLEM STATEMENT: Given coefficient matrix: `A = [[5,1,2],[-3,9,4],[1,2,-7]]` Given Solution Matrix: `X = [[x Read more

WEEK 4: (Effect of Grid-Size on output for the solution of 1D linear wave equation) 1. Problem Setup: Given Partial Differential Equation: `(delu)/(delt) + c(delu)/(delx) = 0` Numerical Discretization:  `u_(i,n+1) = u_(i,n) + (cDeltat)/(Deltax)*(u_(i-1,n+1) Read more


Loading...

The End