3.2.7. Wall Boundary Conditions for Turbulence Models

3.2.7.1. Resolution of Boundary Layer; Momentum

The wall velocity boundary condition is the typical no-slip boundary; a specified value is expected.

3.2.7.2. Resolution of Boundary Layer; Turbulence Quantities

The resolution of the boundary layer is expected when the low Reynolds number or $v^2$ - $f$ model is in use.

For the $v^2$ - $f$ model, the wall turbulent kinetic energy and normal fluctuating stress component are each zero while the dissipation rate is given by

(3.193) $\epsilon_w = 2\nu {\partial k^{1/2} \over \partial x_j }^2.$

For the low Reynolds number, the wall turbulent kinetic energy is again zero while the dissipation rate, here considered to be the isotropic dissipation rate, is given as zero.

3.2.7.3. Resolution of Boundary Layer; Enthalpy

The wall value of enthalpy is computed based on the specified temperature and either reference or local mass fractions. In the case of a heat flux boundary condition, the wall node value is computed based on the control volume balance.

3.2.7.4. Wall Functions for Turbulent Flow Boundary Conditions

Resolution of the near-wall turbulent boundary layer can require extensive mesh points. Adjacent to the wall exists an extremely thin viscous sublayer where these forces dominate and are relatively insensitive to free stream parameters. Following the viscous sublayer is a buffer layer, the so-called “log-layer” and, ultimately, the turbulent core. The Van Driest hypothesis of turbulent flow near solid boundaries can be used to derive the appropriate form of this log-law zone. In general, the use of wall functions eliminates the need to resolve the near wall layers by prescribing the wall shear stress and resulting force based on the law of the wall (Launder and Spalding [35]).

The primary assumptions of the law of the wall are

Local equilibrium of turbulent kinetic energy production and dissipation,

Constant shear stress within the log-law region,

Couette flow (pure shear flow).

3.2.7.5. Wall Functions; Momentum

The wall shear stress enters the discretization of the momentum equations by the term

(3.194) $\int \tau_{ij} n_j dS = -{F_w}_i .$

Wall functions are used to prescribe the value of the wall shear stress rather than resolving the boundary layer within the near-wall domain. The fundamental momentum law of the wall formulation, assuming fully-developed turbulent flow near a no-slip wall, can be written as (Launder and Spalding [35])

(3.195) $u^+ = {u_{\|} \over u_{\tau}} = { 1 \over \kappa } \ln \left(Ey^+\right) ,$

where $u^+$ is defined by the the near-wall parallel velocity, $u_{\|}$ , normalized by the wall friction velocity, $u_{\tau}$ . The wall friction velocity is related to the turbulent kinetic energy by

(3.196) $u_{\tau} = C_\mu^{1/4} k^{1/2}.$

by assuming that the production and dissipation of turbulence is in local equilibrium. Moreover, $y^+$ is defined as the normalized perpendicular distance from the point in question to the wall,

(3.197) $y^+ = {{ \rho Y_p} \over {\mu }}\left(\tau_w \over \rho \right)^{1/2} = {{ \rho Y_p u_{\tau}} \over {\mu }}$

The classical law of the wall is as follows:

(3.198) $u^+ = \frac{1}{\kappa} \ln(y^+) + C$

where $\kappa$ is the von Karman constant and $C$ is the dimensionless integration constant that varies based on authorship and surface roughness. The above expression can be re-written as

(3.199) $u^+ = \frac{1}{\kappa} \ln(y^+) + \frac{1}{\kappa} \ln(\exp(\kappa C))$

or

(3.200) $u^+ = \frac{1}{\kappa} \left(\ln(y^+) + \ln(\exp(\kappa C))\right)$

(3.201) $u^+ = \frac{1}{\kappa} \ln(E y^+)$

where $E$ is referred to in the text as the dimensionless wall roughness parameter and is described by

(3.202) $E = \exp(\kappa C)$

In Fuego, $\kappa$ is set to the value of 0.42 while the value of $E$ is set to 9.8 for smooth walls. White [36] suggests values of $\kappa=0.41$ and $E=7.768$ . The viscous sublayer is assumed to extend to a value of $y^+$ = 11.63.

The wall shear stress, $\tau_w$ , can be expressed as

(3.203) $\tau_w = \rho u_\tau^2 = \rho u_\tau {{u_\|} \over {u^+}} = { {\rho \kappa u_{\tau}} \over {\ln \left(Ey^+\right) } }u_\| = \lambda_w u_\| ,$

where $\lambda_w$ is simply the grouping of the factors from the law of the wall. For values of $y^+$ less than 11.63, the wall shear stress is given by

(3.204) $\tau_w = \mu {u_\| \over Y_p} .$

The force imparted by the wall, for the $i_{th}$ component of velocity, can be written as

(3.205) $F_{w,i}= -\lambda_w A_w u_{i\|} ,$

where $A_w$ is the total area over which the shear stress acts.

The use of a general, non-orthogonal mesh adds a slight complexity to specifying the force imparted on the fluid by the wall. As shown in (3.205), the velocity component parallel to the wall must be determined. Use of the unit normal vector, $n_j$ , provides an easy way to determine the parallel velocity component by the following standard vector projection,

(3.206) $\Pi_{ij} = \left [ \delta_{ij} - n_i n_j \right].$

Carrying out the projection of a general velocity, which is not necessarily parallel to the wall, yields the velocity vector parallel to the wall,

(3.207) $u_{i\|} = \Pi_{ij} {u}_j = u_i\left(1-{n_i}^2\right) -\sum_{j=1;j\neq j}^{n} u_j n_i n_j.$

Note that the component that acts on the particular $i^{th}$ component of velocity,

(3.208) $-\lambda_w A_w \left(1-n_i n_i\right) u_{i\|} ,$

provides a form that can be potentially treated implicitly; i.e., in a way to augment the diagonal dominance of the central coefficient of the $i^{th}$ component of velocity. The use of residual form adds a slight complexity to this implicit formulation only in that appropriate right-hand-side source terms must be added.

3.2.7.6. Wall Functions; Turbulent Kinetic Energy

The near wall turbulent kinetic energy can be obtained by two different procedures. The most common approach is to solve a transport equation for the near wall value of turbulent kinetic energy with a modified production and dissipation term on the right hand side of the turbulent kinetic energy equation, (3.118). As will be shown below, the form of the near wall production and dissipation term are determined based on equilibrium arguments, i.e., $P_k = \rho\epsilon$ .

Another common approach is to assign the value of turbulent kinetic energy that strictly results in the equality $P_k = \rho\epsilon$ . In this formulation, it is assumed that the convection and diffusive flux is zero across the control volume.

Both procedures, which formally do not address the role of buoyancy production, begin with the determination of the near wall value of the production of turbulent kinetic energy. The turbulent kinetic energy production term is consistent with the law of the wall formulation and can be expressed as

(3.209) ${P_k}_w = \tau_w {{\partial u_{\|}} \over {\partial y}}.$

The parallel velocity, $u_{\|}$ , can be related to the wall shear stress by

(3.210) $\tau_w {u^+ \over y^+ } = \mu {u_{\|} \over Y_p }.$

Taking the derivative of both sides of (3.210), and substituting this relationship into (3.209) yields,

(3.211) ${P_k}_w = {\tau_w^2 \over \mu} {{\partial u^+} \over {\partial y^+}}.$

Applying the derivative of the law of the wall formulation, (3.195), provides the functional form of ${\partial u^+ / \partial y^+}$ ,

(3.212) ${\partial u^+ \over \partial y^+} = {\partial \over \partial y^+} \left[{ 1 \over \kappa } \ln \left(Ey^+\right) \right] = {1 \over \kappa y^+}.$

Substituting (3.212) within (3.211) yields a commonly used form of the near wall production term,

(3.213) ${P_k}_w = {{\tau_w}^2 \over \rho\kappa u_{\tau} Y_p}.$

Assuming local equilibrium, $P_k = \rho\epsilon$ , and using (3.213) and (3.196) provides the form of the near wall turbulence dissipation,

(3.214) $\epsilon = {{u_\tau^3} \over {\kappa Y_p}} = { {C_\mu^{3/4} k^{3/2}} \over {\kappa Y_p} } ,$

while the form of the wall shear stress is given by,

(3.215) $\tau_w = \rho C_\mu^{1/2} k$

Under the above assumptions, the near wall value for turbulent kinetic energy, in the absence of convection, diffusion, or accumulation is given by,

(3.216) $k = {{u_\tau^2} \over {C_\mu^{1/2}}}.$

If the second method (Dirichlet condition on near wall turbulent kinetic energy) is to be used, the value of the wall friction velocity, $u_\tau$ , can be obtained in an iterative manner (Sondak and Pletcher [37]) by use of (3.195). This method has been used and shown to be satisfactory (Elkaim [38]) and strictly enforces the assumptions of the law of the wall that have already been outlined.

