<!-- source-page: 391 -->

$$
R e _ {y} = \frac {\rho k ^ {1 / 2} y}{\mu},
$$

where $y$ is the wall-normal distance. The two-layer approach combines the algebraic equation with the transport equation. To provide a smooth transition, the work of Jongen (1998) is followed by implementing a smooth blending function to bridge both formulations:

$$
\lambda_ {\varepsilon} = \frac {1}{2} \left[ 1 + \tanh \left(\frac {\mathrm{Re} _ {\mathrm{y}} - \mathrm{Re} _ {\mathrm{y}} ^ {*}}{\mathrm{A}}\right) \right],
$$

where defines the width of the blending function. We define a width so $\lambda _ { \varepsilon }$ will be within 1% of its far-field value. Given a variation of $\Delta R e _ { y }$ one obtains

$$
A = \frac {| \Delta R e _ {y} |}{\tanh 0 . 9 8}.
$$

The final blending of is done through combining the transport and algebraic equations at the discretized level:

$$
a _ {p} \varepsilon_ {p} + \lambda_ {\varepsilon} \Sigma_ {n b} a _ {n b} \varepsilon_ {n b} = \lambda_ {\varepsilon} S _ {\varepsilon} + (1 - \lambda_ {\varepsilon}) a _ {p} \varepsilon_ {\mathrm{2layer}},
$$

where $a _ { p }$ and $a _ { n b }$ are the diagonal and off-diagonal terms, respectively, and $S _ { \varepsilon }$ is the right-hand side of the discrete transport equation. Finally, the blending of the eddy viscosity is conducted in a straightforward way since it only involves two algebraic equations:

$$
\mu_ {T} ^ {\prime} = \lambda_ {\varepsilon} \mu_ {T} + (1 - \lambda_ {\varepsilon}) \mu_ {\mathrm{T2layer}}.
$$

# Wall functions

The two-layer approach takes care of the near-wall integration of the equation; however, it does not provide any modeling for the turbulent kinetic energy equation. Therefore, to model the budget of the turbulent kinetic energy at the inner layer correctly, an additional near-wall treatment needs to go into the k equation. A hybrid-wall function treatment is implemented to integrate the k equation throughout the inner layer. The hybrid wall function uses the law-of-the-wall and equilibrium assumptions to approximate the sources of k when the mesh resolution is not fine enough to resolve the inner layer. However, if the mesh is fine enough to resolve the entire inner-layer, the hybrid-wall function recovers the viscous sublayer relations for the k sources. The law-of-the-wall and viscous sublayer sources are finally blended to provide a smooth formulation that allows the model to provide accurate solutions irrespective of the mesh’s near-wall resolution.

The-law-of-the-wall is a universal velocity profile that wall-bounded flows develop in the absence of pressure gradients. The velocity profile consists of two well-defined regions:

<!-- source-page: 392 -->

