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

Comparison of simulation results for flow through pipe with different boundary condition using OpenFoam

SSmitesh Badnapurkar26 Jan 2021
CFDMATLAB

Aim:-> To obtain and compare simulation result with "symmetry" boundary condition and "wedge" boundary condition at different wedge angles (10, 25 & 45 degrees).

Theory:-

As we are considering circular pipe which is symmetric about central axis so instead of taking whole pipe cross section let us consider wedge (sector) of pipe only. Taking such sector and simulating it by using symmetrical boundary condition reduces cost and computation time.

Hagen–Poiseuille law:

In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.

Poiseuille's Equation describes the pressure drop due to the viscosity of the fluid. Other types of pressure drops may still occur in a fluid. The assumptions for this concept are- the fluid must be incompressible and Newtonian; the flow should be laminar through a pipe; and pipe should be of constant circular cross-section and there must be no acceleration of fluid in the pipe.

`Deltap=(8muLQ)/(piR^4)`

Where,

`Deltap=` Pressure difference between inlet and outlet,

`L=` Total length of pipe,

`mu=` Dynamic viscosity,

`Q=` Volumetric flow rate,

`R=` Radius of pipe cross-section

As per Darcy-Weisbach relation,

`R_e=(rhovd)/mu`

Where,

`rho=` Density of fluid,

`v=` Velocity of fluid,

`d=` Diameter of pipe.

Let us assume,

As flow is laminar,

Reynolds number`=R_e=2100`

Diameter of pipe`=D=0.01`

Radius of pipe`=R=D/2=0.005`

L=1.5

Consider fluid as water at 25 Degree celsius

Density`=rho=997`

Dynamic viscosity`=mu=0.89xx10^(-3)`

Kinematic viscosity`=nu=0.8297xx10^(-6)`

`theta=3^o`

Hence we can calculate,

Average fluid velocity`=v_(avg)=(mu*R_e)/(rhoD)=0.18746`

Maximum fluid velocity`=v_(max)=2*v_(avg)=0.37492`

Pressure difference`=Deltap = (8muL*v_(max))/(R^2)=160.17`

Kinematic pressure difference `=DeltaP_k = Delta_p/rho=0.16065`

MATLAB code:-

clear all
close all
clc

Re=2100; %Reynolds number
D=0.01; %Diameter of cross section of pipe in m
R=D/2; %Radius of cross section of pipe in m
%Entrance length calculation 
%(Length required to developed velocity profile fully)
L_e=0.06*Re*D; %in m
%Assumed total length of pipe (muxst be > L_e)
L=1.5; %in m
%Properties of water at 25 Degree celcius
rho=997; %Density in Kg/m^3
mu=0.89e-3; %Dynamic viscosity
nu=0.8297e-6; %Kinematic viscosity

%In degrees
printf('Choose angle for sector "theta"');
printf('n');
theta= input('Angle is ');

%Boundary conditions
%Average velocity
v_avg=mu*Re/(D*rho);
%Max velocity
v_max=2*v_avg;
%Pressure difference
del_p = (8*mu*L*v_max)/(R^2);
%Kinematic pressure difference
kinematic_p = del_p/rho;

%Creating vertex
v_0 = [0 0 0];
v_1 = [0 (R*cosd(theta/2)) -(R*sind(theta/2))];
v_2 = [0 (R*cosd(theta/2)) (R*sind(theta/2))];
v_3 = [L 0 0];
v_4 = [L (R*cosd(theta/2)) -(R*sind(theta/2))];
v_5 = [L (R*cosd(theta/2)) (R*sind(theta/2))];;
block = 'hex (0 3 4 1 0 3 5 2) (200 30 1) simpleGrading (1 1 1)';

%Writing blockmesh file
f1 = fopen('blockMeshDict.txt','w');

%Header 
h1 = '/*----------------------------*- C++ -*------------------------------*'; 
h2 = '| =========                 |                                         |';
h3 = '|       /  F ield         | OpenFOAM: The Open Source CFD Toolbox   |';
h4 = '|      /   O peration     | Website:  https://openfoam.org          |';
h5 = '|     /    A nd           | Version:  8                             |';
h6 = '|    /     M anipulation  |                                         |';
h7 = '*-------------------------------------------------------------------*/';

fprintf(f1,'%s n',h1);
fprintf(f1,'%s n',h2);
fprintf(f1,'%s n',h3);
fprintf(f1,'%s n',h4);
fprintf(f1,'%s n',h5);
fprintf(f1,'%s n',h6);
fprintf(f1,'%s n',h7);

fprintf(f1, 'nFoamFilen');
fprintf(f1, '{ntversion           2.0;');
fprintf(f1, 'ntformat            ascii;');
fprintf(f1, 'ntclass       	     dictionary;');
fprintf(f1,	'ntobject            blockMeshDict;n}n');
fprintf(f1, '-------------------------------------------------------------------//');

%To convert values into SI unit (meter)
fprintf(f1, 'nconvertToMeters 1;n');

%Vertices definition of sector
fprintf(f1, 'nverticesn(n');
fprintf(f1, 't (%f  %f  %f)n',v_0(1),v_0(2),v_0(3));
fprintf(f1, 't (%f  %f  %f)n',v_1(1),v_1(2),v_1(3));
fprintf(f1, 't (%f  %f  %f)n',v_2(1),v_2(2),v_2(3));
fprintf(f1, 't (%f  %f  %f)n',v_3(1),v_3(2),v_3(3));
fprintf(f1, 't (%f  %f  %f)n',v_4(1),v_4(2),v_4(3));
fprintf(f1, 't (%f  %f  %f)n);n',v_5(1),v_5(2),v_5(3));

