Week 12 - Symmetry vs Wedge vs HP equation
Objective:
1. Perform the flow through a pipe using a wedge boundary condition.
2. Perform the flow through a pipe using a symmetry boundary condition.
3. Calculate the velocity for the flow using Hagen Poisuelle's equation.
4. Write a Matlab code to create a blockMeshDict file for the following wedge angle (10,25,45).
5. Compare the results obtained.
Introduction:
The flow through pipes is usually involved in heating, cooling, and fluid distribution network applications. The fluid is forced to flow through the pipe with the help of a pump which gives an initial velocity in this case is 0.17m/s. Mostly circular pipes are used for the transport of fluids because circular pipes can withstand larges pressure differences between inside and outside without undergoing significant distortion.
The fluid velocity inside the pipe changes from zero at the wall because of the no-slip condition to the maximum at the center. It is convenient to work with an average velocity which remains constant for the incompressible flow as the cross-sectional area of the pipe is constant. The average velocity changes at heating and cooling applications because of changes with temperature and density. But here the fluid properties are evaluated at some average temperature and treated as constants.
Laminar flow is characterized by smooth streamline and highly ordered motion of the fluid particles.
Turbulent flow is characterized by velocity fluctuations and disordered motion of the fluid particles.
Reynolds Number: It is the ratio of inertia force to viscous force.
When Reynolds number is low which indicates that there is a large viscous force present which suppresses the fluctuations and maintains the ordered motion of the fluid molecules in other words the flow is laminar.
When Reynolds number is high large inertia forces are present which causes disorderedness in the motion of the molecule in other words the flow is turbulent.
`R_e`< = 2300 Laminar flow
`R_e` >= 4000 Turblent Flow.
2300`<=` `R_e``<=` 4000 Transistion flow.
The flow region where the viscous shearing forces due to fluid viscosity is felt is called the velocity boundary layer also called the boundary layer. Due to this hypothetical layer divides the flow region into a boundary layer where the viscous effects are felt and changes in the velocity are significant and an irrotational flow region where the frictional effects are negligible and velocity is constant.
The thickness of the boundary layer increases in the direction of the fluid flow until the reaches the center of the pipe and velocity is fully developed at a distance from the inlet.
The distance from the inlet of the pipe to a point where the velocity is fully developed is called hydrodynamic entrance length.
In the hydrodynamic entrance region, shear wall stress is highest at the inlet of the pipe and where the thickness of the boundary layer is minimum and shear stress decreases along the fluid flow direction to fully developed flow.
Entry length
The hydrodynamic entry length is taken to be the distance from the entrance region to the point where the shear stress is about 2 percent of the fully developed flow.
L/D = 0.005Re.
Velocity profile can be calculated as
`u = 2*V_(avg)*(1-(r)^2/(R)^2)`
where u is velocity profile in m/s.
`V_(avg)` is the average velocity in m/s.
r is the radial distance in m.
R is the radius of the pipe in m.
The maximum velocity is at the center of the pipe which can be obtained by substituting r=0 in the above formula
`V_(max) = 2* V_(avg)`
Hagen Poiseuille Equation:
Hagen Poiseuille equation evaluates the pressure drop in an incompressible and Newtonian fluid in a laminar flow that has a constant cross-sectional area and also the length of the pipe should be larger than its diameter.
Assumptions in this equation are that the fluid is an incompressible and Newtonian fluid whose length is larger than its diameter.
This equation calculates pressure drop due to viscosity in the fluid.
The equation to calculate the pressure drop is given by:
`triangleP = (128*mu*L*Q)/(pi*D^4)`
where Q is the volumetric flow rate in `m^3/s`
`Q=V_(avg)*A`
A is the area of the pipe in m.
D is the diameter of the pipe in m.
L is the length of the pipe in m.
`mu` is the dynamic viscosity in `kg.m^-1.s^-1`
Pressure Drop
Pressure drop is calculated as it is related to the power requirement of the pump to maintain the flow inside the pipe.
Pressur drop can also be calculated using the below formula.
`(dp)/dx = (P_2-P_1)/L`
`triangleP = (8* mu*L*V_(avg))/R^2`
`mu` is a dynamic viscosity in `kg.m^-1.s-1`.
L is the length of the pressure drop region.
`V_(avg)` is the average velocity in m/s.
R is the radius of the pipe in m.
From the above formula, the pressure loss is proportional to viscosity which when its value is zero then there is no pressure drop. Therefore the pressure drop from P1 to P2 is entirely due to viscous for constant diameter and fixed length of the pipe.
Matlab Code to automate the computational mesh for any wedge angle.
% Code to automate the computational mesh for any wedge angle
close all
clear all
clc
%Defining the wedge angle from which the vertices are calculated using
%trignometry fucntions.
theta = input('angle:');%<--User input for the anlge
theta = theta/2;
x = cosd(theta)*0.005;
y = sind(theta)*0.005;
%length of the pipe.
l = 2;
%type of boundary condition to be applied to the wedge face.
str = "symmetry";
%Creating a blockMeshDict file
f = fopen('blockMeshDict.txt','w');
%Below are the codes for blockMeshDict file.
fprintf (f ,'%s \n' ,'/*--------------------------------*- C++ -*----------------------------------*\');
fprintf (f,'%s \n',' ========= |');
fprintf (f,'%s\n',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf (f,'%s\n',' \\ / O peration | Website: https://openfoam.org');
fprintf (f,'%s\n',' \\ / A nd | Version: 8');
fprintf (f,'%s\n',' \\/ M anipulation |');
fprintf (f,'%s\n','\*---------------------------------------------------------------------------*/');
fprintf (f,'%s\n','FoamFile');
fprintf (f,'%s\n','{');
fprintf (f,'%s\n','');
fprintf (f,'%s\n',' version 2.0;');
fprintf (f,'%s\n',' format ascii;');
fprintf (f,'%s\n',' class dictionary;');
fprintf (f,'%s\n',' object blockMeshDict;');
fprintf (f,'%s\n','}');
fprintf (f,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','convertToMeters 1;');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','vertices');
fprintf (f,'%s\n','(');
fprintf (f,'%s\n',' (0 0 0)');
fprintf (f,'%s\n',' (0 0 0)');
fprintf (f,'%s %f %f %f %s \n',' (',y,x,0,')'); %<-- giving the calculated values as vertices for defining the wedge
fprintf (f,'%s %f %f %f %s \n',' (',-y,x,0,')');
fprintf (f,'%s %f %f %f %s \n',' (',0,0,l,')');
fprintf (f,'%s %f %f %f %s \n',' (',0,0,l,')');
fprintf (f,'%s %f %f %f %s \n',' (',y,x,l,')');
fprintf (f,'%s %f %f %f %s \n',' (',-y,x,l,')');
fprintf (f,'%s\n',');');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','blocks');
fprintf (f,'%s\n','(');
fprintf (f,'%s\n',' hex (0 1 2 3 4 5 6 7) (1 50 200) simpleGrading (1 0.1 1)');
fprintf (f,'%s\n',');');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','edges');
fprintf (f,'%s\n','(');
fprintf (f,'%s\n',');');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','boundary');
fprintf (f,'%s\n','(');
fprintf (f,'%s\n',' inlet');
fprintf (f,'%s\n',' {');
fprintf (f,'%s\n',' type patch;');
fprintf (f,'%s\n',' faces');
fprintf (f,'%s\n',' (');
fprintf (f,'%s\n',' (0 3 2 1)');
fprintf (f,'%s\n',' );');
fprintf (f,'%s\n',' }');
fprintf (f,'%s\n','');
fprintf (f,'%s\n',' outlet');
fprintf (f,'%s\n',' {');
fprintf (f,'%s\n',' type patch;');
fprintf (f,'%s\n',' faces');
fprintf (f,'%s\n',' (');
fprintf (f,'%s\n',' (4 7 6 5)');
fprintf (f,'%s\n',' );');
fprintf (f,'%s\n',' }');
fprintf (f,'%s\n','');
fprintf (f,'%s\n',' fixedWalls');
fprintf (f,'%s\n',' {');
fprintf (f,'%s\n',' type wall;');
fprintf (f,'%s\n',' faces');
fprintf (f,'%s\n',' (');
fprintf (f,'%s\n',' (3 2 6 7)');
fprintf (f,'%s\n',' );');
fprintf (f,'%s\n',' }');
fprintf (f,'%s\n',' wedge1');
fprintf (f,'%s\n',' {');
fprintf (f,'%s %s\n',' type',str); %<-- str variable stores the type of Boundary conditions to applied on the wedge patch
fprintf (f,'%s\n',' faces');
fprintf (f,'%s\n',' (');
fprintf (f,'%s\n',' (0 3 7 4)');
fprintf (f,'%s\n',' );');
fprintf (f,'%s\n',' }');
fprintf (f,'%s\n',' wedge2');
fprintf (f,'%s\n',' {');
fprintf (f,'%s %s\n',' type',str);
fprintf (f,'%s\n',' faces');
fprintf (f,'%s\n',' (');
fprintf (f,'%s\n',' (1 5 6 2)');
fprintf (f,'%s\n',' );');
fprintf (f,'%s\n',' }');
fprintf (f,'%s\n',' line1');
fprintf (f,'%s\n',' {');
fprintf (f,'%s\n',' type empty;');
fprintf (f,'%s\n',' faces');
fprintf (f,'%s\n',' (');
fprintf (f,'%s\n',' (0 1 5 4)');
fprintf (f,'%s\n',' );');
fprintf (f,'%s\n',' }');
fprintf (f,'%s\n','');
fprintf (f,'%s\n',');');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','mergePatchPairs');
fprintf (f,'%s\n','(');
fprintf (f,'%s\n',');');
fprintf (f,'%s\n','');
fprintf (f,'%s\n','// ************************************************************************* //');
fprintf (f,'%s\n','');Explanation of the Matlab code:
In the above code, the desired wedge angle is defined for the theta variable which is used to calculate the vertices for the wedge angle using trigonometric functions, and also the length of the pipe is defined in the variable l.
BlockMeshDict.txt file is created to which codes are written which are similar to the original blockMeshDict which includes the vertices for the model, defining of the model satisfying right-hand thumb rule and also defining the type of boundary conditions along with the faces of the wedge.
blockMeshDict for wedge angle 10 degrees for symmetry Boundary Conditions.
/*--------------------------------*- 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0 0)
( 0.000436 0.004981 0.000000 )
( -0.000436 0.004981 0.000000 )
( 0.000000 0.000000 2.000000 )
( 0.000000 0.000000 2.000000 )
( 0.000436 0.004981 2.000000 )
( -0.000436 0.004981 2.000000 )
);
blocks
(
hex (0 1 2 3 4 5 6 7) (1 50 200) simpleGrading (1 0.1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 3 2 1)
);
}
outlet
{
type patch;
faces
(
(4 7 6 5)
);
}
fixedWalls
{
type wall;
faces
(
(3 2 6 7)
);
}
wedge1
{
type symmetry
faces
(
(0 3 7 4)
);
}
wedge2
{
type symmetry
faces
(
(1 5 6 2)
);
}
line1
{
type empty;
faces
(
(0 1 5 4)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Explanation:
In the above code, the first seven lines contain the header information and the code information is in a foam file which is enclosed in '{}'.
After the file information, the dimensions are converted into meters by multiplying the scaling factor by 1.
Then vertices for wedge block are defined using the Right-hand thumb rule in which the z-axis is normal to the surface.
Blocks are defined after the definition of vertices for the blocks where the hex is defining the shape of the block which is hexahedra. The order of the vertices in defining the block should be done carefully because this line also specifies a local axis for the vertices to which the number of cells is defined along with grading. Here simple grading is cell expansion. In the parenthesis, we specify the ratio which is the width of the end cell along one edge of a block to the width of the start cell along that edge.
There are no curved edges specified in the geometry i.e., no points are specified within the parenthesis.
Boundary Conditions:
For the inlet and outlet, the type of boundary conditions specified are:
type: patch type on which particular boundary conditions are applied and faces: mean a list of block faces that make up the patch. The type of Boundary conditions chosen is inlet velocity and outlet pressure.
FixedWalls is of wall type since particular physical properties are defined at these walls, in this case, is the dynamic viscosity in `kg.m^-1s^-1`
Wedge 1 and Wedge2 walls are of type symmetry to simulate the flow in 3 dimensional. The symmetry type boundary conditions allow for the slip for the flow variables at the wedge surface.
Wedge vertices for wedge angle 25 degrees:
vertices
(
(0 0 0)
(0 0 0)
( 0.001082 0.004881 0.000000 )
( -0.001082 0.004881 0.000000 )
( 0.000000 0.000000 2.000000 )
( 0.000000 0.000000 2.000000 )
( 0.001082 0.004881 2.000000 )
( -0.001082 0.004881 2.000000 )
);
Wedge vertices for wedge angle 45 degrees:
vertices
(
(0 0 0)
(0 0 0)
( 0.001913 0.004619 0.000000 )
( -0.001913 0.004619 0.000000 )
( 0.000000 0.000000 2.000000 )
( 0.000000 0.000000 2.000000 )
( 0.001913 0.004619 2.000000 )
( -0.001913 0.004619 2.000000 )
);Initial Conditions for Pressure with Symmetry boundary conditions:
/*--------------------------------*- 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;
}
fixedWalls
{
type zeroGradient;
}
wedge1
{
type symmetry;
}
wedge2
{
type symmetry;
}
line1
{
type empty;
}
}
// ************************************************************************* //Explanation:
In the above code, the first 7 lines represent the header information and the next nine-lines contain file information for which heading is Foam file.
The dimension line represents the units of the physical quantity to which the boundary conditions are applied. Here it is kinematic pressure i.e., `m^2s^-2`
Internal field species the value of the internal field in which a single value is assigned to all elements within the field.
Then Initial Conditions are specified as,
The type of boundary conditions considered is Velocity inlet and pressure outlet boundary conditions since the flow is incompressible.
Hence for the inlet, the pressure is zero gradients which means the normal gradient of P is zero which represents the Neumann boundary conditions and a fixed value of zero pressure is specified which represents the static pressure. This is used when the fluid flows outside the domain at a patch face and a fixed value is used when the fluid flows into the domain.
Since no-slip walls are connected to the patch at which physical condition is defined it is specified as zero gradients and for wedge 1 and wedge 2 symmetry type boundary conditions are given since the flow is to be simulated in 3D.
line1 is the bottom line of the wedge. As it is represented with a face in this problem and no solution is required type is empty.
Initial Conditions for Velocity with Symmetry boundary conditions:
/*--------------------------------*- 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 0 0.17);
}
outlet
{
type zeroGradient;
}
fixedWalls
{
type noSlip;
}
wedge1
{
type symmetry;
}
wedge2
{
type symmetry;
}
line1
{
type empty;
}
}
// ************************************************************************* //Explanation:
The unit for the velocity is m/s which is specified in the dimensions as [0 1 -1 0 0 0 0]
Internal points in the model are specified as 0.
In the velocity boundary conditions since there is an initial velocity a fixed value of 0.17 m/s is defined at the inlet which is calculated from the Reynolds number formula where Re = 2000 since the flow considered is laminar.
At the outlet, the velocity is specified as zero-gradient since the flow attains the steady-state when it reaches the outlet of the model.
No-slip is specified at the fixed wall because of the boundary layer. noSlip fixes the values of the velocity to zero at the walls.
For wedge 1 and wedge 2 it is type symmetry because slip is allowed at the wedge surface.
line1 is the bottom line of the wedge. As it is represented with a face in this problem and no solution is required type is empty.
ControlDict:
/*--------------------------------*- 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 20;
deltaT 0.0085;
writeControl timeStep;
writeInterval 6;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
#includeFunc residuals
}
// ************************************************************************* //Explanation:
This code contains the start time, end time, time step size which can be read from the controlDict dictionary in the system directory.
Start time starts from 0 which means that OpenFoam needs to read field data from a directory named 0. For the fully developed flow, the end time is chosen as 20 seconds.
The next step is to define the time step size which can be calculated from the courant number formula `C= (trianglet*U)/(trianglex)`. where C =1, U =0.17 m/s,
`trianglex = d/n`
where d is the length of the model and n is no. of cells the model is divided into in x-direction.
`trianglex = 0.01`
`trianglet` = 0.01 which when used resulted in the error due to an increase in courant number hence the time step size is decreased further to obtain a stable solution.
The `trianglet` obtained is done for simple grading 0.1. With the change in grading factor`trianglex` changes which in turn changes the courant number.
Write interval 6 in the code specifies the program to store the values after 6-time step iterations.
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.55e-7;
// ************************************************************************* //Explanation:
The only physical property for the solver icoFoam is kinematic viscosity which is stored in the transport properties in the constant directory.
For this case, water properties at `30^o` C is chosen at which
density (`rho`) = 996.6 `(kg)/m^3`.
dynamic viscosity (`mu`) = 0.00085 `(kg)/(ms)`
`nu = (mu)/(rho)`
`nu = 0.000852/996.6`
`nu = 8.55 e^-7`
Steps in the execution of the problem:
1. Choosing the required solver from the tutorial directory. In this case, the given problem is incompressible transient laminar flow hence icoFoam solver is used.
2. Cavity model is used as it is simple and easy to modify.
3. Copy the cavity directory to the Foam run directory because no changes have to be made in the tutorial directories.
4. Inside the run folder the navigate to the cavity directory and to the system directory to modify the BlockmeshDict with the changes mentioned above and similarly navigate to controlDict to make the necessary changes.
5. Apply initial conditions by navigating to the 0 directory and then use the gedit command to edit the Pressure and velocity initial conditions.
6. Once the grading is applied the mesh has to be updated using blockMesh. This command is repeated each time the grading is changed.
7. After execution of blockMesh command if there is no error then the solver is executed with icoFoam command.
8. After successful execution the post-processing is done using the paraFoam command which opens the ParaView window
9. In the ParaView window the flow is analyzed for each flow variable.
Matlab script for comparing the velocity and shear stress for different wedge angles for wedge and symmetry boundary conditions.
close all
clear all
clc
%wedge BC
v10 = readmatrix ('v10.CSV');
u10 = v10(:,1);
r10 = v10(:,2);
v25 = readmatrix ('v25.CSV');
u25 = v25(:,1);
r25 = v25(:,2);
v45 = readmatrix ('v45.CSV');
u45 = v45(:,1);
r45 = v45(:,2);
v45n = readmatrix('v45n.CSV');
u45n = v45n(:,1);
r45n = v45n(:,2);
%Analytical Solution
Vavg=0.17; %<-- Average velocity in m/s
R= 0.005; %<-- Radius of the pipe in m.
u = 2*Vavg*(1-(r10).^2./(R)^2); %<-- Formula to calculate the velocity profile.
figure(1)
plot (u10,r10)
hold on
plot (u25,r25)
plot (u45,r45)
plot (u45n,r45n)
plot (u,r10,'linewidth',2,'color','k')
xlabel ('Velocity in m/s')
ylabel ('Radial distance in m')
title ('Velocity Vs Radial Distance for Wedge BC')
legend ('Wedge angle 10 degree','Wedge angle 25 degree','Wedge angle 45 degree','Improved Velocity wedge angle 45 degree','Analytical solution')
%symmetry BC
v10s = readmatrix ('v10s.CSV');
u10s = v10s(:,1);
r10s = v10s(:,2);
v25s = readmatrix ('v25s.CSV');
u25s = v25s(:,1);
r25s = v25s(:,2);
v45s = readmatrix ('v45s.CSV');
u45s = v45s(:,1);
r45s = v45s(:,2);
v45sn = readmatrix ('v45sn.CSV');
u45sn = v45sn (:,1);
r45sn = v45sn(:,2);
figure(2)
plot (u10s,r10s)
hold on
plot (u25s,r25s)
plot (u45s,r45s)
plot (u45sn,r45sn)
plot (u,r10,'linewidth',2,'color','k')
xlabel ('Velocity in m/s')
ylabel ('Radial distance in m')
title ('Velocity Vs Radial Distance for Symmetry BC')
legend ('Wedge angle 10 degree','Wedge angle 25 degree','Wedge angle 45 degree','Improved Velocity profile','Analytical Solution')
%Shear stress for wedge BC
ss10 = readmatrix('ss10.CSV');
s10 = ss10(:,2);
ss25 = readmatrix ('s25.CSV');
s25 = ss25(:,2);
ss45 = readmatrix('s45.CSV');
s45 = ss45 (:,2);
ss45n = readmatrix('s45n.CSV');
r45n = ss45n(:,1);
s45n = ss45n(:,2);
%Calculation of the shear stress analytically.
dv = 0.000852; %<-- dynamic viscoity in kg/ms.
rd = linspace (0,0.005,1001); %<-- defining the radial distance
r = rd(2)-rd(1); %<-- length between the adjacent points.
s = dv*(diff(u)/r); %<-- Formula to calculate the shear stress analytically
s = flipud(s); %<-- the shear stress obtained is flipped
radial = linspace (0,0.005,1000); %<-- defining radial distance for the plot
figure(3)
plot (s10,rd)
hold on
plot (s25,rd)
plot (s45,rd)
plot (s45n,rd)
plot (s,radial)
xlabel('Shear Stress in N/m^2')
ylabel ('radial distance in m')
title ('Shear Stress Vs radial distance for Wedge BC')
legend ('Wedge angle 10 degree','Wedge angle 25 degree','Wedge angle 45 degree','Improved velocity wedge angle 45 degree','Analytical Solution')
%Shear stress for symmetry BC
ss10s = readmatrix('ss10s.CSV');
s10s = ss10s(:,2);
ss25s = readmatrix ('s25s.CSV');
s25s = ss25s(:,2);
ss45s = readmatrix('s45s.CSV');
s45s = ss45s (:,2);
s45sn = readmatrix ('s45sn.CSV');
ss45sn = s45sn (:,2);
%Calculation of the shear stress analytically.
dv = 0.000852; %<-- dynamic viscoity in kg/ms.
rd = linspace (0,0.005,1001); %<-- defining the radial distance
r = rd(2)-rd(1); %<-- length between the adjacent points.
ss = dv*(diff(u)/r); %<-- Formula to calculate the shear stress analytically
ss = flipud(ss); %<-- the shear stress obtained is flipped
radial = linspace(0,0.005,1000); %<-- defining radial distance for the plot
radial = (radial)';
figure(4)
plot (s10s,rd)
hold on
plot (s25s,rd)
plot (s45s,rd)
plot (ss45sn,rd)
plot (ss,radial)
xlabel('Shear Stress in N/m^2')
ylabel ('radial distance in m')
title ('Shear Stress Vs radial distance for Symmetry BC')
legend ('Wedge angle 10 degree','Wedge angle 25 degree','Wedge angle 45 degree','improved Shear Stress','Analytical Solution')Explanation of the Matlab script:
1. The velocity profile data is loaded into Matlab using the readmatrix command.
2. Assigning the velocity data and radial distance data to the appropriate variables.
3. Calculating the analytical velocity with the formula `u = 2.V_(avg)*(1-(r^2)/R^2)`
4. Plot is created for the different velocity profiles created for different wedge angles for both wedge and symmetry boundary conditions.
5. Inorder to plot the shear stress numerically the following steps are done in post-processing in the ParaView window.
- The magnitude of Velocity is calculated with the calculator in the common toolbar.
- With the magnitude of velocity is known the gradient of the velocity is calculated using the gradients of the unstructured data set filter from the filter search option.
- Use the calculator from the common toolbar to calculate the shear stress using the formula `tau = mu*((du)/(dx))`
- While analyzing the flow choose the shear stress variable in the field variable box.
- Use plot over line command to create a plot of shear stress against the arc length.
- Extract the shear stress using the export scene in the file menu.
6. Load the shear stress into the Matlab and assign the appropriate variables to shear stress values and radial distance values.
7. For the calculation of the analytical values of shear stress, use the same formula used in the calculator in the ParaView window. But the velocity for analytical calculation of shear stress is analytically calculated velocity values.
8. A plot is created to compare the shear stress numerically and analytically for the fully developed flow.
Mesh at the inlet of the pipe:
For 10 degrees wedge angle
For 25 degrees wedge angle:
For 45 degrees wedge angle:
For 45 wedge angle with an increase in the number of cells:
Side view of the model
Flow simulation
At inlet:
At outlet:
Output Graphs:
Velocity profiles using wedge-type boundary conditions:
Velocity profiles for the symmetry type boundary conditions:
Shear Stress profiles for wedge boundary conditions:
Shear Stress profiles for symmetry type boundary conditions:
Explanation of the graphs:
From the above plots, we can see that the numerical solution is close to the analytical solution. There is a small variation in numerical solution in the analytical solution this is because of the discretization error and the round-off error induced by the computer.
From the velocity plot, it can be inferred that the velocity is zero at the wall because of the velocity boundary layer and increases as the flow moves towards the center of the pipe.
From the shear stress plot, we can infer that the shear stress is maximum at the wall due to the viscosity of the fluid and decreases along the radial direction of the pipe, and is zero at the center of the pipe.
With decrease in the wedge angle, the numerical solution is closer to the analytical solution. For the wedge angle 45 degrees, with an increase in the number of cells the numerical solution results in unrealistic solutions which can be compared in the below comparison table.
Comparison Table:
| Wedge Boundary conditions | Symmetry Boundary conditions | Analytical Solution | |||||||
| Particulars | 10-degree wedge angle | 25-degree wedge angle | 45-degree wedge angle | Improved mesh for 45-degree wedge angle | 10-degree wedge angle | 25-degree wedge angle | 45-degree wedge angle | Improved mesh for 45-degree wedge angle | |
| Velocity in m/s | 0.33912 | 0.33917 | 0.33927 | 0.5031 | 0.33912 | 0.33917 | 0.33927 | 0.3382 | 0.34 |
| Shear Stress in N/m^2 | 0.11499 | 0.11745 | 0.12408 | 0.177 | 0.11511 | 0.11745 | 0.12408 | 0.1165 | 0.1149 |
| Pressure in N/m^2 | 0.02348 | 0.02444 | 0.02725 | 0.0370 | 0.02348 | 0.02444 | 0.02726 | 0.0237 | 0.0232 |
Conclusion:
From the above comparison table, we can infer that numerical solutions are closer to the analytical solution for smaller wedge angle and also for the improved mesh for 45-degree wedge angle that is with increase in the number of mesh cells accuracy of the solution increases. As the wedge angle increases the accuracy of the solution decreases.
When the angle greater than 10 degrees is used, to obtain an accurate solution use a finer mesh grid which results in values closer to the analytical solution as in this problem for wedge angle 45 degrees.