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

Simulation of Flow through a pipe in OpenFoam

SShubham Patil19 May 2020
CFD

Question :

  • Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes
    • Hint: wedge angle should be less than 5 degrees
    • Make sure that number of cells along the "axi symmetric"
  • 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 and shear stress.

Answer :

Given :-

  • Working fluid - water
  • Wedge angle should be less than 5 degrees
  • Reynolds number 2100.

1] Assumptions:-

  • D = 0.01m .....Diameter
  • R = 0.005m .....Radius
  • mu = 1.002e-3 .....In SI units (water at 20 degree celsius)
  • rho = 998.21 .....In SI units (water at 20 degree celsius)
  • Re = 2100 .....Reynolds number
  • Assumed Temperature of water = 20 degree C
  • Wedge angle = 4

When the fluid enters the pipe its boundary layer increases till the centre the pipe and fills the pipe,this region is entrance region where the flow is developing gradually and the remaining length is fully developed length.

Entrance length is a length which a flow travels before it becomes fully developed flow. When a boundary layer expands to fill pipe, the flow starts to develop. In the the above figure we can see that the flow is developing in the Entrance region and after the entrance region it gets fully developed and thereafter there are no changes in the flow characteristics. Viscous forces are resposible for the velocity profiles.

Entrance Length of Laminar flow pipe is calculated by (Le/D) = (0.06*Re)

2] Analytical values:-

  • L_e = 0.06*Re*D ..... Entrance Length
  • L_t = L_e + 1.24 ..... Total Length
  • v_avg = (Re*mu)/(rho*D) ..... Velocity
  • v_max = 2*v_avg ...... Max Velocity
  • Pressure Drop = ( 32*mu*v_avg*L_t ) / (D^2)
  • Kinematic pressure = Pressure drop/rho.

3] Steps:-

  • First we write the code which generates the blockmeshDict file for any wedge angle and grading schemes.
  • Inwhich we input the analytical fomulaes and then parse it.
  • Then we copy the cavity directory from FoamTUTORIALS to FOAM_RUN and rename it.
  • after getting the blockmeshDict file we copy it in the system directory.
  • Run the blockMesh command in new renamed directory.
  • Then we apply Initial and Boundary conditions for p and U also modify the transproperties according to working fluid's propeties.
  • Then modify the controlDict file present in system directory, accordingly.
  • Run the IcoFoam command to solve.
  • Then use ParaFoam command for postprocessing for velocity and pressure profiles.

Below are the matlab codes for the blockmeshDict file and also updated p, U, controlDict, Transportproperties files.

4] Code

  • Matlab code for blockMeshDict :
clear all
close all
clc

% Inputs

Re = 2100;
D = 0.01;
R = D/2;
theta = input('Enter value for Wedge angle =');
mu = 1.002e-3;
rho = 998.21;
nu = mu/rho;

% Analytical Values 

L_e = 0.06*Re*D;
L_t = L_e + 1.24;
v_avg = (Re*nu)/D;
v_max = 2*v_avg;
Delta_P = (32*mu*v_avg*L_t)/D^2;
kinematic_P = Delta_P/rho;

% For finding 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_t 0 0];
v4 = [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5 = [L_t R*cosd(theta/2) R*sind(theta/2)];

% Script for blockmeshDict Directory

f1 = fopen('blockMeshDict.txt','w')

fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n=========                 |\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:  7\n');
fprintf(f1,'    \\/     M anipulation  |\n');
fprintf(f1,'%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\nFoamFile\n');

b1 = blanks(4);
b2 = blanks(5);
b3 = blanks(6);
b4 = blanks(7);

fprintf(f1,'{\n');
fprintf(f1,'%s','    version',b2,'2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s','    format',b3,'ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s','    class',b4,'dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s','    object',b3,'blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'}\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

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

% for Vertices

fprintf(f1,'%s','vertices');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,b1,'%s','(');
v_0 = fprintf(f1,'(%f %f %f)\n',v0(1),v0(2),v0(3));
fprintf(f1,b1,'%s','(');
v_1 = fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b1,'%s','(');
v_2 = fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b1,'%s','(');
v_3 = fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b1,'%s','(');
v_4 = fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b1,'%s','(');
v_5 = fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');