%Block definition of sector
fprintf(f1, 'nblocksn(n');
fprintf(f1, 't%sn);',block);

fprintf(f1, 'nedgesn(n');
fprintf(f1, 't arc 1 2 (%d %d %d)n',0,R,0);
fprintf(f1, 't arc 4 5 (%d %d %d)n);',L,R,0);

%Boundary condition part
fprintf(f1, 'nboundarynn(n');

fprintf(f1, 'tinletnnt{n');
fprintf(f1, 'tttype patch;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 2 1 0)ntt);');
fprintf(f1, 'nt}nn');

fprintf(f1, 'toutletnnt{n');
fprintf(f1, 'tttype patch;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (3 4 5 3)ntt);');
fprintf(f1, 'nt}nn');

fprintf(f1, 'tfrontnnt{n');
fprintf(f1, 'tttype symmetry;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 3 5 2)ntt);');
fprintf(f1, 'nt}nn');

fprintf(f1, 'tbacknnt{n');
fprintf(f1, 'tttype symmetry;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 1 4 3)ntt);');
fprintf(f1, 'nt}nn');

fprintf(f1, 'twallnnt{n');
fprintf(f1, 'tttype wall;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (1 2 5 4)ntt);');
fprintf(f1, 'nt}nn');

fprintf(f1, 'taxisnnt{n');
fprintf(f1, 'tttype empty;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 3 3 0)ntt);');
fprintf(f1, 'nt}n');

fprintf(f1, ');nn')

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

fclose(f1);

After running above MATLAB code following blockMeshDict file generates. In above code we have to provide input of "theta" which is wedge (sector) angle.

Blockmeshdict for 10degree:

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

\\-------------------------------------------------------------------//

convertToMeter 1;

vertices
(
	 (0.000000  0.000000  0.000000)
	 (0.000000  0.004981  -0.000436)
	 (0.000000  0.004981  0.000436)
	 (1.500000  0.000000  0.000000)
	 (1.500000  0.004981  -0.000436)
	 (1.500000  0.004981  0.000436)
);



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

edges

(
	 arc 1 2 (0 0.005 0)
	 arc 4 5 (1.5 0.005 0)
);

boundary

(
	inlet

	{
		type patch;
		faces
		(
			 (0 2 1 0)
		);
	}

	outlet

	{
		type patch;
		faces
		(
			 (3 4 5 3)
		);
	}

	front

	{
		type symmetry;
		faces
		(
			 (0 3 5 2)
		);
	}

	back

	{
		type symmetry;
		faces
		(
			 (0 1 4 3)
		);
	}
      

       wall

	{
		type wall;
		faces
		(
			 (1 2 5 4)
		);
	}


	axis

	{
		type empty;
		faces
		(
			 (0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

Blockmeshdict for 25degree:

vertices
(
	 (0.000000  0.000000  0.000000)
	 (0.000000  0.004881  -0.001082)
	 (0.000000  0.004881  0.001082)
	 (1.500000  0.000000  0.000000)
	 (1.500000  0.004881  -0.001082)
	 (1.500000  0.004881  0.001082)
);

Blockmeshdict for 45degree:

vertices
(
	 (0.000000  0.000000  0.000000)
	 (0.000000  0.004619  -0.001913)
	 (0.000000  0.004619  0.001913)
	 (1.500000  0.000000  0.000000)
	 (1.500000  0.004619  -0.001913)
	 (1.500000  0.004619  0.001913)
);

After generating mesh we need to imply initial conditions of pressure and velocity as per problem.

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

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

internalField   uniform (0 0 0);

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

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }
    
    wall
    {
        type            noSlip;
    }

    axis
    {
        type            empty;
    }

}

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

Initial pressure conditions:-

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

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

    wall
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
}

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

Transport properties:-

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


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

After defining all boundary conditions, initial conditions and meshing we need to initiate solver to solve problem. Following controlDict file controls the simulation by timesteps.

Total time for simulation = (Length of pipe)/(Velocity of fluid).

controlDict code:-

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   40;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Meshed geometry:-

Above figure shows meshed geometry for above meshdict code file.

Result:-

Analytical value for velocity = 0.37492 (m/s^2)

Error in velocity value is calculated by formula,

`"Error "=("Theoretical value "-" Numerical value")/"Theoretical value"xx100`

Boundary conditions courant number Execution time (S) velocity (m/s) Error
(%)
mean max
`10^o` symmetry 0.025492 0.05731 1831 0.372544 0.633735
Wedge 0.025492 0.05736 1335 0.372544
`25^o` symmetry 0.025492 0.057936 1714 0.372836 0.555852
Wedge 0.025492 0.057949 1345 0.372836
`45^o` symmetry 0.025492 0.059629 1589 0.373447 0.392884
Wedge 0.025492 0.059536 1323 0.373447

Conclusion:- From above results it is clear that-

  1. For solving with "symmetry" boundary condition required higher computation time than "wedge" boundary condition.
  2. As the wedge angle increases error decreases.
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.