**2. Basic transport equations**

For the description of the transport equations, the *x* and *z* coordinates of the Cartesian system will be taken as the independent variables, to which correspond, for velocity, the *u* and the *w* components. The coordinate transformation for the generalized mesh will be presented later on the present work.

The Navier-Stokes equations describe momentum conservation and, for a 2D situation, may be stated as follows:

$$\frac{\partial \left(\rho u\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left(\rho u^2\right) + \frac{\partial}{\partial \mathbf{z}} \left(\rho uw\right) = \frac{\partial}{\partial \mathbf{x}} \left[\Gamma \left(2\frac{\partial u}{\partial \mathbf{x}} - \frac{2}{3} \text{div}\vec{V}\right)\right] + \frac{\partial}{\partial \mathbf{z}} \left[\Gamma \left(\frac{\partial u}{\partial \mathbf{z}} + \frac{\partial w}{\partial \mathbf{x}}\right)\right] - \frac{\partial p}{\partial \mathbf{x}} \tag{1}$$

$$\frac{\partial \left(\rho u\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left(\rho u^2\right) + \frac{\partial}{\partial z} \left(\rho u \mathbf{w}\right) = \frac{\partial}{\partial \mathbf{x}} \left[\Gamma \left(2\frac{\partial u}{\partial \mathbf{x}} - \frac{2}{3} d \dot{\mathbf{u}} \vec{\mathcal{V}}\right)\right] + \frac{\partial}{\partial z} \left[\Gamma \left(\frac{\partial u}{\partial z} + \frac{\partial \mathbf{w}}{\partial \mathbf{x}}\right)\right] - \frac{\partial p}{\partial \mathbf{x}} \tag{1a}$$

where *p* [N/m*<sup>2</sup>* ] stands for pressure and *I* represents the buoyancy forces. The diffusion coefficient includes the contributions of the dynamic and turbulent viscosity, i.e., *<sup>Γ</sup>* <sup>=</sup>*<sup>μ</sup>* <sup>+</sup> *μt* <sup>N</sup> <sup>s</sup> / <sup>m</sup><sup>2</sup>

In turn, the conservation of mass law, or continuity equation, is

$$\frac{\partial \mathcal{P}}{\partial t} + \frac{\partial}{\partial \mathbf{x}} (\rho u) + \frac{\partial}{\partial z} (\rho w) = 0 \tag{2}$$

The energy conservation equation is obtained considering the transport equation for the enthalpy *ϕ* =*cp T* , where *cp* J/kg K is the specific heat and *T* [K] is its temperature. For a fluid medium, the corresponding equation is

$$\frac{\partial \left(\rho \mathbf{c}\_p T\right)}{\partial \mathbf{t}} + \frac{\partial}{\partial \mathbf{x}} \left(\rho \mathbf{c}\_p \mathbf{u} T\right) + \frac{\partial}{\partial \mathbf{z}} \left(\rho \mathbf{c}\_p \mathbf{w} T\right) = \frac{\partial}{\partial \mathbf{x}} \left[\boldsymbol{\Gamma} \frac{\partial \boldsymbol{T}}{\partial \mathbf{x}}\right] + \frac{\partial}{\partial \mathbf{z}} \left[\boldsymbol{\Gamma} \frac{\partial \boldsymbol{T}}{\partial \mathbf{z}}\right] + \mathbf{S}\_{\Gamma} \tag{3}$$

with the source term *ST* representing the heat generation rate per unit volume [W/m3 ] and *Γ* =(*μ* / Pr + *μt* / Pr*t*)*cp*, with Pr and Prt as the laminar and the turbulent Prandtl number. For solid regions, heat transfer is governed by conduction:

$$\frac{\partial \left(\rho \mathbf{c}\_p T\right)}{\partial t} = \frac{\partial}{\partial \mathbf{x}} \left(\lambda \frac{\partial T}{\partial \mathbf{x}}\right) + \frac{\partial}{\partial z} \left(\lambda \frac{\partial T}{\partial z}\right) + S\_T \tag{4}$$

where *λ* W/mK is the material thermal conductivity. In a multi-component fluid flow situation, different fluids are present, sharing the same velocity, pressure, temperature and turbulence quantities. The concentration of one of the components (the secondary fluid) is computed with a normal transport equation for a scalar:

$$\frac{\partial \left(\rho \,\phi\_{c2}\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left(\rho u \,\phi\_{c2}\right) + \frac{\partial}{\partial z} \left(\rho w \phi\_{c2}\right) = \frac{\partial}{\partial \mathbf{x}} \left(\Gamma \frac{\partial \phi\_{c2}}{\partial \mathbf{x}}\right) + \frac{\partial}{\partial z} \left(\Gamma \frac{\partial \phi\_{c2}}{\partial z}\right) \tag{5}$$

where *ϕc2* is the mass fraction, or concentration, of the secondary fluid. The diffusion coefficient is given by

$$
\Gamma = \rho D\_{\phi\_t} + \frac{\mu\_t}{\mathbf{Sc\_t}} \tag{6}
$$

where *D<sup>ϕ</sup><sup>c</sup>* [m2 /s] is the kinematic diffusivity that characterizes the fluids pair and *Sc*<sup>t</sup> is the turbulent Schmidt number. The concentration of the secondary fluid plays an active role in the flow field since fluid properties depend on the concentration of each component. The mixture properties (such as viscosity, specific heat, etc.) are determined by a weighted average of the properties of the components. Thus, for a generic property *X*, we have

$$X = \phi\_{c1} X\_1 + \phi\_{c2} X\_2 \tag{7}$$

where the constraint *ϕc1* + *ϕc2* =*1* holds. Exception is made for density, where a simple analysis may prove that the mixture density *ρ* is given by

$$\rho = \left(\frac{\phi\_{cJ}}{\rho\_I} + \frac{\phi\_{c\_2}}{\rho\_2}\right)^{-I} \tag{8}$$

### **3. Turbulence modelling**

**2. Basic transport equations**

264 Numerical Simulation - From Brain Imaging to Turbulent Flows

( ) ( )

 ¶

 ¶

rr

medium, the corresponding equation is

¶

 ¶ r

regions, heat transfer is governed by conduction:

¶ r

( *p* )

¶

( )

*p*

*Γ* =(*μ* / Pr + *μt* / Pr*t*)*cp*, with Pr and Prt

¶

¶ r  ¶

rr

 ¶

 ¶

> ¶

 ¶

In turn, the conservation of mass law, or continuity equation, is

¶r ¶

¶¶

( ) ( *w* )

 ¶  r

 ¶

¶¶

> ¶

( ) ( ) ( ) *<sup>2</sup> <sup>u</sup> u 2 uw p u uw 2 divV t x z x x3 z z x x*

+ + = - + +- <sup>ê</sup> ç ÷ ç÷ úê ú <sup>ë</sup> è ø èø ûë û

coefficient includes the contributions of the dynamic and turbulent viscosity, i.e.,

( ) ( )

The energy conservation equation is obtained considering the transport equation for the enthalpy *ϕ* =*cp T* , where *cp* J/kg K is the specific heat and *T* [K] is its temperature. For a fluid

> ¶¶

+ + =++ êúêú

*c T T T c uT c T <sup>S</sup> t x z xxzz*

with the source term *ST* representing the heat generation rate per unit volume [W/m3

*c T T T <sup>S</sup> t xx zz*

æöæö =++ ç÷ç÷ èøèø

 ¶¶

 ¶¶

 l

¶¶

 ¶¶

l

*p p T*

 ¶

++ = *u w0*

 r

r

 ¶

 ¶

 ¶

+ + = - + +- <sup>ê</sup> ç ÷ç÷ úê ú êë è øèø ûë û <sup>r</sup> *<sup>2</sup> <sup>u</sup> u 2 u w <sup>p</sup> u uw 2 divV*

*t x z x x3 z z x x* (1)

 G

 G

on the present work.

be stated as follows:

¶

¶

 ¶

 ¶

( )

¶

¶ r

¶

where *p* [N/m*<sup>2</sup>*

*<sup>Γ</sup>* <sup>=</sup>*<sup>μ</sup>* <sup>+</sup> *μt* <sup>N</sup> <sup>s</sup> / <sup>m</sup><sup>2</sup>

¶ r

For the description of the transport equations, the *x* and *z* coordinates of the Cartesian system will be taken as the independent variables, to which correspond, for velocity, the *u* and the *w* components. The coordinate transformation for the generalized mesh will be presented later

The Navier-Stokes equations describe momentum conservation and, for a 2D situation, may

¶

 ¶

 ¶

é æ ö æö ùé ù

] stands for pressure and *I* represents the buoyancy forces. The diffusion

 ¶

é æ öæö ùé ù

 ¶¶

 ¶¶

 ¶¶

 ¶¶

G

*tx z* (2)

 ¶¶

éùéù

ëûëû

*T*

as the laminar and the turbulent Prandtl number. For solid

GG

¶¶¶¶

G

¶

 ¶

 ¶

<sup>r</sup> (1a)

 ¶

(3)

] and

(4)

Two of the most popular turbulence models are presented next.

#### **3.1. The** *k-ε* **turbulence model**

The standard formulation of this turbulence model is described in, e.g., [2]. The turbulent viscosity is given by

$$
\mu\_t = C\_\mu \frac{\rho k^2}{\varepsilon} \tag{9}
$$

The turbulence kinetic energy, *k* [m2 /s2 ], and its dissipation rate, *ε* [m2 /s3 ], are computed with the following transport equations:

$$\frac{\partial \left(\rho k\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left(\rho uk\right) + \frac{\partial}{\partial z} \left(\rho vk\right) = \frac{\partial}{\partial \mathbf{x}} \left[\left(\mu + \frac{\mu\_t}{\sigma\_k}\right) \frac{\partial k}{\partial \mathbf{x}}\right] + \frac{\partial}{\partial z} \left[\left(\mu + \frac{\mu\_t}{\sigma\_k}\right) \frac{\partial k}{\partial z}\right] + P\_k - \rho \mathbf{c} + G\_T \tag{10}$$

$$\begin{aligned} \frac{\partial \left(\rho \boldsymbol{x}\right)}{\partial t} + \frac{\partial}{\partial \boldsymbol{x}} \left(\rho \boldsymbol{u} \boldsymbol{\varepsilon}\right) + \frac{\partial}{\partial \boldsymbol{z}} \left(\rho \boldsymbol{u} \boldsymbol{\varepsilon}\right) &= \\ \frac{\partial}{\partial \boldsymbol{x}} \left[ \left(\mu + \frac{\mu\_t}{\sigma\_\varepsilon}\right) \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{x}} \right] + \frac{\partial}{\partial \boldsymbol{z}} \left[ \left(\mu + \frac{\mu\_t}{\sigma\_\varepsilon}\right) \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} \right] + \frac{\varepsilon}{k} \left(C\_I P\_k - C\_2 \rho \boldsymbol{z} + C\_3 G\_T\right) \end{aligned} \tag{11}$$

The term *Pk* represents the production rate of *k* as the results of the velocity gradients:

$$P\_k = \mu\_t \left[ 2\left(\frac{\partial u}{\partial \mathbf{x}}\right)^2 + 2\left(\frac{\partial w}{\partial \mathbf{z}}\right)^2 + \left(\frac{\partial u}{\partial \mathbf{z}} + \frac{\partial w}{\partial \mathbf{x}}\right)^2 \right] \tag{12}$$

while the term *G*T accounts for the production or destruction of *k* and *ε* due to the thermal gradients:

$$\mathbf{G}\_T = -\beta \mathbf{g} \frac{\mu\_t}{\text{Pr}\_t} \frac{\partial T}{\partial \mathbf{z}} \tag{13}$$

with *g* [m/s2 ] being the gravity acceleration and *β* [K−1] being the fluid dilatation coefficient. The remaining constants are

$$C\_{\mu} = 0.09; \; \sigma\_{k} = 1.0; \; \sigma\_{\varepsilon} = 1.3; \; C\_{I} = 1.44; \; C\_{2} = 1.92; \; C\_{3} = 1.44 \tag{14}$$

In the proximity of a wall, the previous equations should be modified to account for the viscous effects that become predominant. Wall functions ensure the connection between the viscous sub-layer and the inertia layer, at a location established by the *y*<sup>+</sup> value:

$$\mathbf{y}^+ = \frac{\rho u\_\tau y}{\mu} = \frac{\rho C\_\mu^{0.25} \sqrt{k} y}{\mu} \tag{15}$$

where *u<sup>τ</sup>* represents the friction velocity and *y* [m] is the distance to the wall. The utilization of wall laws does not require mesh nodes in the viscous sub-layer, which is a major advantage from the computational standpoint. The momentum flux per unit area along the direction normal to the wall is given by the wall shear stress, which is computed differently depending on the location of the wall neighbour node relatively to the transition between the viscous and the inertia sub-layers. Denoting by *v* the generic velocity component parallel to the wall (which may, actually, be *u* or *w*) and by *y* the distance to the wall, we have

The turbulence kinetic energy, *k* [m2

 ¶

266 Numerical Simulation - From Brain Imaging to Turbulent Flows

 ¶

¶

¶¶

m

e

m

 s¶

 r

( ) ( ) ( )

*u w tx z*

¶e

re

the following transport equations:

r

¶

¶ re

> ¶

¶

 ¶

¶ r

¶

gradients:

with *g* [m/s2

The remaining constants are

m

s

/s2

*t x z x xz z*

 r e

> m

m

 ¶

 ¶

 ¶

 ¶

++ =

 ¶

é ùé ù æö æö

 ¶

ë ûë û èø èø

2 2 *k t*

m

( ) ( ) ( ) *t t k T*

m

 s¶

*CP C CG x xz z*

 m

ê úê ú ç÷ ç÷ + + + + -+

 s¶

 e

The term *Pk* represents the production rate of *k* as the results of the velocity gradients:

*u w uw <sup>P</sup> x z zx*

*<sup>k</sup> k k uk wk P G*

+ + = + + + +- + ê úê ú ç÷ ç÷

 ¶¶

( ) *t t 1k 2 3 T <sup>k</sup>*

22 2

é ù æö æ öæ ö ¶ ¶ ¶¶ = + ++ ê ú ç÷ ç ÷ç ÷ èø è øè ø ¶ ¶ ¶¶ ë û

while the term *G*T accounts for the production or destruction of *k* and *ε* due to the thermal

Pr *<sup>t</sup> <sup>T</sup>*

*<sup>T</sup> G g <sup>z</sup>* m b

<sup>2</sup> ;; ; ; ; *C 0.09 1.0 1.3 C 1.44 C 1.92 C 1.44*

In the proximity of a wall, the previous equations should be modified to account for the viscous effects that become predominant. Wall functions ensure the connection between the viscous

> *0.25 u y C ky*

> > m

where *u<sup>τ</sup>* represents the friction velocity and *y* [m] is the distance to the wall. The utilization of wall laws does not require mesh nodes in the viscous sub-layer, which is a major advantage

r m

= = == = =

 e

 s

sub-layer and the inertia layer, at a location established by the *y*<sup>+</sup>

*y*

r t

m

*t*

] being the gravity acceleration and *β* [K−1] being the fluid dilatation coefficient.

*k1 3* (14)

¶e

e

], and its dissipation rate, *ε* [m2

*k k*

 ¶

é ùé ù æö æö

/s3

 ¶

ë ûë û èø èø (10)

re

¶ = - ¶ (13)

value:

<sup>+</sup> = = (15)

 m

 s¶

m ], are computed with

(11)

(12)

re

$$
\sigma\_1 \mathbf{y}^+ \le I \\
1.63 \quad \Rightarrow \quad \tau\_0 = \mu \frac{\nu\_0 - \nu}{\mathbf{y}} \tag{16}
$$

$$\text{v.y}^+ \ge II.63 \quad \Rightarrow \quad \tau\_0 = \mu \frac{\rho C\_\mu^{0.25} \text{yr} \sqrt{k}}{\ln \left( E^\ast \text{y}^+ \right)} \text{(v}\_0 - \nu \text{)}\tag{16a}$$

where *v0* is the wall velocity, *E*\* =*9.793* for smooth walls and *χ*=*0.4187* is the von Karman constant. In the wall neighbourhood, the production term given by Eq. (12) is computed assuming a Couette flow:

$$P\_k = \mu\_t \left(\frac{\partial \nu}{\partial \mathbf{y}}\right)^2 \tag{17}$$

where, once again, *v* is the generic velocity component parallel to the wall and *y* is the generic distance to the wall. Due to its significant variations near the wall, *ε* is averaged for the term *ρε* in Eq. (10):

$$\dot{q}\mathbf{y}^+ \le I\mathbf{2} \quad \Rightarrow \quad \dot{q} = \frac{\mu c\_p}{\text{Pr}\,\mathbf{y}} (T\_0 - T) \tag{18}$$

$$
\mu y^\* > 11.63 \quad \Rightarrow \quad \rho \varepsilon = \rho C\_{\mu}^{\ 075} k^{1.5} \frac{\ln \left( E^\* y^\* \right)}{\varkappa y} \tag{18a}
$$

For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. For the dissipation rate, Eq. (11) is not employed in the node adjacent to the wall. Instead, the dissipation rate is given by

$$
\varepsilon\_0 = \frac{C\_{\mu}^{\prime 0.75} k^{\prime 1.5}}{\mathcal{Z}\mathcal{Y}} \tag{19}
$$

As for momentum, the energy flux is computed differently depending on the *y*<sup>+</sup> values. In this case, we have

$$\dot{q}\mathbf{y}^+ \le I\mathbf{2} \quad \Rightarrow \quad \dot{q} = \frac{\mu c\_p}{\text{Pr}\,\mathbf{y}}(T\_0 - T) \tag{20}$$

$$\dot{q}\mathbf{y}^{+} > I\mathbf{2} \quad \Rightarrow \quad \dot{q} = \frac{\rho c c\_{p} C\_{\mu}{}^{0.25} \sqrt{k}}{\text{Pr}\_{\text{f}} \left[ \frac{I}{\mathcal{X}} \ln \left( E^{\star} \mathbf{y}^{+} \right) + P^{\star} \right]} \text{( $T\_{0}$  - T)} \tag{20a}$$

where *T*0 is the wall temperature and *P*\* is the so-called Jayatillaka function:

$$P^\* = 9.24 \left[ \left( \frac{\text{Pr}}{\text{Pr}\_t} \right)^{0.75} - I \right] \left[ I + 0.28 \left( \exp \left( -0.007 \frac{\text{Pr}}{\text{Pr}\_t} \right) \right) \right] \tag{21}$$

#### **3.2. The** *k-ω***SST turbulence model**

The *k-ω* SST model [3] represents a combination of the *k-ε* and the *k*-*ω* models. According to Menter et al. [4], the *k-ω* model is more accurate near the wall but presents a high sensitivity to the *ω* values in the free-stream region, where the *k-ε* model shows a better behaviour. The *k-ω* SST model represents a blend of the two, through a weighting factor computed based on the nearest wall distance. The governing equations are

$$\frac{\partial(\rho k)}{\partial t} + \frac{\partial(\rho uk)}{\partial \mathbf{x}} + \frac{\partial(\rho vk)}{\partial \mathbf{z}} = \overline{P\_k} - \beta^\* \,\rho \alpha k + \frac{\partial}{\partial \mathbf{x}} \Big( (\mu + \sigma\_k \mu\_t) \frac{\partial k}{\partial \mathbf{x}} \Big) + \frac{\partial}{\partial \mathbf{z}} \Big( (\mu + \sigma\_k \mu\_t) \frac{\partial k}{\partial \mathbf{z}} \Big) \tag{22}$$

$$
\frac{
\hat{\sigma}
\left(
\rho \alpha \right)
}{
\hat{\alpha}t
} + \frac{
\hat{\sigma}
\left(
\rho \alpha \alpha \right)
} + \frac{
\hat{\sigma}
\left(
\rho \gamma \alpha \alpha \right)
}{
\hat{\alpha}z
} = \frac{
\hat{\sigma}
}{
\hat{\alpha}
} \left(
\left(
\mu + \sigma\_{\alpha 0} \mu\_t
\right) \frac{
\hat{\alpha}\alpha
}{
\hat{\alpha}z
} \right) + \frac{
\hat{\sigma}
}{
\hat{\alpha}z
} \left(
\left(
\mu + \sigma\_{\alpha 0} \mu\_t
\right) \frac{
\hat{\alpha}\alpha
}{
\hat{\alpha}z
} \right) + 
\frac{
\hat{\sigma}
\left(
\left(
\mu + \sigma\_{\alpha 0} \mu\_t
\right) \frac{
\hat{\alpha}z
}{
\hat{\alpha}z
} \right)}{
\hat{\alpha}z
} \tag{23}
$$

$$
+ \frac{
\alpha \overline{P\_k}}{
\nu\_l} - \beta \rho \alpha \nu^2 + 2 \left(
I - F\_1
\right) \rho \sigma\_{\alpha 2}
\frac{I}{
\hat{\alpha}
} \left(
\frac{
\hat{\alpha}k
}{
\hat{\alpha}z
} \frac{
\hat{\alpha}\alpha
}{
\hat{\alpha}z
} + \frac{
\hat{\alpha}k
}{
\hat{\alpha}z
} \frac{
\hat{\alpha}\alpha
}{
\hat{\alpha}z
} \right) \tag{24}
$$

where *ω* is the frequency of dissipation of turbulent kinetic energy [s-1]. The production of turbulent kinetic energy is limited to prevent the build-up of turbulence in stagnant regions:

$$
\overline{P\_k} = \min\left(P\_k, 10\,\beta^\*\,\rho k\alpha\right) \tag{24}
$$

The weighting function *F*1 is given by

$$F\_I = \tan h \left[ \left\{ \min \left[ \max \left( \frac{\sqrt{k}}{\beta^\* \alpha y}; \frac{\\$00 \,\mathrm{v}}{\mathrm{y}^2 \alpha} \right); \frac{4 \,\rho \sigma\_{oo2} k}{C D\_{ko} y^2} \right] \right\}^4 \right] \tag{25}$$

and

(20)

(21)

As for momentum, the energy flux is computed differently depending on the *y*<sup>+</sup> values. In this

m<sup>+</sup> £ Þ= - & *<sup>p</sup>*

*<sup>c</sup> y 12 q T T*

r

&

> Þ= - é ù

*y 12 q T T*

*p*

c

*t*

( ) Pr

*y*

.

*cC k*

*0 25*

+ ê ú ë û

+

\* \* Pr ln

The *k-ω* SST model [3] represents a combination of the *k-ε* and the *k*-*ω* models. According to Menter et al. [4], the *k-ω* model is more accurate near the wall but presents a high sensitivity to the *ω* values in the free-stream region, where the *k-ε* model shows a better behaviour. The *k-ω* SST model represents a blend of the two, through a weighting factor computed based on

> m s m

> > w

w

æ ö ¶¶ ¶¶

¶¶ ¶ ¶ ¶¶ ¶ æ öæ ö + + =- + + + + ç ÷ç ÷ ¶¶ ¶ ¶ ¶¶ ¶ è øè ø (22)

( ) ( ) ( ) \* () () *<sup>k</sup> k t k t k uk wk k k P k tx z x xz z*

( ) ( ) ( ) () () *t t*

*<sup>P</sup> 1k k 21 F <sup>v</sup> xx zz* w

+ - +- + ç ÷

*P P 10 k k k* = min , \* (

rs

*t x z x xz z* w

*2*

where *ω* is the frequency of dissipation of turbulent kinetic energy [s-1]. The production of turbulent kinetic energy is limited to prevent the build-up of turbulence in stagnant regions:

> b rw

w

¶¶ ¶ ¶ ¶¶ ¶ æ öæ ö + + =+ ++ + ç ÷ç ÷ ¶ ¶ ¶ ¶ ¶¶ ¶ è øè ø

m sm

b rw m

*0*

( ) ( )

is the so-called Jayatillaka function:

*0*

*<sup>1</sup> Ey P* (20a)

 m s m

> w

è ø ¶¶ ¶¶ (23)

) (24)

 m sm

 w  w

case, we have

+

268 Numerical Simulation - From Brain Imaging to Turbulent Flows

where *T*0 is the wall temperature and *P*\*

**3.2. The** *k-ω***SST turbulence model**

rr

rw

the nearest wall distance. The governing equations are

*u w*

brw  r w

( <sup>1</sup>) *<sup>k</sup> <sup>2</sup>*

 r

 r w

*t*

a

$$\hat{I}\_{CD\_{k0}} = \max\left(2\,\rho\,\sigma\_{o\nu2}\frac{l}{o\nu}\frac{\partial k}{\partial \mathbf{x}\_j}\frac{\partial o}{\partial \mathbf{x}\_j}, I0^{-I0}\right) \tag{26}$$

where *y* represents the distance to the neighbour wall and *v* is the laminar dynamic viscosity. *F*<sup>1</sup> is zero away from the wall (*k-ε* model) and changes to unit inside the boundary layer (*k*-*ω* model) with a smooth transition based on *y*. The turbulent viscosity is computed as

$$\text{av}\_t = \frac{a\_I k}{\max\left(a\_I \text{or}; SF\_2\right)}\tag{27}$$

where *S* is the invariant measure of the strain rate given by

$$S = \sqrt{S\_{ij} S\_{ij}} \quad ; \quad S\_{ij} = \frac{I}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{28}$$

and

$$F\_2 = \tanh\left\{ \left[ \max\left( \frac{2\sqrt{k}}{\beta^\* \alpha y}; \frac{500\,\nu}{y^2 \alpha} \right) \right]^2 \right\} \tag{29}$$

The constants are computed as a blend of the *k-ε* and *k-ω* models by the following generic equation:

$$
\alpha = F\_I \alpha\_I + (1 - F\_I)\alpha\_2 \tag{30}
$$

The constants are *α*1= 5/9; *β*1= 3/40; *σ<sup>k</sup>*<sup>1</sup> = 0.85; *σω*<sup>1</sup> = 0.5; *α*2= 0.44; *β*<sup>2</sup> = 0.0828; *σ<sup>k</sup>*<sup>2</sup> = 1; *σω*2 = 0.856; *β\** = 0.09.

The near wall treatment for momentum and turbulence equations follows the proposal described in [5]. The basic principle behind the automatic wall functions is to switch from a low-Reynolds number formulation to a wall function based on the grid nodes proximity to the wall. According to these authors, the automatic wall treatment avoids the deterioration of results typical of low-Reynolds models when applied on too coarse meshes.

The known solutions for *ω* in the viscous (linear) and in the logarithmic near wall region are

$$
\rho\_{\rm vis} = \frac{6\,\text{\nu}}{0.075\,\text{\textnu}\,\text{\textnu}^2} \; ; \quad \rho\_{\rm log} = \frac{u\_{\rm tr}}{0.3\,\text{\textnu}\,\text{\textnu}} \tag{31}
$$

The imposed value for *ω* at the first node close to a wall is

$$
\alpha o\_I = \sqrt{\alpha\_{\rm vis}^2 + \alpha\_{\rm log}^2} \tag{32}
$$

For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. In turn, for the momentum equations, a similar reasoning applies, with expressions for the shear velocity in the viscous and in the logarithmic region:

$$u\_{\tau}^{\text{vis}} = \frac{U\_I}{\mathbf{y}^+}; \quad u\_{\tau}^{\text{log}} = \frac{\mathcal{X}U\_I}{\ln\left(E\,\mathbf{y}^+\right)}\tag{33}$$

with *U*1 being the fluid velocity relative to the wall velocity. The wall shear stress is computed as follows:

$$
\mu\_{\tau} = \sqrt[4]{\left(\mu\_{\tau}^{\mathrm{vis}}\right)^{4} + \left(\mu\_{\tau}^{\mathrm{log}}\right)^{4}} \tag{34}
$$
