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

Simulation of flow through a pipe using OpenFoam - Part 2

SSwaroop B R1 Feb 2021
CFDMATLAB

Simulation of flow through a pipe for different values of theta ( 10, 25, and 45) in degrees

1. Wedge

for angle theta = 10

MATLAB Coding

% Simulation of flow through a pipe using OpenFoam
% Creating blockMeshDict file 

clear all
close all
clc

% Defining variables and assumptions
% Diameter of pipe
D = 0.02;

% Radius 
r =D/2;

% Dynamic viscosity N/m^2
mu = 0.00089;

% Density of fluid kg/m^3
rho = 997;

% Reynold's Number for laminar flow
Re = 2100;

% Entry length of pipe 
L = 0.06*Re*D;

% Total lenghth of pipe
L_p = L+0.3;

% Wedge angle 
disp('Enter the value for wedge angle')
theta = input('theta = ');

% Hagen Poiseuill's law

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

% Maximum velocity
V_max = 2*V_avg;

% Pressure difference 
del_P = (8*mu*L*V_avg)/(r^2);

% Kinematic pressure difference
K_P = del_P/rho;

% Radial distance
R = linspace(0,r,1000);
% Calculating shear stress (wall shear stress)
for i = 1:length(R)
    tau = (2*mu*V_max*R)/(r^2);
end

% Now creating blockMeshDict file for OpenFoam to generate the geometry

% Adding headers
header1 = '/*--------------------------------*- C++ -*----------------------------------*';
header2 = '  =========                 |';
header3 = '  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
header4 = '   \    /   O peration     | Website:  https://openfoam.org';
header5 = '    \  /    A nd           | Version:  8';
header6 = '     \/     M anipulation  |';
header7 = '*---------------------------------------------------------------------------*/';
header8 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
header9 = '// ************************************************************************* //';
% Creating vertices
v0 = [0 0 0]; %point 0
v1 = [ 0 r*cosd(theta/2) -r*sind(theta/2)]; %point 1
v2 = [ 0 r*cosd(theta/2) r*sind(theta/2)]; %point 2
v3 = [L_p 0 0]; %point 3
v4 = [L_p r*cosd(theta/2) -r*sind(theta/2)]; %point 4
v5 = [L_p r*cosd(theta/2) r*sind(theta/2)]; %point 5

% Creating file and writing
f1 = fopen('blockMeshDict','w')

fprintf(f1, '%sn', header1)
fprintf(f1, '%sn', header2)
fprintf(f1, '%sn', header3)
fprintf(f1, '%sn', header4)
fprintf(f1, '%sn', header5)
fprintf(f1, '%sn', header6)
fprintf(f1, '%sn', header7)

fprintf(f1, 'FoamFilen{n')
fprintf(f1, '%12s t 2.0;n','version')
fprintf(f1, '%12s t ascii;n','format')
fprintf(f1, '%12s t dictionary;n','class')
fprintf(f1, '%12s t blockMeshDict','object')
fprintf(f1, '}n')
fprintf(f1, '%snn', header8)

fprintf(f1,'%s 1;nn','covertToMeters')

% Creating vertices
fprintf(f1, '%sn','vertices')
fprintf(f1, '(n')
fprintf(f1, 't(%d %d %d)n',v0)
fprintf(f1, 't(%d %d %d)n',v1)
fprintf(f1, 't(%d %d %d)n',v2)
fprintf(f1, 't(%d %d %d)n',v3)
fprintf(f1, 't(%d %d %d)n',v4)
fprintf(f1, 't(%d %d %d)n',v5)
fprintf(f1, '%2snn',');')

% Creating blocks
fprintf(f1, 'blocksn(n')

% Simple grading factor
grading = 0.2;
% Number of cells
b1 = 50;
b2 = 1;
b3 = 200;
fprintf(f1, 't %s (%d %d %d) %s (%d 0 0)n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading)
fprintf(f1, ');nn')

% Creating edges
fprintf(f1, '%sn','edges')
fprintf(f1, '(n')
fprintf(f1, 't arc %d %d (%d %d %d)n',1,2,0,r,0)
fprintf(f1, 't arc %d %d (%d %d %d)n',4,5,L_p,r,0)
fprintf(f1, ');nn')

% Creating boundary
fprintf(f1, '%sn','boundary')
fprintf(f1, '(n')
% Inlet
fprintf(f1, 't inlet nt {n')
fprintf(f1, 't   type  patch;nt')
fprintf(f1, 't   facesnt')
fprintf(f1, 't   (nt')
fprintf(f1, 'tt (0 1 2 0)nt')
fprintf(f1, 't   );nt }n')
% Outlet
fprintf(f1, 't outlet nt {n')
fprintf(f1, 't   type  patch;nt')
fprintf(f1, 't   facesnt')
fprintf(f1, 't   (nt')
fprintf(f1, 'tt (3 4 5 3)nt')
fprintf(f1, 't   );nt }n')
% Top
fprintf(f1, 't top nt {n')
fprintf(f1, 't   type  patch;nt')
fprintf(f1, 't   facesnt')
fprintf(f1, 't   (nt')
fprintf(f1, 'tt (1 2 5 4)nt')
fprintf(f1, 't   );nt }n')
% Front
fprintf(f1, 't front nt {n')
fprintf(f1, 't   type  symmetry;nt')
fprintf(f1, 't   facesnt')
fprintf(f1, 't   (nt')
fprintf(f1, 'tt (0 3 5 2)nt')
fprintf(f1, 't   );nt }n')
% back
fprintf(f1, 't back nt {n')
fprintf(f1, 't   type  symmetry;nt')
fprintf(f1, 't   facesnt')
fprintf(f1, 't   (nt')
fprintf(f1, 'tt (0 3 4 1)nt')
fprintf(f1, 't   );nt }n')

fprintf(f1, 'n);')

% Mergre patch
fprintf(f1, 'n%s','mergePatchPairs')
fprintf(f1, 'n(')
fprintf(f1, 'n);n')

fprintf(f1, 'n%s', header9)
fclose(f1);

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

covertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.961947e-03 -8.715574e-04)
	(0 9.961947e-03 8.715574e-04)
	(2.820000e+00 0 0)
	(2.820000e+00 9.961947e-03 -8.715574e-04)
	(2.820000e+00 9.961947e-03 8.715574e-04)
);

