CFD simulation of flow through a pipe using icoFoam solver in OpenFoam

Abstract:

The given problem is described to demonstrate a flow simulation over a pipe to observe the boundary layer formation. The velocity profile is captured at three different zones which are at start of entry length, at the end of entry length and after the entry length. Similarly, pressure drop is also observed at the fully developed region. The computational numerical data are validated against the analytical results to support the simulation algorithm, which efficiently predicted the physics of fluid flow along the entire geometrical domain.

Theory:

When a real fluid flows past a stationary solid boundary, a layer of fluid which comes in contact with the boundary surfaces adheres to it (on account of viscosity) and condition of no slip occurs. Thus, the layer of fluid which does not slips away undergoes retardation, this retarded layer further causes retardation in the vicinity of boundary layers in which the velocity of the fluid increases rapidly from zero at the boundary surface to maximum velocity at the main stream. The layer adjacent to the boundary is known as boundary layer. The formation is given below:

Objective:

1. Create a MATLAB code to generate computational mesh automatically for any wedge angle and grading schemes.
2. Predict the velocity profile in the fully developed region.
3. For a base line mesh the velocity profile must match with Hagen-Poiseuilles equation.
4. Post process velocity and shear stress in PARAVIEW.

Problem specification:

geometry

Considering a sector of wedge of the pipe

Initial conditions

Length of the pipe= 1.5m

Density=1000 kg/m³

Diameter=0.01 m

Dynamic viscosity= 8.9*10^-4

Mesh specification

No of cells in x direction = 20

No of cells in y direction = 1

No of cells in z direction = 300

Mesh grading in x direction=0.2

Mesh grading in y direction=1

Mesh grading in z direction=1

blocks

(0 2 1 0 3 5 4 3)

boundary wall conditions

Algorithm:

1. Using file parsing technique, the blockMeshDict file is generated from MATLAB.
2. Flow calculation regarding the average velocity, pressure calculation using Hagen-Poiseulles equation, flow rate is computed.
3. The generated blockMeshDict file is transferred to system run folder.
4. The simulation start and end time is defined in controlDict file.
5. The pressure and velocity boundary conditions are given with respect to wall faces.
6. Run the icofoam solver and converge the solution.
7. Post process the results in PARAVIEW to generate graphical simulation data.
8. Plot shear stress data in the fully developed region. The Hagen-Poiseulles equation for pressure and velocity data are numerically validated with the analytical data.

Workflow diagram:

MATLAB code for blockMeshDict file

clc
clear all
close all

% %------------simulating pipe flow-------------------------------------%
% % specify the condition for the flow as well as the dimension of domain
Re=2100;
rho=1000;
mu=8.9e-4;
d=0.01;
r=d/2;
l=5;
th=3;
area=pi/4*d^2;
l_entry=0.05*d*Re;
v_avg=Re*mu/(rho*d);
v_max=2*v_avg;
Q=v_avg*area;
pressure=32*mu*v_avg*l/d^2;
kinematic_pressure=pressure/rho;
% specifying the number of mesh cells in x,y,z direction 
nx=20;
ny=1;
nz=300;
% specifying the grading 
gx=0.2;
gy=1;
gz=1;

%specifying geometry of the wedges
w_0=[0 0 0];
w_1=[0 r*cosd(th/2) r*sind(th/2)];
w_2=[0 r*cosd(th/2) -r*sind(th/2)];
w_3=[5 0 0];
w_4=[5 r*cosd(th/2) r*sind(th/2)];
w_5=[5 r*cosd(th/2) -r*sind(th/2)];
arc_1=[0.0 r 0.0];
arc_2=[5 r 0.0];

%Using file parsing technique to create blockmeshdict file

line1='/*--------------------------------*- C++ -*------------------------*';
line2=' ==========             	                                         ';
line3='   	 /  F ield     	   | OpenFOAM: The Open Source CFD Toolbox   ';
line4='      /   O peration 	   | Website:  https://openfoam.org          ';
line5='     /	A nd       	   | Version:  7                             ';
line6='    / 	M anipulation  | CodeDeveloper: SOUDIP HAZRA             ';
line7='*-----------------------------------------------------------------*/';
line8='FoamFile ';
line9=' { ';