In the method that elects to solve a near wall turbulent kinetic energy transport equation, the production and dissipation terms in the turbulent kinetic energy transport equation are [potentially] given by (3.213) and

(3.217) $-\rho\epsilon = -\rho{ {C_\mu^{3/4} k^{3/2}} \over {\kappa Y_p} } ,$

Unfortunately, there does not seem to be one universal description of the near wall turbulent kinetic energy production term and dissipation term, (3.213) and (3.217), respectively. For example, in the law of the wall formulation, given by Launder and Spalding [39], the near wall production term is given by,

(3.218) ${P_k}_w = \tau_w {{u_{\|}} \over {y_p}}.$

In this formulation, the wall shear stress is given by the law of the wall formulation, (3.203), providing the value of $y^+$ is greater than 11.63 (otherwise, it is given by the laminar shear stress, (3.204)). The dissipation term, $-\rho \epsilon$ is given by

(3.219) $-\rho\epsilon = -\rho{ {C_\mu^{3/4} k^{3/2}} \over {\kappa Y_p} } \ln{E y^+} .$

Note that in the absence of convection, diffusion or accumulation, the above two forms of the near wall production and dissipation source terms revert to (3.216). Therefore, if the modeled flow is consistent with the law of the wall formulations, all methods should yield similar limiting behavior. Under conditions of non equilibrium, i.e., a separated flow, or values of $y^+$ within the viscous sublayer, some models may perform better. However, it is important to note that if the flow to be simulated includes separation and reattachment, or the computation mesh is such that $y^+$ is within the viscous sublayer, the law of the wall formulation can provide nonsensical results.

In Fuego, there are currently two general supported methods from which to choose when applying the near wall turbulent kinetic energy boundary condition. The first method, which can be activated by the command line omit near wall turbulent ke transport equation, is the form of (3.216) that enforces a Dirichlet condition. The second method is to solve a full control volume balance for the near wall turbulent kinetic, with convection and diffusion terms, with a modified production and dissipation term given by either

Equations (3.213) and (3.214).

Equations (3.218) and (3.219)

The use of (3.213) and (3.214) can be activated by the command line (within the wall boundary condition block) use equilibrium production model which is based on the ability to express the wall shear stress consistent with the assumptions of full equilibrium between production and dissipation, (3.215). In all cases that do not set a Dirichlet condition for the turbulent kinetic energy, the assembled buoyancy source terms are not removed.

3.2.7.7. $k - \omega$ $SST$ Wall Functions; Turbulent Kinetic Energy

When a Dirichlet condition is not set for turbulent kinetic energy, the approach in modifying the near wall production and dissipation terms is followed.

In this approach, the equation for $k$ is solved near the wall to remove the assumptions of log layer flow one level. However, we invoke the log layer assumption to write,

(3.220) $P_{k} = {\tau_{w}^{2} \over \rho \kappa u_{\tau} Y_{p}}.$

Balancing production and dissipation in the $k-\omega$ model allows us to write,

