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

Week 11 - Simulation of Flow through a pipe in OpenFoam

SSandeep Kolaparthi29 Jul 2022
CFDMATLAB

In this challenge, we have to simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.

You will also write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.

Expected Results

  1. Validate Hydro-dynamic length with the numerical result
  2. Validate the fully developed flow velocity profile with its analytical profile
  3. Validate maximum velocity and pressured drop for fully developed flow
  4. Post-process Shear stress and validate wall shear stress for fully developed flow

MATLAB CODE :

Matlab Program

clear all
close all
clc

% The Reynolds number
Re = 2100
% The Diameter of pipe
D = 0.004
% The Entry length of pipe
L =0.06*Re*D
% The total length of pipe
L_p = L+0.7
% Radius of pipe
r = D/2
% Angle of wedge (less than 5)
theta = 4
% Dynamic viscosity
mu = 0.89e-3
% Density 
rho = 997

% Hagen poiseuille flow equations
v_avg = (Re*mu)/(rho*D)
v_max = 2*(v_avg)
del_p = (8*mu*L_p*v_avg)/(r^2)
kinematic_p = del_p/rho

% calculating shear stress
R = linspace(0,(D/2),1000)

for i = 1:length(R)
    tau(i) = 2*mu*v_max*(abs(R(i)))/r^2
end

