Week 11 - Simulation of Flow through a pipe in OpenFoam

FLOW PHYSICS:

Let us consider fluid entering a circular pipe with uniform velocity. the fluid particles or more precisely the fluid layer in contact with the wall surface will come to a complete sto resulting in zero velocity. this is because of no slip condition at the wall surface this layer also causes the fluid particles in thr adjacent layer to show down gradually as result of friction. to make up for this velocity reduction the velocity of fluid at the mid section of pipe has to increase to keep mass flow rate through out the pipe constant. because of this there is formation of velocity gradient along the pipe.

The region of flow in which the effect of viscous shearing forces are caaused by fluid viscosity are felt is called as 'velocity boundary layer' or normally a boundary layer . the hypotheticaly boundary condition surface divides the flow in a pipe into region boundary layer region in which the viscous effects and velocity changes are significant and the irrotational flow region in which the friction effects are negligible and the velocity remains constant in radial direction . the thickness of boundary condition layer increases in the flow direction untill the boundary layer reaches the pipe center as shown in fig . the region from the pipe inlet to the point at which boundary layer merges at the centerline is called as hydrodynamics entrance region and kength of this region is hydrodynamic enterance length . in this region the velocity profile is not fully devrlop but it is developing as a function of pipe length . the region beyond the entrance region where velocity profile is develop is called as hydrodynamically fully develop region the velocity profile in the hydrodynamic fully develop region is parabolic in nature.

IN PUT PARAMETRS AND PIPE GEOMETRY CHARACTERISTICS:

type of fluid = water

temperature of fluid = 20 deg C

dynamic viscosity at 20deg C= 1.002e-3 n.s/m^2

kinematic viscosity at 20deg C = 1.004e-6 m^2/s

shape of pipe = cylindrical

pipe diameter = 4mm or0.004m

pipe length =0.7m

wedge angle= 4deg

scaled model =1:1

as our fluid flow is purely laminar and incompressible, 

1 reynold number 

assume reynold number of Re=2100 for laminar conditions.

2 average velocity.

the average velocity of fluid for flow through pipe under laminar flow condition can be given as

3 hydrodynamic enterance length 

4 total length of pipe =0.7m

here 0.28m is added to pipe length to have a fully develop region as hydrodynamic length is 0.42m which is just an entery region and velocity profile is not fully develop here. so as to have both the hydrodynamic enterance and fully develop region we added 0.28m to the pipe length

tthese above stated are input parameters which will be used to create matlab script and acts as an analytical datum for our CFD results.

MATLAB CODE 

clear all 
close all
clc

%% variable

%reynold number
R=2100;
% wedge angle (in degree)
theta = 3;
% diametr of pipe (m)
D=0.004;
%radius of a pipe(m)
rad= D/2;
% hydrodynamic length (m)
L_h=0.05*R*D;
%total length of pipe 
L= 0.28+L_h;

%% calculating vertices
v_0 = [0 0 0];
vr_1 = [0 rad*cosd(theta/2) -rad*sind(theta/2)];
v_1= [vr_1(1) vr_1(2) vr_1(3)];
vr_2 = [0 rad*cosd(theta/2) rad*sind(theta/2)];
v_2= [vr_2(1) vr_2(2) vr_2(3)];
vr_3=[L 0 0];
v_3=[vr_3(1) vr_3(2) vr_3(3)];
vr_4=[L rad*cosd(theta/2) -rad*sind(theta/2)];
v_4=[vr_4(1) vr_4(2) vr_4(3)];
vr_5=[L rad*cosd(theta/2) rad*sind(theta/2)];
v_5=[vr_5(1) vr_5(2) vr_5(3)];

%% file writing
line1='/*------------------------------*-c++ -*-------------------------*';
line2=' ==========                  |';
line3='           /  F ield         | OpenFOAM:  The Open Source CFD Toolbox';
line4='          /   O peration     | Website: https://openfoam.org';
line5='         /    A nd           | Version:  8';
line6='        /     M anipulation  |';
line7='*---------------------------------------------------------------------------*/';

% file open
block_mesh = fopen('blockMeshDict.txt','w');
% to write header
fprintf(block_mesh,'%sn',line1);
fprintf(block_mesh,'%sn',line2);
fprintf(block_mesh,'%sn',line3);
fprintf(block_mesh,'%sn',line4);
fprintf(block_mesh,'%sn',line5);
fprintf(block_mesh,'%sn',line6);
fprintf(block_mesh,'%sn',line7);

% to write dictionary
fprintf(block_mesh,'%sn','FoamFile');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%sn','  version 2.0;');
fprintf(block_mesh,'%sn','  format  ascii;');
fprintf(block_mesh,'%sn','  class   dictionary;');
fprintf(block_mesh,'%sn','  object  blockMeshDict;');
fprintf(block_mesh,'%sn','}');
fprintf(block_mesh,'%sn','//**********************************//');
fprintf(block_mesh,'n');