fprintf(f1,'\n');
fprintf(f1,'%s','blocks');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,b1,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');

% Mesh 
m1 = input('Enter no.of mesh elements along X-axis =');
m2 = input('Enter no.of mesh elements along Y-axis =');
m3 = input('Enter no.of mesh elements along Z-axis =');

fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%2.0f %1.0f %3.0f)', m1, m2, m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');

% edges
fprintf(f1,'\n');
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s   arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s   arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');

% boundary conditions

% For Inlet
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
fprintf(f1,'%s','   inlet');
fprintf(f1,'\n');
fprintf(f1,'   {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n   }\n');

% for Outlet
fprintf(f1,'%s','   outlet');
fprintf(f1,'\n');
fprintf(f1,'   {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n   }\n');

% for front
fprintf(f1,'%s','   front');
fprintf(f1,'\n');
fprintf(f1,'   {\n');
fprintf(f1,'%s',blanks(10),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(1 0 3 4)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n   }\n');

% for back
fprintf(f1,'%s','   back');
fprintf(f1,'\n');
fprintf(f1,'   {\n');
fprintf(f1,'%s',blanks(10),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n   }\n');

% for fixedWalls
fprintf(f1,'%s','   fixedWalls');
fprintf(f1,'\n');
fprintf(f1,'   {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n   }\n');

% for axis
fprintf(f1,'%s','   axis');
fprintf(f1,'\n');
fprintf(f1,'   {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n   }\n');

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

fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');

fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');


fclose(f1)

After running the above Code in matlab, we get the following result as the BlockMeshDict which we copy in the system directory and replace the old one to get the meshed wedge geometry.

  • blockMeshDict :
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
    (0.000000 0.000000 0.000000)
    (0.000000 0.004997 -0.000174)
    (0.000000 0.004997 0.000174)
    (2.500000 0.000000 0.000000)
    (2.500000 0.004997 -0.000174)
    (2.500000 0.004997 0.000174)
);

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

edges
(
    arc 1 2 (0 0.005000 0)
    arc 4 5 (2.500000 0.005000 0)
);

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

mergePatchPairs
(
);

// ************************************************************************* //
  • ControlDict :
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.5;

deltaT          0.005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //
  • p :
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     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.1693;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
}

// ************************************************************************* //
  • U :
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     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.2108 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

    fixedWalls
    {
        type            noSlip;
    }

    axis
    {
        type            empty;
    }
}

// ************************************************************************* //
  • TransportProperties :
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 1.00379e-6;


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

Below are the Outputs of post processing the solution.

5] OutPuts

  • Mesh Output

This is the mesh Output generated for the geometry with wedge angle of 4 degrees.

  • Velocity Distribution at inlet of pipe
  • Velocity Distribution at outlet of pipe

Below are the Velocity Profile Outputs at various positions of length of pipe.

  • Velocity Profile Outputs

i) for x = 0.01m from inlet

ii) for x = 0.05m from inlet

iii) for x = 0.08m from inlet

vii) for x = 1m from inlet

viii) for x = 1.26m from inlet

ix) for x = 2m from inlet

x) for x = 2.5m from inlet

  • From the above outputs we can see that the profile is gradually developing w.r.t length of pipe.
  • we can see that at 0.01m from inlet the profile is not fully developed as because no slip conditions are imposed on walls but later in increased length the flow gets developed.
  • At 1.26m the profile is fully developed as it is the entrance length distance after which starts the fully developed regions generating fully developed profiles as seen in later output of 2m and 2.5m length.
  • The max Velocity value gets closer to analytical values as the profile develops
  • Henceforth it can be said that Hp equation is applicable to pipe.
  • Pressure Distribution for entire Length of Pipe
  • Pressure Profile
  • Above Output is the pressure Distribution output in which we can see that presuure is decreasing with increase in length of pipe.
  • Also the pressure is non linear at the starting point of pipe, due to higher shear stress at that region.
  • Shear Stress Profile

From the above plot we ca see that the shear stress is zero at the centre , whereas higher at the walls

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.