plot(R,tau)
xlabel(\'Radius\')
ylabel(\'Shear stress\')
grid on
legend (\'Shear stress\')


header1=\'/*----------------------------------*- C++ -*------------------------------------*\\\'
header2=\'  =========                 |                                                    \' 
header3=\'   \\\\      /  F ield                | OpenFOAM: The Open Source CFD Toolbox      \'
header4=\'    \\\\    /   O peration        | Website:  https://openfoam.org                 \'
header5=\'     \\\\  /    A nd                 | Version:  7                                 \'
header6=\'      \\\\/     M anipulation   |                                                  \'
header7=\'\\*----------------------------------------------------------------------------------*/\'

f1 = fopen(\'blockMeshDict\',\'w\')
fprintf(f1, \'%s\\n\', header1)
fprintf(f1, \'%s\\n\', header2)
fprintf(f1, \'%s\\n\', header3)
fprintf(f1, \'%s\\n\', header4)
fprintf(f1, \'%s\\n\', header5)
fprintf(f1, \'%s\\n\', header6)
fprintf(f1, \'%s\\n\', header7)

fprintf(f1, \'FoamFile\\n{\\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;\\n\',\'object\')
fprintf(f1, \'}\\n\')
header8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\'
fprintf(f1, \'%s\\n\\n\',header8)


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

fprintf(f1, \'vertices\\n\')
fprintf(f1, \'(\\n\')
fprintf(f1, \'\\t (%d  %d %d)\\n\',0,0,0)
fprintf(f1, \'\\t (%d  %d %d)\\n\',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, \'\\t (%d  %d %d)\\n\',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, \'\\t (%d  %d %d)\\n\',L_p,0,0)
fprintf(f1, \'\\t (%d  %d %d)\\n\',L_p,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, \'\\t (%d  %d %d)\\n\',L_p,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, \');\\n\\n\')


fprintf(f1, \'blocks\\n\')
fprintf(f1, \'(\\n\')
fprintf(f1, \'\\t hex (0 3 5 2 0 3 4 1) (%d %d %d)\',300,30,1)
fprintf(f1, \' SimpleGrading (1 0.1 1)\\n\\n\')
fprintf(f1, \');\\n\\n\')

fprintf(f1, \'edges\\n\')
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, \');\\n\\n\')

fprintf(f1, \'boundary\\n\')
fprintf(f1, \'(\\n\')
fprintf(f1, \'\\t inlet \\n\')
fprintf(f1, \'\\t {\\n\')
fprintf(f1, \'\\t\\t type patch;\\n\')
fprintf(f1, \'\\t\\t faces\\n\')
fprintf(f1, \'\\t\\t (\\n\')
fprintf(f1, \'\\t\\t\\t (0 1 2 0)\\n\')
fprintf(f1, \'\\t\\t );\\n\')
fprintf(f1, \'\\t }\\n\')


fprintf(f1, \'\\t outlet\\n\')
fprintf(f1, \'\\t {\\n\')
fprintf(f1, \'\\t\\t type patch;\\n\')
fprintf(f1, \'\\t\\t faces\\n\')
fprintf(f1, \'\\t\\t (\\n\')
fprintf(f1, \'\\t\\t\\t (3 5 4 3)\\n\')
fprintf(f1, \'\\t\\t );\\n\')
fprintf(f1, \'\\t }\\n\')


fprintf(f1, \'\\t top\\n\')
fprintf(f1, \'\\t {\\n\')
fprintf(f1, \'\\t\\t type wall;\\n\')
fprintf(f1, \'\\t\\t faces\\n\')
fprintf(f1, \'\\t\\t (\\n\')
fprintf(f1, \'\\t\\t\\t (1 4 5 2)\\n\')
fprintf(f1, \'\\t\\t );\\n\')
fprintf(f1, \'\\t }\\n\')

fprintf(f1, \'\\t front\\n\')
fprintf(f1, \'\\t {\\n\')
fprintf(f1, \'\\t\\t type wedge;\\n\')
fprintf(f1, \'\\t\\t faces\\n\')
fprintf(f1, \'\\t\\t (\\n\')
fprintf(f1, \'\\t\\t\\t (0 3 4 1)\\n\')
fprintf(f1, \'\\t\\t );\\n\')
fprintf(f1, \'\\t }\\n\')

fprintf(f1, \'\\t back\\n\')
fprintf(f1, \'\\t {\\n\')
fprintf(f1, \'\\t\\t type wedge;\\n\')
fprintf(f1, \'\\t\\t faces\\n\')
fprintf(f1, \'\\t\\t (\\n\')
fprintf(f1, \'\\t\\t\\t (0 2 5 3)\\n\')
fprintf(f1, \'\\t\\t );\\n\')
fprintf(f1, \'\\t }\\n\')

fprintf(f1, \'\\t axis\\n\')
fprintf(f1, \'\\t {\\n\')
fprintf(f1, \'\\t\\t type empty;\\n\')
fprintf(f1, \'\\t\\t faces\\n\')
fprintf(f1, \'\\t\\t (\\n\')
fprintf(f1, \'\\t\\t\\t (0 3 3 0)\\n\')
fprintf(f1, \'\\t\\t );\\n\')
fprintf(f1, \'\\t }\\n\')
fprintf(f1, \'); \\n\\n\')

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

header9=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\'
fprintf(f1, \'%s\\n\',header9)

fclose(f1)

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  0 0)
	 (0  1.998782e-03 6.979899e-05)
	 (0  1.998782e-03 -6.979899e-05)
	 (1.204000e+00  0 0)
	 (1.204000e+00  1.998782e-03 6.979899e-05)
	 (1.204000e+00  1.998782e-03 -6.979899e-05)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.000000e-03 0)
	 arc 4 5 (1.204000e+00 2.000000e-03 0)
);

boundary
(
	 inlet 
	 {
		 type patch;
		 faces
		 (
			 (0 1 2 0)
		 );
	 }
	 outlet
	 {
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	 }
	 top
	 {
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	 }
	 front
	 {
		 type wedge;
		 faces
		 (
			 (0 3 4 1)
		 );
	 }
	 back
	 {
		 type wedge;
		 faces
		 (
			 (0 2 5 3)
		 );
	 }
	 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         1;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Velocity Boundary Condition :

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

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

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }

    top
    {
       type             fixedValue;
       value uniform (0 0 0);
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

}

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

Pressure Boundary Condition :

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

    axis
    {
        type            empty;
    }

    top
    {
        type            zeroGradient;
    }

    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:  7
     \\\\/     M anipulation  |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    \"constant\";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 8.9256e-07;


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

Results

Generated Geometry and Mesh

Post Processing

  1. After opening paraview click on generate to see the geometry. View the geometry from different sides across X,Y and Z axes.
  2. Press on the play button to post process to end time.
  3. Click on filters and select plot over line and plot velocity magnitude and pressure with respect to radial distance.
  4. Export the velocity data in .csv format.
  5. For calculating the shear stress import the .csv file into MATLAB/Octave.
  6. Select the suitable column to plot shear stress.
  7. Write a MATLAB/Octave code to calculate shear stress.
  8. The shear stress should be maximum at walls.
  9. Plot the result for shear stress.
  • Calculating Shear Stress
clear all
close all
clc

D = 0.02

R = D/2

r = linspace(-0.01,0.01,1000)

Re = 2100

mu = 8.9e-4

rho = 997.130

V = (Re*mu)/(rho*D)

u_max = 2*V

for i= 1:length(r)
  Tau(i) = (2*mu*u_max*abs(r(i)))/(R^2);
end
figure(1)
plot(r,Tau)
xlabel(\'Radial Distance in m\')
ylabel(\'Shear Stress(Tau) in N/m^2\')
title(sprintf(\'change in shear stress for fully developed flow at x=2.82m with radial distance\'))
grid minor



u=csvread(\'velocity.csv\',\'A2:A1002\');
y = linspace(-0.005,0.005,1000);
dy = y(2)-y(1);

for k=336:length(y);
  du(k)=u(k+1)-u(k);
  tau(k)=mu*(abs(du(k))/dy);
end




figure(2)
plot(r,tau)
xlabel(\'Radial Distance in m\')
ylabel(\'Shear Stress(tau) in N/m^2\')
title(sprintf(\'change in shear stress for fully developed flow at x=2.82m with radial distance\'))
grid minor

Shear stress :

1. When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.

2. For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.

3. Fluids Obeying Newton\'s law where the values of “mu” are constant are known as Newtonian fluids. If ‘mu’ is constant, the shear stress is linearly dependent on the velocity gradient.

By considering the velocity of this relative motion, the Shear rate can be calculated shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layer of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and one is stationary.

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.