Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Aim: To simulate a validation study for Symmetry BC vs Wedge BC in OpenFOAM vs Analytical HP equation with angles to test as below.
Angles to test:
- `10^o`
- `25^o`
- `45^o`
Reference Data: The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. also to validate results with the Hagen- Poiseuille's equations for the fully developed flow. write code in Matlab to automate the generation of blockMeshDict files. The code should take inputs like pipe length, diameter, wedge angle, etc, and output a fully functional blockMeshDict file.
Data Given:
Wedge angle to be less than 5 degrees, Flow is axi-symetrics and incompressible Laminar flow.
Renolds number - 2100
Working fluid: Water.
Diameter of the pipe: 0.02 m.
Dynamic viscosity(μ
)=0.89x10^-3 Pas
Kinematic viscosity(ν
)=0.8927x10^-6 m^2/s
theta=10, 25, 45
Expected Results
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressure drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
- Highlight the differences in both boundary conditions in terms of their usage and application
Procedure
Average velocity
Vavg=ReμρD����=����� m/s
Maximum velocity
Vmax=2⋅Vavg�max=2⋅���� m/s
Pressure drop
ΔP=32μLVD2Δ�=32����2 pa
kinematic pressure drop=ΔPρΔ��
Openfoam Procedure
1. Creating Blockmesh dict file using Matlab code. Here pipe diameter is taken as 10mm and the length 1.2 meters.
2. Mesh is generated.
3. Inputs are added in OpenFoam's Velocity and Pressure files.
4. Run time and Time steps are added in the ControlDict file of OpenFoam.
5. Simulation is run as we have to take IcoFoam as a solver.
6. Results are viewed in ParaFoam Command, Paraview.
7. Plot results with Velocity for shear stress.
Programs for Theta 45
The code for Wedge Boundary Condition at a wedge angle of 45 deg is as follows:
clear all
close all
clc
%Inputs:
l= 1.5; %Length of the pipe.
theta = 45.0; %Wedge angle.
d = 0.01; %Diameter of the pipe.
c = 1.0; %Conversion Factor.
r = d/2; %Radius of the Cylinder.
nx = 500; %number of Mesh points in x-direction on the block.
ny= 25; %Number of Mesh points in the y-direction on the block.
nz= 1; %number of Mesh Points in the z-direction on the block.
gx = 1; %Grading along the x-direction.
gy = 1; %Grading along the y-direction.
gz = 1; %Grading along the z-direction.
f1 = fopen('BlockMeshDict.txt','w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
fprintf(f1,' ========= | n');
fprintf(f1,' / F ield | OpenFOAM: The Open Source CFD Toolbox n');
fprintf(f1,' / O peration | Website: https://openfoam.org n');
fprintf(f1,' / A nd | Version: 8 n');
fprintf(f1,' / M anipulation | n');
fprintf(f1,'/*------------------------------------------------------------------*/ n');
fprintf(f1, 'FoamFile n');
fprintf(f1,'{ n');
fprintf(f1,' version 2.0; n');
fprintf(f1,' format ascii; n' );
fprintf(f1,' class dictionary; n');
fprintf(f1,' object blockMeshDict; n');
fprintf(f1,'} n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // n');
fprintf(f1,'convertToMeters %0.1f; n',c);
fprintf(f1,'n');
fprintf(f1,'vertices n');
fprintf(f1,'( n');
fprintf(f1,'(0 0 0) n');
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1,'(0 0 %0.3f) n',-l);
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'blocks n');
fprintf(f1,'( n');
fprintf(f1,'hex (0 3 5 2 0 3 4 1) (%0.0f %0.0f %0.0f) simpleGrading (%0.0f %0.0f %0.0f) n',nx,ny,nz,gx,gy,gz);
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'edges n');
fprintf(f1,'( n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'boundary n');
fprintf(f1,'( n');
fprintf(f1,'in n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 1 2 0) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'out n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 5 4 3) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'ft n');
fprintf(f1,'{ n');
fprintf(f1,'type wedge; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 3 4 1) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'bk n');
fprintf(f1,'{ n');
fprintf(f1,'type wedge; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 0 2 5) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'top n');
fprintf(f1,'{ n');
fprintf(f1,'type wall; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(5 2 1 4) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'axis n');
fprintf(f1,'{ n');
fprintf(f1,'type empty; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 3 3 0) n');
fprintf(f1,'); n');
fprintf(f1, '} n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'mergePatchPairs n');
fprintf(f1,'( n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'// ************************************************************************* // n');
%END
----------------------------------------------------------------------------------------------------------------
The code for Symmetry Boundary Condition at a wedge angle of 45 deg is as follows:
clear all
close all
clc
%Inputs:
l= 1.5; %Length of the pipe.
theta = 5.0; %Wedge angle.
d = 0.01; %Diameter of the pipe.
c = 1.0; %Conversion Factor.
r = d/2; %Radius of the Cylinder.
nx = 500; %number of Mesh points in x-direction on the block.
ny= 25; %Number of Mesh points in the y-direction on the block.
nz= 3; %number of Mesh Points in the z-direction on the block.
gx = 1; %Grading along the x-direction.
gy = 1; %Grading along the y-direction.
gz = 1; %Grading along the z-direction.
f1 = fopen('BlockMeshDict.txt','w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
fprintf(f1,' ========= | n');
fprintf(f1,' / F ield | OpenFOAM: The Open Source CFD Toolbox n');
fprintf(f1,' / O peration | Website: https://openfoam.org n');
fprintf(f1,' / A nd | Version: 8 n');
fprintf(f1,' / M anipulation | n');
fprintf(f1,'/*------------------------------------------------------------------*/ n');
fprintf(f1, 'FoamFile n');
fprintf(f1,'{ n');
fprintf(f1,' version 2.0; n');
fprintf(f1,' format ascii; n' );
fprintf(f1,' class dictionary; n');
fprintf(f1,' object blockMeshDict; n');
fprintf(f1,'} n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // n');
fprintf(f1,'convertToMeters %0.1f; n',c);
fprintf(f1,'n');
fprintf(f1,'vertices n');
fprintf(f1,'( n');
fprintf(f1,'(0 0 0) n');
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1,'(0 0 %0.3f) n',-l);
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'blocks n');
fprintf(f1,'( n');
fprintf(f1,'hex (0 3 5 2 0 3 4 1) (%0.0f %0.0f %0.0f) simpleGrading (%0.0f %0.0f %0.0f) n',nx,ny,nz,gx,gy,gz);
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'edges n');
fprintf(f1,'( n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'boundary n');
fprintf(f1,'( n');
fprintf(f1,'in n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 1 2 0) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'out n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 5 4 3) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'ft n');
fprintf(f1,'{ n');
fprintf(f1,'type symmetry; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 3 4 1) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'bk n');
fprintf(f1,'{ n');
fprintf(f1,'type symmetry; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 0 2 5) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'top n');
fprintf(f1,'{ n');
fprintf(f1,'type wall; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(5 2 1 4) n');
fprintf(f1,'); n');
fprintf(f1,'} n');
fprintf(f1,'axis n');
fprintf(f1,'{ n');
fprintf(f1,'type empty; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 3 3 0) n');
fprintf(f1,'); n');
fprintf(f1, '} n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'mergePatchPairs n');
fprintf(f1,'( n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'// ************************************************************************* // n');
%END
-------------------------------------------------------------------------------------------
Similar Theta angles for 10 & 25 Deg has to be progammed.
controlDict file
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
Pressure File
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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;
}
wedge_front
{
type wedge;
}
wedge_back
{
type wedge;
}
top
{
type zeroGradient;
}
axis
{
type empty;
}
}
Velocity file
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2006 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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.09373 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
wedge_front
{
type wedge;
}
wedge_back
{
type wedge;
}
axis
{
type empty;
}
}
Results:
Shear Stress Variation with respect to the Radial Distance for 10 Deg
Shear Stress Variation with respect to the Radial Distance for 25 Deg
Shear Stress Variation with respect to the Radial Distance for 45 Deg
Velocity Variation with respect to Radial Distance for 10 deg.
Velocity Variation with respect to Radial Distance for 25 deg.
Velocity Variation with respect to Radial Distance for 45 deg.
Pressure Drop across the axial Distance for 10 deg.
Pressure Drop across the axial Distance for 25 deg.
Pressure Drop across the axial Distance for 45 deg.
Validation
Pressure drop can be calculated using formula as below.
△P=Vavg⋅μ⋅LD2△�=����⋅�⋅��2
Conclusion:
The Courant number is increasing slightly in the case of wedge Boundary Condition when compared with symetry BC. hence solving time increased.
other then that both values are same and by this assesment, we can conclude as changing BC to symmetry or wedge, the results remain same.