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~_{n}is 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~_{B}oru_{iB}which means that

Diffusion term: Following the terminology in Figure3-7
S~_{f}·(∂u_{i}

∂x_{j})_{f} =|S_{f}| u_{iB}−u_{iP}

|d_{P f}|

| {z }

OrthognalContribution

+ k~_{f}·(∂u_{i}

∂x_{j})_{B}

| {z }

Non−orthonalContribution

(3.46)

wherek~_{f} = _{|S}^{S}^{~}^{f}

f|−_{|d}^{d}^{~}^{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 F_{B}=S_{f}·V~_{B}. The convection boundary contribution is obtained by substituting in
Equation (3.19) resulting in:

ρ_{f}F_{f}~V_{f} =ρ_{B}F_{B}~V_{B} (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)

∇p^{m+1}
a_{p}

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 (P_{B}) and it is treated in a similar
manner to Equation (3.46) as follows

∇p^{m+1}
a_{p}

B·S~_{f} = 1
a_{p}

B|S_{f}| p_{B}−p_{P}

|d_{P f}|

| {z }

OrthognalContribution

+ ~k_{f}·(∂p

∂x_{j})_{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~_{B}oru_{iB}which 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 (P_{k}) in Equations (2.10), (2.11), (2.14), (2.15),
(2.25) and (2.26) is modified as

P_{k}=|τ_{w}| u_{τ}

κ|d_{n}| (3.52)

where|d_{n}|is the normal distance to the wall computed from the wall which is compute
from

~d_{n}=|d~_{P f}·S~_{f}
S~_{f}

S~_{f}·S~_{f}| (3.53)

τ_{w}is the wall shear stress

τw=µ|V~_{P}−V~_{B}|

|d_{n}| (3.54)

CHAPTER 3. NUMERICAL METHOD

u_{τ} is the frictional velocity computed from
u_{τ} =C^{1/4}_{µ} p

k_{P} (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/4}k_{P}^{3/2}

κ|d_{n}| (3.56)

ω_{P}=

√k_{P}

κC_{µ}^{1/4}|d_{n}| (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^{+}= ^{ρ|d}^{n}^{|u}^{τ}

µ with u_{τ} computed from Equation (3.55) and the normalized velocity
parallel to the wall u^{+} = ^{|}^{V}^{~}^{P}^{−}_{u}^{V}^{~}^{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} 10^{0} 10^{1} 10^{2}

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}·(∂u_{i}

∂x_{j})_{f} =|S_{f}| V~_{P}

|d_{P f}|
n+1

| {z }

Implict Part

− |S_{f}| V~_{P}

|d_{P f}|
n

+|S_{f}| 2V~_{n}

|d_{P f}|
n

| {z }

Explicit Part

(3.61)

whereV~_{n}is 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