Project

SIMULATING LAMINAR FLOW THROUGH A PIPE USING ICOFOAM SOLVER AND SYMMETRY BOUNDARY CONDITION IN OPENFOAM (PART 2)

RRajamahanti Harshith18 Nov 2020

Aim :

To write a matlab code for symmetry boundary conditions of flow through a pipe
To compare results of wedge boundary conditions with symmetric boundary condition

Theory :

The inputs to the problem are same as the one in wedge boundary condition problem. BlockMeshDict file needs to be updated according to required value of angles i.e 10, 25 and 45

Preprocessing steps :

The initial steps for this challenge remain as the wedge challenge
change the value of theta as per the requirement in each case in matlab and obtain the blockMeshDict file
change the wedge type from wedge to symetry in blockMeshDict file in system folder
similarly make changes in pressure and velocity file in 0 folder and change wedge type
Other inputs remain the same

Matlab code for generating blockMeshDict file and computing Analytical data :

close all
clear all
clc
%The Renoylds number
Re=2100;
%Diameter of the pipe in meters
D=0.01;
%The entry length of the pipe 
Le=0.06*Re*D;
%Total length of the pipe 
L=Le+0.3;
%Radius of the pipe 
R=D/2;
%Angle of the wedge in degrees
Theta=45;
%Dynamic viscosity
mu=0.89e-3;
%Density in kg/m^3
rho=997;
% Haigen poisuelle flow equations
%Avg velocity
Vavg=(Re*mu)/(rho*D);
%Maximum velocity 
Vmax=2*Vavg;
%Piezo metric pressure drop 
Delp=32*mu*Vavg*L/D^2;
%Kinematic pressure 
Kinematic_p=Delp/rho;
%creating set of Diameters for plotting the shear stress
Rvar=linspace(-R,R,1000);
%Calculating shear stress
for i=1:length(Rvar)
    Tau(i)=2*mu*(Rvar(i))*Vmax/R^2;
    v(i)=Vmax*(R^2-Rvar(i)^2)/R^2;
end
%Maximum shear stress value
Tau_max=max(Tau)
%plotting
figure(1)
plot(Tau,Rvar)
hold on 
axis([-0.5 0.5 -R R])
xL = xlim;
yL = ylim;
line(xL,[0 0],'linestyle','--','color','r','linewidth',1);  %x-axis
ylabel('Radius')
xlabel('Shear stress')
title('Analytical variation of shear stress along the radial distance')
grid on
legend('Shear stress','Axis of the pipe')
figure(2)
plot(v,Rvar)
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1);  %x-axis
grid on
xlabel('Velocity')
ylabel('Radius')
legend('Velocity','Axis of the pipe')
title('Analytical variation of velocity along the radial distance')
%creating the blockmesh file
header1='/*---------------------------*- C++ -*--------------------------------*/'
header2=' ==========                 |                                           '
header3='  \\        / F ield          |  OpenFOAM: The open source CFD Toolbox  '
header4='   \\      /  O perration    |  website:  https://openfoam.org          '
header5='    \\    /   A nd                      |  version :  8                 '
header6='     \\  /    M anipulation       |                                     '
header7='\*------------------------------------------------------------------*/'

f1=fopen('blockMeshDict.c','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,0,0)
fprintf(f1,'\t (%d %d %d)\n',L,R*cosd(Theta/2),R*sind(Theta/2))
fprintf(f1,'\t (%d %d %d)\n',L,R*cosd(Theta/2),-R*sind(Theta/2))
fprintf(f1,');\n\n')

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 2 1 0 3 5 4 3) (%d %d %d)',20,1,100)
fprintf(f1, 'simpleGrading (0.2 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,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\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\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\n')

fprintf(f1,'\t front \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type symmetry;\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\n')

fprintf(f1,'\t back \n')
fprintf(f1,'\t {\n')
fprintf(f1,'\t\t type symmetry;\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\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\n')
fprintf(f1,');\n')

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

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

