Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation

Difference Between wedge and symmetry boundary condition

in the last project, we have used the wedge boundary condition, which is used for 2-dimensional axisymmetric cases e.g. a cylinder, the geometry is specified as a wedge of small angle (e.g. < 5^o ) and 1 cell thick running along the plane of symmetry, straddling one of the coordinates planes. The axisymmetric wedge planes must be specified as separate patches of wedge type.

On the contrary, symmetry boundary conditions use for any (non-planar) patch which uses the symmetry plane slip condition. Usually, this boundary condition used on objects which are symmetrical or mirror by its various side. As we are using it for axisymmetrical pipe or object we can use any angle we want, and for a bigger angle, we can also discretize it along that angular portion which gives us finer results. This mostly considers 3D simulation conditions.

 

WEDGE BOUNDARY CONDITION USAGE AND APPLICATION:

• The wedge boundary condition (BC) is used in computational fluid dynamics to define an axisymmetry situation, for eg. Cylinder. The model and flow must be axisymmetry along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.
• This boundary condition is specified by two planes that must be selected on separate sides of the domain (referred to as front and back) on the surface running along the axis. 

 

SYMMETRY BOUNDARY CONDITION USAGE AND APPLICATION:

• This boundary condition is used to apply mirror-symmetry conditions to a structure.
• It can be applied to the faces of a structure and no other user input is needed. If a symmetry plane condition is applied to a face, the displacement of this face is locked in a normal direction but free to slide in tangential directions.

 

Now, let us compare the results obtained from symmetrical boundary condition simulation on the same previous problem.

Problem statement

For flow in pipe of diameter D, experimental observations show that for “fully developed” flow, laminar flow occurs when $R{e}_D$ < 2300. Therefore, let us consider Re for flow is 2000. Assume diameter of pipe (D)=0.02m and length (L) = 2.5 m. also consider working fluid as water at temperature 25^oC so the values of density (rho)=997kg/$m^3$ and Dynamic Viscosity (mu)=8.90X10-4Pa.sec.

                                                  

 

Pre-processing

• For this project, we are going to use the same procedure to simulate our results using symmetry boundary conditions. So, First start terminal in Linux Operating system using ctrl+alt+T buttons on the keyboard or right-click -> Open Terminal.
• Here, exactly as the last project, we have to simulate our problem in OpenFOAM and we are going to solve our problem using tutorials present in OpenFOAM, therefore, just type ‘tut’ to open the OpenFOAM tutorial directory. And type ‘ls’ to identify what kind of tutorials present inside this directory.
• We have to solve incompressible flow simulation for this problem. So select an incompressible directory by typing 'cd incompressible' and then again type 'ls' to explore folders in incompressible folders. And select icoFoam solver which is used for the laminar-viscous flow.
• icoFoam is the solver in OpenFOAM which is used as a transient solver for incompressible, laminar flow of Newtonian fluids.
• (Newtonian fluid is a fluid in which the viscous stress arising from its flow, at every point, are linearly correlated to the local strain rate - the rate of change of its deformation over time. Newtonian fluids are the simplest mathematical models of fluids that account for viscosity.)
• Therefore, we will select the icoFoam folder by typing 'cd icoFoam' and by typing 'ls' we explore this folder in the terminal.
• Here, we have, again to option, cavity and elbow. They both are part of the icoFoam solver and cavity section used to create a required domain and black mesh so select the 'cavity' folder.
• In the cavity, we are, again, going to see the cavity folder, where we have to define our fluid domain and generate block mesh.
• We are going to copy this cavity folder using the 'cp -r cavity' command. Next, rename it as ‘Pipe(45)’, as we did for the ‘Pipe’ folder in our last project, and paste it in OpenFOAM run folder by typing 'cp -r cavity $FOAM_RUN/Pipe(45)'.
• And then open the run folder by typing 'cd $FOAM_RUN'. where we can see the 'Pipe(45)' file.
• Now, we can make any changes to the Pipe file, which will never affect the OpenFOAM official files. Open the Pipe folder. Here, we can provide specifications for our block mesh in the system folder.
• There are four script file system folders and by opening blockMeshDict, we can create a domain and provide the required specification for block mesh. 
• Here have to create a blockMeshDict file with Matlab and replace that with the original one. But we already have the program file for wedge boundary condition so we just edit that program for symmetry boundary condition.

 

