Week 10 - Simulation of Flow through a pipe in OpenFoam

Aim

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 have to validate results with the Hagen- Poiseuille's equations for the fully developed flow.

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.

The Hagen-Poiseuille equation, also known as Poiseuille's law, describes the flow of a viscous fluid through a cylindrical pipe or tube. It relates the flow rate of the fluid to various parameters such as the pressure difference across the pipe, the pipe length, the pipe radius, and the viscosity of the fluid. The equation is named after Gotthilf Hagen and Jean Léonard Marie Poiseuille, who independently derived it.

The Hagen-Poiseuille equation is given as follows:

Q = (π * ΔP * r^4) / (8 * η * L),

where: Q is the volumetric flow rate of the fluid, ΔP is the pressure difference across the pipe, r is the radius of the pipe, η is the dynamic viscosity of the fluid, and L is the length of the pipe.

The equation shows that the flow rate (Q) is directly proportional to the pressure difference (ΔP) and the fourth power of the pipe radius (r^4), while it is inversely proportional to the viscosity (η) and the length of the pipe (L).

The Hagen-Poiseuille equation assumes laminar flow, where the fluid flows smoothly in parallel layers with no turbulence. It also assumes that the fluid is incompressible and the pipe has a constant circular cross-section.

This equation finds applications in various fields such as fluid dynamics, engineering, and medicine, where the flow of fluids through pipes or blood vessels needs to be analyzed or controlled.

Procedure

for creating blockmesh we need to write code in matlab .Here for solving the problem numerically we will use the icoFoam solver with the cavity tutorial make the geometry along with the mesh and make the blockmeshDict file with the help of matlab as shown below:

matlab code

 matlab code

close all
clc
%% Input %%
Re= 2100; % Reynold's Number
D= 0.02; % Diameter of pipe
R=D/2; % Radius
theta= input('Enter Wedge angle ='); % Wedge angle
mu=0.89e-3; %Dynamic viscocity @ 25 degree Celsius
rho=997; %Density of water @ 25 degree Celsius
nu= mu/rho; %Kinematic viscocity @ 25 degree Celsius
% Analytical Velocity, pressure
L_e= 0.06*Re*D; % Entrance length
L_t= L_e + 0.3; % Total length
v_avg= (Re*mu)/(rho*D); % Average velocity
v_max= 2*v_avg; % Maximum velocity
Delta_P= (32*mu*v_avg*L_t)/(D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding vertices
v0= [0 0 0];
v1= [0 R*cosd(theta/2) -R*sind(theta/2)];
v2= [0 R*cosd(theta/2) R*sind(theta/2)];
v3= [L_t 0 0];
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];
%blockMeshDict file.
f1= fopen('blockMeshDict.txt','w'); % create blockmesh file
fprintf(f1,'%s \n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'%s \n',' ========= |');
fprintf(f1,'%s \n',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'%s \n',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'%s \n',' \\ / A nd | Version: 10');
fprintf(f1,'%s \n',' \\/ M anipulation |');
fprintf(f1,'%s \n','\*---------------------------------------------------------------------------*/');
fprintf(f1,'%s \n','FoamFile');
fprintf(f1,'%s \n','{');
fprintf(f1,'%s \n',' version 2.0;');
fprintf(f1,'%s \n',' format ascii;');
fprintf(f1,'%s \n',' class dictionary;');
fprintf(f1,'%s \n',' object blockMeshDict;');
fprintf(f1,'%s \n','}');
fprintf(f1,'%s \n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%s \n','convertToMeters 1;');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
% Finding vertices
fprintf(f1,'vertices \n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s \n','(');
v_0=fprintf(f1,'(%f %f %f) \n',v0(1),v0(2),v0(3));
fprintf(f1,b4,'%s','(');
v_1=fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','(');
v_2=fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b4,'%s','(');
v_3=fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b4,'%s','(');
v_4=fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b4,'%s','(');
v_5=fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');
% blocks creation
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% mesh generation
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:');
m3= input('Enter no. of mesh element along Z-axis:');
fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%2.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% edges creation
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');
% Apply boundary conditions
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%s',' axis');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% fixedwalls
fprintf(f1,'%s',' fixedWalls');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% front
fprintf(f1,'%s',' front');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(1 0 3 4)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% back
fprintf(f1,'%s',' back');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% outlet
fprintf(f1,'%s',' outlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
fprintf(f1,'\n);\n');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
 

This is the blockmesh file created by the matlab code above

boundary condition

matlab code 
 
clear all
close all
clc
τ
=
μ
(
d
u
d
y
)
% Initial Parameters
Re=2100;
D= 0.02;
R=D/2;
mu=0.00089; %Dynamic viscocity
rho=997;
nu=mu/rho;
%analytical calculation
L_e=0.06*Re*D;
L_t=L_e+0.3;
v_avg=(Re*nu)/D;
v_max=2*v_avg;
Delta_P=(8*mu*v_avg*L_t)/R^2;
kin_P=Delta_P/rho;
j= linspace(0,R,1000);
for i=1:length(j)
tau(i)=2*mu*v_max*(abs(j(i)))/R^2;
end
figure(1)
plot(j,tau)
grid on
hold on
title("Shear Stress")
xlabel("Radius")
ylabel('Shear Stress')