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

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 different points of the pipe for fully developed flow.
3. Comparison of Wedge Boundary Condition and Symmetry Boundary Condition against Hagen Poiseuille Equation.

_______________________________________________________________________

PART [2/2]:

[I] INTRODUCTION:

In the previous part we looked into the simulation results for a laminar, incompressible flow through a pipe in OpenFOAM. Now in this section, we would attempt the same problem with the symmetry boundary conditions.

[II] MATLAB Program:

The required MATLAB Program to automatically generate the blockMeshDict file for a given wedge angle is given below. The syntax is same as before with the exception of the boundary conditions.

%% 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 symmetry;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 symmetry;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);



The generated blockMeshDict file fro wedge angle = 25 degree is given 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.009763 -0.002164)
    (0 0.009763 0.002164)
    (2.820000 0 0)
    (2.820000 0.009763 -0.002164)
    (2.820000 0.009763 0.002164)
);

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 symmetry;
        faces
        (
            (0 3 5 2)
        );
    }
    wedgeBack
    {
        type symmetry;
        faces
        (
            (0 1 4 3)
        );
    }
);

mergePatchPairs
(
);

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

The boundary conditions are same as before however the wedge faces have a symmetry condition imposed.

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            symmetry;
    }

    wedgeBack
    {
        type            symmetry;
    }
}

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

Velocity Boundary Conditions:

/*--------------------------------*- 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		symmetry;
    }
  
    wedgeBack
    {
 	type		symmetry;
    }
}

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

[III] MESH PROFILE:

The mesh geometry for wedge angle = 25 degree is given below:

wedge 25 mesh

Symmetry Boundary Condition: The symmetry boundary condition defines a mirror face/surface. The symmetry condition is used if the physical object or geometry and the expected flow field pattern of the developed solution is mirrored along that surface. This boundary condition helps to reduces the computational domain in size and enables the modeling of a sub-domain of the complete setup.

[IV] SIMULATION RESULTS:

1. WEDGE ANGLE: `25^o`

A. Geometry near the entrance:

wedge 25 ent

The simulation shows a developing hydrodynamic boundary layer at the entrrance. The free stream velocity keeps increases as we move farther from the entrance of the pipe.

B. Geometry near the exit (Fully Developed flow):

wedge 25 exit

C. Velocity Profile curves:

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

wedge 0.02 25deg

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 from the entrance:

25 deg 0.1m

The above velocity profile shows a shrinking turbulent core.

iii. Velocity profile at x = 0.25 m from the entrance:

25 deg 0.25m

iv. Velocity profile at x = 0.5m from the entrance:

25 deg 0.5m

v. Velocity profile at x = 1.0m from the entrance:

25 deg 1.0m

vi. Velocity profile at entrance length `L_(entrance) = 2.52m` (Fully developed flow):

25 deg 2.52m

(uni-directional flow along the x axis indicating fully laminar velocity profile)

Shear Stress Profile:

shear stress sym

Almost linear shear stress profile obtained by the simulation. Shear stress is zero at the walls and maximum at the centerline.

Linear fit of Shear stress data:

shear stress fit sym

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

2. SIMULATION RESULTS FOR OTHER ANGLES:

A. WEDGE ANGLE: `10^o`

i. Mesh Profile:

10 deg mesh

ii. Geometry near the entrance:

wedge 10 ent

iii. Geometry at the fully developed flow region:

10 deg exit

 Cross sectional velocity profiles at different distances from entrance:

iv. Velocity Profile at x = 0.2m from the entrance:

10 deg 0.2m

v. Velocity Profile at x = 2.52m from the entrance showing fully developed flow:

10 deg 2.52m

Exactly identical as before with minor variations in magnitude of `V_(max)`

B. WEDGE ANGLE `45^o`

i. Mesh Profile:

45 deg mesh

ii. Geometry near the entrance region:

45 deg ent

iii. Geometry at the fully developed region:

45 deg exit

iv. Velocity Profile at x = 0.2m from the entrance:

45 deg 0.2m

v. Velocity Profile Depicting the fully formed region:

45 deg 2.52m

 [V] CONCLUSION AND REMARKS:

Comparison of maximum velocities `V_(max)`

symm cop

Comparison of courant numbers and execution times for various wedge angles:

eec symm

 CONCLUSIONS:

1. SYMMETRY BOUNDARY CONDITION IS MORE ACCURATE THAN WEDGE BOUNDARY CONDITION.

2. EXECUTION TIME IS GREATER FOR SYMMETRY BOUNDARY CONDITION THAN WEDGE BOUNDARY CONDITION FOR THE SAME COMPUTATIONAL MESH.

3. COURANT NUMBER FOR SIMULATION INCREASES WITH INCREASING WEDGE ANGLES

4. TIME STEP OF THIS PROBLEM IS 0.01s WITH A WRITE TIME OF 1000 TIMESTEPS.


Projects by Priyotosh Bairagya

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 bounda 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