%note as all calculated variable are in keeping  scale to be as one (1)
sacle=1;  % scale input

fprintf(block_mesh,'%s%.4f%sn','convertToMeters',sacle,';');
fprintf(block_mesh,'n');
% vertices
fprintf(block_mesh,'%sn','vertices');
fprintf(block_mesh,'%s','(');
fprintf(block_mesh,'n');
fprintf(block_mesh, '(%d %d %d)n', v_0);
fprintf(block_mesh, '(%d %d %d)n', v_1);
fprintf(block_mesh, '(%d %d %d)n', v_2);
fprintf(block_mesh, '(%d %d %d)n', v_3);
fprintf(block_mesh, '(%d %d %d)n', v_4);
fprintf(block_mesh, '(%d %d %d)n', v_5);
fprintf(block_mesh, '%snn',');');

% block parameters
fprintf(block_mesh, 'n');
mesh ='hex';
vertice = [0 1 2 3 4 5];
vertice_seq=[0 1 2 0 3 4 5 3];
nx = 30;
ny = 1;
nz= 200;
xgf = 0.2
ygf = 1;
zgf = 1;
arc_1 = [1 2];
arc_2 = [4 5];
% origin and extension along x
x0 = 0;
y0 = 0;
z0= 0;
x1 = L;
 
% write block
fprintf(block_mesh,'%sn','blocks');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%s (%d %d %d %d %d %d %d %d) (%d %d %d) %s (%d %d %d)',mesh,vertice_seq,nx,ny,nz,'simpleGrading',xgf,ygf,zgf);
fprintf(block_mesh,'n');
fprintf(block_mesh,'%snn',');');
% edges
fprintf(block_mesh,'%sn','edges');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%s %d %d (%d %d %d)n','arc',arc_1(1),arc_1(2),x0,rad,z0);
fprintf(block_mesh,'%s %d %d (%d %d %d)n','arc',arc_2(1),arc_2(2),x0,rad,z0);
fprintf(block_mesh,'%snn',');');
% boundary
fprintf(block_mesh,'%sn','boundary');
boundary_type ={'patch','empty','wall','symmteryPlane','wedge'}; %types of wall boundary condition
% inlet
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%sn','inlet');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{1},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(3),vertice(2),vertice(1));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% outlet
fprintf(block_mesh,'%sn','outlet');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{1},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(4),vertice(5),vertice(6),vertice(4));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% top wall
fprintf(block_mesh,'%sn','topWall');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{3},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(2),vertice(3),vertice(6),vertice(5));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% bottom wall
fprintf(block_mesh,'%sn','axis');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{3},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(4),vertice(4),vertice(1));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
%front wedge
fprintf(block_mesh,'%sn','frontWedge');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{5},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(4),vertice(6),vertice(3));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% back wedge
fprintf(block_mesh,'%sn','backWedge');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{5},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(2),vertice(5),vertice(4));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
fprintf(block_mesh,'%sn',');');

fprintf(block_mesh,'n');
fprintf(block_mesh,'%sn','mergePatchPairs');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','//****************************************************************//');
fclose(block_mesh);

BLOCKMESHDICT FILE

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

convertToMeters 1.0000;

vertices
(
    (0 0 0)
    (0 1.998782e-03 -6.979899e-05)
    (0 1.998782e-03 6.979899e-05)
    (7.000000e-01 0 0)
    (7.000000e-01 1.998782e-03 -6.979899e-05)
    (7.000000e-01 1.998782e-03 6.979899e-05)
);

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

edges
(
    arc 1 2 (0 2.000000e-03 0)
    arc 4 5 (7.000000e-01 2.000000e-03 0)
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 2 1 0)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (3 4 5 3)
        );
    }
    topwall
    {
        type wall;
        faces
        (
            (1 2 5 4)
        );
    }
    axis
    {
        type emty;
        faces
        (
            (0 3 3 0)
        );
    }
    frontWedge
    {
        type wedge;
        faces
        (
            (0 3 5 2)
        );
    }     
    backkWedge
    {
        type wedge;
        faces
        (
            (0 1 4 3)
        );
    }                  
);

mergePatchPairs
(
);

inlet boundary condition for pressure

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

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    inlet
    {
        type            zeroGradient;
    }

    frontWedge
    {
        type            wedge;
    }
    backWedge
    {
        type            wedge;
    }
    topWall
    {
        type            zeroGradient;
    }
    axis
    {
        type            empty;
    }           
}

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

inlet boundary condition for velocity

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

    outlet
    {
        type            zeroGradient;
    }

    frontWedge
    {
        type            wedge;
    }
    backWedge
    {
        type            wedge;
    }
    topWall
    {
        type            noslip;
    }
    axis
    { 
        type            empty;
    }            
}

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

CONTROLDICT

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.008;

writeControl    runTime;

writeInterval   5;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

VIEW IN PARAFOAM velocity distribution

inition boundary layer formation

velocity distribution along radial distance x=0.01m

velocity distribution along radial distance x=0.35

velocity distribution along radial distance x=0.65

 pressure across the pipe length

pressure plot across pipe length

MATLAB CODE FOR SHEAR STRESS FULLY DEVELOP CONDITION

clear all
close all
clc
% required variable
% renold number
R= 2100;
mu=1.002e-3; %dynamic viscosity of water at t=20deg c
% kinematic viscosity
k_v = 1.004e-6;
% diam of pipe
D=0.004;
 
% read the text file extracted from simukltion
f=fopen('arc_length_tau.txt','r');
g=fopen('tau.txt','r');
h=fopen('velocity.txt','r');
% skipping header line
a= fgetl(f);
b=fgetl(g);
c=fgetl(h);

% initialize variables
U= ones(1,1001);
r_d=zeros(1,1001);
U(1)=str2double(fgetl(h));
tau_an=zeros(1,1001);

% calculating shear stress newton second law of viscosity
for i=2:1002
    r_d(i)=str2double(fgetl(f));
    U(i)=str2double(fgetl(h));
    tau_an(i)=mu*(U(i-1)-U(i))/(r_d(i)-r_d(i-1));
    U_min = min(U);
U_max=max(U);
end

U_avg=(U_max-U_min)/2;
u_avg_analytical = R*k_v/(d);
U_max_analytical=2*U_avg_analytical;

error_avg=U_avg_analytical-U_avg;
error_max=U_max_analytical-U_max;

% PLOTTING SHEAR STRESS
figure(1)
plot(r_d,tau_an,'r')
hold on
plot(-r_d,tau_an,'r')
grid on
xlabel('Radius (m)')
ylabel('shear stress (N/m^2) (pa)')
title('shear stress variation along radial distance')

% plotting velocity
figure(2)
plot(r_d,U,'b')
hold on
plot(-r_d,U,'b')
grid on
hold off
xlabel('Radius (m)')
ylabel('velocity (m/s)')
title('Velocity  variation along radial distance')

SHEAR STRESS PLOT

above graph shear stress plot shows v shaped shear stress distribution as per our reduced model(wedge shape).

velocity plot:

the above plot shows perfect parabolic velocity profile in hydrodynamically fully develop region.  

* kinematic pressure gradient:

matlab code

clear all
close all
clc

R=2100; %reynolds number
theta=4; %wedge angle
D=0.004; %diam of pipe
rad=D/2; %radius of pipe
L_h=0.05*R*D; %hydrodynamic of length
L=0.28+L_h; %totla length of pipe
l=L-L_h; %length of fully develop
n=1250;%totalnum of time step
charac_1=0.18; %characetristic number
rho=1000;%density of number
k_v=1.004e-6;%kinematic viscosity
mu=1.008e-3;%dynamic viscosity

v_avg=R*k_v/(D);

press_drop_simple=(8*mu*v_avg*charac_1)/(rad^2);

f=fopen('arc_length_pressure.txt','r');
g=fopen('press.txt','r');

skip_arc=fgetl(f);
skip_press=fgetl(g);

for i=2:1002
    a=fgetl(g);
    kinep(i)=str2double(a);
    c=fgetl(f);
    arc(i)=str2num(c);

end
%calculating kinematic pressure gradient
error_0=0;
for m=2:length(kinep)-1
    kinepgrad_num(m)=abs(kinep(m+1)-kinep(m));
    for x=2:length(arc)-1
        press_grad(x)=(8*mu*v_avg)/(rad^2)*(arc(x+1)-arc(x));
        kinepgrad_analytical(x)=press_grad(x)/(rho);
        error(m)=abs(kinepgrad_analytical(x)-kinepgrad_num(m));
    end
end

arc_length=arc(2:end);


mean_err=mean(error);
mea_pr_grad=mean(press_grad);
mean_num=mean(kinepgrad_num);
mean_an=mean(kinepgrad_analytical);

figure(1)
plot(arc_length,kinepgrand_analytical,'r')
hold on
plot(arc_length,kinepgrad_num,'b')
title('kinematic pressure gradient variation along x-distance')
xlabel('X-distsnce (fully develop region)')
ylabel('kinematic pressure gradient')
legend('kinematic pressure gradient (analytical)','kinematic pressure gradient (numerical)')

kinematic pressure plot:

*disscusion:

in this plot pressure gradient in hydrodynamically fully develop region is calculated in the x-direction in this region the velocity gradient in uniform along the radial distance which means the pressure at discrete radial distance is constant . so the drop of pressure in the x direction will also be constant as it is hydrodynamically fully develop flow as we are dealling with the uncompressible flows the pressure is normalized with density of fluid which results in kinematic pressure so in this case our numerically kinematic pressure is closely matching to analytical kinematic pressure.

conclusion:

we have successfully executed a set up for numerical case for flow through pipe .there were a lot of parameter in case to setup ,flow physics,many others.the reduced model of the cyclinder pipe gave congruent results as that of the analytical solution.