Simulation of flow through a pipe in OpenFOAM Part 1
Simulation of flow through a pipe in OpenFOAM
Part-I
-By Mandeep Shetty
This project is about the simulation laminar flow of water along a smooth pipe, using the transient incompressible solver, ‘icofoam’, in OpenFOAM. It is assumed that the water is flowing with a Reynolds number of 2100, and a temperature of 25o C. As the flow is considered to be taking place in a cylindrical pipe, I am simulating the flow as a 2D using an axisymmetric model, with the application of a wedge-shaped block. The wedge angle is considered to be 4o. The aim of the project is to compare the computed values with the analytical values calculated using Hagen Poiseuille’s equation and to show that the analytically calculated entry length is sufficient to produce a fully developed flow in the pipe.
Data:
Fluid: Water at 25oC
Reynolds number = 2100
Pipe diameter: 1Cm or 0.01m
Density: 997 Kg/m3
Dynamic Viscosity: 8.9e-4
Wedge angle: 4o
Number of cells in x-direction: 500
Number of cells in y-direction: 20
Mesh grading factor in x-direction: 5
Mesh grading factor in y-direction: 0.5
Calculations:
1.Entry length for laminar flow
Lent = 0.05*Re*D
Lent = 1.05m
The pipe length is assumed to be 0.5m more than the entry length
Therefore, L = 1.55m
2.Average inlet flow velocity
`V = (Re*μ)/(ρ*D)`
V = 0.1875 m/s
3.Pressure Drop
`∆P= (32*μ*L*V)/D^2`
`∆P=82.77 Pa`
Since icofoam is an incompressible flow solver, it uses the simplified incompressible flow equation, which uses kinematic viscosity and calculates kinematic pressure
`nu= μ/ρ= 8.92e-7`m2/s
Kinematic pressure drop = `∆P/ρ= 0.0830`m2/s2
Even though the flow is 2D, the velocity in y-direction is considered to be zero to calculate the time step using the Courants criteria, C=1
`∆t= V_x/(C*∆x)+V_y/(C*∆y)`
The maximum x-component of velocity, Vx, is 0.368m/s and minimum mesh dimension in x-direction, , is 0.0012m.
Therefore,
The current simulation is running with a timestep of 0.003s
A Matlab program, pipe_prob.m, is coded to automate the process of creating the mesh file, blockMeshDict.
pipe_prob.m
clc
clear all
close all
Re = 2100;
D = 0.01; %input(\"Enter the diameter of the pipe, in m: \\n\");
r = D/2;
fprintf(\"The fluid flowing in the pipe is water at a temperature of 25 deg and at Re = 2100\");
Mu = 8.90e-4;
rho = 997;
Nu = Mu/rho;
%Ref: https://www.engineersedge.com/physics/water__density_viscosity_specific_weight_13146.htm
v = (Re*Mu)/(rho*D); %average velocity calculated using the input Re
%Laminar flow entrance length
Lent = 0.05*Re*D ;
L = Lent + 0.5; %A little extra length to get a fully developed flow and avoid backflow, if any
dp = (32*Mu*L*v)/(D^2); %Pressure drop
dpk = dp/rho; %Kinematic pressure drop
%Wedge angle
a = 4 ; %In degree %input(\"Enter the wedge angle, in degrees: \\n\");
xcells = 500; %input(\"Enter the number of cells along x-direction: \\n\");
ycells = 20; %input(\"Enter the number of cells along y-direction: \\n\");
xgrad = 5; %input(\"Enter the (simple)grading factor along x-axis \\n\");
ygrad = 0.5; %input(\"Enter the (simple)grading factor along y-axis \\n\");
%The 6 vertices for the wedge
v1 = [0 0 0];
v2 = [L 0 0];
v3 = [L r*cosd(a/2) -r*sind(a/2)];
v4 = [0 r*cosd(a/2) -r*sind(a/2)];
v5 = [0 r*cosd(a/2) r*sind(a/2)];
v6 = [L r*cosd(a/2) r*sind(a/2)];
fid1 = fopen(\'blockMeshDict\',\'w\');
fprintf(fid1,[\' /*--------------------------------*- C++ -*----------------------------------*\\\\ \\n\' ...
\' ========= | \\n\' ...
\' \\\\\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\n\' ...
\' \\\\\\\\ / O peration | Website: https://openfoam.org \\n\' ...
\' \\\\\\\\ / A nd | Version: 6 \\n\' ...
\' \\\\\\\\/ M anipulation | \\n\' ...
\'\\\\*---------------------------------------------------------------------------*/ \\n\']);
fprintf(fid1,[\'FoamFile \\n\'...
\'{ \\n\' ...
\' version 2.0; \\n\' ...
\' format ascii; \\n\' ...
\' class dictionary; \\n\' ...
\' object blockMeshDict; \\n\' ...
\'} \\n\' ...
\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \\n\' ]);
fprintf(fid1,[\'\\nconvertToMeters 1; \\n\']);
fprintf(fid1,[\'\\nvertices\\n\'...
\'( \\n\'...
\' (%f %f %f) \\n\'...
\' (%f %f %f) \\n\'...
\' (%f %f %f) \\n\'...
\' (%f %f %f) \\n\'...
\' (%f %f %f) \\n\'...
\' (%f %f %f) \\n\'...
\'); \\n\'],v1,v2,v3,v4,v5,v6);
fprintf(fid1,[\'\\nblocks \\n\'...
\'( \\n\'...
\' hex (0 1 2 3 0 1 5 4) (%d %d 1) simpleGrading (%f %f 1) \\n\'...
\'); \\n\'],xcells,ycells,xgrad, ygrad);
fprintf(fid1,[\'\\n edges \\n\'...
\'( \\n\'...
\' arc 3 4 (0 %f 0) \\n\' ...
\' arc 2 5 (%f %f 0) \\n\' ...
\'); \\n\' ],r, L, r);
fprintf(fid1,[\'\\n boundary \\n\'...
\'( \\n\' ...
\' pipewall \\n\'...
\' { \\n\'...
\' type wall; \\n\'...
\' faces \\n\'...
\' ( \\n\'...
\' (4 5 2 3) \\n\'...
\' ); \\n\'...
\' } \\n\'...
\' backwedge \\n\'... %For axi-symmentric case, the 2 wedge planes should be specified as seperate patches of type \'wedge\'.Hence the frontwedge and backwedge.
\' { \\n\'...
\' type wedge; \\n\'...
\' faces \\n\'...
\' ( \\n\'...
\' (0 3 2 1) \\n\'...
\' ); \\n\'...
\' } \\n\'...
\' frontwedge \\n\'...
\' { \\n\'...
\' type wedge; \\n\'...
\' faces \\n\'...
\' ( \\n\'...
\' (0 1 5 4) \\n\'...
\' ); \\n\'...
\' } \\n\'...
\' inlet \\n\'...
\' { \\n\'...
\' type patch; \\n\'...
\' faces \\n\'...
\' ( \\n\'...
\' (0 4 3 0) \\n\'...
\' ); \\n\'...
\' } \\n\'...
\' outlet \\n\'...
\' { \\n\'...
\' type patch; \\n\'...
\' faces \\n\'...
\' ( \\n\'...
\' (1 2 5 1) \\n\'...
\' ); \\n\'...
\' } \\n\'...
\' axis \\n\'...
\' { \\n\'...
\' type empty; \\n\'...
\' faces \\n\'...
\' ( \\n\'...
\' (0 1 1 0) \\n\'...
\' ); \\n\'...
\' } \\n\'...
\' ); \\n\' ]);
fprintf(fid1,[\'\\nmergePatchPairs \\n\'...
\'( \\n\'...
\'); \\n\'...
\'\\n// ************************************************************************* // \\n\']);
Results:
Once the simulation is done, post-process the result in Parafoam and the required values of various variables at all the cell vertices can be exported in the CSV format by using the File -> Save data option. Save the ‘point data’ using the name ‘data/csv’. Use this CSV file, containing the different components of velocity, pressure and the coordinates of the points, as input to the Matlab code,shear_stress_calculator.m, to calculate the shear stress and plot the velocity and pressure along the length of the pipe and the variation of shear stress at the entrance length. Please not to run shear_stress_calculator.m after running pipe_prob.m as shear_stress_calculator utilizes certain variables that get stored in workspace on running pipe_prob script.
shear_stress_calculator.m
%Run this after runnin the pipe_prob program. This uses variables stored in
%the workspace from running pipe_prob program.
clc
close all
data = readmatrix(\'data.csv\');
S_L= 1; %Sampling lenght, (Choose a value along the lenght of pipe)
[~,inx] = (min(abs(data(:,5)-S_L)));
test_len = data(inx,5);
count = 1;
for i = 1:length(data(:,5))
if(data(i,5) == test_len && data(i,7)>0)
test_index(count) = i;
count = count+1;
end
end
test_vel = data(test_index,1);
test_rad = data(test_index,6);
for j = 2:length(test_vel)
dv = test_vel(j) - test_vel(j-1);
dy = test_rad(j) - test_rad(j-1);
tau(j-1) = -8.9e-4*(dv/dy);
end
figure(1)
plot(test_rad(1:end-1),tau(1:end)) %Wall shear stress
xlabel(\'Radial distance,m\');
ylabel(\'Wall shear stress,N/m^2\');
title(\'Shear stress vs Raidal distance\');
figure(2)
plot(data(1:xcells+1,5),data(1:xcells+1,4)*rho) %Pressure drop along the length
hold on
plot([data(1,5),data(xcells+1,5)],[dp,0],\'r\');
xlabel(\'Pipe Length,m\');
ylabel(\'Pressure,Pa\');
title(\'Pressure vs Pipe length\');
legend(\'CFD\',\'Analytical\');
hold off
figure(3)
plot(data(1:xcells+1,5),data(1:xcells+1,1),data(1:xcells+1,5),2*v*ones(xcells+1),\'r\') %Velocity variation along the length.Analytical velocity used is the max velocity,2*v, and not average velocity,v
xlabel(\'Pipe length,m\');
ylabel(\'Velocity,m/s\');
title(\'Velocity vs Pipe length\');
legend(\'CFD\',\'Analytical\');Conclusion:
The development of the velocity profile along the length of the pipe can be seen in Figures 2 to 8. The maximum velocity at the analytically calculated entry length, 1.05m, resulting from the simulation is 0.368 m/s which is 98.133% identical to the analytically calculated max velocity this point, 0.375 m/s, for the fully developed flow. Figure 9 shows that the velocity from the simulation asymptotically approaches the analytically calculated max velocity. There is a large pressure drop a the entrance of the pipe which the analytical calculation fails to capture. The shear stress value is maximum at the pipe wall and minimum at the pipe center, which is in accordance with the experimental findings.