(3.221) $P_{k} = \rho {u_{\tau}^{3} \over \kappa Y_{p}} = \rho {(\beta ')^{3/4} k^{3/2} \over \kappa Y_{p}}.$

The dissipation rate is also modified accordingly such that the production equality with dissipation is retained. An alternative method is to use the approximation of of Launder and Spaulding which prescribes production as,

(3.222) $P_{k} = \tau_{w} {u_{||} \over Y_{p}}.$

In practice, this formulation seems to be less stable since the production and dissipation terms are now in near-equilibrium.

3.2.7.8. Wall Functions; Turbulence Dissipation Transport

Consistently within the literature, the near wall turbulence dissipation is assigned the Dirichlet value given by (3.214). Frequently, this expression is lagged by one sub-iteration in an effort to maintain consistency between the Dirichlet wall condition and the freezing of the $\epsilon / k$ ratio of the turbulence dissipation equation, (3.119).

3.2.7.9. Wall Functions; Turbulent Frequency Transport

3.2.7.9.1. Low Reynolds Number Treatment

The low Reynolds approach for $k-\omega$ uses a sequence of Dirichlet conditions similar to what is used for $k-\epsilon$ . However, unlike the latter, $k-\omega$ requires no extra damping terms near the wall. When the wall is resolved, exact Dirichlet conditions are known for both the velocity and $k$ :

(3.223) $\vec{u} = 0, \qquad k = 0.$

A Dirichlet condition is also used on $\omega$ . While the $k-\epsilon$ model is rendered less stable because $k$ appears in this boundary condition, the $\omega$ equation depends only on the near-wall grid spacing. The boundary condition is

(3.224) $\omega = {6 \nu \over \beta y^{2}},$

which is valid for $y^{+} < 3$ . Above, $\beta$ depends on the model type. If SST is in use, $\beta = \beta_1$ while if the Wilcox model is in use, $\beta = \beta_0$ .

3.2.7.9.2. High Reynolds Number Treatment

The high Reynolds approach is also quite similar to the $k-\epsilon$ model except $\omega$ is handled differently.

3.2.7.9.3. Automatic Wall Functions

Because $\omega$ has analytic solutions in both the log layer and viscous sub-layer, an automatic treatment is developed that blends those two solutions to provide Dirichlet conditions for all $y$ . Let $\omega_{h}$ be the high Reynolds number formulation and $\omega_{l}$ be the low Reynolds version. Then the Dirichlet condition on $\omega$ is

(3.225) $\omega = \omega_{l} \sqrt{ 1 + \left({\omega_{h} \over \omega_{l} }\right)^{2}}.$

However, $u_{\tau}$ for the high Reynolds $\omega$ value is computed based on the parallel velocity: The velocity equation is augmented by a traction force based on the friction velocity $u_{\tau}$ . This quantity may be solved for iteratively using the law of the wall. A Dirichlet condition is also used for $k$ , assuming it is in the log region, which is similar to the $k-\epsilon$ model:

(3.226) $k = {u_{\tau}^{2} \over \sqrt{\beta^*}}.$

In the case of $\omega$ , an analytic expression is known in the log layer:

(3.227) $\omega = {u_{\tau} \over \sqrt{\beta^*} \kappa y},$

which is independent of $k$ . Note that some implementations use a predefined constant instead of $\sqrt{\beta'}$ , although the standard values are consistent with these expressions. Because all these expressions require $y$ to be in the log layer, they should absolutely not be used unless it can be guaranteed that $y^{+} > 10$ , and $y^{+} > 25$ is preferable.

(3.228) $u_{\tau} = \sqrt{ \nu \bigg| {u_{||} \over y} \bigg|}.$

The automatic wall function approach is obtained by removing the omit near wall turbulent ke equation line command and activating either the SST or KW turbulence models.

3.2.7.10. Wall Functions; Enthalpy Transport

For non-adiabatic boundaries, heat loss to the wall must be considered. The use of the Reynolds analogy provides a functional form of the energy transport similar to the that of the logarithmic law-of-the-wall momentum formulation. The thermal boundary layer is modeled either as a linear profile ( $y^+ < 11.63$ ) where the thermal boundary layer is dominated by conduction or a logarithmic profile where the effects of turbulence dominate over thermal conduction, Versteeg and Malalasekera [40].

The law-of-the-wall used in Fuego has the following form,

(3.229) $q_w = {{\rho C_p \left(T_w - T_p\right) u_{\tau}} \over {T^+}} ,$

where

(3.230) $T^+ = \sigma_T \left[u^+ + P\right].$

The role of $T^+$ is to account for the fact that the thickness of the thermal conduction layer is [practically] of a different size than that of the viscous sublayer (momentum).

In the above equation, $P$ is the universal “P function” (Jayatilleke [41]) and can be expressed as a function of the molecular and turbulent Prandtl number,

(3.231) $P = {9.24{\left[ \left({{\sigma} \over {\sigma_T}} \right)^{0.75} - 1 \right] }} \left(1+0.28exp\left[-0.007 {\sigma \over \sigma_T} \right]\right) ,$

where $\sigma_T$ and $\sigma$ represent the turbulent and molecular Prandtl number, respectively.

Therefore, it is seen that the so-called “P function” is the parameter that functionally changes the thickness of the thermal conduction layer from that of the viscous sublayer. For example, if one were to model a high-Prandtl number fluid such as common vegetable oil, one would note that the thickness of the viscous sublayer is far greater than that of the thermal sublayer. However, for low-Prandtl number fluids, the opposite is true. The subsequent value of $T^+$ ensures this functionality.

In the case of a user defined heat flux at a wall boundary condition, the full quantity is assembled as a right-hand-side source term. As a post processing step, (3.229) (in temperature form) is rearranged to provide the wall temperature. In practice, the heat flux boundary condition block is to be defined on an already defined wall boundary condition block (without temperature specification). In this manner, multiple boundary conditions are “painted” on a particular sideset.

3.2.7.11. Wall Functions; Scalar Transport

Wall functions for use in a convective diffusive problem, e.g., diffusional transport of fuel (through multicomponent evaporation) from a jet fuel pool, are not currently supported.