kEpsilonOperator (1d, 2d, 3d)

Estimates the turbulent viscosity using k-epsilon model (

The kEpsilonOperator implements the right-hand side of the “Standard” Menter SST Two-Equation Model:

\[\notag \begin{align} \frac{\partial \left( \rho k \right)}{\partial t} + \nabla \cdot \left[ \rho \mathbf{u} k \right] = P - \rho \epsilon + \nabla \cdot \left[ \left( \mu + \frac{\mu_\mathrm{turb}}{\sigma_{k}} \right) \nabla k \right] + \rho L_{k} \\ \frac{\partial \left( \rho \epsilon \right)}{\partial t} + \nabla \cdot \left[ \rho \mathbf{u} \epsilon \right] = \frac{C_{\epsilon_{1}} f_{1} \epsilon}{k} P - \frac{C_{\epsilon_{2}} f_{2} \rho \epsilon^{2}}{k} + \nabla \cdot \left[ \left( \mu + \frac{\mu_\mathrm{turb}}{\sigma_{\epsilon}} \right) \nabla \epsilon \right] + \rho L_{\epsilon} \end{align}\]

The full details of this model, including the definition of the various constants, etc. can be found at (

The kEpsilonOperator operator computes the right-hand side of this model:

\[\notag \begin{align} \mathcal{S}_{\rho k} = P - \rho \epsilon + \nabla \cdot \left[ \left( \mu + \frac{\mu_\mathrm{turb}}{\sigma_{k}} \right) \nabla k \right] + \rho L_{k} \\ \mathcal{S}_{\rho \omega} = \frac{C_{\epsilon_{1}} f_{1} \epsilon}{k} P - \frac{C_{\epsilon_{2}} f_{2} \rho \epsilon^{2}}{k} + \nabla \cdot \left[ \left( \mu + \frac{\mu_\mathrm{turb}}{\sigma_{\epsilon}} \right) \nabla \epsilon \right] + \rho L_{\epsilon} \end{align}\]

The advective terms, \(\nabla \cdot \left[ \rho \mathbf{u} k \right]\) and \(\nabla \cdot \left[ \rho \mathbf{u} \omega \right]\) can be computed using classicMusclUpdater (1d, 2d, 3d) combined with multiSpeciesSingleVelocityEqn.


in (string vector of 7, required)
Fluid Model (nodalArray, 5-components, required)

The vector of conserved quantities for the fluid model, \(\mathbf{q}\) has 5 entries:

  1. \(\rho\): mass density
  2. \(\rho\,u_{\hat{\mathbf{i}}} = \rho \mathbf{u} \cdot \hat{\mathbf{i}}\): momentum density in the \(\hat{\mathbf{i}}\) direction
  3. \(\rho\,u_{\hat{\mathbf{j}}} = \rho \mathbf{u} \cdot \hat{\mathbf{j}}\): momentum density in the \(\hat{\mathbf{j}}\) direction
  4. \(\rho\,u_{\hat{\mathbf{k}}} = \rho \mathbf{u} \cdot \hat{\mathbf{k}}\): momentum density in the \(\hat{\mathbf{k}}\) direction
  5. \(E = \frac{P}{\gamma -1} + \tfrac{1}{2}\rho|\mathbf{u}|^2\): total energy density
Turbulence model (nodalArray, 2-components, required)

The vector of conserved quantities for the turbulence model:

  1. \(\rho k\)
  2. \(\rho \epsilon\)
Fluid velocity (nodalArray, 3-components, required)

Vector of fluid velocities, required if enableViscous = true:

  1. \(u_{\hat{\mathbf{i}}} = \mathbf{u} \cdot \hat{\mathbf{i}}\): velocity in the \(\hat{\mathbf{i}}\) direction
  2. \(u_{\hat{\mathbf{j}}} = \mathbf{u} \cdot \hat{\mathbf{j}}\): velocity in the \(\hat{\mathbf{j}}\) direction
  3. \(u_{\hat{\mathbf{k}}} = \mathbf{u} \cdot \hat{\mathbf{k}}\): velocity in the \(\hat{\mathbf{k}}\) direction

Fluid Temperature (nodalArray, 1-components, required)

Dynamic Viscosity (nodalArray, 1-components, required)

Thermal Conductivity (nodalArray, 1-components, required)

Distance from Wall (nodalArray, 1-components, required)

out (string vector of 4, required)

Vector of Fluid Model Source terms (nodalArray, 5-components, required)

  1. \(\mathcal{S}\left( \rho \right)\): mass source
  2. \(\mathcal{S}\left( \rho\,u_{\hat{\mathbf{i}}} \right)\): \(\hat{\mathbf{i}}\) momentum source
  3. \(\mathcal{S}\left( \rho\,u_{\hat{\mathbf{j}}} \right)\): \(\hat{\mathbf{j}}\) momentum source
  4. \(\mathcal{S}\left( \rho\,u_{\hat{\mathbf{k}}} \right)\): \(\hat{\mathbf{k}}\) momentum source
  5. \(\mathcal{S}\left( E \right)\): total energy source

Vector of Turbulence Model Source terms (nodalArray, 2-components, required)

  1. \(\mathcal{S}\left( \rho k \right)\)
  2. \(\mathcal{S}\left( \rho \epsilon \right)\)

Turbulent viscosity (nodalArray, 1-component, required)

Maximum turbulent diffusion (nodalArray, 1-component, required)


The kEpsilonOperator updater accepts the parameters below, in addition to those required by Updater:

numberOfInterpolationPoints (integer, required)

Number of points to be considerd for the least squares fit. This parameter varies from mesh to mesh and should be determined by computing a known function on the mesh.

The numberOfInterpolationPoints must be greater than (or equal to) the number of coefficients in the polynomial approximation. This means that in 1d the value is 4, in 2D the value is at least 6 and in 3D the value is at least 10.

These choices do not guarantee that a matrix inverse will be found. The following values though appear to be adequate in general: in 1D 4; in 2D 8 and in 3D 20.

orderAccuracy (integer, option)
Order of the polynomial that is used to form the operator. Choice of 1, 2 or 3 corresponding, respectively to first, second and third order accuracy. The appropriate choice of order varies on the problem type and the mesh used. Defaults to 2.
turbulentPrandtlNumber (float, required)
Prandtl number for turbulent flows, which is the ratio of eddy diffusivities of momentum and heat transfer
Cp (float, required)
Specific heat at constant pressure
minDt (float, required)
The dissipation time step is proportion to \(k / \epsilon\), which can go to infinity. We limit this ratio to the value specified by minDt.
computeViscosity (bool, required)
Whether to compute the turbulent viscosity
computeSource (bool, required)
Whether to compute source terms for the turbulence model.


<Updater computeRansSource>
  kind = kEpsilonOperator2d
  onGrid = domain
  coefficient = 1.0
  numberOfInterpolationPoints = 16
  turbulentPrandtlNumber = 0.85
  minDt = MINDT
  computeViscosity = false
  computeSource = true
  Cp = CP
  in = [q, kEpsilon, velocity, temperature, visc, cond, distance]
  out = [dummySource, kEpsilonSource, turbulentViscosity, maxTurbulentDiffusion]