Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Symmetry Boundary Condition:
Generally we will encounter an analysis where you need more computer memory and/or CPU speed that you can obtain. Therefore you should get familiar with methods for using efficient modeling techniques. One method of efficiently using finite volume method(FVM) modeling is to exploit the planes of symmetry method. When we take advantage of symmetry, we need to model only a portion of the actual structure to reduce the analysis run time and memory required .The lines or planes of symmetry in a FVM model can be simulated by providing proper restraints to the symmetrical faces or edges. To model symmetry, the geometry and the restraints must be symmetric about a plane.
Wedge Boundary Condition:
The wedge boundary condition is used to define an axisymmetric situation, for eg. 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.
Here wedge angle(`theta`) was restricted to `5^o` only.In case if wedge angle(`theta`) was assumed greater than `5^o`(`theta>5^o`),then incorrect solutions will occur.
In such conditions,symmetric boundary conditions will be useful because symmetric boundary doesn't have angle restrictions.
CONCEPT:
- write Matlab code to generate automatic blockMeshDict file for different angles and also to calculate theoritical values of Hagen-poiseuille's equation.
- Check once created geometry and mesh condition using "blockMesh" command.
- Update initial condition and constant properties
- Update controlDict file
- Run icoFoam command
- Run paraFoam command to navigate to Paraview.
(Note:Hence this challenge is a continuation of the week 11 challenge.,All initial and constant properties are same)
MATLAB CODE:
for `theta=10`
clear all
close all
clc
%constant properties
Re=2100; %Reynolds number
D=0.02; %Diameter of pipe
R=D/2; %Radius of pipe
L=3; %length of pipe
mu=1.002*10^(-3); %Dynamic viscosity
rho = 998; %Density
theta=10; %angle of wedge
nu=mu/rho; %Kinematic viscosity
% Mesh details
nx=10;ny=1;nz=100;%number of grid points along respective axis
gx=1;gy=1;gz=1; % grading factor along respective axis
%analytical solution
L_e=0.05*Re*D
vz_avg=(Re*mu)/(rho*D) % average velocity
vz_max=2*vz_avg %maximum velocity
dp= (-8*mu*vz_avg/(R^2))*L;
r=linspace(0,R,nx);
vz=vz_max*(1-(r.^2./R^2))
T_wall=4*mu*vz_avg/R
plot(r,vz,'linewidth',2)
xlabel('radius')
ylabel('velocity')
title('velocity profile at fully developed region')
%vertices
v0=[0 0 0];
v1=[R*cosd((theta/2)) -R*sind((theta/2)) 0];
v2=[R*cosd((theta/2)) R*sind((theta/2)) 0];
v3=[0 0 L];
v4=[R*cosd((theta/2)) -R*sind((theta/2)) L];
v5=[R*cosd((theta/2)) R*sind((theta/2)) L];
f1=fopen('blockMeshDict','w');
fprintf(f1,'%s\n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'%s\n',' ========= |');
fprintf(f1,'%s\n',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'%s\n',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'%s\n',' \\ / A nd | Version: 8');
fprintf(f1,'%s\n',' \\/ M anipulation |');
fprintf(f1,'%s\n','\*---------------------------------------------------------------------------*/');
fprintf(f1,'%s\n','FoamFile');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n',' version 2.0;');
fprintf(f1,'%s\n','format ascii;');
fprintf(f1,'%s\n','class dictionary;');
fprintf(f1,'%s\n','object blockMeshDict;');
fprintf(f1,'%s\n','}');
fprintf(f1,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%s \n\n','convertToMeters 1 ;');
fprintf(f1,'%s \n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,' (%d %d %d) \n', [v0(1) v0(2) v0(3)]);
fprintf(f1,' (%d %d %d) \n', [v1(1) v1(2) v1(3)]);
fprintf(f1,' (%d %d %d) \n', [v2(1) v2(2) v2(3)]);
fprintf(f1,' (%d %d %d) \n', [v3(1) v3(2) v3(3)]);
fprintf(f1,' (%d %d %d) \n', [v4(1) v4(2) v4(3)]);
fprintf(f1,' (%d %d %d) \n', [v5(1) v5(2) v5(3)]);
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s',' hex (0 1 2 0 3 4 5 3)');
fprintf(f1,' (%d %d %d)',[nx ny nz] );
fprintf(f1,' %s (%d %d %d)','simpleGrading', [gx gy gz]);
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',' inlet');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type patch;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 1 2 0)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' outlet');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type patch;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (3 5 4 3)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' wedgeside1');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type symmetry;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 3 4 1)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' wedgeside2');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type symmetry;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 2 5 3)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' pipewall');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type wall;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (2 1 4 5)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',' axis');
fprintf(f1,'%s\n',' {');
fprintf(f1,'%s\n',' type empty;');
fprintf(f1,'%s\n',' faces');
fprintf(f1,'%s\n',' (');
fprintf(f1,'%s\n',' (0 0 3 3)');
fprintf(f1,'%s\n',' );');
fprintf(f1,'%s\n',' }');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');
fclose(f1);Above code results,below blockMeshDict file:
/*--------------------------------*- 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)
(9.961947e-03 -8.715574e-04 0)
(9.961947e-03 8.715574e-04 0)
(0 0 3)
(9.961947e-03 -8.715574e-04 3)
(9.961947e-03 8.715574e-04 3)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (10 1 100) simpleGrading (1 1 1));
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wedgeside1
{
type symmetry;
faces
(
(0 3 4 1)
);
}
wedgeside2
{
type symmetry;
faces
(
(0 2 5 3)
);
}
pipewall
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 0 3 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //INITIAL CONDITION:
Update Initial 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 0 0.11);
}
outlet
{
type zeroGradient;
}
wedgeside1
{
type symmetry;
}
wedgeside2
{
type symmetry;
}
pipewall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //Update Initial 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 1;
}
wedgeside1
{
type symmetry;
}
wedgeside2
{
type symmetry;
}
pipewall
{
type zeroGradient;
}
axis
{
type 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 30;
deltaT 0.005;
writeControl timeStep;
writeInterval 30;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Results:
Analytical Solution:
Analytical solution was from Matlab code by substituting values in formulae.
Numerical Solution:
1. For `theta=10^o`
As we move way from centerline,its velocity was decreasing and its shear stress was increasing.
After length of pipe=2.4,its velocity profile was unchanged.Hence it is its Hydrodynamic length
2.For `theta=25`
As we move way from centerline,its velocity was decreasing and its shear stress was increasing.
After length of pipe=2.4,its velocity profile was unchanged.Hence it is its Hydrodynamic length
Here only vertices are updated for `theta` value.
3. For `theta=45`
As we move way from centerline,its velocity was decreasing and its shear stress was increasing.
After length of pipe=2.4,its velocity profile was unchanged.Hence it is its Hydrodynamic length
Here only vertices are updated for `theta` value.
Validation:
| parameters | Analytical solution | `theta=10` | `theta=25` | `theta=45` |
| Time taken | 43 | 40 | 35 | |
| Hydro-dynamic length | 2.1 | 2.4 | 2.4 | 2.4 |
| velocity profile | Parabolic | Parabolic | Parabolic | Parabolic |
| maximum velocity | 0.21 | 0.218 | 0.217 | 0.2175 |
| wall shear stress | 0.0423 | 0.03 | 0.032 | 0.035 |
Difference Between Wedge and Symmetric Boundary Condition:
- A wedge is used for axial and spherical symmetry and needs to perform a coordinate transformation for the wedge angle and not a simple mirroring.
- Symmetry was used only when all physical parameters are mirror to plane.
- Wedge geometry was constrained to angle(`theta<5^o`).So such suitations symmetry can be used,because symmetry boundary doesn't have any angle restriction.
- wedge boundary condition is applicable only for Wedgeshape geometry.But symmetric can used of any geometry.