k-ε-φ-f

Properties

• Three transport-equation linear-eddy-viscosity RANS turbulence closure model with an elliptic relaxation function:
  • Turbulent kinetic energy, $k$,
  • Turbulent kinetic energy dissipation rate, $\epsilon$,
  • Normalised wall-normal fluctuating velocity scale, $\phi$,

whose role is to predict *viscous damping*,

• Elliptic relaxation factor, $f$, whose role is to predict effects of

pressure-rate-of-strain, i.e. to account for the non-local, wall-echo and
blocking effects of solid walls.

• Based on: Laurence, Uribe, and Utyuzhnikov (2005) 38
• The name of the original variable replacing $\texttt{v2}$ is $\texttt{phi}$ in the original research paper (38, Eq. 14). However, the name $\texttt{phi}$ preexisted in OpenFOAM; therefore, this name was replaced by $\texttt{phit}$.

Model equations

The turbulent kinetic energy transport equation, $k$, 38, Eq. 3]:

\partial_t k + u_i \partial_{x_i} k = P - \epsilon + \partial_{x_j} \left{ \left( \nu + \frac{\nu_t}{\sigma_k} \right) \partial_{x_j} k \right}

Where:

$k$ : Turbulent kinetic energy [$\text{m}^2 \text{s}^{-3}$]

$P$ : Turbulent kinetic energy production rate [$\text{m}^2 \text{s}^{-3}$]

$\epsilon$ : Turbulent kinetic energy dissipation rate [$\text{m}^2 \text{s}^{-3}$]

$\nu$ : Kinematic viscosity of fluid [$\text{m}^2 \text{s}^{-1}$]

$\nu_t$ : Turbulent viscosity [$\text{m}^2 \text{s}^{-1}$]

$\sigma_k$ : Turbulent Prandtl number for $k$ [-]

The turbulent kinetic energy dissipation rate transport equation, $\epsilon$, 38, Eq. 4]:

\partial_t \epsilon + u_i \partial_{x_i} \epsilon = \frac{C_{\epsilon_1} P}{T} - \frac{C_{\epsilon_2} \epsilon}{T} + \partial_{x_j} \left{ \left( \nu + \frac{\nu_t}{\sigma_\epsilon} \right) \partial_{x_j} \epsilon \right}

Where:

$T$ : Modelled turbulent time scale [$s$]

$C_{\epsilon_*}$ : Empirical model constants [-]

$\sigma_\epsilon$ : Turbulent Prandtl number for $\epsilon$ [-]

The normalised wall-normal fluctuating velocity scale transport equation, $\phi$, 38, Eq. 17]:

\partial_t \phi + u_i \partial_{x_i} \phi = f - P \frac{\phi}{k} + \frac{2 \nu_t}{k \sigma_k} \partial_{x_j} \phi \partial_{x_j} k + \partial_{x_j} \left{ \left( \frac{\nu_t}{\sigma_k} \right) \partial_{x_j} \phi \right}

where:

$f$ : Elliptic relaxation factor [$s^{-1}$]

$\phi$ : Normalised wall-normal fluctuating velocity scale [-]

$\sigma_\phi$ : Turbulent Prandtl number for $\phi$ [-]

The elliptic relaxation factor equation, $f$, 38, Eq. 18]:

$$L^2 \partial^2_{x_j} f - f = \frac{1}{T} (C_1 - 1) \left[ \phi - \frac{2}{3} \right] - C_2 \frac{P}{k} - 2 \frac{\nu }{k} \partial_{x_j} \phi \partial_{x_j} k - \nu \partial^2_{x_j} \phi$$

Where:

$L$ : Modelled turbulent length scale [$m$]

$C_{f_*}$ : Empirical model constants [-]

The turbulent time scale equation, $T$, 38, Eq. 7]:

$$T = \max \left[ \frac{k}{\epsilon}, 6.0 \sqrt{\frac{\nu}{\epsilon}} \right]$$

The turbulent length scale equation, $L$, 38, Eq. 7]:

$$L = C_L \max \left[ \frac{k^{1.5}}{\epsilon}, C_\eta \left( \frac{\nu^3}{\epsilon} \right)^{0.25} \right]$$

Where:

$C_L$ : Empirical model constant [-]

$C_\eta$ : Empirical model constant [-]

The turbulent viscosity equation, $\nu_t$, 38, p. 173]:

$$\nu_t = C_\mu \phi k T$$

Where:

$C_{\mu}$ : Model coefficient for the turbulent viscosity [-]

$\nu_t$ : Turbulent viscosity [$\text{m}^2 \text{s}^{-1}$]

OpenFOAM implementation

Equations

The main differences between the implemented and original governing equations are as follows:

• Change of variables in the $\epsilon$-equation: $k/T \rightarrow \epsilon$ (when the term $P$ was expanded into its constituents),
• Change of variables in the $k$-equation: $\epsilon/k \rightarrow 1/T$,
• Inclusion of nu into DphitEff even though it is not present in 38, Eq. 17] (This change has provided higher level of resemblance to benchmarks for the tests considered, particularly for the peak skin friction, yet pressure-related predictions were unaffected. Users can switch off nu in DphitEff by using includeNu entry in kEpsilonPhitFCoeffs as shown below in order to follow the reference paper thereat. includeNu is left true by default.)
• The implementation is compressible and multiphase even though the derivation of the models in the original research paper did not consider compressibility or multiphase effects.

The implementation of the $\epsilon$ transport equation:

\ddt{\alpha \rho \epsilon} + \div \left(\alpha \rho \u \epsilon \right) - \laplacian \left(\alpha \rho D_\epsilon \epsilon \right) = \alpha \rho C_{\epsilon_1} \frac{G}{T} - \left( \left{ \frac{2}{3} C_{\epsilon_1} \right} \left{ \alpha \rho \div \u \right} \epsilon \right) - \left( \alpha \rho \frac{C_{\epsilon_2}}{T} \epsilon \right) + S_{\text{fvOptions}}

