Week 11 - Simulation of Flow through a pipe in OpenFoam
%calculating shear stress
clear all
close all
clc
D = 0.02;
R = D/2;
Re = 2100; % Reynolds number for laminar flow
mu = 8.9e-4; % dynamic viscosjity of water
rho = 997.139; % denisty of water
V = (Re*mu)/(rho*D); %average velocity
U_max = 2*V; %maximum velocity
r = linspace(-0.01,0.01,100);
for i =1:length(r)
tau (i) = (2*mu*U_max*abs(r(i)))/(R^2)
end
figure(1)
plot(r,tau)
xlabel('Radial distance in m')
ylabel('Shear stress(tau) in N/m^2')
% BlockMesh file Created in MATLAB
clear all
close all
clc
L = 1;
D = 0.02; % diameter of pipe
R = D/2; % radius of pipe
theta= 4; % wedge angle
nx = 500;
ny = 50;
nz = 10;
u = 0.0937;
Re = 2100;
L_c = 0.06*Re*D; % hydrodynamic entry length
L_p = L_c+0.7; % characterstics length
mu = 0.000889;
rho = 997;
%nu = mu/rho; % kinematic viscosity
nu = 8.9238e-9;
%finding vertices
%creating blockMeshDict file:
v0 = [0 0 0];
v1 = [0 R*cosd(0.5*theta) -R*sind(0.5*theta)];
v2 = [0 R*cosd(0.5*theta) R*sind(0.5*theta)];
v3 = [L_p 0 0];
v4 = [L_p R*cosd(0.5*theta) -R*sind(0.5*theta)];
v5 = [L_p R*cosd(0.5*theta) R*sind(0.5*theta)];
f1 = fopen('blockMeshDict.txt','w'); % creating a text file
fprintf(f1,'%s', '/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'\n');
fprintf(f1,'%s', '========== ');
fprintf(f1,'\n');
fprintf(f1,'%s', '\\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / O peration | Version: 2.2.1 ');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / A nd | Web: www.OpenFOAM.com ');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\/ M anipulation | ');
fprintf(f1,'\n');
fprintf(f1,'%s','*---------------------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'%s', 'FoamFile');
fprintf(f1,'\n');
fprintf(f1,'%s','{');
fprintf(f1,'\n');
fprintf(f1,'%s', 'version 2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s', 'format ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s', 'class dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s', 'object blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'%s','}');
fprintf(f1,'\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s', 'convertToMeters 1');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'vertices');
fprintf(f1,'\n');
fprintf(f1,'(');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(7),'(');
fprintf(f1,'%f %f %f' ,v0(1),v0(2),v0(3));
fprintf(f1,')');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(7),'(');
fprintf(f1,'%f %f %f' ,v1(1),v1(2),v1(3));
fprintf(f1,')');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(7),'(');
fprintf(f1,'%f %f %f' ,v2(1),v2(2),v2(3));
fprintf(f1,')');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(7),'(');
fprintf(f1,'%f %f %f' ,v3(1),v3(2),v3(3));
fprintf(f1,')');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(7),'(');
fprintf(f1,'%f %f %f' ,v4(1),v4(2),v4(3));
fprintf(f1,')');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(7),'(');
fprintf(f1,'%f %f %f' ,v5(1),v5(2),v5(3));
fprintf(f1,')');
fprintf(f1,'\n');
fprintf(f1,');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'blocks');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(3));
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
fprintf(f1,'%s', blanks(3));
fprintf(f1,'%s',' (50 1 200)');
fprintf(f1,'%s', blanks(3));
fprintf(f1,'%s','simpleGrading (1 0.1 1)');
fprintf(f1,'\n');
fprintf(f1,');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'edges');
fprintf(f1,'\n');
fprintf(f1,'%s', '(');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(4));
fprintf(f1,'arc 1 2 (');
fprintf(f1,'%f %f %f' ,0,R,0);
fprintf(f1,'%s', ')');
fprintf(f1,'\n');
fprintf(f1,'%s', blanks(4));
fprintf(f1,'arc 4 5 (');
fprintf(f1,'%f %f %f' ,L,R,0);
fprintf(f1,'%s', ')');
fprintf(f1,'\n');
fprintf(f1,'%s', ');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'boundary');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'axis');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'inlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'pipe_surface');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(1 4 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'outlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(3 4 5 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'front');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'back');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 1 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','mergePatchPairs');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'%s','//**************************************************************//');
fclose(f1);
# Pressure(P)
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
pipe_surface
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0169;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Objectives of the Project :
1. Calculation of flow variables.
2. Generations of blockMeshDict using the wedge boundary condition. create the blockMesgDict file in MATLAB such that the mesh is automatically created for less than 5degree wedge angle .
3. Using the pipeflow in OpenFoam using ICOFOAM solver.
4. Post process the Velocity profile and shear stress and compare the simulation results with analytical results obtained using the Hagen-Poiseuille equation.
Entrance length :
The length of the entrance is the distance that fluid travels after entering a pipe before the full development of the flow.
Hagen poiseuille's Law : In non linear fluid dynamics, Hagen-Poiseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a long cylinder pipe of constant cross section.
Average velocity Vavg = Re.mu/(rho*D)
Maximum velocity Vmax = 2*Vavg.
Summary :
Consider a flow entering a pipe. Let us think of the entering flow being uniform . As soon as the flow 'hits' the pipe many changes take place. The most important of these is that viscosity imposes itself on the flow and the 'No slip' condition at the wall of the pipe comes into effect. Consequently the velocity components are each zero on the wall u = v = 0. The flow adjacent to the wall decelearates continuously.The length of the pipe between the start and the point where the fully developed flow begins is called the entrance length.
Methodology -
The icoFoam solver available in OpenFoam has been used. Due to unavailability of a steady state solver for laminar incompressible flow in OpenFoam, icoFoam solver has been used for a time of 50s and a time step size of 0.01.
Assumptions -
1. The fluid is Newtonian.
2. No slip at the pipe wall.
3. There is no velocity in the radial direction.
4. Reynolds number for the laminar flow is 2100.
Procedure -
Firstly, i need to calculate the average velocity of the given fluid. for that i need to give proper inputs of reynolds number, Diameter,Viscosity and Density values.
1. Be carefull with parenthesis like commas,brackets and slash symbols which may leads us to the error.
2. Angles for the wedg has been given in the program.
3. Using all these values we calculate average velocity, maximum velocity, pressure difference and kinematic pressure difference. The reason to calculate the kinematic pressure is in the OpenFoam pressure file.
4. The angle theta taken less than 3 degree is because only one line has to pass through the center of the axis of symmetry.
5. By using the file parsing technique, we need to write the matlab code that can be useful to perform the simulations using different angles.
6. After getting the BlackMeshDict file,copy it to the case setup file which is located in the run folder.
7. Like this, change other parameters like velocity,Pressure,and Transport properties file which is obtained from the analytical solution means values obtained from the Matlab Code.
8. Now set the start and end time in the control-dict file, and Mention the delta T value also.
9. Now place the command directory file in the cavity, run the file using command "block-Mesh".
10. The Above command gives us the mesh file.
11. Now check the mesh to see whether any error in mesh, Non-alignment and some other using 'checkMesh' command.
12. Onces the mesh check is over and success. Now use 'icoFoam' solver to run the case setup.
13. After running successfully, Use 'paraFoam' command to post-process the velocity and shear stress of the fluid.
Shear Stress
1. when the fluid is in motion , shear stress are developed due to the particles in the fluid moving relative to one another.
2. For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.
3. Fluids obeying Newton law where the value of 'mu' are constant are known as Newtonian fluids. if 'mu'(dynamic viscosity) is constant, the shear stress is linearly dependent on the velocity gradient .
here,
tau = mu(du/dy)
Conclusion :
- The velocity profile is plotted at varying distances from the inlet of the pipe.
- The velocity at a distance x = 0.01 from the inlet seems to be slow and stable around the center of the pipe which shows that the velocity is gradually developing. It suddenly becomes zero at the wall because the non-slip condition.
- The velocity profile at a distance of x = 0.5 , x = 1.26, x = 2.52 is parabolic in nature because the velocity profile moves to the hydrodynamically fully developed region.
- The shear stress plot is plotted as given below.
Plot of shear Stress
Plotting shear stress is not straight forward as plotting the velocity profile. Because in paraView tau is calculated directly,so it becomes relevent that we first calculate tau and then use plotcoverline command to plot the shear stress.
we have Newton law for viscous flow , tau = mu*(du/dy) thi scan be used to calculate the shesr stress. Steps of calculations are given below:
1. firstly, we need to calculate magnitude of vector U, this can be achieved using the mag option present in the calculator. If we provide mag(U) in the calculator we will get the magnitude of velocity.
2. Once magnitude of velocity is computed next step is calculated its gradient , this can be achieved by inserting gradient, this is computational option present in the paraview. Once the option is inserted from filters we can provide the calculated velocity gradient as input and hence the gradient of velocity can be calculated.
3. Once the gradient is calculated,dynamic viscosity of water at 25degree C is multiplied by 0.000889 so hence tau is calculated.
# Velocity (U)
/*--------------------------------*- 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.0937 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0937 0 0);
}
pipe_surface
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //# blockMeshDict
/*--------------------------------*- C++ -*----------------------------------*
==========
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Version: 2.2.1
\\ / A nd | Web: www.OpenFOAM.com
\\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009962 -0.000872)
(0.000000 0.009962 0.000872)
(2.520000 0.000000 0.000000)
(2.520000 0.009962 -0.000872)
(2.520000 0.009962 0.000872)
);
blocks
(
hex (0 3 4 1 0 3 5 2 ) (200 50 1) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
pipe_surface
{
type wall;
faces
(
(1 4 5 2)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//**************************************************************//# 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 27;
deltaT 0.001;
writeControl timeStep;
writeInterval 50;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //