This section presents the discrete numerical implementation of the boundary conditions employed in this work, as previously listed in Section2.5.2.Figure3-7gives a drawing of a near boundary cell (i.e. cell with boundary face).
P
f
Boundary face(B)
Figure 3-7: Drawing of a near boundary cell
The near boundary cell and boundary face will be given the subscripts (P) and (B) respectively. The surface normal vector (S~f) is taken such that it points outside the computational domain and d~P f is the vector from the cell-center to the center of the boundary face. The vectord~nis the vector from the cell-center normal to the boundary face. The detailed numerical implementation is as follows:
3.6.1 Inlet Boundary
At the inlet boundary, all variables are known which implies a Dirichlet boundary except the pressure which is unknown and thus treated as Neumann boundary.
43 3.6. BOUNDARY CONDITIONS Momentum equation: The velocity is prescribed asV~BoruiBwhich means that
Diffusion term: Following the terminology in Figure3-7 S~f·(∂ui
∂xj)f =|Sf| uiB−uiP
|dP f|
| {z }
OrthognalContribution
+ k~f·(∂ui
∂xj)B
| {z }
Non−orthonalContribution
(3.46)
wherek~f = |SS~f
f|−|dd~P f
P f|, and the velocity gradient is obtained by extrapolation from the domain interior.
Convection term: The convection term is computed based on the first order upwind only, thus simplifying its implementation.
The volumetric flux (F) is obtained based on the specified velocity(V~B) at the bound-ary FB=Sf·V~B. The convection boundary contribution is obtained by substituting in Equation (3.19) resulting in:
ρfFf~Vf =ρBFB~VB (3.47) whereρf is the density value at the boundary face computed based on the volume frac-tion (αB) , which must be specified on the boundary, using Equation (3.20). The con-vection contribution is constant and thus it is explicitly.
Pressure equation: The pressure gradient is set to zero and its contribution vanishes from the pressure equation Equation (3.42)
∇pm+1 ap
B·S~f =0 (3.48)
Turbulence: The distribution of the turbulence field variables is prescribed and its values are typically computed based on a specific turbulence intensity and length scale ( or turbulence viscosity ratio ) using Equations (2.29) to (2.32).
CHAPTER 3. NUMERICAL METHOD
3.6.2 Open Boundary
Open boundary is treated as either an outlet boundary or inlet boundary based on the sign of the volumetric flux computed such that mass conservation is satisfied. If volumetric flux is negative, the boundary is treated as an inlet following the previously outlined procedure. If volumetric flux is positive, the boundary is treated as an outlet as follows
Momentum equation: The velocity is computed from the domain interior as
V~B=V~P (3.49)
Diffusion term: Substituting Equation (3.49) in Equation (3.46) results in the van-ishing of the orthogonal contribution in the term and only the non-orthogonal contribu-tion remains which is often neglected.
Convection term: When first order upwind is used for treating this term, its contri-bution vanishes.
Pressure equation: The pressure value is prescribed (PB) and it is treated in a similar manner to Equation (3.46) as follows
∇pm+1 ap
B·S~f = 1 ap
B|Sf| pB−pP
|dP f|
| {z }
OrthognalContribution
+ ~kf·(∂p
∂xj)B
| {z }
Non−orthonalContribution
m+1
(3.50)
where 1
ap
Bis computed by extrapolation from the domain interior.
Turbulence: Zero normal gradient is applied for all turbulence variables (k,ε,ω, and µt ).
∂k
∂n =0 ∂ ε
∂n =0 ∂ ω
∂n =0 ∂ µt
∂n =0 (3.51)
45 3.6. BOUNDARY CONDITIONS
3.6.3 Rigid Wall Boundary
Wall boundary can either be stationary or it could have specified velocity which is used here for generality.
Momentum equation: The velocity is prescribed asV~BoruiBwhich means that Diffusion term: Equation (3.46) can be used.
Convection term: Owing to no penetration condition, the volumetric flux is set to zero thus this term vanishes.
Pressure equation: The pressure gradient is set to zero and its contribution vanishes from the pressure equation Equation (3.42)
Turbulence: In this work, walls are assumed smooth. For RANS models, the law of the wall[47] is used to model the field in the cells adjacent to the wall(share a boundary face with the wall).
Firstly, the turbulence production term (Pk) in Equations (2.10), (2.11), (2.14), (2.15), (2.25) and (2.26) is modified as
Pk=|τw| uτ
κ|dn| (3.52)
where|dn|is the normal distance to the wall computed from the wall which is compute from
~dn=|d~P f·S~f S~f
S~f·S~f| (3.53)
τwis the wall shear stress
τw=µ|V~P−V~B|
|dn| (3.54)
CHAPTER 3. NUMERICAL METHOD
uτ is the frictional velocity computed from uτ =C1/4µ p
kP (3.55)
andκ=0.41 andCµ =0.09
The equations for the turbulence dissipation (ε) and specific turbulence dissipation (ω) are not solved for the near boundary cell (P) and its value is set as
εP=Cµ3/4kP3/2
κ|dn| (3.56)
ωP=
√kP
κCµ1/4|dn| (3.57)
The turbulence viscosity is computed by starting from the definition of the wall shear stress (τw) in Equation (3.54) and using the definition of the normalized distance to the wally+= ρ|dn|uτ
µ with uτ computed from Equation (3.55) and the normalized velocity parallel to the wall u+ = |V~P−uV~B|k
τ where |V~P−V~B|k is the magnitude of the velocity parallel to the wall.
µt =µ y+
u+ −1
(3.58) Depending on the y+ value (which positions the cell inside the turbulent boundary layer), the value ofu+ is set as
u+ =
y+ ify+≤y+cut off log(E+y+) ify+>y+cut off
(3.59)
whereE+ is a constant approximately taken around 10 . Different values fory+cut off are reported in various source but all them lie between 11 and 12. In this work,y+cut off=11 is employed throughout the numerical calculation. This method yields good results
47 3.6. BOUNDARY CONDITIONS when used with a mesh having 30<y+ <300. Depending on the required level of accuracy, the physics of the simulation, the Reynolds number and mesh to be used, other types of wall functions can be used[62].
In the case of LES models, a zero gradient boundary condition is often applied in
10−2 10−1 100 101 102
0 5 10 15
20 Viscous Buffer Turbulent
y+
u+
Figure 3-8: Relation betweenu+ andy+ in the turbulent boundary layer
conjunction with an analytical viscosity decay function such as van Driest[23] but it was not implemented in this work.
3.6.4 Shell Boundary
A shell boundary is an interior face having cells on both its sides as shown in Figure3-9 where it is treated as wall boundary to each side. Consequently, its treatment is the same as the previously outlined wall boundary for each side.
CHAPTER 3. NUMERICAL METHOD
zero thickness shell face
B
Fluid Fluid
zero thickness shell face
B
Fluid Fluid
zero thickness shell face
B
Fluid Fluid
Figure 3-9: Discrete shell boundary condition
3.6.5 Symmetric-Plane Boundary
Momentum equation: The velocity component normal to the plane is set to zero but its normal gradient is not zero. Additionally, and the normal gradient of the tangential components is set to zero but their values are not zero. This would result in a different system matrix for each velocity component which in turn would increase the memory requirements of the method and significantly complicate the PISO algorithm. To avoid that, the deferred correction algorithm is employed. First the velocity component at the boundary are extrapolated from the domain interior using Equation (3.60)
V~B=V~P+d~P f·
∇~V
P (3.60)
Diffusion term: Substituting in Equation (3.46) results in S~f·(∂ui
∂xj)f =|Sf| V~P
|dP f| n+1
| {z }
Implict Part
− |Sf| V~P
|dP f| n
+|Sf| 2V~n
|dP f| n
| {z }
Explicit Part
(3.61)
whereV~nis the normal velocity vector defined as V~n=
V~B·~n
~n (3.62)
Convection term: Since no flow can cross the boundary owing to symmetry, the
49 3.7. SOLUTION PROCEDURE
