Most of nowadays CFD simulations are conducted with the Reynolds Averaging approach. Reynolds Averaged Navier-Stokes (RANS) simulation, the “working horse” of industrial CFD is based on the Reynolds decomposition according to which a flow variable is decomposed into mean and fluctuating quantities. When the decomposition is applied to Navier-Stokes equation an extra term known as the Reynolds Stress Tensor arises and a modelling methodology is needed to close the equations. The “closure problem” is apparent as higher and higher moments of the set of equations may be taken, more unknown terms arise and the number of equations never suffices. This is of course an obvious consequence to the fact that taking these higher moments is simply a mathematical endeavor and has no physical contribution what so ever.
Levels of modeling are related to the number of differential equations added to Reynolds Averaged Navier-Stokes equations in order to “close” them.
0-equation (algebraic) models are the simplest form of turbulence models, a turbulence length scale is specified in advance through experimenting. 0-equations models are very limited in applications as they fail to take into account history effects, assuming turbulence is dissipated where it’s generated, a direct consequence of their algebraic nature.
1-equation and 2-equations models, incorporate a differential transport equation for the turbulent velocity scale (or the related the turbulent kinetic energy) and in the case of 2-equation models another transport equation for the length scale (or time scale), subsequently invoking the “Boussinesq Hypothesis” relating an eddy-viscosity analog to its kinetic gasses theory derived counterpart (albeit flow dependent and not a flow property) and relating it to the Reynolds stress through the mean strain.
In this sense 2-equation models can be viewed as “closed” because unlike 0-equation and 1-equation models (with exception of 1-equations transport for the eddy viscosity itself as described in a former post: Understanding The Spalart-Allmaras Turbulence Model) these models possess sufficient equations for constructing the eddy viscosity with no direct use for experimental results.
2-equations models do however contain many assumptions along the way for achieving the final form of the transport equations and as such are calibrated to work well only according to well-known features of the applications they are designed to solve. Nonetheless although their inherent limitations, today industry need for rapid answers dictates CFD simulations to be mainly conducted by 2-equations models whose strength has proven itself for wall bounded attached flows at high Reynolds number (thin boundary layers) due to calibration according to the law-of-the-wall.
A drawback evident in almost all eddy-viscosity models is the inability to inherently account for rotation and curvature. This drawback is resulted from relating the Reynolds stress to the mean flow strain and in fact is the major difference between such a modeling approach and a full Reynolds-stress model (RSM). The RSM approach accounts for the important effect of the transport of the principal turbulent shear-stress. On the other hand, RSM simulations are not computationally cost-effective, in as much that one does get an improved physical fidelity that is worth the time and computational resources consumed, not only that, they often do not converge.
Standard k-ε Turbulence Model
The k-ε turbulence model still remains among the most popular, most known is the Jones-Launder k-ε turbulence model.
To kick-off the brief description I shall remind the above paragraphs explanation, meaning 2-equation closure models carry two additional transport equations, in the case k-ε turbulence model they are the turbulence kinetic energy k and the turbulence dissipation ε: This eddy viscosity dimensional grounds based relation relating the Reynolds Stresses to the mean strain rate:
Realizable k-ε Turbulence Model
Although being perhaps the most popular and known turbulence model, the standard k-ε turbulence model carries along some harsh shortcomings which are important to acknowledge.
First, it is important to note that the model is essentially a high Reynolds model, meaning the law of the wall must be employed and provide velocity “boundary conditions” away from solid boundaries (what is termed “wall-functions”). From a mathematical standpoint, even if one could impose Dirichlet conditions for ε on solid boundary, after meshing it would still be difficult to numerically approach the problem due to what is termed in numerical analysis as stiffness of the numerical problem, partially related to the high gradients.
In order to integrate the equations through the viscous/laminar sublayer a “Low Reynolds” approach must be employed. This is achieved as additional highly non-linear damping functions are needed to be added to low-Reynolds formulations (low as in entering the viscous/laminar sublayer) to be able to integrate through the laminar sublayer (y+<5). This again produces numerical stiffness and in case is problematic to handle in view of linear numerical algorithms.
Furthermore the low Reynolds methodology should not be confused with transitional Reynolds modeling, as sometimes many practitioners rely on low Reynolds models to achieve a measure of transition from laminar to turbulence prediction. As the low Reynolds methodology is devised to handle the near wall viscous/laminar sublayer, there is no reason to expect it also satisfactory transition predictions, especially noting that there are many mechanisms for transition onset. Perhaps the only predictions which shall be close to satisfactory by low Reynolds model (maybe “pseudo-transition” behavior) are the ones related to bypass transition of which high levels of turbulence in the free stream occur and transition is dominated by diffusion effects (a brief description of the mechanism in my former post: A Forest of Hairpins – on the quest for turbulence coherent structures ).
Another major drawback is the model lack of sensitivity to adverse pressure-gradient. It was observed that under such conditions it overestimates the shear stress and by that delays separation. Menter’s k-ω SST alleviates this drawback through the Shear Stress Transport (SST) concept. This drawback is evident in almost all eddy-viscosity models, relating the Reynolds stress to the mean flow strain and in fact is the major difference between such a modeling approach and a full Reynolds-stress model (RSM). The RSM approach accounts for the important effect of the transport of the principal turbulent shear-stress. The ingenious idea of Menter to include it in the revised k-ω model (termed the Baseline (BSL) model) is related observed success in implementing what is termed as the Bradshaw’s assumption, that the shear-stress in the boundary layer is proportional to the turbulent kinetic energy (a brief evaluation of the model in former post: Understanding The k-ω SST Model ).
A step to circumvent some of the above presented drawbacks whilst still in the framework of the ε transport equation is made by invoking realizability constraints.
There are a number of such constraints, the usual ones are that all normal stresses should remain positive and the correlation coefficients for the shear stress should not exceed one:
Although the criteria is rarely used to close the RANS equations with two additional equations, satisfying the condition is of much importance for eddy viscosity models (EDMs), especially for stagnation flows.
If it is assumed that the flow that on one axis the flow approaches a wall, the Boussinesq Hypothesis for the normal stress becomes:
It could be seen that if s11 is too large then
which is nonphysical, i.e. non-realizable.
It is customary at this stage to introduce the concept of “invariant”, meaning something that is independent on a coordinate system. For the above case this relates to rotation.
Now, for a symmetric tensor the eigenvalues are known to be real, which means that one may rotate the coordinate system so that the off-diagonal components vanish and get the principal coordinate directions.
Symmetric (S) and anti-symmetric (Ω) parts of the velocity gradient tensor J
For the strain tensor (which is derived out of the velocity gradient tensor as explained in the above figure) this means that the off-diagonal components vanish but also that this is the direction of which the diagonal components become the largest. A direct consequence by which this coordinate system most susceptible for:
Now for the task of computing the invariants…
To solve for the eigenvalues of a tensor (given in 2-D):
Which in 2-D is:
Upon solving the equation we get:
Remembering that this set of equations is the same for any orientation of the coordinate system, we may deduce that its coefficients I1 andf I2 are invariants.
(We never forget that turbulence is 3D and rotational as the displaying the concept for a 2-D incompressible flow is the same…).
Continuing with the invariant:
The fact that it is zero is due to continuity.
Dealing with the second invariant we get:
Independent of orientation of the coordinate system:
The eigenvalues of S correspond to the strains in the principal axis, since we have applied the equation on the principle axis, s11 is replaced by the largest eigenvalue such that:
This simple modification to an eddy viscosity model ensures that the normal stresses stay positive.
Another realizability constraint appears when we require that if:
It shall be done smoothly.
A way to ensure such a behavior is to require that:
(Here d/dt is the Lagrangian derivative in Lagrangian coordinates, meaning following a fluid parcel as it approaches the wall).
We also impose a requirement that when approaches zero, the transport equation for it (left side/right side) shall do so to.
The main attention shifts to the pressure-strain term, and it is interesting to show how it behaves in the near wall region for the case of:
This is one of these cases where the normal stress goes to zero faster (O)(x4) than the parallel one (O(x2)) and creates the state of turbulence called the two component limit.
There are many models to justify such behavior, like the well known CL96, its interesting feature is that it does not need a prescribed wall distance parameter which makes it favorable for complex geometries.