Where:

$\alpha$ : Phase fraction of the given phase [-]

$\rho$ : Density of the fluid [$\text{kg} m^{-3}$]

$S_{\text{fvOptions}}$ : Source terms introduced by fvOptions dictionary

The implementation of the $k$ transport equation:

$$\ddt{ \alpha \rho k } + \div \left( \alpha \rho \u k \right) - \laplacian \left( \alpha \rho D_k k \right) = \alpha \rho G - \left( \left[ \frac{2}{3} \alpha \rho \div \u k \right] \right) - \left( \frac{\alpha \rho}{T} k \right) + S_{\text{fvOptions}}$$

Where

$G$ : Anisotropic contribution part of the turbulent kinetic energy production rate, i.e. $P = G - 2/3 k \div \u$

The implementation of the $f$ elliptic relaxation equation:

$$- \laplacian f = - \left( \frac{1}{L^2} f \right) - \left( (C_{f_1} - 1) \frac{\phi - \frac{2}{3}}{T} - \frac{C_{f_2} G}{k} + C_{f_2} \frac{2}{3} \div \u - \frac{2 \nu (\nabla \phi \cdot \nabla k )}{k} - \nu \laplacian \phi \right)\frac{1}{L^2}$$

The implementation of the $\phi$ transport equation:

\ddt{\alpha \rho \phi} + \div \left(\alpha \rho \u \phi \right) - \laplacian \left(\alpha \rho D_\phi \, \phi \right) = \alpha \rho f - \left( \alpha \rho \left{ \frac{G}{k} - \frac{2}{3} \div \u - \frac{2 \nu (\nabla \phi \cdot \nabla k )}{k \sigma_k \phi } \right} \phi \right) + S_{\text{fvOptions}}

The implementation of the $T$ equation:

$$T = \max \left( \frac{k}{\epsilon}, C_T \frac{\sqrt{\max{\nu, 0}}}{\epsilon} \right)$$

The implementation of the $L$ equation:

$$L = C_L \max \left( \frac{k^{1.5}}{\epsilon}, C_\eta \left[ \frac{(\max{\nu, 0})^3}{\epsilon} \right]^{0.25} \right)$$

Default model coefficients

The model coefficients are 38, Eqs. 19-20]:

$$\begin{aligned} C_\mu = 0.22; \quad C_{\epsilon_1} = 1.4(1.0 + 0.05\sqrt{1.0/\phi} ); \quad C_{\epsilon_2} = 1.9; \quad C_T = 6.0; \quad C_L = 0.25;\\ C_{f_1} = 1.4; \quad C_{f_2} = 0.3; \quad C_\eta = 110.0; \quad \sigma_k = 1.0; \quad \sigma_\epsilon = 1.3; \quad \sigma_\phi = 1.0. \end{aligned}$$

Please note that 38, p. 176] stated that the model constants above have been calibrated to match the flow statistics of the smooth-wall plane channel flow at $\text{Re}_\tau = 395$.

Usage

The model can be enabled by using constant/turbulenceProperties dictionary:

simulationType    RAS;
RAS
{
  // Mandatory entries
  RASModel    kEpsilonPhitF;

  // Optional entries
  turbulence      on;
  printCoeffs     on;
  kEpsilonPhitFCoeffs
  {
      includeNu   true;    // include nu in DphitEff()
      Cmu         0.22;    // Turbulent viscosity constant
      Ceps1a      1.4;     // Model constant for epsilon
      Ceps1b      1.0;     // Model constant for epsilon
      Ceps1c      0.05;    // Model constant for epsilon
      Ceps2       1.9;     // Model constant for epsilon
      Cf1         1.4;     // Model constant for f
      Cf2         0.3;     // Model constant for f
      CL          0.25;    // Model constant for L
      Ceta        110.0;   // Model constant for L
      CT          6.0;     // Model constant for T
      sigmaK      1.0;     // Turbulent Prandtl number for k
      sigmaEps    1.3;     // Turbulent Prandtl number for epsilon
      sigmaPhit   1.0;     // Turbulent Prandtl number for phit = sigmaK
  }

Please note that $\phi$ is a symmetric, and $f$ is an asymmetric matrix.

Initial conditions

For $k$ and $\epsilon$ fields, the initial conditions can be estimated by using the recommendations made for the kEpsilon model.

For $\phi$:

object                phit;
dimensions            [0 0 0 0 0 0 0];
internalField         uniform 0.66;

For $f$:

object                f;
dimensions            [0 0 -1 0 0 0 0];
internalField         uniform 1.0;

Boundary conditions

For $k$ and $\epsilon$ fields, please use the boundary conditions that are available for the kEpsilon model.

For $\phi$ and $f$:

Inlet/Outlet:

• Zero gradient

Walls:

• No wall function is required for $\phi$ or $f$ fields,
• fixedValue with $\phi$ and $f$ zero at the wall can be used.

For example, $\phi$:

object          phit;
dimensions      [0 0 0 0 0 0 0];
wall
{
    type      fixedValue;
    value     uniform 1e-10;
}
"(inlet|outlet)"
{
    type      zeroGradient;
}

For example, $f$:

object          f;
dimensions      [0 0 -1 0 0 0 0];
wall
{
    type      fixedValue;
    value     uniform 0;
}
"(inlet|outlet)"
{
    type      zeroGradient;
}

Further information

Source code:

• Foam::RASModels::kEpsilonPhitF

References:

• : Laurence, Uribe, and Utyuzhnikov (2005) 38