blocks
(
	 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 0 0)
);

edges
(
	 arc 1 2 (0 1.000000e-02 0)
	 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

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

);
mergePatchPairs
(
);

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

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
{
    axis
    {
    	type		empty;
    }
    inlet
    {
        type            zeroGradient;
    }

    outlet 
    {
        type            fixedValue;
        value           uniform 0.018;
    }
    
    top
    {
    	type		zeroGradient;
    }
    
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
}

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

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.0937312 0 0);

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

    outlet
    {
        type            zeroGradient;
    }
    
    top
    {
    	type		fixedValue;
    	value		uniform (0 0 0);
    }
    
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
}

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

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] 1e-06;


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

controlDict

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   200;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

output:

Profile and meshing

Pressure

Velocity

For L=0.05m For L=2.52m

Shear Stress Courant number And run time

For theta = 25

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

covertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.762960e-03 -2.164396e-03)
	(0 9.762960e-03 2.164396e-03)
	(2.820000e+00 0 0)
	(2.820000e+00 9.762960e-03 -2.164396e-03)
	(2.820000e+00 9.762960e-03 2.164396e-03)
);

blocks
(
	 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 0 0)
);

edges
(
	 arc 1 2 (0 1.000000e-02 0)
	 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

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

);
mergePatchPairs
(
);

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

Output:

Profile and mesh

Pressure

Velocity

for L=0.05m for L=2.52m

Shear stress Courant number and run time

For theta = 45

Profile and mesh

Pressure

Velocity

for L=0.05m for L=2.52m

2. Symmetry

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

covertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.961947e-03 -8.715574e-04)
	(0 9.961947e-03 8.715574e-04)
	(2.820000e+00 0 0)
	(2.820000e+00 9.961947e-03 -8.715574e-04)
	(2.820000e+00 9.961947e-03 8.715574e-04)
);

blocks
(
	 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 0 0)
);

edges
(
	 arc 1 2 (0 1.000000e-02 0)
	 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
	 inlet 
	 {
	   type  patch;
		   faces
		   (
			 (0 1 2 0)
		   );
	 }
	 outlet 
	 {
	   type  patch;
		   faces
		   (
			 (3 4 5 3)
		   );
	 }
	 top 
	 {
	   type  patch;
		   faces
		   (
			 (1 2 5 4)
		   );
	 }
	 front 
	 {
	   type  symmetry;
		   faces
		   (
			 (0 3 5 2)
		   );
	 }
	 back 
	 {
	   type  symmetry;
		   faces
		   (
			 (0 3 4 1)
		   );
	 }

);
mergePatchPairs
(
);

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

Profiles

theta = 10 theta = 25

theta = 45

Pressure

theta = 10

theta = 25

theta = 45

Vleocity:

theta = 10

theta = 25

theta = 45

Shear stress and Courant number

theat = 10

theta = 25

theta = 45

Results:

Solutions

Theta = 10

Theta = 25

Theta = 45

Symmetry

0.186385

0.186522

0.186809

Wedge

0.186384

0.186522

0.186809

Numerical solution

0.1875

0.1875

0.1875

Error (NS ~ Symmetry)

0.001115

0.000978

0.000691

Error (NS ~ Wedge)

0.001116

0.000978

0.000691




Symmetry

Theta = 10

Theta = 25

Theta = 45

Courant Number Mean

0.0689128

0.0689136

0.689152

Courant Number Max

0.341419

0.34465

0.353631

Execution Time (in sec)

390.32

380.71

375.88

Clock Time (in sec)

392

394

388




Wedge

Theta = 10

Theta = 25

Theta = 45

Courant Number Mean

0.0689128

0.0689136

0.689152

Courant Number Max

0.341313

0.344619

0.353614

Execution Time (in sec)

349.51

350.56

341.57

Clock Time (in sec)

450

458

453
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.