f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%sn',line1);
fprintf(f1,'%sn',line2);
fprintf(f1,'%sn',line3);
fprintf(f1,'%sn',line4);
fprintf(f1,'%sn',line5);
fprintf(f1,'%sn',line6);
fprintf(f1,'%sn',line7);
fprintf(f1,'n');
fprintf(f1,'%sn',line8);
fprintf(f1,'%snn',line9);

fprintf(f1,'%16s %10s n','version','2.0;          ');
fprintf(f1,'%15s %12s n','format' ,'ascii;        ');
fprintf(f1,'%14s %11s n','class'  ,'dictionary;   ');
fprintf(f1,'%15s %11s n','object' ,'blockMeshDict;');
fprintf(f1,'n } nn');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //n');
fprintf(f1,' convertToMeters %d;n n',1);

% defining geometry%
fprintf(f1,'tvertices n');
fprintf(f1,'tt( n');
fprintf(f1,'ttt(%f %f %f)n',w_0);
fprintf(f1,'ttt(%f %f %f)n',w_1);
fprintf(f1,'ttt(%f %f %f)n',w_2);
fprintf(f1,'ttt(%f %f %f)n',w_3);
fprintf(f1,'ttt(%f %f %f)n',w_4);
fprintf(f1,'ttt(%f %f %f)n',w_5);
fprintf(f1,'tt);n n');
fprintf(f1,'ttblocks n');
fprintf(f1,'tt( n');
fprintf(f1,'ttt hex (0 2 1 0 3 4 5 3) (%d %d %d) simplegrading (%d %d %d)ntt); n n',nx,ny,nz,gx,gy,gz);

% defining sides of geometry%
% edge of pipe%
fprintf(f1,'ttedges n');
fprintf(f1,'tt( n');
fprintf(f1,'tttarc 1 2 (%f %f %f) n',arc_1);
fprintf(f1,'tttarc 4 5 (%f %f %f) n',arc_2);
fprintf(f1,'tt);n n');

% defining boundary conditions%
fprintf(f1,'ttboundaryntt(n');
% inlet section of pipe
fprintf(f1,'%16sn','inletWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type patch');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(0 1 2 0)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
% outlet section of pipe
fprintf(f1,'%16sn','outletWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type patch');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(3 5 4 3)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');

%front side of pipe
fprintf(f1,'%16sn','frontWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type wedge');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(0 3 4 1)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');

%back side of pipe
fprintf(f1,'%16sn','backWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type wedge');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(2 5 3 0)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');

%top side of wall
fprintf(f1,'%16sn','topWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type wall');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(1 4 5 2)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');

%bottom of pipe
fprintf(f1,'%16sn','bottomWall');
fprintf(f1,'%10sn','{');
fprintf(f1,'%25sn','type empty');
fprintf(f1,'%19sn','faces');
fprintf(f1,'%15sn','(');
fprintf(f1,'%28sn','(0 3 3 0)');
fprintf(f1,'%16sn',')');
fprintf(f1,'%10sn','}');
fprintf(f1,'n');
fprintf(f1,');n n');

%merge patch pair
fprintf(f1,'%20sn','mergePatchPairs');
fprintf(f1,'%6sn','(');
fprintf(f1,'%5sn',');');
fprintf(f1,'%sn','//* * * * * * * * * * * * * * * * * * * * * * * * * *  * * * * * * * * * * * * *//');
fclose(f1)

ControlDict file code

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

deltaT      	1e-2;

writeControl	timeStep;

writeInterval   100;

purgeWrite  	0;

writeFormat 	ascii;

writePrecision  6;

writeCompression off;

timeFormat  	general;

timePrecision   6;

runTimeModifiable true;


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

Pressure boundary condition code