FILE WRITING (in Matlab) for symmetry  boundary condition

But, this time, we are going to create blockMeshDict file using Matlab coding by file parsing techniques. Prior, lest define all the given data first in our main program

Condition of working fluid according to problem statement.

Reynolds number ($R_e$) = 2000

Diameter of the pipe (D) = 0.02m

Length of the pipe (L) = 2.5m

Density of working fluid ($\rho$) = 997 kg/$m^3$

Dynamic viscosity of working fluid ($\mu$) = 8.9X10^-4 Pa.sec

%% Input parameters
% Reynold's Number
Re = 2000;
% Total Length of pipe(m)
L = 2.5;
% Diameter of pipe in meters
D = 0.02;
% Wedge Angle in degrees
th = input('<0> Wedge Angle in degrees - ');
% Density (kg/m^3)
rho = 997;
% Dynamic Viscosity (Pa.sec)
mu = 8.90e-4;
% Radius of pipe(m)
R = D/2;

 

Now, let us create separate function file to write blockMeshDict file in Matlab program. Create blockMeshDict file and start writing it using following ‘fopen’ command.

function[f1] = blockMesh(R,th,L,nx,ny,nz)

%% Writing BlockMeshDict file for openFoam
f1 = fopen('blockMeshDict','w');

 

Before defining any data we are going to state with printing header information in the form of comment which is not executable in OpenFoam

fprintf(f1,"/*--------------------------------*- C++ -*----------------------------------*\n");
fprintf(f1,'  =========                 | \n');
fprintf(f1,'  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox \n');
fprintf(f1,'   \\    /   O peration     | Website:  https://openfoam.org \n');
fprintf(f1,'    \\  /    A nd           | Version:  8 \n');
fprintf(f1,'     \\/     M anipulation  |\n');
fprintf(f1,"\*---------------------------------------------------------------------------*/\n");

 

Next, we set units in meter,  by converting all type of length in to meter.

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

 

Vertices:-

Afterward, in codes, we provide list of vertices by providing x,y,and z coordinates locations; its always start from vertex (0,0,0). so we have to create following geometry.

                                                                                   

 

Mesh generation

The mesh of the Pipe case includes a new feature: a block with fewer than 8 vertices will be used. There is an axi-symmetry in the case so it is not necessary to mesh the whole pipe but only a wedge of it. Since $\frac{\partial }{\partial \theta }=0$ for the whole domain, the results of the simulation will be the same while the number of elements of the mesh will be lower. Following figure represents a schematic view of the geometry  that it is going to be created:

                                                  

 

First of all, it is important to comprehend the collapsing of the vertices. It can be seen that vertices 0 and 3 are repeated in an adequate position within the parentheses to indicate that the two vertices that would be occupying this position have been collapsed.

Secondly, in the patches definition, the same vertices 0 and 3 are joined creating an empty patch forming the axis of the pipe. The patch type wedge is used in two faces to indicate that an axi-symmetry exists and there are no physical walls.

fprintf(f1,'vertices(\n');
fprintf(f1,'    (0 0 0)\n');
fprintf(f1,('    (%f %f 0)\n'),R*cosd(th/2),R*sind(th/2));
fprintf(f1,('    (%f %f %f)\n'),R*cosd(th/2),R*sind(th/2),L);
fprintf(f1,('    (%f 0 0)\n'),L);
fprintf(f1,('    (%f %f 0)\n'),R*cosd(th/2),-R*sind(th/2));
fprintf(f1,('    (%f %f %f)\n);\n'),R*cosd(th/2),-R*sind(th/2),L);

 

Next, we are going to define block for meshing. If we look at he mesh format first letter, in block parenthesis, is ‘hex’ which define type of mesh which is hexahedral mesh.

Then, again in parenthesis, we provide vertices number according to order which follow right hand thumb rule, meaning, thumb is the z axis and order of vertices follow rotational order of our other finger.

Then comes, number of mesh grid that we want to divide along x, y, and z length. We define it as (nx,ny,nz) format in main program.

Given mesh specification:

1. Number of cells along the z-direction(longer dimension) = 500
2. Number of cells along with the longer x-direction = 20
3. Number of the cell for y-direction =theta/5 (vary according to theta angle)

 

% Number of cells in x,y,z Directions
nx = 20;
ny = th/5;
nz = 500;

Lastly, we give simple grading which is cell expansion ratios

 

Grading factor or expansion ratios:-

The expansion ratios enables the mesh to be graded, or refined, in specified directions. The ratio is that of the width of the end cell ${\delta }_e$ along one edge of a block to the width of the start cell ${\delta }_s$ along the edge as shown in below figure.

                                                                                   

In simple format,

                                       Expansion ratio = (size of end cell) / (size of starting cell)

The simple description specifies uniform expansion in the local ,${)}_2$and $x_3$directions respectively with only 3 expansion ratio. So here, we have 1 wedge shape block so we have to define hex mesh for 1 block.

There is a grading towards the wall to compute a more exact value of wall shear stress and towards the inlet where the flow is not fully developed and therefore the velocity profile changes.

fprintf(f1,('hex (0 4 1 0 3 5 2 3) (%d %d %d) simpleGrading (0.2 1 1)\n);\n'),nx,ny,nz);

 

In this case, for edge 4-1 and 5-2 we have radius curve so we will define this in edges section following way.

fprintf(f1,'edges(\n');
fprintf(f1,'    arc 4 1 (0 %f 0)\n',R);
fprintf(f1,'    arc 5 2 (%f %f 0)\n);\n',L,R);

 

Next, define boundary of the problem, in this we have to define faces of the our 3D block this should be in order according to right hand thumb rule. So set faces for this problem following way.

fprintf(f1,'boundary(\n');
fprintf(f1,'     inlet\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type patch;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 4 1 0)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     outlet\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type patch;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (3 5 2 3)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     front\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type symmetry;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 1 2 3)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     back\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type symmetry;\n');
fprintf(f1,'           faces \n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 3 5 4)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     pipewall\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type wall;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (1 4 5 2)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     axis\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type empty;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 3 3 0)\n');
fprintf(f1,'           );\n');
fprintf(f1,'      }\n');
fprintf(f1,');\n');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
fprintf(f1,'// ************************************************************************* // \n');

 After running this MATLAB code for different angles 

Copy blockMeshDict file and replace it in the system directory 

                                     

Furthermore, we are editing the ‘controlDict’ file to provide time limits - start time and end time. Also to represent $\delta t$. While the selection of $\delta t$ CFL value should be less than 1. otherwise, the solution could be stopped forcibly or it may diverge at the time of converging. So make the following changes in the ‘controlDict’ file.

 

Boundary and initial conditions

Then, we have to define the initial value for pressure and velocity. In order to do that, we have to go back to Pipe(45) directory using the ‘cd ..’ command and select a folder named ‘0’, which stands for initial values. In this case initial value for pressure and velocity.

When solving circulating internal flow, the most common boundary conditions are to fix an inlet velocity and an outlet pressure.

And for symmetry boundary condition just define to side faces to symmetry.

Initial setting for pressure

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

    back
    {
        type            symmetry;
    }

    axis
    {
        type            empty;
    }

    pipewall
    {
        type            zeroGradient;
    }
}

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

 

In the U dictionary, one can find out a new type of boundary condition: inlet and outlet. It switches U and p between fixed value and zero gradients: the first when the flow is in-going and the second when it is outgoing.

For example: in a circular pipe case, the outlet velocity will always be adapted to the conditions according to zero gradients except if there is an inflow at the outlet patch. If this happens, to this inflow velocity it is assigned the value specified under the inlet-outlet instruction and so the flow will always be outgoing.

Lastly, define side faces to symmetry instead of wedge boundary condition.

Initial setting for Velocity

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

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.09373);
    }

    outlet
    {
        type            zeroGradient;
    }
    
    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

    axis
    {
        type            empty;
    }

    pipewall
    {
        type            noSlip;
    }
}

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

 

Next, we are going to keep the same dynamic viscosity as we kept in the previous project of wedge boundary condition.

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


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

 

After deploying mesh, now we are ready to run the simulation. Like last time, this time also we are running the simulation for viscus-laminar flow that's why we will choose our solver as icoFoam. Hence, just type ‘icoFoam’ and simulation will start. 

