TurbulenceMethods: complete list of utilities

This document focuses exclusively on the utilities implemented in the code and explains, for each function, exactly which mathematical expressions are implemented and how they are built from the input quantities.

1. computeStrainRotationInvariants

This function computes velocity-gradient-based invariants used in turbulence modeling: the strain-rate invariant , the rotation-rate invariant , and the divergence .

  1. It first obtains the gradient of the streamwise component : Analogously, when provided, it obtains and for the other velocity components.

  2. The code then defines the velocity-gradient components as

  3. The divergence of the velocity is computed as with terms omitted in lower dimensions. In code this is div_u = dux + dvy + dwz where the absent components default to zero.

  4. The symmetric strain-rate tensor is constructed as

  5. The antisymmetric rotation tensor is

  6. Finally the invariants returned are

The struct StrainRotationInvariants stores S2 = S^2, W2 = W^2, and div_u for use by other turbulence utilities.

2. cl_from_Cmu

This helper computes a near-wall length-scale constant as a function of :

This constant is used by the two-layer models (twoLayerWolfstein, twoLayerNorrisReynolds).

3. twoLayerWolfstein

This function implements a Wolfstein-type two-layer near-wall model. It returns a TwoLayerLengths struct with

  • l_eps: dissipation length scale ,

  • mu_ratio: viscosity ratio .

Using from cl_from_Cmu, where

  • is the wall distance,

  • is a Reynolds number based on wall distance.

4. twoLayerNorrisReynolds

This function implements a Norris–Reynolds-type two-layer model. It returns

with as before.

5. twoLayerXu

This function implements an Xu-type two-layer correlation using a modified wall coordinate (denoted yv_star in the code). It returns

Note that Cmu and Re_d are not used in this formulation.

6. f2_SKE_LRe

This function provides a low-Reynolds-number damping function for the standard model: where is a model constant and is a turbulent Reynolds number (typically ).

7. fmu_SKE_LRe

This function computes a low-Re damping factor for the eddy viscosity:

This is used to modify as for the first cell near the wall.

8. f2_RKE

This is another damping function (used in realizable or low-Re variants):

It approaches 0 near the wall (small ), and 1 away from the wall.

9. Cmu_realizable

This function computes a realizable eddy-viscosity constant based on the invariants and :

  1. First compute the ratio and then

  2. The realizable is then

with the constants pre-defined realizability constants , , , and .

This quantity is used in realizable models to define .

10. Cmu_nonlinear

This function provides a nonlinear-model-compatible . If or , it returns 0. Otherwise:

  1. Use inv.S2 and inv.W2 from computeStrainRotationInvariants and clamp them to non-negative values.

  2. Compute

  3. Then

This has the same functional form as Cmu_realizable, but uses the invariants packaged in StrainRotationInvariants and is used with nonlinear constitutive relations.

11. computeGk

This function implements the shear-production term for the -equation.

  • Always includes the basic term

  • If include_compressibility_terms == true, it further subtracts the compressibility corrections:

So the full expression is

Note, however, that for incompressible flows . So, include_compressibility_terms should have no effect on the shear-production term.

12. computeGb

This function computes the buoyancy production term in the -equation as where

  • is the thermal expansion coefficient,

  • is the turbulent Prandtl number,

  • is the gravity vector,

  • is the temperature gradient.

13. computeGammaM

This function returns the compressibility correction added to the -equation:

where

  • is a model constant,

  • is the speed of sound.

14. computeGammaY

This function implements the Yap correction for the equation.

  1. If , return 0.

  2. Otherwise compute where is the bulk eddy lenght scale, then take

  3. The final correction is

So injects extra dissipation whenever .

15. computeGprime

This function computes an additional low-Re production term for the -equation.

  1. First build where

    • is the shear-production term,

    • is wall distance.

  2. Then where and are constants, is typically a low-Re damping function, and is a Reynolds number based on wall distance.

16. computeVelocityGradient

This function constructs the full velocity gradient tensor as a RankTwoTensor.

  • It queries the gradients of , and, when present, and and fills:

This tensor is then used by the nonlinear constitutive relation and curvature correction.

17. Nonlinear constitutive relation utilities

The remaining utilities implement quadratic and cubic nonlinear eddy-viscosity models and associated terms.

There are fixed model coefficients in an anonymous namespace:

  • Quadratic/cubic model coefficients: CNL1–CNL7,

  • Ca0–Ca3 for the nonlinear expression. These are used internally by computeTRANS_NL and computeGnl.

17.1 computeTRANS_NL

This function builds the nonlinear part of the Reynolds-stress tensor associated with a quadratic or cubic constitutive relation.

  1. If model == None or any of are non-positive, it returns the zero tensor.

  2. Construct the strain and rotation tensors from the velocity gradient:

  3. Use the invariants in inv to compute

  4. Compute the nonlinear-compatible and then form the denominator

  5. From this denominator define

  6. Build the identity tensor and the quadratic tensor products:

  7. Construct the quadratic bracket:

  8. The quadratic nonlinear Reynolds-stress contribution is

    • If model == Quadratic, this is the final tensor returned.

  9. For the cubic model, additional terms are included:

    • Compute

    • Define and

    • Let and assemble the cubic bracket

    • With we obtain the cubic contribution

  10. In the cubic case, the total nonlinear stress tensor returned is .

17.2 computeGnl

This function returns the nonlinear production term associated with the nonlinear Reynolds stress:

  1. If model == None, it returns 0.

  2. Otherwise it first computes via computeTRANS_NL.

  3. The production term is then i.e. the double contraction between the nonlinear stress tensor and the velocity-gradient tensor.

This contribution is added to the standard production when nonlinear models are enabled.

18. computeCurvatureFactor

This function computes a curvature/rotation correction factor (of Spalart–Shur type) that multiplies the shear-production term.

  1. If model == None, it returns .

  2. Otherwise it starts from the invariants (clamped to non-negative):

  3. If , it again returns 1. Otherwise define

  4. Construct dimensionless measures of rotation vs strain:

  5. With internal constants it forms a Spalart–Shur-like correction

  6. The function then clamps to be at least and finally defines

The resulting is used to scale production in curvature/rotation-corrected models, enhancing or damping turbulence depending on the local ratio of rotation to strain.