/*--------------------------------*- 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
{
	inletWall
	{
    	type        	zeroGradient;
	}
    
	outletWall
	{
    	type        	fixedValue;
    	value   	  uniform 0.2661;   	 
	}

	topWall
	{
    	type        	zeroGradient;
	}
    
	bottomWall
	{
   	 type   	  empty;
	}
    
	frontWall
	{
   	 type   	  wedge;
	}

	backWall
	{
    	type        	wedge;
	}
}

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

Velocity boundary condition code

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

boundaryField
{
	inletWall
	{
    	type        	fixedValue;
    	value       	uniform (0.1869 0 0);
	}

	outletWall
	{
    	type        	zeroGradient;
   	 
	}

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

	bottomWall
	{
    	type        	empty;
	}
    
	frontWall
	{
    	type        	wedge;
	}

	backWall
	{
    	type        	wedge;
	}
}

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

MATLAB code for validation of numerical solution with analytical solution

%  numerical vs analytical solution at fully developed region
clc
clear all
T=readtable('velocity.csv');
radial=T{:,1}';% radial distance
numerical_velocity=T{:,2}';% numerical velocity data
R=0.005;% radius of pipe
v_max=0.3738;% anlytical max velocity
M=readtable('pressure.csv');
pressure=M{:,2}';
max_pressure=max(pressure);
min_pressure=min(pressure);
numerical_pressure_difference=max_pressure-min_pressure;%numerical hagen poiseulli pressure
for i=1:length(radial)
    analytical_velocity=v_max.*(1-(radial./R).^2);%analytical velocity data
    error=abs(analytical_velocity-numerical_velocity);
end
figure(1)
plot(radial,analytical_velocity,'-*','color','b','linewidth',1.5)
hold on
plot(radial,numerical_velocity,'-s','color','g','linewidth',1.5)
legend('analytical velocity','CFD numerical velocity')
xlabel('Radial distance');
ylabel('velocity');
title('Analytical velocity vs Numerical velocity at fully developed region');
figure(2)
bar(error,'facecolor',[0.12 0.4 0.8]);
title('Error plot')
xlabel('velocity')
ylabel('error')

Outputs

1. Velocity data

Radial distance	velocity
0	0.37126
4.9983e-05	0.370947
9.9966e-05	0.370635
0.000149949	0.370323
0.000199932	0.370011
0.000249915	0.369699
0.000299898	0.369387
0.000349881	0.369075
0.000399864	0.368763
0.000449847	0.368451
0.00049983	0.368104
0.000549813	0.367099
0.000599796	0.366094
0.000649779	0.365089
0.000699762	0.364084
0.000749745	0.363079
0.000799728	0.362074
0.000849711	0.361069
0.000899694	0.360064
0.000949677	0.359059
0.00099966	0.357426
0.00104964	0.355729
0.00109963	0.354032
0.00114961	0.352336
0.00119959	0.350639
0.00124958	0.348942
0.00129956	0.347246
0.00134954	0.345549
0.00139952	0.343541
0.00144951	0.341231
0.00149949	0.338921
0.00154947	0.336611
0.00159946	0.334301
0.00164944	0.33199
0.00169942	0.32968
0.00174941	0.32737
0.00179939	0.324618
0.00184937	0.321748
0.00189935	0.318878
0.00194934	0.316008
0.00199932	0.313139
0.0020493	0.310269
0.00209929	0.307399
0.00214927	0.30417
0.00219925	0.300787
0.00224924	0.297405
0.00229922	0.294022
0.0023492	0.290639
0.00239918	0.287256
0.00244917	0.283786
0.00249915	0.279933
0.00254913	0.276079
0.00259912	0.272225
0.0026491	0.268372
0.00269908	0.264518
0.00274907	0.260579
0.00279905	0.256293
0.00284903	0.252007
0.00289901	0.247721
0.002949	0.243436
0.00299898	0.23915
0.00304896	0.234586
0.00309895	0.229904
0.00314893	0.225221
0.00319891	0.220539
0.0032489	0.215856
0.00329888	0.210938
0.00334886	0.20589
0.00339884	0.200843
0.00344883	0.195796
0.00349881	0.190747
0.00354879	0.185365
0.00359878	0.179984
0.00364876	0.174602
0.00369874	0.16922
0.00374873	0.16361
0.00379871	0.157921
0.00384869	0.152232
0.00389867	0.146542
0.00394866	0.140621
0.00399864	0.134649
0.00404862	0.128677
0.00409861	0.122648
0.00414859	0.116416
0.00419857	0.110185
0.00424856	0.103953
0.00429854	0.0975045
0.00434852	0.0910346
0.0043985	0.0845647
0.00444849	0.0779038
0.00449847	0.071215
0.00454845	0.0645101
0.00459844	0.05762
0.00464842	0.0507299
0.0046984	0.0437449
0.00474839	0.0366699
0.00479837	0.0295685
0.00484835	0.0223237
0.00489833	0.0150789
0.00494832	0.00938317
0.0049983	0

2.Pressure data

length of fully developed region	pressure
0	0.425712
0.0100001	0.42518
0.0200002	0.424647
0.0300003	0.424115
0.0400005	0.423582
0.0500006	0.423049
0.0600007	0.422516
0.0700008	0.421984
0.0800009	0.421452
0.090001	0.420919
0.100001	0.420386
0.110001	0.419854
0.120001	0.419322
0.130002	0.418789
0.140002	0.418257
0.150002	0.417724
0.160002	0.417192
0.170002	0.41666
0.180002	0.416127
0.190002	0.415595
0.200003	0.415062
0.210003	0.41453
0.220003	0.413998
0.230003	0.413465
0.240003	0.412933
0.250003	0.412401
0.260003	0.411868
0.270003	0.411336
0.280003	0.410804
0.290004	0.410272
0.300004	0.409739
0.310004	0.409207
0.320004	0.408675
0.330004	0.408143
0.340004	0.407611
0.350004	0.407079
0.360004	0.406546
0.370005	0.406014
0.380005	0.405482
0.390005	0.40495
0.400005	0.404418
0.410005	0.403885
0.420005	0.403353
0.430005	0.402821
0.440006	0.402289
0.450006	0.401756
0.460006	0.401224
0.470006	0.400692
0.480006	0.40016
0.490006	0.399628
0.500006	0.399096
0.510006	0.398564
0.520006	0.398032
0.530007	0.3975
0.540007	0.396968
0.550007	0.396435
0.560007	0.395903
0.570007	0.395371
0.580007	0.394839
0.590007	0.394307
0.600007	0.393775
0.610008	0.393243
0.620008	0.392711
0.630008	0.392179
0.640008	0.391647
0.650008	0.391115
0.660008	0.390582
0.670008	0.39005
0.680009	0.389518
0.690009	0.388986
0.700009	0.388454
0.710009	0.387922
0.720009	0.38739
0.730009	0.386858
0.740009	0.386326
0.750009	0.385794
0.760009	0.385262
0.77001	0.38473
0.78001	0.384198
0.79001	0.383666
0.80001	0.383134
0.81001	0.382602
0.82001	0.38207
0.83001	0.381538
0.84001	0.381006
0.850011	0.380473
0.860011	0.379942
0.870011	0.37941
0.880011	0.378878
0.890011	0.378346
0.900011	0.377813
0.910011	0.377282
0.920012	0.37675
0.930012	0.376218
0.940012	0.375686
0.950012	0.375154
0.960012	0.374622
0.970012	0.37409
0.980012	0.373558
0.990012	0.373026
1.00001	0.372494

Physics behind simulation:

1. The wedge angle is kept fixed at a certain angle to capture physics of flow from the entry region to the exit region. This allowed to visualize the boundary flow formation.
2. The flow dynamics is captured at different zone to predict the numerical maximum velocity.
3. The pressure gradient flow and shear stress at the fully developed region is numerically predicted to validate against the analytical solutions.

Result and analysis:

creating mesh pattern

velocity at x=0.1m (start of entry length)

The maximum velocity at the start of entry length is almost 0.265m/s.

velocity at x=1m (at end of entry length)

The maximum velocity is almost 0.37m/s which exactly matches the analytical max velocity.

velocity at x=2m (after the entry length)

The maximum velocity is almost 0.37m/s which exactly matches the analytical max velocity.

velocity at x=3m (after the entry length)

The maximum velocity is almost 0.37m/s which exactly matches the analytical max velocity.

pressure distribution

The pressure gradient is measured at two different regions of fully developed phase flow

point1(2,0,0)

point2(3,0.0049,0)

shear stress

The shear stress is also measured at two different regions of fully developed phase flow

point1(2,0,0)

point2(3,0.0049,0)

comparison between the analytical and numerical parameters in fully developed region

From the graph it is seen that the numerical data matches exactly with the analytical data at the fully developed region.

The error plot gets minimised as it moves forward from entry length region to the fully developed region.

Conclusions:

1. The boundary layer formation is captured exactly for a small wedge angle, having theta 3 degree.
2. The velocity data matches approximately equal to the numerical simulation.
3. From Hagen-Poiseulles equation for pressure at the fully developed region, it is seen that the numerical pressure drop is close to the analytical solution.
4. The error plot shows the minimum deviation from the exact results, which essentially captures the physics of the flow behaviour.