fclose(f1)

blockMeshDict file :

/*---------------------------*- C++ -*--------------------------------*/
 ==========                 |                                           
  \\        / F ield          |  OpenFOAM: The open source CFD Toolbox  
   \\      /  O perration    |  website:  https://openfoam.org          
    \\    /   A nd                      |  version :  8                 
     \\  /    M anipulation       |                                     
\*------------------------------------------------------------------*/
FoamFile
{
     version 	 2.0;
      format 	 ascii;
       class 	 dictionary;
      object 	 blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
	 (0 0 0)
	 (0 4.619398e-03 1.913417e-03)
	 (0 4.619398e-03 -1.913417e-03)
	 (1.560000e+00 0 0)
	 (1.560000e+00 4.619398e-03 1.913417e-03)
	 (1.560000e+00 4.619398e-03 -1.913417e-03)
);

blocks
(
	 hex (0 2 1 0 3 5 4 3) (20 1 100)simpleGrading (0.2 1 1)

);

edges
(
	 arc 1 2 (0 5.000000e-03 0)
	 arc 4 5 (1.560000e+00 5.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 symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	 }

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

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

);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

controlDict 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;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.05;

writeControl    timeStep;

writeInterval   10;

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:  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.1875 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }
   
    axis
    {
        type            empty;
    }
   
    top
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    front
    {
        type            symmetry;
    }
   
    back
    {
        type            symmetry;
    }
}

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

Pressure boundary 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       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    axis
    {
        type            empty;
    }
   
    top
    {
        type            zeroGradient;
    }
   
    front
    {
        type            symmetry;
    }
   
    back
    {
        type            symmetry;
    }
   
}

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

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.9268e-7;

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

Results :

10 degree :

Mesh :

Inlet velocity :

25 degree :

Mesh:

Inlet velocity :

45 degree :

Mesh :

Inlet velocity :

Velocity variation for 45 deg wedge angle :

at length=0.01

At length = 0.05

At length =0.1

At length =0.5

At length = 1

At length = 1.5

Calculation of shear stress at 10, 25 and 45 degree wedge angle by importing the .csv files from openfoam

For 45 deg :

clear all
close all
clc
%Dynamic viscosity
mu=0.89e-3;
%The Renoylds number
Re=2100;
%Diameter of the pipe in meters
D=0.01;
%The entry length of the pipe 
Le=0.06*Re*D;
%Total length of the pipe 
L=Le+0.3;
%Radius of the pipe 
R=D/2;

data1=csvread('45deatx=1.5.csv',1,0);
u1=data1(:,2);
y1=data1(:,1);
for i=1:length(y1)
    y1_n(i)=-y1(length(y1)+1-i);
    y1_neg=y1_n';
     u1_n(i)=u1(length(u1)+1-i);
    u1_neg=u1_n';
end
y1_complete=[y1_neg;y1(2:end)];
u1_complete=[u1_neg;u1(2:end)];
plot(u1_complete,y1_complete,'r','linewidth',2)
u1_maximum=max(u1_complete)
u1_avg=u1_maximum/2;

data2=csvread('45degatx=1.csv',1,0);
u2=data2(:,2);
y2=data2(:,1);
for i=1:length(y2)
    y2_n(i)=-y2(length(y2)+1-i);
    y2_neg=y2_n';
     u2_n(i)=u2(length(u2)+1-i);
    u2_neg=u2_n';
end
y2_complete=[y2_neg;y2(2:end)];
u2_complete=[u2_neg;u2(2:end)];
u2_maximum=max(u2_complete)
hold on
plot(u2_complete,y2_complete,'g','linewidth',2)

data3=csvread('45degatx=0.5.csv',1,0);
u3=data3(:,2);
y3=data3(:,1);
for i=1:length(y3)
    y3_n(i)=-y3(length(y3)+1-i);
    y3_neg=y3_n';
     u3_n(i)=u3(length(u3)+1-i);
    u3_neg=u3_n';