results are same to all different angles 

 

                                                          Hydro-dynamic length

                                                      PRESSURE ALONG THE PIPE 

                                                                   SHEAR STRESS NEAR THE WALL

                                        velocity profile at length 2.5m 

FILE WRITING (in Matlab) for wedge boundary condition

function[f1] = blockMesh(R,th,L,nx,ny,nz)

%% Writing BlockMeshDict file for openFoam
f1 = fopen('blockMeshDict','w');

%% Input parameters
% Reynold's Number
Re = 2100;
% Total Length of pipe(m)
L = 3.1;
% Diameter of pipe in meters
D = 0.02;
% Wedge Angle in degrees
th = input('<0> Wedge Angle in degrees - ');
% Density (kg/m^3)
rho = 997;
% Dynamic Viscosity (Pa.sec)
mu = 8.90e-4;
% Radius of pipe(m)
R = D/2;
% Number of cells in x,y,z Directions
nx = 20;
ny = 1;
nz = 500;

fprintf(f1,"/*--------------------------------*- C++ -*----------------------------------*\n");
fprintf(f1,'  =========                 | \n');
fprintf(f1,'  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox \n');
fprintf(f1,'   \\    /   O peration     | Website:  https://openfoam.org \n');
fprintf(f1,'    \\  /    A nd           | Version:  8 \n');
fprintf(f1,'     \\/     M anipulation  |\n');
fprintf(f1,"/*---------------------------------------------------------------------------*/\n");
fprintf(f1,'FoamFile \n');
fprintf(f1,'{ \n');
fprintf(f1,'    version     2.0;\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');
fprintf(f1,'vertices\n');
fprintf(f1,'( \n');
fprintf(f1,'    (0 0 0)\n');
fprintf(f1,('    (%f %f 0)\n'),R*cosd(th/2),R*sind(th/2));
fprintf(f1,('    (%f %f %f)\n'),R*cosd(th/2),R*sind(th/2),L);
fprintf(f1,('    (0 0 %f)\n'),L);
fprintf(f1,('    (%f %f 0)\n'),R*cosd(th/2),-R*sind(th/2));
fprintf(f1,('    (%f %f %f)\n);\n'),R*cosd(th/2),-R*sind(th/2),L);
fprintf(f1,'blocks\n');
fprintf(f1,'( \n');
fprintf(f1,('hex (0 4 1 0 3 5 2 3) (%d %d %d) simpleGrading (0.2 1 1)\n);\n'),nx,ny,nz);
fprintf(f1,'edges\n');
fprintf(f1,'( \n');
fprintf(f1,'    arc 4 1 (0 %f 0)\n',R);
fprintf(f1,'    arc 5 2 (%f %f 0)\n);\n',L,R);
fprintf(f1,'boundary\n');
fprintf(f1,'( \n');
fprintf(f1,'     inlet\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type patch;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 4 1 0)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     outlet\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type patch;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (3 5 2 3)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     front\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type wedge;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 1 2 3)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     back\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type wedge;\n');
fprintf(f1,'           faces \n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 3 5 4)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     pipewall\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type wall;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (1 4 5 2)\n');
fprintf(f1,'           );\n');
fprintf(f1,'     }\n');
fprintf(f1,'     axis\n');
fprintf(f1,'     {\n');
fprintf(f1,'           type empty;\n');
fprintf(f1,'           faces\n');
fprintf(f1,'           (\n');
fprintf(f1,'               (0 3 3 0)\n');
fprintf(f1,'           );\n');
fprintf(f1,'      }\n');
fprintf(f1,');\n');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
fprintf(f1,'// ************************************************************************* // \n');
fclose(f1);


end

 

Furthermore, we are going to edit ‘controlDict’ file to provides time limits - start time and end time. Also to represent  deltat. While selection of deltat CFL value shoud be less than 1. otherwise, the solution could stopped forcibly or it may diverge at the time of converging. So make following changes in ‘controlDict’ file.

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

writeControl    timeStep;

writeInterval   10;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

 

Boundary and initial conditions

Then, we have to define initial value for pressure and velocity. In order to do that, we have to go back in cavity_copy directory using ‘cd ..’ command and select folder named ‘0’, which stand for initial values. In this case initial value for pressure and velocity.

