Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Aim: To obtain results of the simulation of flow through circular pipes by applying symmetry and wedge boundary condition and validation of the results with the analytical solution.
Objectives:
1. Validation of hydrodynamic length with the numerical result.
2. Validation of the obtained velocity profile.
3. Validation of maximum velocity and pressure drop with the numerical result.
4. To validate wall shear stress for fully developed flow.
5. Highlight the differences in both boundary condition in terms of their usage and application.
Introduction:
There are various types of boundary conditions available in OpenFoam but as per the problem specificity the boundary conditions are filtered. In this project we will be using "symmetry" and "wedge" type of boundary conditions to obtain their significance on our flow problem.
Wedge boundary condition:
- The wedge boundary condition in CFD is used to define an axisymmetric situation for e.g: cylinder.
- The model and flow must be axisymmetric along a center line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.
- It works for 2-dimensions only.
- One big advantage of "wedge" boundary condition 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 to as zero.
Symmetry boundary condition:
- The symmetry boundary condition defines a mirror face/surface.
- It should only be used if the physical object or geomtery and the expected flow field pattern of the developed solution are mirrored along that surface.
- By using this boundary condition, the domain can essentially be halved, reducing the time to achieve a solution.
- The symmetry condition applies the following constraints on the flow variables,
a) The fluxes acrosss the symmetry are zero.
b) The normal components of all variables are set to zero.
Project methodology:
In this project, we will simulate the flow through circular pipes by applying symmetry and wedge boundary condition. Moreover, we will test various wedge angles to study the influence on the flow field variables. This project is purely meant to study the effect of boundary condition on the given problem. Initially, the blockMeshDict will automated through a MATLAB script. Then the resulted blockMeshDict will imported to OPENFOAM, with appropriate initial and boundary condition. Each wedge angle (`10^0, 25^0 and 45^0`) with respective boundary condition will simulated. The results will be partially post-processed in PARAVIEW and remaining validation will be executed in MATLAB.
As per the previous project, a MATLAB script is used to automate the geometry and mesh characteristics of the given problem. For various wedge angles input prompt is created, which will automate the blockMeshDict as per the user input.
Following is the MATLAB script for blockMeshDict automation for symmetry type boundary condition.
MATLAB Code:
clear all
close all
clc
%% Variables
% Reynold's number
R = 2100;
% wedge angle (In degree)
theta = input('Enter the wedge angle = ');
% Diameter of the pipe (m)
D = 0.004;
% Radius of the pipe (m)
rad = D/2;
% Hydrodynamic length (m)
L_h = 0.05*R*D;
% Total length of the pipe (m)
L = 0.28 + L_h;
% Calculating vertices
v_0 = [0 0 0];
vr_1 = [0 rad*cosd(theta/2) -rad*sind(theta/2)];
v_1 = [vr_1(1) vr_1(2) vr_1(3)];
vr_2 = [0 rad*cosd(theta/2) rad*sind(theta/2)];
v_2 = [vr_2(1) vr_2(2) vr_2(3)];
vr_3 = [L 0 0];
v_3 = [vr_3(1) vr_3(2) vr_3(3)];
vr_4 = [L rad*cosd(theta/2) -rad*sind(theta/2)];
v_4 = [vr_4(1) vr_4(2) vr_4(3)];
vr_5 = [L rad*cosd(theta/2) rad*sind(theta/2)];
v_5 = [vr_5(1) vr_5(2) vr_5(3)];
%% File writing
line1 = '/*--------------------------------*- C++ -*----------------------------------*\';
line2 = ' ========= |';
line3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4 = ' \\ / O peration | Website: https://openfoam.org';
line5 = ' \\ / A nd | Version: 8';
line6 = ' \\/ M anipulation |';
line7 = '\*---------------------------------------------------------------------------*/';
% File open
block_mesh = fopen('blockMeshDict_10.txt','w');
% To write header
fprintf(block_mesh,'%s\n',line1);
fprintf(block_mesh,'%s\n',line2);
fprintf(block_mesh,'%s\n',line3);
fprintf(block_mesh,'%s\n',line4);
fprintf(block_mesh,'%s\n',line5);
fprintf(block_mesh,'%s\n',line6);
fprintf(block_mesh,'%s\n',line7);
% To write the dictionary
fprintf(block_mesh,'%s\n','FoamFile');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s\n',' version 2.0;');
fprintf(block_mesh,'%s\n',' format ascii;');
fprintf(block_mesh,'%s\n',' class dictionary;');
fprintf(block_mesh,'%s\n',' object blockMeshDict;');
fprintf(block_mesh,'%s\n','}');
fprintf(block_mesh,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(block_mesh,'\n');
% Note: As all the calculated variables are in meter(m), keeping scale to
% be as one (1).
scale = 1; % SCALE input
fprintf(block_mesh,'%s%.4f%s\n','convertToMeters ',scale,';');
fprintf(block_mesh,'\n');
% Vertices
fprintf(block_mesh,'%s\n','vertices');
fprintf(block_mesh,'%s','(');
fprintf(block_mesh,'\n');
fprintf(block_mesh,'(%d %d %d)\n',v_0);
fprintf(block_mesh,'(%d %d %d)\n',v_1);
fprintf(block_mesh,'(%d %d %d)\n',v_2);
fprintf(block_mesh,'(%d %d %d)\n',v_3);
fprintf(block_mesh,'(%d %d %d)\n',v_4);
fprintf(block_mesh,'(%d %d %d)\n',v_5);
fprintf(block_mesh,'%s\n\n',');');
% Block parameters
fprintf(block_mesh,'\n');
mesh = 'hex';
vertice = [0 1 2 3 4 5];
vertice_seq = [0 1 2 0 3 4 5 3];
nx = 30; % Cells along x
ny = 1; % Cells along y
nz = 200; % Cells along z
xgf = 0.2; % Grading along x
ygf = 1; % Grading along y
zgf = 1; % Grading along z
arc_1 = [1 2]; % Arc 1 input
arc_2 = [4 5]; % Arc 2 input
% Origin and extension along x
x0 = 0;
y0 = 0;
z0 = 0;
x1 = L;
% Write block
fprintf(block_mesh,'%s\n','blocks');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s (%d %d %d %d %d %d %d %d) (%d %d %d) %s (%d %d %d)',mesh,vertice_seq,nx,ny,nz,'simpleGrading',xgf,ygf,zgf);
fprintf(block_mesh,'\n');
fprintf(block_mesh,'%s\n\n',');');
% Edges
fprintf(block_mesh,'%s\n','edges');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s %d %d (%d %d %d)\n','arc',arc_1(1),arc_1(2),x0,rad,z0);
fprintf(block_mesh,'%s %d %d (%d %d %d)\n','arc',arc_2(1),arc_2(2),x1,rad,z0);
fprintf(block_mesh,'%s\n\n',');');
% Boundary
fprintf(block_mesh,'%s\n','boundary');
boundary_type = {'patch','empty','wall','symmtery','wedge'}; % Types of wall boundary condition
% Inlet
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s\n','inlet');
fprintf(block_mesh,'%s\n','{')
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{1},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(3),vertice(2),vertice(1));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Outlet
fprintf(block_mesh,'%s\n','outlet');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{1},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(4),vertice(5),vertice(6),vertice(4));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Top wall
fprintf(block_mesh,'%s\n','topWall');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{3},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(2),vertice(3),vertice(6),vertice(5));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Bottom wall
fprintf(block_mesh,'%s\n','axis');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{2},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(4),vertice(4),vertice(1));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Front wedge
fprintf(block_mesh,'%s\n','frontWedge');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{4},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(4),vertice(6),vertice(3));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Back wedge
fprintf(block_mesh,'%s\n','backWedge');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{4},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(2),vertice(5),vertice(4));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'\n');
fprintf(block_mesh,'%s\n','mergePatchPairs');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s','// ************************************************************************* //');
fclose(block_mesh);Note: The above MATLAB code with different wedge angles results in generation of blockMeshDict for different cases. So separate MATLAB code is not attached for all the three cases (`10^0, 25^0 and 45^0`).The resulted blockMeshDict for all the three cases with symmetry boundary condition is attached below,
blockMeshDict (Wedge angle: `10^0`) (BC : symmetry):
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.992389e-03 -1.743115e-04)
(0 1.992389e-03 1.743115e-04)
(7.000000e-01 0 0)
(7.000000e-01 1.992389e-03 -1.743115e-04)
(7.000000e-01 1.992389e-03 1.743115e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topWall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type symmetry;
faces
(
(0 3 5 2)
);
}
backWedge
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //blockMeshDict (Wedge angle: `25^0`) (BC : symmetry):
/*--------------------------------*- 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.0000; vertices ( (0 0 0) (0 1.952592e-03 -4.328792e-04) (0 1.952592e-03 4.328792e-04) (7.000000e-01 0 0) (7.000000e-01 1.952592e-03 -4.328792e-04) (7.000000e-01 1.952592e-03 4.328792e-04) ); blocks ( hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1) ); edges ( arc 1 2 (0 2.000000e-03 0) arc 4 5 (7.000000e-01 2.000000e-03 0) ); boundary ( inlet { type patch; faces ( (0 2 1 0) ); } outlet { type patch; faces ( (3 4 5 3) ); } topWall { type wall; faces ( (1 2 5 4) ); } axis { type empty; faces ( (0 3 3 0) ); } frontWedge { type symmetry; faces ( (0 3 5 2) ); } backWedge { type symmetry; faces ( (0 1 4 3) ); } ); mergePatchPairs ( ); // ****************************** ****************************** ************* //
blockMeshDict (Wedge angle: `45^0`) (BC : symmetry):
*----------------------------Assigning Initial and boundary conditions:
The initial and boundary conditions are assigned to the given problem, as velocity inlet and pressure outlet same as previous study for wedge boundary condition study. To have a common datum for comparison no significant changes are made to the given parameters. The initial and boundary conditions will be same for all three cases (`10^0, 25^0 and 45^0`) with symmetry boundary condition at the front and back wedge.
Pressure:
/*--------------------------------*- 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 { outlet { type fixedValue; value uniform 0; } inlet { type zeroGradient; } frontWedge { type symmetry; } backWedge { type symmetry; } topWall { type zeroGradient; } axis { type empty; } } // ****************************** ****************************** ************* //
Velocity:
/*--------------------------------*- 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.5271 0 0); } outlet { type zeroGradient; } frontWedge { type symmetry; } backWedge { type symmetry; } topWall { type noSlip; } axis { type empty; } } // ****************************** ****************************** ************* //
Time control:
The time control characteristics are kept same to see the variation in CFL number for different cases.
controlDict:
/*----------------------------With this input parameters, case is executed with the solver "icoFoam", which deals with incompressible and laminar flows. Finally, after executing the case, we post-process the results in PARAVIEW.
A] Post-Processing in paraView:
Mesh Visualization:
a) Wedge angle (`10^0`)
b) Wedge angle (`25^0`)
c) Wedge angle (`45^0`)
1] Shear stress:
a) (Wedge angle: `10^0`) (BC : symmetry)
Shear stress:
b) (Wedge angle: `25^0`) (BC : symmetry)
Shear stress:
c) (Wedge angle: `45^0`) (BC : symmetry)
Discussion:
As previously discussed, the shear stress as a function of radial distance will increase with the radial distance. The above plots resulted with different wedge angles does not show any disparity. However, we can observe that, at the hydrodynamic entrance length and in fully developed region (at x = 0.65m), there is negligible shear stress variation. This is because of the uniform flow at the end of hydrodynamic entrance length and to the fully developed region. So here velocity is varying with constant velocity gradient along the radial distance.
2] Velocity:
a) (Wedge angle: `10^0`) (BC : symmetry)
b) (Wedge angle: `25^0`) (BC : symmetry)
c) (Wedge angle: `45^0`) (BC : symmetry)
Discussion:
The velocity profile is plotted for three different wedge angles at three different locations along the pipe length. With this plots, we can observe that the velocity profile at x = 0.42m matches the velocity profile at x = 0.65m with very negligible differences. At the hydrodynamic length, we obtain uniform flow with developed boundary layer resulting in parabolic velocity profile. So is the numerical result at x = 0.42m. As per our analytical hydrodynamic length (`L_h = 0.42m`), the numerical results are closely in agreement with the analytical solution.
3] Kinematic pressure:
a) (Wedge angle: `10^0`) (BC : symmetry)
b) (Wedge angle: `25^0`) (BC : symmetry)
c) (Wedge angle: `45^0`) (BC : symmetry)
Discussion:
The resulted kinematic pressure along the pipe length (x) doesn't show any variation at the macroscopic level. The pressure is non-linear at the hydrodynamic entrance region but from x = 0.4m (start of fully developed region) the pressure variation is linear with constant pressure gradient.
By comparing the plots resulted from these studies does not give a clear insight about the variation of the flow variables. However, we need to post-process the results with their numerical value and need to calculate associated errors. This will help to analyze data in more meaningful manner.
B] Post-Processing in MATLAB:
The purpose of MATLAB post-processing is to validate the kinematic pressure drop and maximum velocity with the analytical solution.
MATLAB Code:
clear all
close all
clc
%% Kinematic pressure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Calculate kinematic pressure gardient %%%%%%%%%%%%%%%%%%%%%%%%%%%5
% Reynold's number
R = 2100;
% Diameter of the pipe (m)
D = 0.004;
% Radius of the pipe (m)
rad = D/2;
% Hydrodynamic length (m)
L_h = 0.05*R*D;
% Total length of the pipe (m)
L = 0.28 + L_h;
% Length for fully developed flow (characteristic length)
l = L - L_h;
% Total number of time steps
n = 1250;
% Characteristic length (m)
charac_l = 0.18;
% Density of water (kg/m^3)
rho = 1000;
% Kinematic Viscosity (N.s/m^2)
k_v = 1.004e-6;
% Dynamic viscosity (N.s/m^2)
mu = 1.004e-3;
%To calculate average velocity (m/s)
v_avg_an = R*k_v/(D);
% To calculate maximum velocity
v_max_an = 2*v_avg_an;
% Analytical pressure drop calc using HAGEN-POISEUILLES EQN
press_drop = (8*mu*v_avg_an*charac_l)/(rad^2);
% Analytical kinematic pressure
kin_p_analytical = press_drop/rho;
%% To validate kinematic pressure drop in fully developed region %%%%%
a = csvread('kinep_xsymm_10.csv',1,1);
del_p_10 = abs(a(673)-a(930));
b = csvread('kinep_xsymm_25.csv',1,1);
del_p_25 = abs(b(673)-a(930));
c = csvread('kinep_xsymm_45.csv',1,1);
del_p_45 = abs(c(673)-c(930));
d = csvread('pressure.csv',1,1);
del_p_wedge = abs(d(1)-d(1001));
% Errors
error_kp_10 = abs(kin_p_analytical - del_p_10);
error_kp_25 = abs(kin_p_analytical - del_p_25);
error_kp_45 = abs(kin_p_analytical - del_p_45);
error_kp_wedge = abs(kin_p_analytical - del_p_wedge);
%% To validate maximum velocity
x = csvread('vel_radsymm_10.csv',1,1);
vel_maxnum_10 = max(x(:,1)); % First column belongs to fully developed region
y = csvread('vel_radsymm_25.csv',1,1);
vel_maxnum_25 = max(y(:,1));
z = csvread('vel_radsymm_45.csv',1,1);
vel_maxnum_45 = max(z(:,1));
% Errors
error_max_vel_10 = abs(v_max_an - vel_maxnum_10);
error_max_vel_25 = abs(v_max_an - vel_maxnum_25);
error_max_vel_45 = abs(v_max_an - vel_maxnum_45);Observation table:
| Boundary conditions | Execution time (s) | CFL | Kinematic pressure drop (`m^2/s`) | Maximum velocity (`m/s`) | Error when compared with analytical solution (Kinematic pressure drop) | Error when compared with analytical solution (Maximum velocity) |
| Wedge (`4^0`) | 50.73 | 4.225.1 | 0.1916 | 1.0482 | 0.0011 | 0.006 |
| Symmetric (`10^0`) | 51.83 | 4.23805 | 0.1927 | 1.0483 | 0.0022 | 0.0059 |
| Symmetric (`25^0`) | 56.03 | 4.28229 | 0.2025 | 1.0488 | 0.0120 | 0.0054 |
| Symmetric (`45^0`) | 55.49 | 4.38226 | 0.2231 | 1.0497 | 0.0325 | 0.0045 |
| Analytical HP solution | NA | NA | 0.1905 | 1.0542 | NA | NA |
In this observation table, we can see the magnitude of flow variables and the numerical errors associated with the boundary conditions. The wedge type boundary condition is way better to execute the results in less computation time as compared to symmetry boundary condition. The resulting CFL of wedge boundary condition is also less because of the effective radial distance. While automating blockMeshDict, we observed that more is the wedge angle, less is the effective radial distance (i.e. `rcostheta`). Because of this less radial distance and grading, the CFL number of wedge angle `45^0` is more and it is decreasing as wedge angle is reduced. This is the case with our numerical set-up, however CFL number will vary with the other characteristics too.
In addition to the results obtained, we can see that wedge boundary condition tends to be efficient in calculating the kinematic pressure drop as compared to symmetry. It has less numerical errors and the errors are increasing with the wedge angle in symmetry boundary condition. In contrast to it, the symmetry boundary condition with `45^0` approximates the maximum velocity of the fluid more accurately as compared to others. Certainly it possess less numerical errors while calculating maximum velocity.
Conclusion:
The application study of symmetry and wedge boundary condition was successfully conducted while keeping the same input parameters to build a common datum. The results obtained are interesting and unbiased to a specific boundary condition. The numerical hydrodynamic length is successfully validated with analytical solution with very negligible differences. Both the wedge and symmetry boundary condition gave similar results for the case of hydrodynamic length. However, while validating kinematic pressure drop and maximum velocity of the flow, no such one boundary condition gave accurate results in both the cases. The wedge boundary condition captured kinematic pressure drop more accurately, whereas the maximum velocity was more approximate in case of symmetry boundary condition with `45^0` wedge angle. Either any of the boundary condition defined superiority over other in all aspects. As the significance of boundary condition is problem specific and should be carefully opt to capture the concerned flow variables.
References:
2. https://www.simscale.com/
3. https://www.simscale.com/