Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
AIM
To simulate laminar flow through pipe having constant cross-section using symmetry boundary condition by taking the wedge angle of 10o,25o and 45 and also comparing it with wedge boundary condition.
OBJECTIVE
- To Generate BlockMesh file using the MATLAB program.
- Simulation of flow through wedge with a given wedge angle to determine the velocity profile, shear stress distribution, maximum velocity and pressure drop using openfoam with symmetry boundary condition and wedge boundary condition.
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- validate wall shear stress for fully developed flow
INTRODUCTION
Liquid or gas flow through pipes or ducts is commonly used in heating and cooling applications and fluid distribution networks. The fluid in such applications is usually forced to flow by a fan or pump through a flow section. We pay particular attention to friction, which is directly related to the pressure drop and head loss during flow through pipes and ducts. The pressure drop is then used to determine the pumping power requirement. A typical piping system involves pipes of different diameters connected to each other by various fittings or elbows to route the fluid, valves to control the flow rate, and pumps to pressurize the fluid.
The fluid velocity in a pipe changes from zero at the wall because of the no-slip condition to a maximum at the pipe centre. In fluid flow, it is convenient to work with an average velocity `V_(avg)`, which remains constant in incompressible flow when the cross-sectional area of the pipe is constant. The average velocity in heating and cooling applications may change somewhat because of changes in density with temperature. But, in practice, we evaluate the fluid properties at some average temperature and treat them as constants.
THEORY
Symmetry Boundary Condition
This boundary condition is used to apply mirror-symmetry conditions to a structure.
It can be applied to the faces of a structure and no other user input is needed. If a symmetry plane condition is applied to a face, the displacement of this face is locked in a normal direction but free to slide in tangential directions.
During the time of numerical analysis, it calculates the flow field variables for the complete cell volume rather than taking a single plane as in Wedge BC.
The symmetry BC allows us to use multiple meshes keeping in mind the Courant number(more the meshes, lesser is the dx value), to be within the stability conditions.
For the Symmetry BC, the angular changes that are considered don't affect or create an unsteady pattern for the path of attaining the convergence.
Wedge Boundary Condition
Wedge means ' taking small portion of entire object.
The wedge boundary condition is used in computational fluid dynamics to define an axisymmetric situation, for eg. a cylinder. The model and flow must be axisymmetric along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.
This boundary condition is specified by two planes that must be selected on separate sides of the domain (referred to as front and back) on surfaces running along the axis. It works for 2D simulations only.
A Wedge boundary is different from the symmetry boundary condition in the sense that the wedge boundary condition applies the Navier-Stokes (NS) equations in cylindrical coordinates whereas the symmetric applies the NS equations in cartesian coordinates. One big advantage of wedge boundary is that it allows for the rotational flow around the axis and this is possible only because the fluxes and normal components are not assumed as zero.
The wedge BC considers only a single straight mid-plane for a given wedge-shaped geometry(the triangular kind of looking shape) and represents the simulated values of the flow field variables only through that plane. The plane looks like when it is rotated 360 degrees it forms a full cylinder.
Usage of Wedge BC limits our input value of a number of meshes along the wedge. Usage of many meshes in Wedge BC leads to errors in the result.
Laminar Flow
Laminar or streamline flow in pipes or tubes occurs when a fluid flows in parallel layers, with no disruption between the layers. At low velocities, the fluid tends to flow without lateral mixing, and adjacent layers slide one another. There are no cross-currents perpendicular to the direction of flow, no swirls of fluids. In laminar flow, the motion of the particles of the fluid is very orderly with all particles moving in straight lines parallel to the pipe walls.
Turbulent Flow:
Turbulent flow is a flow regime characterized by chaotic property changes. This includes a rapid variation of pressure and flows velocity in space and time. Turbulent flow is a type of fluid flow in which the fluid undergoes irregular fluctuations, or mixing in contrast to laminar flow, in which the fluid moves in smooth paths or layers. In turbulent flow, the speed of the fluid at a point is continuously undergoing changes in both magnitude and direction.
If we have been around smokers, we probably noticed that the cigarette smoke rises in a smooth plume for the first few centimetres and then starts fluctuating randomly in all directions as it continues its rise. So, that is one of the example of laminar and turbulent flow.
Reynolds number
The Reynolds number represents the ratio of inertial forces to viscous forces and is a convenient parameter for predicting if a flow condition will be laminar or turbulent. It is defined as a characteristic length multiplied by a characteristic velocity and divided by the kinematic viscosity.
`Re_D="Inertial Forces" /"Viscous Forces"`
`=(ρVD)/μ=(VD)/ν`
Where,V = Velocity of the Fluid (m/s)
ρ = Density of the Fluid (kg/m3)
D= Characteristic Length (m)
μ = Dynamic Viscosity (Pa·s or N·s/m2 or kg/(m·s))
ν = Kinematic Viscosity (m2/s).
The Reynolds number at which the flow becomes turbulent is called the critical Reynolds number, Recr. The value of the critical Reynolds number is different for different geometries and flow conditions. For internal flow in a circular pipe, the generally accepted value of the critical Reynolds number is Recr=2300. Under most practical conditions, the flow in a circular pipe is laminar for Re ≲ 2300, turbulent for Re ≳ 4000, and transitional in between. That is,
(Re ≲ 2300) - laminar flow
(2300 ≲ Re ≲ 4000) - transitional flow
(Re ≳ 4000) - turbulent flow
Hydro-dynamic length
In fluid dynamics, the entrance length is the distance a flow travels after entering a pipe before the flow becomes fully developed.
Entrance length refers to the length of the entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer. When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer change with increased distance along the pipe. Many different entrance lengths exist to describe a variety of flow conditions.
Hydrodynamic entrance length describes the formation of a velocity profile caused by viscous forces propagating from the pipe wall. Thermal entrance length describes the formation of a temperature profile. Awareness of entrance length may be necessary for the effective placement of instrumentation, such as fluid flow meters.
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. This region is characterized by a non-uniform flow. The fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no-slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile. For the conservation of mass to hold true, the velocity of layers of the fluid in the center of the pipe increases to compensate for the reduced velocities of the layers of fluid near the pipe surface. This develops a velocity gradient across the cross-section of the pipe.
Shear Stress
In the hydrodynamic entrance region, the wall shear stress (τw ) is highest at the pipe inlet, where the boundary layer thickness is the smallest. Shear stress decreases along the flow direction. That is why the pressure drop is highest in the entrance region of a pipe, which increases the average friction factor for the whole pipe. This increase in the friction factor is negligible for long pipes.
In a fully developed region, the pressure gradient and the shear stress in flow are in balance.
Now consider flow through circular pipe with uniform velocity,
As a result of the no-slip conditions along the walls of the pipe, the fluid particles that come in contact with the wall come to a complete stop. As a result of friction, this causes the adjacent layers of the fluid to slow down as well. To make up for the reduction in velocity, the fluid velocity at the mid-section of the pipe has to increase to maintain a constant mass flow rate through the pipe.
The various regions of the pipe as shown in the above diagram are as follows,
- The irrotational (core) flow region is one in which the frictional effects are negligible and the velocity remains essentially constant in the radial direction.
- The region of flow in which the effects of viscous shearing forces caused by fluid viscosity are felt and the velocity changes are significant, is referred to as the velocity boundary layer.
- The region from the pipe inlet to the point at which the velocity profile is fully developed is called the hydrodynamic entrance region and the length of this region is referred to as the hydrodynamic entry length.
- Flow in the entrance region is called hydrodynamically developing flow since it is where the velocity profile develops. The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fully developed region.
Hagen- Poiseuille's Theory
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter; and there is no acceleration of fluid in the pipe.
The pressure drop in an incompressible and Newtonian fluid in laminar flow through a long cylindrical pipe of constant cross – section can be derived as,
`∆P=P_2-P_1=(8μLV_(avg))/R^2`
`=(32μLV_(avg))/D^2`
Where μ is the fluid viscosity, L is the length of the pipe(> entry length), Vavg is the average velocity, R is the radius of the pipe and D is the diameter of the pipe. The Hagen–Poiseuille equation describes the relationship between pressure, fluidic resistance, and flow rate, analogous to voltage, resistance, and current, respectively, in Ohm’s law for electrical circuits (V=RI).
The ratio of length to radius of a pipe should be greater than one forty-eighth of the Reynolds number for the Hagen–Poiseuille law to be valid. If the pipe is too short, the Hagen–Poiseuille equation may result in unphysically high flow rates.
A few examples in which the Hagen–Poiseuille law is practically applied is the air flow in lung alveoli, or the flow through a drinking straw or through a hypodermic needle. It is also applied for the process of blood flow through a narrow constriction.
METHODOLOGY
In this project, we will be simulating the same case with different boundary conditions of Wedge and Symmetry for angles of 10, 25 and 45 respectively.
The Wedge boundary condition is used for 2-dimensional axisymmetric simulations, where the Navier-Stokes equations are taken in the Cylindrical co-ordinates. This boundary conditions yields an advantage of simulating a part of geometry with respect to its axis.
While the Symmetry boundary condition is used for the convenience of mirroring of the surface of the geometry which yields a decrease in the computational cost. The Symmetry boundary condition takes the Navier-Stokes equations in Cartesian co-ordinates.
The workflow goes as follows:
- We will first generate our BlockMeshDict files with the Matlab code
- Then, run the simulations in OpenFoam
- Post-processing the data through analysing the resulting plots using Paraview
- And finally, validating the Numerical solutions with the Analytical solutions with the help of Matlab.
PROCEDURE
- As this is a problem of incompressible laminar viscous flow we use cavity folder in tutorial.
- Create a copy of cavity folder from tutorial in the $FOAM_RUN
- Then open folder which has 3 folders 0 ,constant and system.
- Open system folder to open notepad.exe blockmeshdict text file then specify the vertices as per the given dimensions of pipe.
- Specify the hex block mesh with given mesh specification with grading factor of 1 0.2 1 in x y and z directions respectively.
- Specify the inlet and outlet with type patch, front and back walls with type wedge and noslipwall with type wall and axis with type empty
- Save the blockMeshdict text file.
- Type blockMesh in the windows terminal to create meshing.
- Once you create meshing, open controldict text file and specify the start time ,end time,time step, dt and write interval.
- Go to folder 0
- Open U textfile and specify the velocity at the inlet as 0.2113m/s,outlet as zerogradient ,front and back as wedge and noslipwalls as noslip and axis as empty and save.
- Open p textfile and specify the pressure at the inlet as zerogradient,outlet as fixed and uniform 0.17 N/m2 ,front and back as empty and noslipwalls as zerogradient,and axis as empty and then save.
- After specifying the velocity and presssure as boundary conditions open folder constant and open textfile transport properties
- Specify the value of kinematic viscosity then save.
- Then come back to main run folder and type blockMesh and then check the quality of mesh.
- Run icofoam application for simulation.
- open paraview and plot the graph for velocity profile ,shear stress ,pressure drop and hydrodynamic length for the comparsion with numerical result using matlab.
Repeat the above procedure for θ=25o&45o with wedge boundary condition
follow the same procedure for θ=10o,25o&45o with Symmetry boundary condition
Calculations for Analytical Solution
Assumption
- Flow is steady, incompressible (density is constant) and laminar.
- Constant fluid properties like Viscosity and density.
- Working temperature of water 25oC
- Reynolds Number = 2100 (Given)
- Diameter = 2cm = 0.02m (Assumed)
Analytical Solution
- Diameter of the pipe `D=0.02m → R=D/2=0.01m`
- Density of working fluid (water) at 25oC `ρ=997 (Kg)/m^3`
- Wedge angle `θ=10^o ,25^o 45^o`
- Dynamic Viscosity of water at 25oC `μ=8.92×10^(-4) Pa.s`
- Dynamic Viscosity of water at 25oC `ν=μ/ρ=(8.90×10^(-4))/997=8.92×10^(-7) m^2/sec`
- Hydrodynamic Entrance Length `L_h=0.06⋅R_e⋅D= =0.06×2100×0.02=2.52m`
- Total Length of the Pipe `L=2.8m`
- Average Velocity `V_avg=(μ⋅R_e)/(ρ⋅D)==(8.90×10^(-4) × 2100)/(997×0.02)=0.0937 m/s`
- Maximum Velocity `Δp = (32⋅μ⋅V_(avg)⋅L)/D^2 ==(32×8.90×10^(-4)×0.0937×2.8)/(0.02)^2 =18.680 Pa`
- Velocity Profile `u(r)=2⋅V_(avg)⋅(1–r^2/R^2 )`
- Pressure Drop `Δp = (32⋅μ⋅V_(avg)⋅L)/D^2 =(32×8.90×10^(-4)×0.0937×2.8)/(0.02)^2 =18.680 Pa`
- Kinematic Pressure Drop `Δp=(Δp)/ρ=18.680/997=0.0187 Pa`
- Shear Stress `τ=(2⋅μ⋅V_max⋅R_e)/R^2 =(2×8.90×10^(-4)×0.1874×2100)/0.01^2 =7005.012 Pa`
Case 1: Wedge boundary condition with angle 10⁰
Matlab code to create blockMeshDict for Wedge boundary condition with angle 10⁰
Clear all
Close all
clc
%% INPUTS
% Diameter of Pipe (m)
D = 0.02;
% Radius of Pipe (m)
R = D/2;
% Reynold's Number
Re = 2100;
% Wedge Angle (degrees)
theta = 10;
% Total Length of the Pipe (m)
L = 2.8;
% Assuming working fluid is water at 25° C
% Density of water (kg/m³)
rho = 997;
% Dynamic Viscosity of Water (Pa.s)
mu = 8.90 * 10^-4;
% Kinematic Viscosity of Water (m^2/s)
nu = mu/rho;
% Hydrodynamic Entrance Length (m)
L_h = 0.06 * Re * D;
%% HAGEN - POISEUILLE'S FLOW EQUATIONS
% Average Velocity (m/s)
V_avg = (mu * Re)/(rho * D);
% Maximum Velocity (m/s)
V_max = 2 * V_avg;
% Pressure Drop (Pa)
delta_p = (32*mu*V_avg*L)/(D^2);
% Kinematic Pressure Drop (Pa)
kinematic_delta_p = delta_p/rho;
%% Calculating shear stress
r = linspace(0,(D/2),500);
for i = 1:length(r)
tau(i) = 2 * mu * V_max * abs(r(i))/R^2;
end
figure(1)
plot(r,tau)
xlabel('Radius of pipe')
ylabel('shear stress')
title('shear stress')
legend('shear stress')
%% DEFINING THE GEOMETRY
% 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 0 0];
V4 = [L R*cosd(theta/2) -R*sind(theta/2)];
V5 = [L R*cosd(theta/2) R*sind(theta/2)];
%% BLOCKMESHDICT FILE
% Creating the file
% Header
h1 = '/*--------------------------------*- C++ -*----------------------------------*\';
h2 = ' ========= |';
h3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4 = ' \\ / O peration | Website: https://openfoam.org';
h5 = ' \\ / A nd | Version: 9';
h6 = ' \\/ M anipulation |';
h7 = '\*---------------------------------------------------------------------------*/';
h8 = 'FoamFile';
h9 = '{';
h10 = ' format ascii;';
h11 = ' class dictionary;';
h12 = ' object blockMeshDict;';
h13 = '}';
h14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%% Creating BlockMesh File
f1 = fopen('blockMeshDict.txt','w');
% Header
fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'%s\n',h12);
fprintf(f1,'%s\n',h13);
fprintf(f1,'%s\n\n',h14);
% Conversion Factor
fprintf(f1,'%s\n\n','convertToMeters 1;');
% Vertices
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t(%d %d %d)\n',V0);
fprintf(f1,'\t(%d %f %f)\n',V1);
fprintf(f1,'\t(%d %f %f)\n',V2);
fprintf(f1,'\t(%.2f %d %d)\n',V3);
fprintf(f1,'\t(%.2f %f %f)\n',V4);
fprintf(f1,'\t(%.2f %f %f)\n',V5);
fprintf(f1,'%s\n\n',');');
% Blocks
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (400 30 1) simpleGrading (1 0.1 1)\n');
fprintf(f1,'%s\n\n',');');
% Edges
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\tarc 2 1 (%d %.3f %d)\n',0,R,0);
fprintf(f1,'\tarc 5 4 (%.3f %.3f %d)\n',L,R,0);
fprintf(f1,'%s\n\n',');');
% Boundary
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
% Axis
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Inlet
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 1 2 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Outlet
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Front Face
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 4 1)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Back Face
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 2 5 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Top Wall
fprintf(f1,'\tpipeWall\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(1 4 5 2)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');
Output of MATLAB code :-
BlockMeshDict file for Wedge boundary condition with angle 10⁰:-
The blockMeshDict file has been written automatically by the MatLab program as shown above, which we use to create the geometry of the wedge and run the simulation for the same.
ControlDict file
In controlDict file we will change End Time, Delta T and Write Interval values.
- End Time :- To decide suitable end time we will first decide minimum time required to run simulation by formula Velocity = Dist / Time
Time required = Total length / Avg velocity
= 2.8/0.0937 =29.88 Sec
So, we will put value of end time 30 Sec
- Delta T :- We will check for 0.01 if found error we will change value according to required number of Iterations to solve.
- In this case for Delta T = 0.01 it took so much time to run so we will change Delta T = 0.05
- Write Interval :- This defines time after which new interval decides and new folder created for that each interval. So if interval value is less, solver take so many intervals so we will increase value accordingly.
- In this case for Write Interval = 20, it took so many interval to run so we will change Write Interval = 200
We did not make any other changes in ControlDict file and save it.
We will not change in any remaining part in ControlDict file
U file (Velocity) :-
We will make changes in velocity file located in
“0” folder --> U file
We will make changes in boundary field and use same names for boundary fields as used in BlockMeshDict file.
We also define type of boundary field as per geometry.
|
Boundary |
Type |
Description |
|
Axis |
empty |
Empty Patch field |
|
Inlet |
fixedValue |
Predefining values of variables |
|
Outlet |
zeroGradient |
All variable are set to constant |
|
Front |
wedge |
define an axisymmetric situation |
|
Back |
wedge |
define an axisymmetric situation |
|
PipeWall |
noSlip |
No slip in velocity at wall |
After changing all boundary fields and type we will change value of Inlet for velocity. We are considering average velocity of fluid at inlet so change type of inlet to fixedValue and change value of inlet to 0.0937 in X-direction.
P file (Pressure) :-
Pressure file located in “0” folder --> P file
We will change in boundary fields as done in Velocity file only difference is that we are not using predefined value for pressure at inlet. We are keeping pressure constant at inlet and changing value of pressure to zero at outlet as shown below.
Physical Properties file :-
Located in “Constant” folder --> physicalProperties file
As we are considering water as fluid at 25⁰ C so We will change value of kinematic viscosity to 8.92e-7 i.e. (8.92 x 10^-7)
- Now after saving all these changes in editor, we will create Mesh in OpenFoam by running blockMesh command in
- After running blockMesh successfully we will run icoFoam command in terminal.
- After completing all iterations we will run touch command to create and save Paraview model file with “name.foam” extension.
- For this case : “touch wedge_10.foam”
Now we will run paraFoam command to see ParaView of model
ParaView of model :
Detailed View :-
Velocity Distribution :-
Pressure Distribution :-
Pressure Profile :-
Graph of Shear Stress distribution from MatLab:-
Shear Stress distribution in ParaView :-
Wedge Boundary condition with angle 25⁰ :-
For 25⁰ angle, we will copy entire folder created for wedge boundary condition with 10⁰ angle including Matlab code.
We will make changes in Matlab code created for 10⁰ and change value of theta to 25 from 10, as shown below.
Clear all
Close all
clc
%% INPUTS
% Diameter of Pipe (m)
D = 0.02;
% Radius of Pipe (m)
R = D/2;
% Reynold's Number
Re = 2100;
% Wedge Angle (degrees)
theta = 25;
% Total Length of the Pipe (m)
L = 2.8;
We will not make any other changes in remaining code.
Due to changing in theta value, vertices will change in output BlockMeshDict file.
For wedge boundary condition with 25⁰ , we will change only vertices and no other change is necessary.
No need to make any changes in ControlDict, pressure, velocity and physical property file for same boundary condition.
Similarly, we follow same procedure for wedge boundary condition with angle 45⁰
ParaView Model for Wedge Bondary condition :-
Case 2
Symmetry Boundary condition with angle 10⁰ :-
For symmetry boundary condition with 10⁰ angle, we will copy folder created for wedge boundary with 10⁰ angle. We will edit Matlab code and change boundary conditions for front face and back face to Symmetry as shown below
Clear all
Close all
clc
%% INPUTS
% Diameter of Pipe (m)
D = 0.02;
% Radius of Pipe (m)
R = D/2;
% Reynold's Number
Re = 2100;
% Wedge Angle (degrees)
theta = 10;
% Total Length of the Pipe (m)
L = 2.8;
% Assuming working fluid is water at 25° C
% Density of water (kg/m³)
rho = 997;
% Dynamic Viscosity of Water (Pa.s)
mu = 8.90 * 10^-4;
% Kinematic Viscosity of Water (m^2/s)
nu = mu/rho;
% Hydrodynamic Entrance Length (m)
L_h = 0.06 * Re * D;
%% HAGEN - POISEUILLE'S FLOW EQUATIONS
% Average Velocity (m/s)
V_avg = (mu * Re)/(rho * D);
% Maximum Velocity (m/s)
V_max = 2 * V_avg;
% Pressure Drop (Pa)
delta_p = (32*mu*V_avg*L)/(D^2);
% Kinematic Pressure Drop (Pa)
kinematic_delta_p = delta_p/rho;
%% Calculating shear stress
r = linspace(0,(D/2),500);
for i = 1:length(r)
tau(i) = 2 * mu * V_max * abs(r(i))/R^2;
end
figure(1)
plot(r,tau)
xlabel('Radius of pipe')
ylabel('shear stress')
title('Shear stress for Symmetry Boundary condition with angle 10?')
legend('shear stress')
%% DEFINING THE GEOMETRY
% 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 0 0];
V4 = [L R*cosd(theta/2) -R*sind(theta/2)];
V5 = [L R*cosd(theta/2) R*sind(theta/2)];
%% BLOCKMESHDICT FILE
% Creating the file
% Header
h1 = '/*--------------------------------*- C++ -*----------------------------------*\';
h2 = ' ========= |';
h3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4 = ' \\ / O peration | Website: https://openfoam.org';
h5 = ' \\ / A nd | Version: 9';
h6 = ' \\/ M anipulation |';
h7 = '\*---------------------------------------------------------------------------*/';
h8 = 'FoamFile';
h9 = '{';
h10 = ' format ascii;';
h11 = ' class dictionary;';
h12 = ' object blockMeshDict;';
h13 = '}';
h14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%% Creating BlockMesh File
f1 = fopen('blockMeshDictS10.txt','w');
% Header
fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'%s\n',h12);
fprintf(f1,'%s\n',h13);
fprintf(f1,'%s\n\n',h14);
% Conversion Factor
fprintf(f1,'%s\n\n','convertToMeters 1;');
% Vertices
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t(%d %d %d)\n',V0);
fprintf(f1,'\t(%d %f %f)\n',V1);
fprintf(f1,'\t(%d %f %f)\n',V2);
fprintf(f1,'\t(%.2f %d %d)\n',V3);
fprintf(f1,'\t(%.2f %f %f)\n',V4);
fprintf(f1,'\t(%.2f %f %f)\n',V5);
fprintf(f1,'%s\n\n',');');
% Blocks
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (400 30 1) simpleGrading (1 0.1 1)\n');
fprintf(f1,'%s\n\n',');');
% Edges
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\tarc 2 1 (%d %.3f %d)\n',0,R,0);
fprintf(f1,'\tarc 5 4 (%.3f %.3f %d)\n',L,R,0);
fprintf(f1,'%s\n\n',');');
% Boundary
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
% Axis
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Inlet
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 1 2 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Outlet
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Front Face
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype symmetry;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 4 1)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Back Face
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype symmetry;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 2 5 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
% Top Wall
fprintf(f1,'\tpipeWall\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(1 4 5 2)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');
BlockMeshDict file for Symmetry boundary condition with angle 10⁰:-
ControlDict file for Symmetry BC with 10⁰:-
In controlDict file we will kept End Time, Delta T and Write Interval values same as used for wedge boundary conditions. We will not make any changes.
U file (Velocity) and p file (Pressure) :-
We will make changes in velocity and pressure file located in
“0” folder à U file and P file
We will make changes in boundary field and use same names for boundary fields as used in BlockMeshDict file.
We also define type of boundary field for front and back faces type to “symmetry” from “wedge” as shown below.
No need to change anything else except front and back face type.
We will kept physicalProperties file same as used for previous case.
Symmetry Boundary condition with angle 25⁰ :-
For 25⁰ angle, we will copy entire folder created for symmetry boundary condition with 10⁰ angle including Matlab code.
We will make changes in Matlab code created for 10⁰ and change value of theta to 25 from 10, as shown below.
Clear all
Close all
clc
%% INPUTS
% Diameter of Pipe (m)
D = 0.02;
% Radius of Pipe (m)
R = D/2;
% Reynold's Number
Re = 2100;
% Wedge Angle (degrees)
theta = 25;
% Total Length of the Pipe (m)
L = 2.8;
We will not make any other changes in remaining code.
Due to changing in theta value, vertices will change in output BlockMeshDict file.
For symmetry boundary condition with 25⁰ , we will change only vertices and no other change is necessary.
No need to make any changes in ControlDict, pressure, velocity and physical property file for same boundary condition.
Similarly, we follow same procedure for symmetry boundary condition with angle 45⁰
ParaView Model for Symmetry Bondary condition :-
Velocity profile for Wedge BC along Y axis :- (Vmax)
Velocity profile for Symmetry BC along Y axis :- (Vmax)
Comparison of Analytical values and CFD simulation values :-
Average Velocity = 2 x Max. Velocity
Pressure drop in CFD simulation for wedge BC :-
`ΔP=(32⋅μ⋅V_(avg)⋅L)/D^2`
`ΔP=(32 xx 8.90 xx 10^(-4) xx 0.0928 xx 2.8)/(0.02)^2`
`ΔP=18.500 Pa`
|
Term |
Analytical Value |
Simulation Value |
|
Max Velocity |
0.1874 m/s |
0.1856 m/s |
|
Avg Velocity |
0.0937 m/s |
0.0928 m/s |
|
Hydro dynamic length |
2.52 m |
2.3 m |
|
Pressure Drop |
18.680 Pa |
18.5 Pa |
Comparison of Velocity profile for Wedge Boundary conditions at different angles :-
Comparison of Velocity profile for Symmetry Boundary conditions at different angles :-
Comparison of Shear stress distribution for Wedge Boundary conditions at different angles :-
Comparison of Shear stress distribution for Symmetry Boundary conditions at different angles :-
Usage and application of symmetry and wedge boundary condition
- Symmetry boundary condition is used for 3D objects under consideration, whereas the wedge boundary condition for axis- symmetric and 2D objects.
- symmetry is used in order to reduce the computation time and also the number of cell and complexity.
- The conditions for symmetric boundary are no flow across boundary , no scalar flux across boundary.
Conclusion:
- The Symmetry BC yields an accurate result for all angles whereas, the Wedge BC yields an accurate result pathway only for the lower order of wedge angles.
- It is observed that there isn't any large variation of Kinematic Pressure drop for Symmetry BC, whereas for Wedge BC it shows an oscillating value, and this time the value for pressure for the highest angle is not that accurate when compared to what we inferred in velocity part.