end
y3_complete=[y3_neg;y3(2:end)];
u3_complete=[u3_neg;u3(2:end)];
u3_maximum=max(u3_complete)
hold on
plot(u3_complete,y3_complete,'b','linewidth',2)

data4=csvread('45degatx=0.1.csv',1,0);
u4=data4(:,2);
y4=data4(:,1);
for i=1:length(y4)
    y4_n(i)=-y4(length(y4)+1-i);
    y4_neg=y4_n';
     u4_n(i)=u4(length(u4)+1-i);
    u4_neg=u4_n';
end
y4_complete=[y4_neg;y4(2:end)];
u4_complete=[u4_neg;u4(2:end)];
u4_maximum=max(u4_complete)
hold on
plot(u4_complete,y4_complete,'y','linewidth',2)

data5=csvread('45degatx=0.05.csv',1,0);
u5=data5(:,2);
y5=data5(:,1);
for i=1:length(y5)
    y5_n(i)=-y5(length(y5)+1-i);
    y5_neg=y5_n';
     u5_n(i)=u5(length(u5)+1-i);
    u5_neg=u5_n';
end
y5_complete=[y5_neg;y5(2:end)];
u5_complete=[u5_neg;u5(2:end)];
u5_maximum=max(u5_complete)
hold on
plot(u5_complete,y5_complete,'k','linewidth',2)

data6=csvread('45degatx=0.01.csv',1,0);
u6=data6(:,2);
y6=data6(:,1);
for i=1:length(y6)
    y6_n(i)=-y6(length(y6)+1-i);
    y6_neg=y6_n';
     u6_n(i)=u6(length(u6)+1-i);
    u6_neg=u6_n';
end
y6_complete=[y6_neg;y6(2:end)];
u6_complete=[u6_neg;u6(2:end)];
u6_maximum=max(u6_complete)
hold on
plot(u6_complete,y6_complete,'c','linewidth',2)
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1);  %x-axis
xlabel('velocity')
ylabel('Radial distance')
legend('At length=1.5','At length=1','At length=0.5','At length=0.1','At length=0.05','At length=0.01','axis')
title('velocity variation along the radial distance for 45 degree wedge angle')

%calculation of shear stress

for i=1:length(y1_complete)
    Tau(i)=2*mu*(y1_complete(i))*u1_maximum/R^2;
end
Tau_max=max(Tau)
%calculation of pressure difference 
Delp=32*mu*u1_avg*L/D^2;
figure(3)
plot(Tau,y1_complete,'k','linewidth',2)
xlabel('shear stress')
ylabel('Radial distance')
axis([-0.15 0.15 -R R])
xL = xlim;
line(xL, [0 0],'linestyle','--','color','r','linewidth',1);  %x-axis
legend('shear stress','axis')
title('shear stress variation along the radius for 45 degree wedge angle')

similarly it can be done for 10 degree and 25 degree

Comparison :

Variables

Analytical

(HP)

2 degree wedge

3 degree wedge

4 degree wedge

10 degree

symmetrical

25 degree

symmetrical

45 degree

symmetrical

Maximum velocity in m/sec

0.3749

0.3714

0.3712

0.3715

0.371

Wall shear stress in N/m^2

0.1335

0.1322

0.1321

0.1316

0.1291

0.1223

Conclusion :

While the angle increases the error in maximum velocity is decreased and courants number is slightly increasing
There is no significant difference between analytical and simulated results
Symmetric plane boundary condition plays a crucial role in terms of accuracy in result
In case of symmetry boundary condition , it mirrors all the values as it reduces computational efforts
Symetric boundary conditions mirrors faces/surfaces which enable us to go for a higher wedge angle where as wedge boundary condition is constrained to < 5 degree as it is a axi-symetry boundary condition.

SIMULATING LAMINAR FLOW THROUGH A PIPE USING ICOFOAM SOLVER AND SYMMETRY BOUNDARY CONDITION IN OPENFOAM (PART 2)

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.