Turbulence equations (compressible flow)

Here, the two-parameter turbulence models BSL (baseline) and SST (shear stress transport) are treated [51]. The k-\epsilon [37] and the k-\omega -model [90] are special cases of the BSL-model. The two parameters are the turbulent kinetic energy k and the turbulence frequency \omega. The equation for k reads [51]:

\frac{\partial \rho k}{\partial t} + \nabla \cdot (\rho k \boldsymbol{v}) = \nabla \cdot \left [ ( \mu + \sigma_k \mu_t) \nabla k \right ] + t_{ij}^t \frac{\partial v_i}{\partial x_j} - \beta^* \rho \omega k. (713)

Since

t_{ij}^t = 2 \mu_t (d_{ij}-\frac{1}{3} d_{kk} \delta_{ij}) - \frac{2}{3} \rho k \delta_{ij} (714)

and

t_{ij}^l = 2 \mu(d_{ij}-\frac{1}{3} d_{kk} \delta_{ij}), (715)

one can write

t_{ij}^t = \frac{\mu_t}{\mu } t_{ij}^l - \frac{2}{3} \rho k \delta_{ij}. (716)

This leads to:

\frac{\partial \rho k}{\partial t} + \nabla \cdot (\rho k \boldsymbol{v}) = \nabla \cdot \left [ ( \mu + \sigma_k \mu_t) \nabla k \right ] + \frac{\mu_t}{\mu } t_{ij}^l d_{ij} - \frac{2}{3} \rho k d_{kk} - \beta^* \rho \omega k. (717)

Notice that \mu_t t_{ij}^l d_{ij}/\mu can be written as:

\displaystyle \frac{\mu_t}{\mu} t_{ij}^l d_{ij} \displaystyle = \mu_t \left ( 2 \left [ \left ( \frac{\partial v_1}{\partial x_... ... ) ^2 + \left ( \frac{\partial v_3}{\partial x_3} \right ) ^2 \right ] \right . \displaystyle + \left [ \left ( \frac{\partial v_1}{\partial x_2} + \frac{\part... ...ial v_2}{\partial x_3} + \frac{\partial v_3}{\partial x_2} \right ) ^2 \right ]

\displaystyle - \left . \frac{2}{3} (\nabla \cdot \boldsymbol{v})^2 \right ) (718)

In the above conservation equation \sigma_k is a function of the flow characteristics through the blending factor F1 [51] and the dynamic turbulent viscosity satisfies:

\mu_t=\rho \frac{k}{\omega }. (719)

In fact, the BSL model is a linear combination of the k-\omega model and the k-\epsilon-model with coefficients F_1 and 1-F_1, respectively. Near a wall F_1 tends to 1, thus favoring the k-\omega -model which is particularly good in the near-wall region, far away from a wall F_1 tends to zero, leading to a pure k-\epsilon-model.

The blending factor used in iteration (m) is F_1 ^{(m-1)}, i.e. its calculation is mainly based on the results from the previous iteration and satisfies the following equations:

\displaystyle (F_1)_P ^{(m-1)} \displaystyle = \tanh \left [ \left ( {{\text{arg}}_1}_P ^{(m-1)} \right ) ^4 \right ] (720)

\displaystyle {{\text{arg}}_1}_P ^{(m-1)} \displaystyle = \min \left [ \max \left ( \frac{\sqrt{k_P ^{(m-1)} }}{0.09 \ome... ...a_\omega }_2 k_P ^{(m-1)} }{{(\text{CD}_{k \omega}})_P ^{(m-1)} y_P^2} \right ] (721)

\displaystyle {(\text{CD}_{k \omega}})_P ^{(m-1)} \displaystyle = \max \left ( 2 \rho_P ^{(m)} {\sigma_\omega }_2 \frac{1}{\omega... ...1)} } \nabla k_P ^{(m-1)} \cdot \nabla \omega _P ^{(m-1)} , 10^{-20} \right ) , (722)

where y_P is the distance from the element center P to the next solid surface.

Using the blending factor one obtains the correct parameter values, e.g.:

\sigma_k = F_1 {\sigma_k}_1 + (1-F_1) {\sigma_k}_2. (723)

For completeness the parameters are listed here:

\displaystyle \gamma_1 \displaystyle = \frac{\beta_1}{\beta^*} - \frac{{\sigma_\omega}_1 \kappa^2}{\sq... ...2 = \frac{\beta_2}{\beta^*} - \frac{{\sigma_\omega}_2 \kappa^2}{\sqrt{\beta^*}} (724)

\displaystyle {\sigma_k}_1 \displaystyle = 0.5 \;\;\;\;\; {\sigma_\omega }_1 = 0.5 \;\;\;\;\; \beta_1 = 0.075 (725)

\displaystyle {\sigma_k}_2 \displaystyle = 1.0 \;\;\;\;\; {\sigma_\omega }_2 = 0.856 \;\;\;\;\; \beta_2 = 0.0828 (726)

\displaystyle \beta^* \displaystyle = 0.09 \;\;\;\;\; \kappa=0.41 (727)

In the conservation equation for k one very easily identifies the time-dependent, convective, diffusive and body terms. They are treated in a completely analogous way to the energy equation. The convective boundary conditions amount to:

• Inlet: \overrightarrow{k}_f ^{(m)} = 10^{-(2 \rightarrow 5)} \nu_{\infty} (1 \rightarrow 10) \frac{v_{\infty}}{L} is given
• Outlet: \overrightarrow{k}_f ^{(m)} = \overrightarrow{k}_f ^{UD(m)}= k_P ^{(m)}
• Wall and sliding conditions: \dot{m}_f ^{(m)} =0: no contribution

For the convective interpolation of k the modified smart algorithm has not shown any advantages, therefore, the upwind difference scheme is always used.

The diffusion boundary conditions are:

• Inlet:

  (\mu + \sigma_k \mu_t)_f ^{(m-1)} \left [ \frac{\overrightarrow{k} _f ^{(m)} - k_P ^{(m)} }{l_{Pf}}+ (\nabla k)_f ^{(m-1)} (\boldsymbol{n}_f - \boldsymbol{j}_f) \right ]. (728)

  k_f ^{(m)} is known (cf. convective inlet boundary conditions), and

  (\mu + \sigma_k \mu_t)_f ^{(m-1)} \approx \mu_f ^{(m-1)} + (\sigma_k)_P ^{(m-1)} \frac{\overrightarrow{\rho}_f ^{(m)} \overrightarrow{k}_f ^{(m-1)} }{\overrightarrow{\omega}_f ^{(m-1)} } . (729)

  Notice that no facial values are calculated for the turbulent parameters. Therefore, (\sigma_k)_f in the above equation is approximated by (\sigma_k)_P.

• Outlet: \nabla k = 0 in \boldsymbol{ j}-direction; the term reduces to:

  (\mu + \sigma_k \mu_t)_f ^{(m-1)} \left [ (\nabla k)_f ^{(m-1)} \boldsymbol{n}_f \right ]. (730)

• Wall: k=0 is known (therefore also (\mu_t)_f=0), treatment similar as for inlet. One obtains:

  \mu_f ^{(m-1)} \left [ \frac{ - k_P ^{(m)} }{l_{Pf}}+ (\nabla k)_f ^{(m-1)} (\boldsymbol{n}_f - \boldsymbol{j}_f) \right ]. (731)

The source terms are treated in the following way:

• \int_{V}^{} \left[ \frac{\mu_t}{\mu } t_{ij}^l d_{ij} \right ] dV \approx V_P \left [ \frac{\mu_t}{\mu } t_{ij}^l d_{ij} \right ]_P ^{(m-1)}, (732)

  in which the term in brackets depends on \mu_t ^{(m-1)} = \rho ^{(m)} k ^{(m-1)} / \omega ^{(m-1)} and \nabla \boldsymbol{v} ^{(m)} only, cf. Equation (). It can be split in a sum of

  V_P {\mu_t}_P ^{(m-1)} \left [ \frac{t_{ij}^l d_{ij}}{\mu } + \frac{2}{3} d_{kk}^2 \right ]_P ^{(m)}, (733)

  which is positive and corresponds to a source (treated explicitly, i.e. on the right hand side) and

  -V_P {\mu_t}_P ^{(m-1)} \left [ \frac{2}{3} d_{kk}^2 \right ]_P ^{(m)}, (734)

  which is negative (sink) and is treated implicitly (left hand side).

• Depending on whether d_{kk} > 0 or d_{kk}<0 the remaining production term is treated as

  -\frac{2}{3} V_P \rho ^{(m)} k ^{(m)} d_{kk} ^{(m)} (735)

  or

  -\frac{2}{3} V_P \rho ^{(m)} k ^{(m-1)} d_{kk} ^{(m)}, (736)

  respectively.

• The k \omega - term is approximated as follows:

  - \int_{V}^{} \beta^* \rho \omega k dv \approx - \beta^* \rho ^{(m)} \omega ^{(m-1)} k ^{(m)} V_P. (737)

  This term is treated implicitly (K is evaluated at iteration (m) and the term ends up on the left hand side) since it is a negative source.

The equation for the turbulence frequency \omega runs [51]:

\frac{\partial \rho \omega }{\partial t} + \nabla \cdot (\rho \omega \boldsymbol{v}) = \nabla \cdot \left [ ( \mu + \sigma_\omega \mu_t) \nabla \omega \right ] + \frac{\gamma }{\nu_t} t_{ij}^t d_{ij} - \beta \rho \omega ^2 + 2 \rho (1-F_1) {\sigma_\omega} _2 \frac{1}{\omega } \nabla k \cdot \nabla \omega , (738)

which can be rewritten as:

\displaystyle \frac{\partial \rho \omega }{\partial t} + \nabla \cdot (\rho \omega \boldsymbol{v}) \displaystyle = \nabla \cdot \left [ ( \mu + \sigma_\omega \mu_t) \nabla \omega... ...t ] + \frac{\gamma }{\nu_t} \left [ \frac{\mu_t}{\mu } t_{ij}^l d_{ij} \right ]

\displaystyle - \frac{2}{3} \gamma \rho \omega d_{kk} - \beta \rho \omega ^2 + 2 \rho (1-F_1) {\sigma_\omega} _2 \frac{1}{\omega } \nabla k \cdot \nabla \omega . (739)

One easily recognizes the time-dependent, convective, diffusive and source terms. The convective boundary conditions amount to:

• Inlet: \overrightarrow{\omega }_f ^{(m)} = (1 \rightarrow 10) \frac{v_\infty}{L} is given
• Outlet: \overrightarrow{\omega }_f ^{(m)} = \overrightarrow{\omega }_f ^{UD(m)}
• Wall and sliding conditions: \dot{m}_f ^{(m)} =0: no contribution

For the convective interpolation of \omega the modified smart algorithm has not shown any advantages, therefore, the upwind difference scheme is always used.

The diffusion boundary conditions are:

• Inlet:

  (\mu + \sigma_\omega \mu_t)_f ^{(m-1)} \left [ \frac{\omega _f ^{(m)} - \omega _P ^{(m)} }{l_{Pf}}+ (\nabla \omega )_f ^{(m-1)} (\boldsymbol{n}_f - \boldsymbol{j}_f) \right ]. (740)

  \omega _f ^{(m)} is known (cf. convective inlet boundary conditions), and

  (\mu + \sigma_\omega \mu_t)_f ^{(m-1)} \approx \mu_f ^{(m-1)} + (\sigma_\omega)_P ^{(m-1)} \frac{\overrightarrow{\rho} ^{(m)} \overrightarrow{k} ^{(m)} }{\overrightarrow{\omega} ^{(m-1)} } . (741)

  Notice that no facial values are calculated for the turbulent parameters. Therefore, (\sigma_\omega )_f in the above equation is approximated by (\sigma_\omega )_P.

• Outlet: \nabla \omega = 0 in \boldsymbol{ j}-direction; the term reduces to:

  (\mu + \sigma_\omega \mu_t)_f ^{(m-1)} \left [ (\nabla \omega )_f ^{(m-1)} \boldsymbol{n}_f \right ]. (742)

• Wall: \omega is known ([51]):

  \omega = 10 \frac{6 \mu ^{(m-1)} }{\beta_1 (\Delta y_1)^2 \rho ^{(m)} }, (743)

  where \beta_1=0.075 and \Delta y_1 is the distance to the next point away from the wall; same treatment as for inlet (notice, however, hat \mu_t=0 at the wall).

The source terms are treated as follows:

• The first part of the turbulent production term can be written as:

  \int_{v}^{} \frac{\gamma }{\nu_t} \left [ \frac{\mu_t}{\mu } t_{ij}^l d_{ij} \right ] dv \approx \gamma _P ^{(m-1)} \rho_P ^{(m)} f(\nabla \boldsymbol{v}_P ^{(m)} V_P (744)

  where f(\nabla \boldsymbol{v}_P ^{(m)}) is the term in the outer brackets on the right hand side of Equation (). Part of this term goes to the right hand side (sources) and part to the left hand side (sinks) as discussed extensively for the energy equation.

• The second part of the turbulent production is approximated by

  -\frac{2}{3} \gamma_P ^{(m-1)} \rho_P ^{(m)} \omega_P ^{(m)} (d_{kk})_P ^{(m)} V_P (745)

  if (d_{kk})_P^{(m)} > 0 (sink) and by

  -\frac{2}{3} \gamma_P ^{(m-1)} \rho_P ^{(m)} \omega_P ^{(m-1)} (d_{kk})_P ^{(m)} V_P (746)

  if (d_{kk})_P^{(m)} < 0 (source).

• The \omega^2-term is a sink and is treated implicitly:

  - \beta_P ^{(m-1)} \rho_P ^{(m)} \omega_P ^{(m-1)} \omega_P ^{(m)} (747)

• Finally, the cross-diffusion term is approximated by:

  2 \rho_P ^{(m)} (1-F_1 ^{(m-1)}) {\sigma_\omega }_2 \frac{\nabla k_P ^{(m)} \nabla \omega_P ^{(m-1)} }{\omega_P ^{(m-1)} } V_P (748)

  if the term is positive (source), and by:

  2 \rho_P ^{(m)} (1-F_1 ^{(m-1)}) {\sigma_\omega }_2 \frac{\nabla k_P ^{(m)} \nabla \omega_P ^{(m-1)} }{(\omega_P ^{(m-1)})^2 } \omega_P ^{(m)} V_P, (749)

  if the term is negative (sink).