When solving circulating internal flow, the most common boundary conditions are to fix an inlet velocity and an outlet pressure.

Initial setting for pressure

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

    axis
    {
        type            empty;
    }

    pipewall
    {
        type            zeroGradient;
    }
}

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

 

In the U dictionary one can find out a new type of boundary condition: inlet and outlet. It switches U and p between fixedValue and zeroGradient: the first when the flow is in-going and the second when it is outgoing.

Example: in circular pipe case, the outlet velocity will always be adapted to the conditions according to zeroGradient except if there is an inflow at the outlet patch. If this happens, to this inflow velocity it is assigned the value specified under the inletOutlet instruction and so the flow will always be outgoing.

Initial setting for Velocity

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

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.1);
    }

    outlet
    {
        type            zeroGradient;
    }
    
    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

    axis
    {
        type            empty;
    }

    pipewall
    {
        type            noSlip;
    }
}

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

 

Next, we have to mention dynamic viscosity in transposeProperty in constant directory of the Pipe.

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


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

 

 

                                                          Hydro-dynamic length

                                                      PRESSURE ALONG THE PIPE 

                                                                   SHEAR STRESS NEAR THE WALL

                                                             velocity profile

Comparing the both wedge and symmetry we can see that there is no chnage in the results so that we can compare the any of one to the analyatical values 

ANALYTICAL CALCULATIONS

Average Velocity:-

 

% Average Velcity
V_avg = (Re*mu)/(rho*D);

 

 

Maximum velocity:-

$v_{max}=2\cdot v_{avg}$

$v_{max}=2\cdot 0.0937$

$v_{max}=0.1874$

 

Pressure drop across the length of pipe:-

$\Delta P=\frac{32\mu L{V}_{avg}}{D^2}$

length L = 0.5

$\Delta p=3.5$

 

Wall Shear stress

${\tau }_w=\frac{4\cdot \mu \cdot v_{avg}}{R}$

${\tau }_w=0.0356$

 

Velocity at various radius(r) near fully developed flow region

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

% ANALYTICAL
R = 0.01;
r=linspace(-R,R,1000);

V_avg = (Re*mu)/(rho*D);

% Velocity_analytical
V_A=2*V_avg*(1-(r.^2/R^2));

Analytical Hydro-dynamic length

$\frac{l_e}{D}=0.06\cdot R_e$

 

where 

   entrance length'

D  Diameter

 renolds number

 

 

Comparing the NUMERICAL AND ANALYTICAL CALCULATIONS

comparing the velocity profile using MATLAB

clear all
close all
clc


% Reynold's Number
Re = 2100;
% Total Length of pipe(m)
L = 3.1;
% Diameter of pipe in meters
D = 0.02;
% Density (kg/m^3)
rho = 997;
% Dynamic Viscosity (Pa.sec)
mu = 8.90e-4;
% Radius of pipe(m)
R = D/2;

% ANALYTICAL
R = 0.01;
r=linspace(-R,R,1000);

V_avg = (Re*mu)/(rho*D);

% Velocity_analytical
V_A=2*V_avg*(1-(r.^2/R^2));


%NUMERICAL

v_r = csvread("v_p.csv",1,1);
c = csvread("v_p.csv",1);


for i = 1:1000
    
    %velocity_numerical
    V_n(i) = v_r(i)* (1-(c(i)^2));
    
end


figure(1)
plot(r,V_A,'r','linewidth',1)
hold on
plot(r,V_n,'--','color','k','linewidth',2.5)
xlabel('<---- Radius of pipe in m ---->','fontweight','bold','fontsize',14)
ylabel(' Velocity in m/s ------>','fontweight','bold','fontsize',14)
legend(' Analytical method',' Numerical method')
grid on

Below graph shows the velocity profile comparison between  NUMERICAL AND ANALYTICAL

Validate Hydro-dynamic length with the numerical result

ANALYTICAL value is 2.52 m

NUMERICAL value is 0.08 m

	ANALYTICAL	NUMERICAL
	0.1875	0.177416
for length 0.5m	3.5	3.8185
	0.0356	0.0301

In symmetry BC if we change the cells along the Y axis to 3 cells there will be change in the results along the y direction and there will be change in the velocity profile