$$
V ^ {+} = \left\{ \begin{array}{l l} y ^ {+} & \text {if} y ^ {+} \leq y _ {c} ^ {+} \quad \text {Viscous sublayer} \\ \frac {1}{\kappa} \ln (E y ^ {+}) & \text {if} y ^ {+} > y _ {c} ^ {+} \quad \text {Logarithmic layer}, \end{array} \right.
$$

where

$$
V ^ {+} = \frac {V}{v _ {\tau}},
$$

$$
y ^ {+} = \frac {y v _ {\tau}}{\nu},
$$

$$
v _ {\tau} = \sqrt {\frac {\tau_ {w a l l}}{\rho}}.
$$

is the wall tangent velocity, $\nu$ is the kinematic viscosity, $\rho$ is the density, $\tau _ { w a l l }$ is the shear stress at the wall, $y _ { c } ^ { + }$ is the intersection point of the viscous sublayer and the logarithmic layer, $\kappa = 0 . 4 1$ is the Von Karman constant, and $E = 8 . 4$ is the law-of-the-wall’s logarithmic constant. The viscous sublayer is the region where momentum transfer is dominated by molecular diffusion; turbulence here is virtually absent. The logarithmic layer or fully turbulent regime establishes the onset of the fully turbulent regime; here momentum transport is fully dominated by turbulent transport relegating molecular diffusion to a second-order effect.

The standard law of the wall profile is limited in its usage. For example, in recirculating flows the turbulent kinetic energy, $k ,$ becomes zero at separation and reattachment points, where, by definition, $v _ { \tau }$ is zero. This singular behavior causes the predicted results to be erroneous. To overcome this, the standard law of the wall is modified based on a new scale for the friction velocity following the method proposed by Launder and Spalding (1974). The modified friction velocity is given by

$$
v ^ {*} = C _ {\mu} ^ {1 / 4} k ^ {1 / 2},
$$

which does not suffer from a singular behavior at flow reattachment, separation, and at points of flow impingement. Correspondingly, the wall distances are re-scaled as follows:

$$
y ^ {*} = \frac {y v ^ {*}}{\nu} = \frac {y C _ {\mu} ^ {1 / 4} k ^ {1 / 2}}{\nu}.
$$

The modified law of the wall reduces to the standard law of the wall under the conditions of uniform wall-shear stress and when the generation and dissipation of turbulent kinetic energy are in balance $( \mathrm { i . e . }$ , when the turbulence structure is in equilibrium). Under such conditions, $v ^ { * } \approx v ^ { + }$ and, thus, $y ^ { * } \approx y ^ { + }$ .

The wall-shear stress for the modified law of the wall can be evaluated as (Albets-Chico, et al., 2008)

$$
\tau_ {w a l l} = \rho u ^ {*} \frac {V _ {p}}{V ^ {+} (y _ {p} ^ {*})},
$$

<!-- source-page: 393 -->

where the subscript p denotes the wall element center at which all the quantities of interest are evaluated. Using the wall function the k sources are derived to correct the turbulent kinetic energy budget at the inner layer as follows:

$$
\operatorname{Production} (k) = \left\{ \begin{array}{l l} 0. 0 & \text {if} y ^ {*} \leq y _ {c} ^ {+} \\ \frac {\tau_ {w a l l} ^ {2}}{\rho \kappa u ^ {*} y _ {n}} & \text {if} y ^ {*} > y _ {c} ^ {+}. \end{array} \right.
$$

Following a similar procedure the value of the dissipation rate for k is given by

$$
\operatorname{Dissipation} (k) = \left\{ \begin{array}{l l} \frac {2 \mu k}{y ^ {2}} & \mathrm{if} y ^ {*} \leq y _ {c} ^ {+} \\ \frac {\rho u ^ {* 3}}{\kappa y} & \mathrm{if} y ^ {*} > y _ {c} ^ {+}. \end{array} \right.
$$

Having defined the sources in the viscous sublayer and logarithmic layer, a blending approach is implemented to provide a smooth transition that allows us to operate on meshes with arbitrary near-wall resolution:

$$
\mathrm{Production} (k) = (1 - f) \frac {\tau_ {w a l l} ^ {2}}{\rho \kappa u ^ {*} y _ {n}},
$$

$$
\mathrm{Dissipation} (k) = \rho \bigg (f \frac {2 \nu k}{y ^ {2}} + (1 - f) \frac {u ^ {* 3}}{\kappa y} + f (1 - f) \Big (\frac {2 \nu k}{y ^ {2}} + \frac {u ^ {* 3}}{\kappa y} \Big) \bigg),
$$

$$
u ^ {*} = \sqrt {f \frac {\nu V _ {p}}{y _ {p}} + (1 - f) C _ {\mu} ^ {1 / 2} k},
$$

$$
f = e ^ {\frac {- R _ {e y}}{y _ {c} ^ {+}}},
$$

where

$$
R _ {e y} = \frac {\sqrt {k} y}{\nu}.
$$

The boundary conditions enforced at the wall for both the k and transport equations is a zero flux (i.e., homogeneous Neumann boundary conditions) at the walls.

# Implementation in the momentum equation

For cases where the mesh resolution is not enough to capture the near-wall gradients, a near-wall model is required to provide the correct wall-shear stress in coarse meshes. The wall shear is obtained from the wall-function approach through an effective edge viscosity:

$$
\mu_ {\mathrm{eff}} = \frac {\tau_ {w a l l} y}{V} = \left(\frac {\rho v _ {\tau} V}{V ^ {+}}\right) \frac {y}{V}.
$$

<!-- source-page: 394 -->

# Energy wall functions

The wall-function approach can be extended to the energy equation by using the temperature law-ofthe-wall, which is a semi-empirical universal temperature profile obtained in equilibrium wall-bounded flows, when the temperature, T, and the wall-normal distance, y, are normalized with wall units. The standard temperature wall function is defined as

$$
T ^ {+} (y ^ {+}) = \frac {(T _ {w} - T) \rho c _ {p} u _ {\tau}}{\dot {q}} = \left\{ \begin{array}{l l} P _ {r} y ^ {+} & \mathrm{if} y ^ {+} \leq y _ {T c} ^ {+} \\ P _ {r T} \Big (\frac {1}{\kappa} \ln (E y ^ {+}) + P (P _ {r}, P _ {r T}) \Big) & \mathrm{if} y ^ {+} > y _ {T c} ^ {+}, \end{array} \right.
$$

where $y _ { T c } ^ { + }$ is the intersection point of the viscous-sublayer and the logarithmic layer in the temperature wall-function, $T _ { w }$ is the wall temperature, $P _ { r }$ is the Prandtl number, $P _ { r T }$ is the turbulent Prandtl number, $\dot { q }$ is the wall-heat flux, $c _ { p }$ is the specific heat coefficient at constant pressure, and $P$ is computed using the Jayatilleke (1969) expression:

$$
P \left(P _ {r}, P _ {r T}\right) = 9. 2 4 \left[ \left(\frac {P _ {r}}{P _ {r T}}\right) ^ {3 / 4} - 1 \right] \left(1 + 0. 2 8 e ^ {- 0. 0 0 7 P _ {r} / P _ {r T}}\right).
$$

For the hybrid wall-function approach a continuous temperature wall function proposed by Kader (1981) is implemented

$$
T ^ {+} (y ^ {+}) = e ^ {- \Gamma} P _ {r} y ^ {+} + e ^ {- 1 / \Gamma} P _ {r T} \bigg (\frac {1}{\kappa} \mathrm{ln} (E y ^ {+}) + P (P _ {r}, P _ {r T}) \bigg),
$$

with the blending function defined as

$$
\Gamma = \frac {0 . 0 1 (P _ {r} y ^ {+}) ^ {4}}{1 + 5 P _ {r} ^ {3} y ^ {+}}.
$$

Finally, the heat flux is obtained from the precomputed $v _ { \tau }$ flow properties temperature field and the continuous temperature wall function:

$$
\dot {q} = \frac {(T _ {w} - T) \rho c _ {p} v _ {\tau}}{T ^ {+}}.
$$

# Implementation in the energy equation

For cases where the mesh resolution is not enough to capture the near-wall gradients, a near-wall model is required to provide the correct wall-heat flux in coarse meshes. The wall-heat flux is obtained from the wall-function approach through an effective edge heat conductivity:

$$
\kappa_ {\mathrm{eff}} = \dot {q} \frac {y}{(T _ {w} - T)}.
$$

<!-- source-page: 395 -->

# k–omega SST turbulence model

The $\pmb { k } \mathrm { - } \omega$ SST turbulence model is a two-equation turbulence model that uses a different dissipation parameter than the $\pmb { k } { \mathcal { E } }$ turbulence model; namely, the specific energy dissipation rate, $\omega ~ \approx ~ \varepsilon / k$ introduced by Menter in 1994. Although $\omega$ varies even more rapidly than $\varepsilon$ and takes on very large values near the wall, the model is relatively insensitive to the wall values. One of the most attractive properties of $\pmb { k } \mathrm { - } \omega$ models is that they can be employed throughout the viscous sublayer without further modification (as would required by $\pmb { k } \widetilde { \varrho }$ models). However, a weakness of the standard $\pmb { k } \mathrm { - } \omega$ model is that it is much more sensitive than $\pmb { k } \mathrm { - } \varepsilon$ to the values of turbulence parameters in the freestream and, hence, to inflow turbulence values.

Menter (1994) introduced a hybrid model that blends the standard $\pmb { k } \mathrm { - } \omega$ model near to the wall with a transformed version of the standard $\pmb { k } \mathrm { - } \varepsilon$ model (into $\pmb { k } \mathrm { - } \omega$ form) far from the wall. This hybrid model adds an additional cross-diffusion term, in regions away from the wall, to the transport equation for a specific dissipation rate. It is this term that reduces the sensitivity of the model to freestream turbulence values. To complete the $\pmb { k } \mathrm { - } \omega$ SST model, Menter further added a limiter for turbulent shear stress that prevents excessive shear stress levels in boundary layers. The $\pmb { k } \mathrm { - } \omega$ SST model essentially contains two sets of coefficients, one for the $k { \mathrm { - } } \omega$ part and one for the (transformed) $\pmb { k } \widetilde { \varrho }$ part. The blending is achieved by functions of two different turbulent Reynolds numbers, which, in turn, require the normal distance from the wall. Abaqus/CFD automatically computes this distance function. The transport equations for the model are

$$
\begin{array}{l} \frac {d}{d t} \int_ {V} \rho k d V + \int_ {S} \rho k (\mathbf {v} - \mathbf {v} _ {m}) \cdot \mathbf {n} d S = \int_ {S} (\mu + \sigma_ {k} \mu_ {T}) \nabla k \cdot \mathbf {n} d S \\ + \int_ {V} \left(\tau_ {i j} S _ {i j} + G _ {b}\right) d V - \int_ {V} \beta^ {*} \rho \omega k d V, \\ \end{array}
$$

$$
\frac {d}{d t} \int_ {V} \rho \omega d V + \int_ {S} \rho \omega (\mathbf {v} - \mathbf {v} _ {m}) \cdot \mathbf {n} d S = \int_ {S} (\mu + \sigma_ {\omega} \mu_ {T}) \nabla \omega \cdot \mathbf {n} d S
$$

$$
+ \int_ {V} \frac {\rho \gamma}{\mu_ {T}} (\tau_ {i j} S _ {i j} + C _ {\omega 3} G _ {b}) d V - \int_ {V} \beta \rho \omega^ {2} d V + \int_ {V} 2 (1 - F _ {1}) \rho \sigma_ {\omega 2} \frac {1}{\omega} \nabla k \cdot \nabla \omega d V.
$$

The $\pmb { k } \mathrm { - } \omega$ SST model coefficients are computed by blending the original $\pmb { k } \mathrm { - } \omega$ model by Wilcox (1988) and the standard $\pmb { k } { - } \varepsilon$ model coefficients:

$$
\phi = F _ {1} \phi_ {1} + (1 - F _ {1}) \phi_ {2}.
$$

The blending function is defined as

$$
F _ {1} = \tanh (\arg_ {1} ^ {4}),
$$

<!-- source-page: 396 -->

$$
\arg_ {1} = \min \biggl (\max \biggl (\frac {\sqrt {k}}{0 . 0 9 \omega y}, \frac {5 0 0 \nu}{y ^ {2} \omega} \biggr), \frac {4 \rho \sigma_ {\omega 2} k}{C D _ {k \omega} y ^ {2}} \biggr),
$$

$$
C D _ {k \omega} = \max \biggl (2 \rho \sigma_ {\omega 2} \frac {1}{\omega} \frac {\partial k}{\partial x _ {j}} \frac {\partial \omega}{\partial x _ {j}}, 1 0 ^ {- 1 0} \biggr).
$$

In the incompressible formulation the Reynolds stresses are closed as

$$
\tau_ {i j} = 2 \mu_ {T} S _ {i j},
$$

$$
S _ {i j} = \frac {\partial u _ {i}}{\partial x _ {j}} + \frac {\partial u _ {j}}{\partial x _ {i}},
$$

with the turbulent eddy viscosity computed using the shear-stress transport approach:

$$
\mu_ {T} = \frac {\rho a _ {1} k}{\max (a _ {1} \omega , S F _ {2})},
$$

$$
S = \sqrt {2 S _ {i j} S _ {i j}},
$$

$$
F _ {2} = \tanh (\arg_ {2} ^ {2}),
$$

$$
\arg_ {2} = \max \biggl (\frac {2 \sqrt {k}}{0 . 0 9 \omega y}, \frac {5 0 0 \nu}{y ^ {2} \omega} \biggr).
$$

The $\pmb { k } \mathrm { - } \omega$ SST model coefficients are shown in Table 6.6.2–4 and Table 6.6.2–5. In addition, a turbulent Prandtl number, $P r _ { t }$ , can be specified.

Table 6.6.2–4 $k { - } \omega$ SST model coefficients (Wilcox $k { - } \omega$ model). 
<table><tr><td> $\sigma_{k1}$ </td><td> $\sigma_{\omega 1}$ </td><td> $\beta_1$ </td><td> $\gamma_1$ </td><td> $\kappa$ </td></tr><tr><td>0.5</td><td>0.5</td><td>0.075</td><td> $\beta_1/\beta^* - \sigma_{\omega 1}\kappa^2/\sqrt{\beta^*}$ </td><td>0.41</td></tr></table>

<!-- source-page: 397 -->

Table 6.6.2–5 k– SST model coefficients (standard k– ). 

<table><tr><td> $\sigma_{k2}$ </td><td> $\sigma_{\omega2}$ </td><td> $\beta_2$ </td><td> $\gamma_2$ </td><td> $\beta^*$ </td></tr><tr><td>1.0</td><td>0.856</td><td>0.0828</td><td> $\beta_2/\beta^* - \sigma_{\omega2}\kappa^2/\sqrt{\beta^*}$ </td><td>0.09</td></tr></table>

Input File Usage: Use both of the following options:

\*CFD

\*TURBULENCE MODEL, TYPE=KOMEGA SST

Abaqus/CAE Usage: The k– SST turbulence model is not supported in Abaqus/CAE.

# Wall functions

The k– SST model can be integrated throughout the inner-layer of the turbulent boundary layer. However, this requires fine near-wall resolutions. The same model can also be applied using a wall-function approach, in which the near-wall element centroid is located in the logarithmic (fully turbulent) part of the boundary layer. This greatly reduces the mesh refinement requirements and, hence, computational costs, albeit with a potential loss of accuracy.

The conventional wall-function approach is based on the law-of-the-wall, which is a semi-empirical universal velocity profile obtained in equilibrium wall-bounded flows when the flow velocity, , and wall normal distance, , are normalized with the kinematic viscosity, , and friction velocity $v _ { \tau }$ (known as viscous units or wall units):

$$
V _ {(y ^ {+})} ^ {+} = \left\{ \begin{array}{l l} y ^ {+} & \mathrm{if} y ^ {+} \leq y _ {c} ^ {+} \\ \frac {1}{\kappa} \ln (E y ^ {+}) & \mathrm{if} y ^ {+} > y _ {c} ^ {+}, \end{array} \right.
$$

$$
V ^ {+} = \frac {V}{v _ {\tau}},
$$

$$
y ^ {+} = \frac {y v _ {\tau}}{\nu},
$$

$$
v _ {\tau} = \sqrt {\frac {\tau_ {w a l l}}{\rho}},
$$

where is the kinematic viscosity, $\rho$ is the density, $\tau _ { w a l l }$ is the shear stress at the wall, $y _ { c } ^ { + }$ is the intersection point of the linear and logarithmic velocity profile, $\kappa = 0 . 4 1$ is the Von Karman constant, and is the law of the wall constant.

In equilibrium conditions the wall-shear stress is approximately equal to the Reynolds stress:

<!-- source-page: 398 -->

$$
\begin{array}{l} \tau_ {w a l l} = \rho v _ {\tau} ^ {2} \approx \mu_ {T} \frac {\partial v}{\partial y} = \rho \beta^ {*} \frac {k ^ {2}}{\varepsilon} \frac {\partial v}{\partial y} \\ = \rho^ {2} \beta^ {*} \frac {k ^ {2}}{\tau_ {w a l l} \frac {\partial v}{\partial y}} \frac {\partial v}{\partial y} = \frac {\rho^ {2} \beta^ {*} k ^ {2}}{\tau_ {w a l l}}, \\ \tau_ {w a l l} ^ {2} = (\rho v _ {\tau} ^ {2}) ^ {2} = \rho^ {2} \beta^ {*} k ^ {2}. \\ \end{array}
$$

Thus, the friction velocity in equilibrium conditions is shown to be

$$
v _ {\tau} = \beta^ {* 1 / 4} k ^ {1 / 2}.
$$

The wall-shear stress can be linearized with the wall function as

$$
\tau_ {w a l l} = \rho v _ {\tau} \frac {V}{V ^ {+}},
$$

such that the production of the turbulent kinetic energy

$$
P _ {k} = \tau_ {i j} S _ {i j} = 2 \mu_ {T} S _ {i j} S _ {i j},
$$

can be simplified using the wall-function relations

$$
P _ {k} = \tau_ {w a l l} \frac {\partial v}{\partial y} = \frac {\tau_ {w a l l} ^ {2}}{\mu} \frac {d v ^ {+}}{d y ^ {+}},
$$

where the gradient of the law-of-the-wall is evaluated in the logarithmic region. Therefore, the final equation is

$$
P _ {k} = \frac {1}{\mu \kappa y ^ {+}} \left(\rho v _ {\tau} \frac {V}{V ^ {+}}\right) ^ {2}.
$$

The specific energy dissipation rate is defined as

$$
\omega = \frac {\varepsilon}{\beta^ {*} \kappa}.
$$

Using the wall-function relations and equilibrium assumption, the value of can be found as

$$
\rho \varepsilon = P _ {k} = \tau_ {w a l l} \frac {\partial v}{\partial y} = \rho^ {2} v _ {\tau} ^ {4} \frac {1}{\mu} \frac {d v ^ {+}}{d y ^ {+}}.
$$

Thus, the value of is obtained by evaluating the gradient of the law-of-the-wall in the logarithmic region:

<!-- source-page: 399 -->

$$
\omega = \frac {v _ {\tau}}{\beta^ {* 1 / 2} \kappa y}.
$$

The previous wall-function relationships are valid only in the logarithmic layer, where the near-wall resolution is $y ^ { + } > 3 0$ . However, since the $\pmb { k } \mathrm { - } \omega$ model can be implemented throughout the inner layer, it is important to properly handle the fine resolution region $y ^ { + } < 3 0$ . In Abaqus/CFD a hybrid wall-function approach is employed in which the near-wall conditions adjust to their appropriate asymptotic limits for either very fine meshes $( y ^ { + } \sim 1 )$ or for wall-function-type meshes $( y ^ { + } \geq 3 0 )$ . Appropriate blending functions are used such that for intermediate meshes (the near-wall element centroid is located in the buffer region) the accuracy is not impaired significantly. For completeness, we first state the appropriate relationships for the viscous sublayer:

$$
v ^ {+} = y ^ {+},
$$

$$
\frac {V}{v ^ {*}} = \frac {y v ^ {*}}{\nu},
$$

$$
v ^ {*} = \sqrt {\frac {\nu V}{y}},
$$

$$
P _ {k} = 0,
$$

$$
\omega = \frac {6 \nu}{\beta_ {1} y ^ {2}}.
$$

Having defined the viscous-sublayer and the fully turbulent relations, the blending is accomplished in the following form:

$$
v _ {\tau} = \sqrt {g \frac {\mu V}{\rho y} + (1 - g) \beta^ {* 1 / 4} k ^ {1 / 2}},
$$

$$
P _ {k} = (1 - g) \frac {1}{\mu \kappa y ^ {+}} \tau_ {w a l l} ^ {2},
$$

$$
\omega = g \frac {6 \nu}{\beta_ {1} y ^ {2}} + (1 - g) \frac {v _ {\tau}}{\sqrt {\beta^ {*}} \kappa y} + g (1 - g) \bigg (\frac {6 \nu}{\beta_ {1} y ^ {2}} + \frac {v _ {\tau}}{\sqrt {\beta^ {*}} \kappa y} \bigg),
$$

where the blending function is defined as

$$
g = e ^ {- R _ {e y} / y _ {c} ^ {+}},
$$

<!-- source-page: 400 -->

$$
R _ {e y} = \frac {\sqrt {k} y}{\nu}.
$$

Here, $R _ { e y }$ is a local Reynolds number, and $y _ { c } ^ { + }$ is the intersection point of the linear and logarithmic velocity profiles of the law-of-the-wall.

The hybrid wall-function approach is independent of the near-wall resolution; therefore, the law-ofthe-wall $V ^ { + } ( y ^ { + } )$ implemented needs to accurately predict the viscous-sublayer, the logarithmic-layer, and the buffer layer (region that connects the viscous and logarithmic zones) since the cell-center adjacent to the wall can be located anywhere within the inner layer. Therefore, a single smooth correlation that reproduces the entire law-of-the-wall proposed by Reichardt (1951) is implemented:

$$
V ^ {+} (y ^ {+}) = \frac {1}{\kappa} \mathrm{ln} (1 + \kappa y ^ {+}) + C \Big [ 1 - e ^ {- \frac {y ^ {+}}{y _ {c} ^ {+}}} - \frac {y ^ {+}}{y _ {c} ^ {+}} e ^ {- b y ^ {+}} \Big ],
$$

where

$$
C = \frac {1}{\kappa} \mathrm{ln} \left(\frac {E}{\kappa}\right),
$$

$$
b = \frac {1}{2} \left(\frac {y _ {c} ^ {+} \kappa}{C} + \frac {1}{D}\right).
$$

Implementation in the momentum equation

For cases where the mesh resolution is not enough to capture the near-wall gradients, a near-wall model is required to provide the correct wall-shear stress in coarse meshes. The wall shear is obtained from the wall-function approach through an effective edge viscosity:

$$
\mu_ {\mathrm{eff}} = \frac {\tau_ {w a l l} y}{V} = \left(\frac {\rho v _ {\tau} V}{V ^ {+}}\right) \frac {y}{V}.
$$

Energy wall functions

The wall-function approach can be extended to the energy equation by using the temperature law-ofthe-wall which is a semi-empirical universal temperature profile obtained in equilibrium wall-bounded flows, when the temperature, T, and the wall-normal distance, y, are normalized with wall units. The standard temperature wall function is defined as

$$
T ^ {+} (y ^ {+}) = \frac {(T _ {w} - T) \rho c _ {p} u _ {\tau}}{\dot {q}} = \left\{ \begin{array}{l l} P _ {r} y ^ {+} & \mathrm{if} y ^ {+} \leq y _ {T c} ^ {+} \\ P _ {r T} \Big (\frac {1}{\kappa} \ln (E y ^ {+}) + P (P _ {r}, P _ {r T}) \Big) & \mathrm{if} y ^ {+} > y _ {T c} ^ {+}, \end{array} \right.
$$

where $y _ { T c } ^ { + }$ is the intersection point of the viscous-sublayer and logarithmic layer in the temperature wallfunction, $T _ { w }$ is the wall temperature, $P _ { r }$ is the Prandtl number, $P _ { r T }$ is the turbulent Prandtl number,