2.1. Reynolds stress and realisability

The Reynolds stress tensor can be decomposed into isotropic and anisotropic components according to

(2.1) \begin{equation} \tau _{\textit{ij}} = \overline {u^{\prime}_i u^{\prime}_{\!j}} = a_{\textit{ij}} + \frac {2}{3}k\delta _{\textit{ij}}. \end{equation}

Here, $\tau _{\textit{ij}}=\overline {u^{\prime}_i u^{\prime}_{\!j}}$ is the Reynolds stress tensor, where $u^{\prime}_i$ and $u^{\prime}_{\!j}$ represent the fluctuating velocity components. The tensor $a_{\textit{ij}}$ is the anisotropy tensor, $k$ is the turbulent kinetic energy and $ ({2}/{3})k\delta _{\textit{ij}}$ is the isotropic stress (Lumley Reference Lumley1978). The operator, $\delta _{\textit{ij}}$ , is the Kronecker delta.

The anisotropy tensor is a second-order symmetric tensor with five independent components and zero trace. It can be normalised by $2k$ to obtain the normalised anisotropy tensor,

(2.2) \begin{equation} b_{\textit{ij}} = \frac {a_{\textit{ij}}}{2k} = \frac {\overline {u^{\prime}_i u^{\prime}_{\!j}}}{2k} - \frac {1}{3}\delta _{\textit{ij}}. \end{equation}

The normalised anisotropy tensor, $b_{\textit{ij}}$ , is constrained by the requirement that the Reynolds stress tensor is positive semi-definite or, equivalently, that the eigenvalues are non-negative (Schumann Reference Schumann1977). The realisability constraints may be physically interpreted as a condition that the turbulent kinetic energy be non-negative.

The non-negativity of the eigenvalues implies that the minimum eigenvalue of the Reynolds stress tensor is 0. The maximum occurs when only one eigenvalue is non-zero, as that eigenvalue will be equal to $2k$ by definition of the turbulent kinetic energy. Substituting the maximum and minimum eigenvalues into (2.2) allows realisability to be written in terms of the eigenvalues of the normalised anisotropy tensor (Lumley Reference Lumley1978),

(2.3) \begin{equation} -\frac {1}{3} \leqslant \lambda _{3} \leqslant \lambda _2 \leqslant \lambda _1 \leqslant \frac {2}{3}. \end{equation}

Here, $\lambda _3 \leqslant \lambda _2 \leqslant \lambda _1$ are the eigenvalues of $b_{\textit{ij}}$ . Note that it is unnecessary to consider $\lambda _2$ when determining realisability because it will always be between $\lambda _1$ and $\lambda _3$ (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007). Additionally, as a consequence of $b_{\textit{ij}}$ being traceless, the eigenvalues cannot all have the same sign. From this, it is clear that $\lambda _1 \geqslant 0$ and $\lambda _3 \leqslant 0$ .

Because of the bounds on $b_{\textit{ij}}$ , it is possible to represent all turbulent states in a finite region using a coordinate system constructed from the invariants of the normalised anisotropy tensor.

2.2. Current invariant maps

Invariants are quantities associated with a tensor that do not change when the coordinate system is rotated. The principal invariants for an arbitrary second-order symmetric tensor, $\sigma _{\textit{ij}}$ , are

(2.4) \begin{equation} I_1=\sigma _{ii}, \quad I_2 = \tfrac {1}{2}\left(\sigma _{ii}^2 - \sigma _{\textit{ij}} \sigma _{\textit{ij}}\right), \quad I_3=\textrm {det}(\sigma _{\textit{ij}}). \end{equation}

These three invariants are referred to as the first, second and third invariants (Irgens Reference Irgens2019). For the anisotropic part of a second-order symmetric tensor, the first invariant is zero and the other two principal invariants become (Chakrabarty Reference Chakrabarty2006)

(2.5) \begin{equation} J_2=\frac {1}{2}\sigma _{\textit{ij}}^{\textit{dev}}\sigma _{\textit{ij}}^{\textit{dev}}=\frac {1}{3}I_1^2-I_2, \quad J_3 = \textrm {det}\left(\sigma _{\textit{ij}}^{\textit{dev}}\right)=\frac {2I_1^3-9I_1I_2+27I_3}{27}. \end{equation}

Here, $J_2$ and $J_3$ are the principal invariants of $\sigma _{\textit{ij}}^{\textit{dev}}$ , which is the deviatoric part of $\sigma _{\textit{ij}}$ .

Three invariant maps have been commonly used in the study of turbulence and turbulence modelling. The first was introduced by Lumley & Newman (Reference Lumley and Newman1977) using the second and third invariants of the anisotropy tensor,

(2.6) \begin{equation} J_{2,b} = \tfrac {1}{2}b_{\textit{ij}} b_{\textit{ij}}, \quad J_{3,b} = \det {(b_{\textit{ij}})}. \end{equation}

The subscript $b$ indicates that these invariants are associated with $b_{\textit{ij}}$ . Lumley used these invariants to form a region with two curved sides and one straight side that contained all realisable states of turbulence.

The second invariant map was introduced by Choi & Lumley (Reference Choi and Lumley2001) by defining two new invariants derived from the principal invariants,

(2.7) \begin{equation} \xi ^3 = J_{3,b}/2, \quad \eta ^2 = J_{2,b}/3. \end{equation}

The invariants, $\xi$ and $\eta$ , form a mapping with two straight sides and one curved side. This modification was motivated by the study of the return to isotropy of homogeneous turbulence. By using $\xi$ and $\eta$ , a linear return to isotropy as proposed by Rotta (Reference Rotta1951) would show as a straight line through the origin, allowing any nonlinearities to be clearly visualised. Both invariant maps have been used in the development of return-to-isotropy models for Reynolds Stress Transport Models (Lumley & Newman Reference Lumley and Newman1977; Choi & Lumley Reference Choi and Lumley2001).

While these mappings were highly nonlinear, which can result in misleading visualisations of the turbulent state (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007). For example, in an early DNS study of fully developed turbulent channel flow, the nonlinear mapping made it appear that the turbulence was nearly isotropic at the centreline, when in fact, it was far from isotropic (Antonia, Kim & Browne Reference Antonia, Kim and Browne1991). A similar conclusion was made in a different study of turbulent pipe flow (Krogstad & Torbergsen Reference Krogstad and Torbergsen2000).

The most recent invariant map was proposed by Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007) to provide unbiased visualisations of the turbulent state and as a formulation for exact interpolation between states. Banerjee et al. used the eigenvalues to construct a barycentric coordinate system where the turbulent states were convex combinations of the limiting states of turbulence, defined by

(2.8) \begin{equation} \begin{bmatrix} x_B \\ y_B \end{bmatrix} = \begin{bmatrix} x_{1c} & x_{2c} & x_{3c} \\ y_{1c} & y_{2c} & y_{3c} \end{bmatrix} \begin{bmatrix} C_1 \\ C_2 \\ C_3 \end{bmatrix}. \end{equation}

Here, $x_{1c}$ , $y_{1c}$ and so on are the Cartesian coordinates of the corners of the map. Additionally, $C_1$ , $C_2$ and $C_3$ are the barycentric coordinates for the first, second and third corners, respectively. The coordinates $x_B$ and $y_B$ are the Cartesian coordinates of the barycentric map. The barycentric coordinates are defined by

(2.9) \begin{equation} \begin{matrix} C_1 = \lambda _1 - \lambda _2, \quad C_2 = 2(\lambda _2 - \lambda _3), \quad C_3 = 1 + 3\lambda _3, \\[4pt] \textrm {where} \quad C_1 + C_2 + C_3 = 1, \quad C_1, C_2, C_3 \geqslant 0. \end{matrix} \end{equation}

All of these invariant maps are depicted in figure 1. In any invariant map, the corners represent the three limiting states of turbulence: one-component (1C), two-component axisymmetric (2C) and isotropic (3C). The 1C–2C line indicates states of two-component turbulence, but not necessarily axisymmetric. The 1C–3C line indicates axisymmetric expansion (AE), and the 2C–3C line indicates axisymmetric contraction (AC) (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007). It is important to note that early researchers have used these terms to refer to both the magnitude of the principal stresses and the physical shape of the eddies. This was a source of confusion as the identification of the sides corresponding to axisymmetric expansion and contraction would flip depending on which interpretation was used (Simonsen & Krogstad Reference Simonsen and Krogstad2005).

Here, we use the convention used by Lumley & Newman (Reference Lumley and Newman1977), Choi & Lumley (Reference Choi and Lumley2001) and Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007), which refer to the relative magnitudes of principal Reynolds stresses. One-component means one principal stress is non-zero; two-component, axisymmetric means two principal stresses are non-zero and equal in magnitude; and isotropic means all principal stresses are non-zero and equal. The terms axisymmetric expansion and contraction in this context refer to the behaviour of the principal Reynolds stresses as the anisotropy magnitude shrinks. If the two equal principal stresses grow, it is in a state of expansion, and if they shrink, it is in a state of contraction.

Figure 1.

A comparison of the $J_{2,b}$ – $J_{3,b}$ , $\xi$ – $\eta$ and barycentric maps. 1C indicates one-component turbulence, 2C-axi indicates two-component axisymmetric turbulence, and 3C indicates isotropic turbulence. The dashed line indicates the plane-strain limit.

3.1. Analytical form of eigenvalues

It is informative to begin from the eigenvalues of the normalised anisotropy tensor, $b_{\textit{ij}}$ . The eigenvalues of the arbitrary second-order, symmetric tensor, $\sigma _{\textit{ij}}$ , are the solutions to the cubic equation (Pope Reference Pope2000),

(3.1) \begin{equation} \lambda ^3 - I_1 \lambda ^2 + I_2 \lambda - I_3 = 0. \end{equation}

Here, $I_1$ , $I_2$ and $I_3$ are the principal invariants of $\sigma _{\textit{ij}}$ (Irgens Reference Irgens2019). The solution to this equation can be obtained through Vieta’s trigonometric formula for the roots of a cubic polynomial. This gives the following solution for the eigenvalues in terms of the invariants (Zucker Reference Zucker2008):

(3.2a) \begin{align} \lambda _n &= \frac {I_1}{3} + \sqrt {\frac {4J_2}{3}}\cos {\left [\theta - \frac {2\pi (n-1)}{3}\right ]}, \quad n = 1, 2, 3, \end{align}

(3.2b) \begin{align} \theta &= \frac {1}{3}\arccos \left (\frac {3\sqrt {3}}{2} \frac {J_{3}}{J_{2}^{3/2}}\right ). \end{align}

In (3.2), $\theta$ is an invariant that exists for all non-zero, symmetric, 3 × 3 tensors (Chakrabarty Reference Chakrabarty2006). Although we are unaware of any examples of its usage in turbulence modelling, it is an established invariant in solid mechanics that is used in the Haigh–Westergaard coordinate system for modelling material behaviour and visualising stress states (Kassir Reference Kassir2017).

Substituting in the definitions of $I_1$ , $J_2$ and $J_3$ , (3.2a ) may be rewritten in terms of $\sigma _{\textit{ij}}$ and $\theta$ as

(3.3) \begin{equation} \lambda _n = \frac {\sigma _{ii}}{3} + \sqrt {\frac {2}{3}\sigma _{\textit{ij}}^{\textit{dev}}\sigma _{\textit{ij}}^{\textit{dev}}}\cos {\left [\theta - \frac {2\pi (n-1)}{3}\right ]}, \quad n = 1, 2, 3. \end{equation}

An additional simplification can be made by recognising $\sigma _{\textit{ij}}^{\textit{dev}}\sigma _{\textit{ij}}^{\textit{dev}}$ as the square of Frobenius norm, $||\sigma _{\textit{ij}}^{\textit{dev}}||^2$ . The general solution simplifies to

(3.4) \begin{equation} \lambda _n = \frac {\sigma _{ii}}{3} + \frac {\sqrt {6}}{3}||\sigma _{\textit{ij}}^{\textit{dev}}||\cos {\left [\theta - \frac {2\pi (n-1)}{3}\right ]}, \quad n = 1, 2, 3. \end{equation}

Using the same simplification for (3.2b ) yields

(3.5) \begin{equation} \theta = \frac {1}{3}\arccos \left (\frac {3\sqrt {3}}{2(1/2)^{3/2}} \frac {\textrm {det}(\sigma _{\textit{ij}}^{\textit{dev}})}{||\sigma _{\textit{ij}}^{\textit{dev}}||^3}\right ) = \frac {1}{3}\arccos \left (3\sqrt {6} \ \textrm {det}\left (\frac {\sigma _{\textit{ij}}^{\textit{dev}}}{||\sigma _{\textit{ij}}^{\textit{dev}}||}\right )\right ). \end{equation}

In (3.5), bringing the norm of $\sigma _{\textit{ij}}^{\textit{dev}}$ into the determinant follows from the property that $\textrm {det}(c \boldsymbol{A})=c^n \det (\boldsymbol{A})$ for any scalar $c$ and matrix $\boldsymbol{A}$ of rank $n$ (Baker & Porteous Reference Baker and Porteous1990).

Equations (3.4) and (3.5) can now be applied to the normalised anisotropy tensor by substituting $b_{\textit{ij}}$ for $\sigma _{\textit{ij}}$ . Because $b_{\textit{ij}}$ has zero trace, the first term on the right-hand side of (3.4) is eliminated and $b_{\textit{ij}}=\sigma _{\textit{ij}}^{\textit{dev}}$ . Making these substitutions gives

(3.6a) \begin{align} \lambda _n &= \frac {\sqrt {6}}{3}||b_{\textit{ij}}||\cos {\left [\theta - \frac {2\pi (n-1)}{3}\right ]}, \quad n = 1, 2, 3. \end{align}

(3.6b) \begin{align} \theta &= \frac {1}{3}\arccos \left (3\sqrt {6} \ \textrm {det}\left (\frac {b_{\textit{ij}}}{||b_{\textit{ij}}||}\right )\right ) \in \left [0,\frac {\pi }{3}\right ]. \end{align}

The bounds on $\theta$ are a result of the limited range of the arccosine function. Neglecting cases in which $||b_{\textit{ij}}|| = 0$ , the restrictions on the domain of arccosine are satisfied by all symmetric, deviatoric tensors (Chakrabarty Reference Chakrabarty2006). A proof of this may be found in Appendix A.

Evaluating (3.6a ) for $\lambda _1$ , $\lambda _2$ and $\lambda _3$ finally gives explicit forms for the eigenvalues of $b_{\textit{ij}}$ ,

(3.7a) \begin{align} \lambda _1 &= \frac {\sqrt {6}}{3}||b_{\textit{ij}}||\cos {\theta }, \end{align}

(3.7b) \begin{align} \lambda _2 &= \frac {\sqrt {6}}{3}||b_{\textit{ij}}||\cos {\left [\theta - \frac {2\pi }{3}\right ]}, \end{align}

(3.7c) \begin{align} \lambda _3 &= \frac {\sqrt {6}}{3}||b_{\textit{ij}}||\cos {\left [\theta + \frac {2\pi }{3}\right ]}, \end{align}

It should be noted that $\theta$ is undefined when $||b_{\textit{ij}}|| = 0$ due to $b_{\textit{ij}}/||b_{\textit{ij}}||$ becoming indeterminate. The coefficient $||b_{\textit{ij}}||$ in (3.7a )–(3.7c ) drives the eigenvalues to zero, so the value of $\theta$ does not affect the output. This occurs for purely isotropic flow, where the normalised anisotropy tensor becomes a zero matrix. Such a case can be handled numerically by explicitly setting the eigenvalues to zero when $||b_{\textit{ij}}||=0$

3.2. Normalisation of anisotropy magnitude

In (3.7), $||b_{\textit{ij}}||$ and $\theta$ form a set of polar coordinates where $||b_{\textit{ij}}||$ describes the anisotropy magnitude and $\theta$ describes the relationship of the eigenvalues or principal stresses. This is most easily interpreted by viewing the eigenvalues as forming a vector in the principal axes, where $\theta$ parametrises the direction of the vector. Figure 2 depicts this scenario with the unit vector, $\boldsymbol{\hat {\lambda }}$ . In figure 2, the axes are the eigenvectors $\boldsymbol{e}_1$ , $\boldsymbol{e}_2$ and $\boldsymbol{e}_3$ , associated with the principle stresses, $\lambda _1$ , $\lambda _2$ and $\lambda _3$ . The circle represents the plane where $\lambda _1 + \lambda _2 + \lambda _3 = 0$ , and its edge is defined by the unit vectors of the principal anisotropic stresses, denoted as $\boldsymbol{\hat {\lambda }}$ . Equation (3.7a ) constrains the possible combinations of the principal stresses to the sector bounded by the blue lines.

Figure 2.

A geometric interpretation of $\theta$ . The axes are the eigenvectors of $b_{\textit{ij}}$ , the blue lines indicate the bounds on $\theta$ and the red line is a unit vector.

Converting from standard Cartesian invariant maps to a polar invariant map provides an alternative perspective on the realisability bounds as a limit on the normalised anisotropy magnitude. Using this concept of realisability as a limit on the normalised anisotropy magnitude, we can normalise $||b_{\textit{ij}}||$ according to its maximum realisable magnitude at the same $\theta$ as

(3.8) \begin{equation} \alpha \equiv \frac {||b_{\textit{ij}}||}{||b_{\textit{ij}}||_{\textit{max}}}. \end{equation}

Here, $\alpha$ is the fraction of the maximum anisotropy. It can be interpreted as a normalised distance from isotropy in that a value of zero implies perfectly isotropic turbulence, while a value of one implies maximally anisotropic turbulence. In this context, $\theta$ takes on the role of determining the limiting anisotropy magnitude.

The maximum value of $||b_{\textit{ij}}||$ can be determined as a function of $\theta$ by substituting (3.7a ) and (3.7c ) into the inequality in (2.3). This gives the set of inequalities,

(3.9a) \begin{align} \lambda _1 &= \frac {\sqrt {6}}{3}||b_{\textit{ij}}||\cos {\theta } \leqslant 2/3, \end{align}

(3.9b) \begin{align} \lambda _3 &= \frac {\sqrt {6}}{3}||b_{\textit{ij}}||\cos {\left [\theta + \frac {2\pi }{3}\right ]} \geqslant -1/3. \end{align}

A negative can be factored out of $\lambda _3$ by introducing a phase shift of $\pi$ radians. Then, the constraints on $||b_{\textit{ij}}||$ may be expressed as

(3.10a) \begin{align} ||b_{\textit{ij}}|| &\leqslant \frac {\sqrt {6}}{3\cos {\theta }}, \end{align}

(3.10b) \begin{align} ||b_{\textit{ij}}|| &\leqslant \frac {\sqrt {6}}{6\cos {[\theta - \frac {\pi }{3}}]}. \end{align}

It is simple to verify that the first expression always exceeds or equals the second for $\theta \in [0,\pi /3]$ . Thus, the limiting anisotropy magnitude is determined by the more restrictive inequality, (3.10b ), to be

(3.11) \begin{equation} ||b_{\textit{ij}}||_{\textit{max}} = \frac {\sqrt {6}}{6\cos {\left[\theta - \dfrac {\pi }{3}\right]}}. \end{equation}

Using $||b_{\textit{ij}}||_{\textit{max}}$ , an explicit form for $\alpha$ may be obtained as

(3.12) \begin{equation} \alpha = \frac {||b_{\textit{ij}}||}{||b_{\textit{ij}}||_{\textit{max}}} = \sqrt {6}||b_{\textit{ij}}||\cos {\left (\theta - \frac {\pi }{3}\right )}. \end{equation}

3.3. Polar invariant map properties

With $\alpha$ and $\theta$ known, the polar invariant map can be plotted as

(3.13) \begin{equation} x_p = \alpha \cos {\theta }, \quad y_p = \alpha \sin {\theta }, \end{equation}

where $x_p$ and $y_p$ are the Cartesian x and y coordinates, respectively. An annotated diagram of the polar invariant map is depicted in figure 3. Like all invariant maps, the three corners correspond to the limiting states of turbulence: one-component; two-component, axisymmetric; and isotropic. The line $\theta =0$ represents states of axisymmetric expansion and the line $\theta = \pi /3$ represents states of axisymmetric contraction. Here, the terms axisymmetric expansion and contraction refer to the behaviour of the principal Reynolds stresses as the anisotropy magnitude shrinks. If the two equal principal stresses grow, it is in a state of expansion, and if they shrink, it is in a state of contraction (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007).

Figure 3.

An annotated diagram of the polar $\alpha$ – $\theta$ invariant map.

Splitting the invariant map into two halves is the plane-strain line at $\theta = \pi /6$ . Below the plane-strain line, $\lambda _2$ is less than 0, implying $\lambda _1$ has a greater magnitude than $\lambda _3$ . Conversely, above the plane-strain line, $\lambda _2$ is greater than 0, implying $\lambda _1$ has a smaller magnitude than $\lambda _3$ . At the plane-strain line, $\lambda _2 = 0$ , so $\lambda _1$ and $\lambda _3$ have equal magnitude.

Lines of constant $\theta$ on the map imply a fixed ratio of the principal stresses, while lines of constant $\alpha$ imply a fixed ratio of the actual and maximum anisotropy magnitudes. Additional insight can be gained by substituting the definition of $\alpha$ into (3.7) to obtain eigenvalue formulae in terms of $\alpha$ and $\theta$ ,

(3.14a) \begin{align} \lambda _1 &= \frac {\alpha }{3}\cos {\theta }\sec {\left [\theta - \frac {\pi }{3}\right ]}, \end{align}

(3.14b) \begin{align} \lambda _2 &= \frac {\alpha }{3}\cos {\left [\theta - \frac {2\pi }{3}\right ]}\sec {\left [\theta - \frac {\pi }{3}\right ]}, \end{align}

(3.14c) \begin{align} \lambda _3 &= -\frac {\alpha }{3}. \end{align}

Two important observations may be made from (3.14). First, $\lambda _3$ depends solely on $\alpha$ . Additionally, when $\alpha = 1$ , $\lambda _3$ reaches its minimum value of –1/3. Thus, the maximum anisotropy magnitude always occurs when $\lambda _3 = -1/3$ . This corresponds to the two-component turbulence arc on the polar invariant map.

Second, because $\alpha$ is positive by definition and $\theta$ is defined on $[0,\pi /3]$ for all symmetric, deviatoric tensors, all non-realisable states must have $\alpha \gt 1$ . Consequently, knowing $\alpha$ , or equivalently $\lambda _3$ , is sufficient to determine whether a point in the flow field violates realisability. Furthermore, all combinations of $\alpha \in [0,1]$ and $\theta \in [0,\pi /3]$ are realisable states.

The preceding observations lead to two convenient properties. First, the nearest realisable point to any non-realisable point is one with the same $\theta$ , but at $\alpha =1$ . Second, because of the form of the eigenvalues, any non-realisable state, $\alpha \gt 1$ , can be brought to the nearest realisable state, $\alpha =1$ , by simply scaling the original tensor by the multiplicative inverse, $1/\alpha$ . More generally, a tensor can be shifted in the polar invariant map along a line of constant $\theta$ using a scaling factor of

(3.15) \begin{equation} \beta = \frac {\alpha _1}{\alpha _0}, \end{equation}

where $\alpha _0$ is the original $\alpha$ value and $\alpha _1$ is the target $\alpha$ . This provides a direct link to turbulence models that can be written as a nonlinear eddy viscosity model by representing realisability limits in terms of a scaling factor analogous to the eddy viscosity. Such an approach would be especially useful for the myriad of machine learning-based turbulence modelling concepts, which are generally challenging to enforce realizability with. The potential applications will be elaborated upon in § 4.

4.1. Mathematical relationship

To more precisely show the relationship between the polar invariant map and the other invariant maps, the equations for their coordinates can be rewritten in terms of $\alpha$ and $\theta$ . As a first step, (2.6) may be rewritten in terms of $\lambda _1$ and $\lambda _3$ as

(4.1) \begin{equation} J_{2,b} = \lambda _1^2 +\lambda _1 \lambda _3 + \lambda _3^2, \quad J_{3,b} = -\lambda _1^2 \lambda _3 - \lambda _1 \lambda _3^2. \end{equation}

Next, define $\zeta = \cos {(\theta )}\sec {(\theta -\pi /3)}$ . Substituting (3.14a ) and (3.14c ) into (4.1) allows the invariants to be written in terms of $\alpha$ and $\theta$ ,

(4.2) \begin{equation} J_{2,b} = \left (\frac {\alpha }{3}\right )^2 [\zeta ^2 - \zeta + 1], \quad J_{3,b}=\left (\frac {\alpha }{3}\right )^3 [\zeta ^2 - \zeta ]. \end{equation}

Note that $\zeta$ is equivalent to $\lambda _1/\lambda _3$ . Using the invariants as given in (4.2), $\xi$ and $\eta$ become

(4.3) \begin{equation} \eta = \frac {\alpha }{3} \left [\frac {\zeta ^2 - \zeta + 1}{3}\right ]^{1/2}, \quad \xi = \frac {\alpha }{3} \left [\frac {\zeta ^2 - \zeta }{2}\right ]^{1/3}. \end{equation}

It is clear from (4.3) that the nonlinearity in the $\xi$ – $\eta$ invariant map comes exclusively from the $\theta$ dependence. This indicates that the $\xi$ – $\eta$ invariant map represents the relationship between eigenvalues in a nonlinear manner, and the relationship between the actual and maximum anisotropy magnitude in a linear manner.

To develop mathematical relationships between the polar and barycentric maps, it is helpful to begin by using the property of the normalised anisotropy tensor having zero trace to eliminate $\lambda _2$ from the barycentric coordinates. This gives

(4.4) \begin{equation} C_1 = 2\lambda _1 + \lambda _3, \quad C_2 = -2(\lambda _1 + 2\lambda _3), \quad C_3 = 3\lambda _3 + 1. \end{equation}

Equations (3.14a ) and (3.14c ) can be substituted into (4.4) to obtain the barycentric coordinates in terms of $\alpha$ and $\theta$ ,

(4.5) \begin{equation} C_1 = \frac {\alpha }{3}(2\zeta - 1), \quad C_2 = \frac {2\alpha }{3}(2 - \zeta ), \quad C_3 = 1 - \alpha . \end{equation}

Note that $C_3$ is a linear function of only $\alpha$ . This is unsurprising as the $C_3$ axis in the barycentric map describes the level of isotropy. In contrast, $C_1$ and $C_2$ depend on both $\alpha$ and $\theta$ . This also makes sense as the $C_1$ and $C_2$ axes have both vertical and horizontal components.

Banerjee et al. (Reference Banerjee, Krahl, Durst and Zenger2007) defined the turbulence triangle’s Cartesian coordinates using $(x_{1c},y_{1c}) = (1,0)$ , $(x_{2c},y_{2c}) = (0,0)$ and $(x_{3c},y_{3c}) = (1/2, \sqrt {3}/2)$ . Using these coordinates for the corners of the barycentric map in (2.8) simplifies the transformation of barycentric to Cartesian coordinates to

(4.6) \begin{equation} x_B = C_1 + \frac {1}{2}C_3, \quad y_B = \frac {\sqrt {3}}{2}C_3. \end{equation}

Substituting (4.5) into (4.6) gives the Cartesian coordinates in terms of $\alpha$ and $\theta$ ,

(4.7) \begin{equation} x_B = \frac {1}{2} + \frac {\alpha }{3}\left (2\zeta -\frac {5}{2}\right ), \quad y_B = \frac {\sqrt {3}}{2}(1 - \alpha ). \end{equation}

In (4.7), similar to the barycentric coordinates, the Cartesian coordinates show that the vertical axis is solely dependent on $\alpha$ , while the horizontal axes are dependent on both. The dependence on $\alpha$ and $\theta$ is necessary because the range of $x_B$ contracts as $y_B$ increases.

Using (4.3) and (4.7), $\theta$ and $\alpha$ isolines can be plotted on the barycentric and $\xi$ – $\eta$ invariant maps to visualise the distortion relative to the polar map. This is depicted in figure 5.

Figure 5.

The $\xi{ -}\eta$ map and barycentric map overlaid with lines of $\alpha$ and $\theta$ . Red lines indicate constant $\alpha$ and blue lines indicate constant $\theta$ . The black dashed line indicates the plane-strain limit.

The most significant difference between both the $\xi$ – $\eta$ and barycentric maps compared with the polar invariant map is the spacing of the $\theta$ isolines. In the barycentric map, the $\theta$ isolines tend to compress when moving from the one-component to two-component corners. This is consistent with the bias towards the two-component corner observed for the polar invariant map in figure 4. In comparison, the $\xi$ – $\eta$ map shows almost the opposite distortion in the $\theta$ isolines as the barycentric map. Here, the isolines are strongly pushed apart near the plane-strain limit at $\theta = \pi /6$ and are strongly compressed near the straight edges.

Interestingly, despite the differences in the shape of the $\alpha$ isolines, both maps preserve the linearity of $\alpha$ as the normalised distance from the isotropic corner. In the barycentric map, this manifests as $\alpha$ consisting of equally spaced straight lines. The $\xi$ – $\eta$ map instead displays $\alpha$ isolines that are equally spaced along any given $\theta$ isoline, but with a different spacing for each isoline.

Additionally, the polar invariant map sheds light on a disagreement between Mandler & Weigand (Reference Mandler and Weigand2022) and Blauw (Reference Blauw2019) over whether all non-realisable points fall below the 1C–2C line on the barycentric map. The polar invariant map clarifies that any point with a valid $\theta$ can only be non-realisable by exceeding $\alpha =1$ , corresponding to points that fall below the 1C–2C line on the barycentric map. Points that fall outside the 1C–3C and 2C–3C lines correspond to points outside of $\theta =[0,\pi /3]$ . This would require a violation of the range of (3.6b ), which is valid for all symmetric, deviatoric tensors. It is clear then that the points Mandler & Weigand (Reference Mandler and Weigand2022) noted to fall outside the 1C–3C and 2C–3C lines during the solution convergence process must occur due to the violation of either the zero-trace or symmetric properties of the anisotropy tensor. This is unsurprising for a solution in the process of converging, as intermediate iterations frequently show violations of continuity, which result in a non-deviatoric anisotropy tensor.

4.2. Comparison with sample data

To illustrate the differences in the invariant maps on real data, we will compare DNS data of a fully developed channel flow and a cross-section of a square duct in each mapping.

4.2.1. Fully developed channel flow

First, consider a fully developed plane channel flow. A time-averaged DNS dataset was obtained from Hoyas & Jiménez (Reference Hoyas and Jiménez2006). The flow field is statistically two-dimensional, and the velocity and Reynolds stress profiles are self-similar (Pope Reference Pope2000). Figure 6 depicts the DNS Reynolds stress profile for ${Re}_\tau =2003$ visualised in each invariant map.

Figure 6.

A comparison of fully developed channel flow at ${Re}_\tau =2003$ in each invariant map. The data are coloured by $y^+$ . Lines interpolated from the DNS data are underlaid for clarity.

All three invariant maps show that the turbulence starts near the wall close to the plane-strain limit on the 1C–2C line, migrates towards the one-component corner and then moves closer to the centreline, before following a straight path to terminate near the 1C–3C side. Both the polar and barycentric maps show that the turbulence consistently moves closer to isotropy, though it remains significantly anisotropic throughout the entire profile. In contrast, the $\xi$ – $\eta$ map obscures this trend. Due to the expansion of the near plane-strain region, the profile appears significantly closer to the axisymmetric expansion line than it truly is. Similarly, the trajectory in the near-wall region appears to start much closer to the one-component corner in the $\xi$ – $\eta$ map than in any other mapping. It also shows a kink in the profile at approximately $y^+$ = 250–750, which causes it to appear to be moving away from the isotropic corner, even though the barycentric and polar maps clearly show the profile continually moves towards that corner. These observations demonstrate how strongly nonlinear mappings like the $\xi$ – $\eta$ map can produce misleading visualisations and that the polar invariant map avoids this disadvantage.

While all of the maps provide information about the componentiality, the polar invariant map complements this interpretation by also making clear the relationship between eigenvalues and their relative magnitudes. The polar invariant map shows that the observation that the anisotropy approaches the isotropic corner is mathematically equivalent to stating the anisotropy magnitude is decreasing relative to its maximum value. Furthermore, the radial axis quantifies this change, indicating that the anisotropy magnitude at the centre of the channel is approximately 20 % of its maximum value for the given $\theta$ . It can also be seen that the relationship of the eigenvalues changes much more slowly past approximately y+ = 250 based on the low variation in $\theta$ . Even the barycentric map, which provides very similar visualisations to the polar invariant map, appears to show a much greater rate of change in the relationship of the eigenvalues due to the expanded spacing of $\theta$ noted in figure 5. Thus, although the information is contained in other invariant maps, it is not readily accessible or quantifiable like in the polar invariant map.

4.2.2. Square duct flow

Next, we will compare the invariant maps for fully developed square duct flow. The DNS dataset was obtained from the turbulence modelling database of McConkey, Yee & Lien (Reference McConkey, Yee and Lien2021), using the DNS dataset generated by Pinelli et al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010). A diagram of the square duct is displayed in figure 7. The flow domain has a height and width of $2H$ , and a length of $L$ . All of the square duct data are provided as time-averaged quantities.

Figure 7.

A diagram of the square duct flow. The flow field is symmetric across the dashed diagonals.

Square ducts are one of the simplest geometries that exhibit secondary flows. Here, the secondary flows are pairs of counter-rotating vortices located in each corner of the duct. The DNS data from Pinelli et al. (Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010) were collected for sixteen different Reynolds numbers; however, only the maximum Reynolds number, ${Re} = 3500$ , will be shown here.

Figure 8 compares the Reynolds stress in a cross-stream plane in the $\xi{ -}\eta$ , barycentric and polar invariant maps. The points are coloured by the distance to the nearest wall, $d/H$ . Note that $d/H = 1$ corresponds to the centreline of the duct. Similar to the channel flow case, each invariant map shows that the turbulence starts out on the 1C–2C line near the wall, then moves towards the centre before terminating near the isotropic-one-component line. The movement towards the plane-strain limit, as well as the distance between the turbulent trajectories and the axisymmetric expansion line, are clearly shown in the barycentric and polar maps but are nearly imperceptible in the $\xi$ – $\eta$ mapping. This again shows the limitations of the $\xi$ – $\eta$ invariant map for visualisation purposes.

Figure 8.

A comparison of the square duct Reynolds stress visualised in each invariant map. The points are coloured by the wall distance, $d/H$ .

In contrast to the channel flow case, all of the turbulent trajectories start considerably closer to the one-component corner. This is highlighted in the polar invariant map where the channel flow turbulent trajectory initiates near the wall at approximately $\theta = \pi /7$ , whereas most points in the square duct never exceed $\theta =\pi /12$ . Additionally, the trajectories fan out in the range of approximately 0.1 to $0.4d/H$ , likely due to the profiles passing through different areas of the counter-rotating vortices.

To better visualise how figure 8 relates to physical space, contour maps of the invariant coordinates can be plotted. This is displayed in figure 9, overlaid with streamlines. Here, $C_1$ and $C_3$ from the barycentric map were plotted separately to make a fair comparison.

Figure 9.

Contour maps of invariant coordinates overlaid with streamlines for the square duct case. The scales for $\theta$ and $\xi$ have been reduced to half-range for clarity. (a) Polar map. left: $\alpha$ , right: $\theta$ , (b) Barycentric map. left: $C_3$ , right: $C_1$ , (c) $\xi$ – $\eta$ map. Left: $\eta$ , Right: $\xi$ .

As expected from (4.5), the $\alpha$ and $C_3$ contours are mirror images of each other, the only difference being the reversal of the colour map. The $C_1$ , $\eta$ and $\xi$ contours also show several similarities to the $\alpha$ contours. In each case, the variable is largest near the wall and decreases towards the centre. In addition, at any given distance from the wall, the points near the diagonal axes of symmetry show a lower magnitude. The invariants $\eta$ and $\xi$ display nearly identical contours as expected from the points forming a nearly straight line of unit slope in figure 8.

The area of lowest magnitude at the centre corresponds to the core flow, which is expected to be the region closest to isotropy, although the polar invariant map indicates the anisotropy magnitude remains approximately 15 % of the maximum value. The region of low magnitude along the symmetry planes is associated with the space between the corner vortices. This is an indication that the corner vortices are increasing the relative anisotropy magnitude.

The only variable that shows a significantly different distribution than the others is $\theta$ . In contrast to the cross-shape observed in the other contour maps, $\theta$ is only smaller along the axes of symmetry near the wall. Beyond the near-wall region, it increases nearly uniformly in a circular pattern, then decreases again in the core flow region. This reflects the movement of the anisotropy towards the plane strain line and then back towards the axisymmetric expansion line as seen in figure 8.

The reason for the difference between the contours of $\theta$ and the other invariant map coordinates is due to the interdependence of the bounds on the coordinates in other invariant maps. For example, as $\eta$ decreases, the range of $\xi$ also decreases. Likewise, as $C_3$ increases, the upper limit on $C_1$ and $C_2$ decreases to satisfy the additivity property. Consequently, a correlation between the contours of each invariant coordinate is expected. In contrast, the range of acceptable values of $\alpha$ and $\theta$ never change, as any combination of $\alpha \in [0,1]$ and $\theta \in [0,\pi /3]$ is allowed. As a result, the $\theta$ contours show the similarity of the anisotropy along the symmetry lines, which is lost in the other representations.

5.1. Current approaches to realisability

Most current approaches to realisability in ML turbulence modelling rely on the barycentric map, which is useful for interpolation between turbulent states due to its linearity (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007). The most direct usage of the barycentric map is by applying a correction to a baseline RANS model in barycentric coordinates, then transforming the new normalised anisotropy tensor back to the original coordinate system, as done by Tracey et al. (Reference Tracey, Duraisamy and Alonso2013) and Wang et al. (Reference Wang, Wu and Xiao2017).

A comprehensive discussion of the workflow for constructing and implementing models based on the barycentric map was given by Wang et al. (Reference Wang, Wu and Xiao2017) and Wu et al. (Reference Wu, Xiao and Paterson2018). The prediction process can be divided into several steps:

1. (i) run a baseline RANS simulation to convergence;

2. (ii) calculate eigenvalues and eigenvectors of the baseline solution;

3. (iii) convert eigenvalues to the Cartesian coordinates on the barycentric mapping;

4. (iv) use an ML model to predict the discrepancies in the eigenvalues, eigenvectors and turbulent kinetic energy, then add the discrepancies to the baseline solution;

5. (v) apply the eigenvectors as transformation matrices to revert to standard coordinates.

Outside of running the baseline simulation, steps (ii) and (v) are the most expensive. The cost in step (ii) comes from performing the eigendecomposition of the flow field, while the cost in step (v) comes from the use of transformation matrices.

As a method for ensuring an ML model makes realisable predictions, the barycentric map is highly effective, but it has limited applicability to realisability enforcement outside of eigendecomposition-based ML turbulence modelling approaches. Another common approach to ML turbulence modelling is based on Pope’s tensor basis formulation (Pope Reference Pope1975). Using the Cayley–Hamilton theorem, it can be shown that all eddy viscosity and algebraic Reynolds stress models can be represented as the weighted sum of a finite number of tensors. The tensor basis is given by

(5.1) \begin{equation} b_{\textit{ij}} = \sum _{i=1}^{10}{G^{(n)}T^{(n)}}, \end{equation}

where $T^{(n)}$ are basis tensors formed from the strain-rate and rotation-rate tensors, $S_{\textit{ij}}$ and $W_{\textit{ij}}$ , defined respectively as

(5.2) \begin{equation} S_{\textit{ij}} = \frac {1}{2}\left (\frac {\partial u_i}{\partial x_j} + \frac {\partial u_j}{\partial x_i} \right ), \quad W_{\textit{ij}} = \frac {1}{2}\left (\frac {\partial u_i}{\partial x_j} - \frac {\partial u_j}{\partial x_i} \right ). \end{equation}

Here, $G^{(n)}$ are scalar coefficients formed from the invariants of the basis tensors. The full list of definitions of tensors and invariants can be found from Pope (Reference Pope1975). Ling et al. (Reference Ling, Kurzawski and Templeton2016b ) was the first to propose using machine learning to predict the coefficients of the tensor basis as a way to embed invariance properties. The approach, which they refer to as the tensor basis neural network (TBNN), and variants using other ML algorithms have been widely used in ML turbulence modelling research (Kaandorp & Dwight Reference Kaandorp and Dwight2020; Mandler & Weigand Reference Mandler and Weigand2022; Cai et al. Reference Cai, Angeli, Martinez, Damblin and Lucor2024). In this case, there is no existing framework for creating models that inherently satisfy the realisability constraints, as purely invariant map-based approaches can. For tensor basis models, the need for baseline eigenvalues and eigenvectors limits current invariant maps to being used in post-processing corrections and training.

Blauw (Reference Blauw2019) and Mandler & Weigand (Reference Mandler and Weigand2022) used the barycentric map to shift non-realisable predictions to the boundaries of the barycentric map in post-processing. Both works used the same process for points below the 1C–2C line:

1. (i) calculate the eigendecomposition and convert to barycentric coordinates;

2. (ii) if the point is non-realisable, shift it to the 1C–2C boundary;

3. (iii) convert the new barycentric coordinates to eigenvalues.

4. (iv) transform back to standard coordinates using the eigenvectors.

Mandler & Weigand (Reference Mandler and Weigand2022) also included additional corrections for points above the 1C–2C line on the barycentric map, which they found were necessary during the convergence process. As discussed in § 3.1, this only occurs when the predicted normalised anisotropy tensor is not traceless.

Blauw (Reference Blauw2019) gives the new Cartesian coordinates on the barycentric map as

(5.3) \begin{align} x_{B,c} &= x_B + \frac {x_{3c}-x_{B}}{y_{3c}-y_B}(y_{1c}-y_B),\nonumber\\ y_{B,c} &= y_{1c}. \end{align}

Here, $x_{B,c}$ and $y_{B,c}$ are the corrected Cartesian coordinates on the barycentric map, and the original coordinates are $x_B$ and $y_B$ . The quantities $x_{1c}$ , $x_{3c}$ , $y_{1c}$ and $y_{3c}$ are the Cartesian coordinates of the corners of the barycentric map. This correction shifts the turbulent state along a line passing through the isotropic corner of the barycentric map.

Once the corrected coordinates are available, they are transformed to the corrected eigenvalues, $\lambda _{1,c} \geqslant \lambda _{2,c} \geqslant \lambda _{3,c}$ , by solving the system of equations,

(5.4) \begin{equation} \begin{bmatrix} x_{1c} & (-x_{1c}+2 x_{2c}) & (-2 x_{2c}+3 x_{3c}) \\ y_{1c} & (-y_{1c}+2 y_{2c}) & (-2 y_{2c}+3 y_{3c}) \\ 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} \lambda _{1,c} \\ \lambda _{2,c} \\ \lambda _{3,c} \end{bmatrix} = \begin{bmatrix} x_{B,c} - x_{3c} \\ y_{B,c} - y_{3c} \\ 0 \end{bmatrix}. \end{equation}

Finally, the matrix of eigenvalues is transformed back to standard coordinates using

(5.5) \begin{equation} b_{ij,c}=\boldsymbol{V}\boldsymbol{\varLambda }_c \boldsymbol{V}^{-1}, \end{equation}

where $b_{ij,c}$ is the corrected anisotropy tensor, $\boldsymbol{V}$ is the matrix of eigenvectors and $\boldsymbol{\varLambda }_c$ is the diagonal matrix of eigenvalues. This approach guarantees realisability, but does not affect the underlying model.

More recently, Riccius, Agrawal & Koutsourelakis (Reference Riccius, Agrawal and Koutsourelakis2023) and McConkey, Yee & Lien (Reference McConkey, Yee and Lien2024) have proposed realisability-informed training, inspired by PINNs. Both added an additional term in the loss function based on realisability constraints to penalise non-realisable predictions. This approach explicitly ‘teaches’ the model realisability constraints. As a representative example, McConkey et al. (Reference McConkey, Yee and Lien2024) used the equations,

(5.6) \begin{align} f_1 &= \max \left [(b_{\textit{ij}} - 2/3), -(b_{\textit{ij}} + 1/3), 0\right ], \quad i=j, \nonumber\\ f_2 &= \max \left [(b_{\textit{ij}} - 1/2), -(b_{\textit{ij}} + 1/2), 0\right ], \quad i \neq j, \nonumber\\ f_3 &= \max \left [(3|\lambda _2| - \lambda _2)/2 - \lambda _1, 0 \right ], \nonumber\\ f_4 &= \max \left [(\lambda _1 - (1/3 - \lambda _2)), 0 \right ]. \end{align}

The first two equations give the violation of the realisability conditions on each component, while the second two give the violation of the eigenvalue conditions described by Schumann (Reference Schumann1977). If a point is realisable, then the value is zero. These are combined in the penalty function,

(5.7) \begin{equation} \mathcal{L}_M=\gamma \Bigg \{\frac {1}{6}\Bigg [\sum _{i=j}{f_1^2} + \sum _{i \neq j}{f_2^2}\Bigg ] + \frac {1}{2}\left [f_3^2 + f_4^2 \right] \Bigg \}. \end{equation}

Here, $\mathcal{L}_M$ is the penalty function and $\gamma$ is a tuning parameter controlling the strength of the penalty. Although it does not guarantee realisable predictions, McConkey et al. (Reference McConkey, Yee and Lien2024) reported a substantial reduction in the number of non-realisable predictions. For example, they found that including the realisability-informed loss function reduced the number of non-realisable predictions in the square duct case from 9.2 % to only 0.4 %.

5.2. Polar invariant map strategies

The polar invariant map provides a framework that can simplify post-processing and realisability-informed training methods. Moreover, it can potentially enable the embedding of realisability constraints in tensor basis models. It allows realisability to be enforced through the use of $\alpha$ as a scaling factor, opening the door to new approaches. This is possible because of the representation of realisability as a limiting magnitude, which does not require an explicit eigendecomposition to calculate.

The simplest application for the polar invariant map is to shift non-realisable predictions to the nearest realisable state. In terms of the polar invariant map, the current method of shifting points to the 1C–2C boundary along a line connecting the isotropic corner is simply scaling along a line of constant $\theta$ . This is exactly the transformation described by

(5.8) \begin{equation} b_{ij,c} = \frac {1}{\alpha _0}b_{\textit{ij}}, \end{equation}

where $\alpha _0$ is the $\alpha$ value of $b_{\textit{ij}}$ . The only matrix operations required are a norm and a determinant for calculating $\alpha _0$ . Removing the need for a cumbersome process of eigendecomposition, solving a linear system and applying rotation matrices can be expected to yield a substantial reduction in computation time. A comparison of the calculation time for the two methods is described in the following subsection.

The application of the polar invariant map for realisability-informed training is also straightforward. In current approaches, the penalty term is formulated using the component-wise constraints and the boundaries on the eigenvalues. The eigenvalue constraints require two equations to account for the upper and lower bounds. Using the polar invariant map, the eigenvalue constraints can be simplified to a single equation penalising $\alpha \gt 1$ . Removing the need for an explicit eigenvalue calculation may yield substantial savings in training time over the thousands of iterations needed to train an ML model.

A more complex application is in turbulence models with adaptive coefficients. This avenue of research focuses on improving current turbulence models by using machine learning to tune the coefficients based on the local flow physics. Most of these are based on a linear eddy viscosity model using the Boussinesq hypothesis, which posits a linear relationship to the velocity gradients of the form (Pope Reference Pope2000)

(5.9) \begin{equation} \tau _{\textit{ij}} = -2\nu _t S_{\textit{ij}} + \frac {2}{3}k\delta _{\textit{ij}} = a_{\textit{ij}} + \frac {2}{3}k\delta _{\textit{ij}}. \end{equation}

In this case, the anisotropy is simply a scaling of the strain-rate tensor by the eddy viscosity, $\nu _t$ . Because $\theta$ is independent of magnitude, the maximum anisotropy magnitude can be determined from the strain-rate tensor alone at every point in the flow field. Then, the limit on $\nu _t$ can be easily determined as the value that produces $\alpha = 1$ ,

(5.10) \begin{equation} \left | \frac {\nu _{t,max}}{k} \right | = \frac {1}{\alpha _S}. \end{equation}

Here, $\alpha _S$ is the value of $\alpha$ when $S_{\textit{ij}}$ is substituted for $b_{\textit{ij}}$ in (3.12). Having an inexpensive way to calculate the local maximum $\nu _t$ may be useful for actively constraining the model coefficients of adaptive turbulence models.

The knowledge of the local maximum $\nu _t$ and model definition could enable the development of adaptive models that can actively constrain the predicted coefficients based on local conditions.

The polar invariant map’s potential applications are not limited to just simplifying existing approaches, but may also allow for incorporating realisability constraints into the structure of ML turbulence models. To see this, consider the output of an arbitrary turbulence model, $\mathcal{M}(\boldsymbol{x})$ . Using the polar invariant map, the output can be separated into a tensor, $\tilde {\mathcal{M}}(\boldsymbol{x})$ , with arbitrary magnitude, and a scaling in terms of $\alpha$ ,

(5.11) \begin{equation} b_{\textit{ij}} = \mathcal{M}(\boldsymbol{x})=\frac {\alpha _1}{\alpha _0}\mathcal{\tilde {M}}(\boldsymbol{x}). \end{equation}

Here, $\alpha _0$ is the $\alpha$ value of $\mathcal{M}(\boldsymbol{x})$ and $\alpha _1$ is the $\alpha$ value of $b_{\textit{ij}}$ . Equation (5.11) can be expanded to

(5.12) \begin{equation} b_{\textit{ij}} =\frac {\alpha _1}{\sqrt {6}||\mathcal{\tilde {M}}(\boldsymbol{x})||\cos \left(\theta -\dfrac {\pi }{3}\right)}\mathcal{\tilde {M}}(\boldsymbol{x})=\frac {\alpha _1}{\sqrt {6}\cos \left(\theta -\dfrac {\pi }{3}\right)}\frac {\mathcal{\tilde {M}}(\boldsymbol{x})}{||\mathcal{\tilde {M}}(\boldsymbol{x})||}. \end{equation}

Finally, (5.12) may be simplified by recognising the denominator as $1/||b_{\textit{ij}}||_{\textit{max}}$ from (3.11) to the form

(5.13) \begin{equation} b_{\textit{ij}} =\alpha _1||b_{\textit{ij}}||_{\textit{max}}\frac {\mathcal{\tilde {M}}(\boldsymbol{x})}{||\mathcal{\tilde {M}}(\boldsymbol{x})||}. \end{equation}

Equation (5.13) shows that the output is guaranteed to be realisable if and only if $\alpha _1\in [0,1]$ . The practical consequence is that realisability can theoretically be incorporated into the structure of an arbitrary model by predicting $\alpha _1\in [0,1]$ in a second prediction step.

A general way to accomplish this is to extract the coefficient of (5.12) as a scaling factor of the form

(5.14) \begin{equation} F(\mathcal{M(\boldsymbol{x})},\boldsymbol{y}) = \frac {\alpha _1(\boldsymbol{y})}{\sqrt {6}||\mathcal{\tilde {M}}(\boldsymbol{x})||\cos \left(\theta -\dfrac {\pi }{3}\right)}. \end{equation}

Here, $F(\mathcal{M(\boldsymbol{x})},\boldsymbol{y})$ is a function of the initial prediction $\mathcal{M(\boldsymbol{x})}$ and additional input variables $\boldsymbol{y}$ . The inputs, $\boldsymbol{y}$ , are used to predict $\alpha _1$ , and $M(\boldsymbol{x})$ is required to calculate $\theta$ and $||M(\boldsymbol{x})||$ .

Hypothetically, the above-mentioned strategy could be applied to any model architecture that can provide an initial prediction of the normalised anisotropy tensor with the correct $\theta$ . The resulting output will be realisable as long as $\alpha _1(\boldsymbol{y})$ remains between 0 and 1. Enforcing that range for $\alpha _1(\boldsymbol{y})$ can be easily accomplished by mapping the output to a sigmoid function – a common technique used in machine learning to predict classification probabilities. The usefulness of this approach would depend on the properties of $\mathcal{M}(\boldsymbol{x})$ ; for example, there would be no benefit to barycentric map-based models, which already enforce realisability. In contrast, it could be highly useful as a way to incorporate realisability constraints in the structure of a broad range of tensor basis models. The full range of benefits, drawbacks and details of implementation for such an approach remains a subject for further research.

5.2.1. Numerical considerations

All of the previously discussed potential applications for the invariant map are based on the assumption that $\theta$ is calculated from a symmetric, traceless tensor. These assumptions are based on the properties of the normalised anisotropy tensor, which, by definition, is symmetric and traceless. However, in practice, the traceless property of $b_{\textit{ij}}$ may not be satisfied while a CFD solution is converging, leading to violations of the domain of the arccosine in calculating $\theta$ . It is therefore important to consider how to handle these cases.

The reason for the violation of the traceless property of $b_{\textit{ij}}$ can be easily seen by considering a linear eddy viscosity model,

(5.15) \begin{equation} b_{\textit{ij}}=\frac {a_{\textit{ij}}}{2k}=-\frac {\nu _t}{k}S_{\textit{ij}}. \end{equation}

It is clear from (5.15) that an eddy viscosity model will only predict a traceless $b_{\textit{ij}}$ if $S_{\textit{ij}}$ is also traceless. Physically, this is guaranteed in an incompressible flow because the trace of $S_{\textit{ij}}$ is identically zero as a consequence of the conservation of mass (White & Majdalani Reference White and Majdalani2021). In a solver, however, there is no guarantee that this will be the case before a solution converges. In fact, violations of the conservation of mass are frequently used as a convergence criterion for this reason. If conservation of mass is not satisfied at a point, then $b_{\textit{ij}}$ will incorrectly have a trace. This reasoning also extends to the full tensor basis, which includes $S_{\textit{ij}}$ as the first basis tensor.

The simplest method of handling non-traceless input tensors with the polar invariant map is to use a modified strain-rate tensor formed by removing the trace,

(5.16) \begin{equation} S_{\textit{ij}}^* = S_{\textit{ij}}-\frac {1}{3}S_{kk}\delta _{\textit{ij}}. \end{equation}

Equation (5.16) only impacts the solution when a converging solution violates the conservation of mass. As the solution converges, $S_{\textit{ij}}^*$ converges to $S_{\textit{ij}}$ . By using $S_{\textit{ij}}^*$ to calculate the basis tensors, the zero trace property of $b_{\textit{ij}}$ can be ensured throughout the entire convergence process.

It is interesting to note that $S_{\textit{ij}}^*$ is mathematically identical to the traceless strain-rate tensor used by eddy viscosity models in compressible flows (Wilcox Reference Wilcox2006). Although $S_{\textit{ij}}^*$ is motivated by numerical artefacts, unlike the compressible flow case, where it is necessary for physical reasons, the same constitutive relation is successfully used by general-purpose CFD codes throughout both regimes (Pletcher, Tannehill & Anderson Reference Pletcher, Tannehill and Anderson2016). This may be taken as an indication that $S_{\textit{ij}}^*$ is a numerically feasible solution to the problem of non-traceless input tensors.

5.3. Post-processing corrections: sample calculations

We will now demonstrate the potential of the polar invariant for simplifying and reducing the time for calculations by making a comparison between the correction method of (5.8) and the formulation of Blauw (Reference Blauw2019). The comparison was conducted in Python on randomly generated, non-realisable, normalised anisotropy tensors. The methodology can be divided into four steps:

1. (i) randomly generate a set of $N$ symmetric, deviatoric matrices;

2. (ii) apply both correction methods sequentially to each point;

3. (iii) record the total calculation time;

4. (iv) repeat the calculations on the same set $M$ times and average the calculation times.

To calculate the percent reduction in calculation time, the following formula was used:

(5.17) \begin{equation} \textrm {Percent Reduction in Time = }100\boldsymbol{\cdot }\left (\frac {T_{\textit{bary}} - T_{\textit{polar}}}{T_{\textit{bary}}}\right ), \end{equation}

where $T_{\textit{bary}}$ is the time to correct $N$ tensors using the barycentric formulation and $T_{\textit{polar}}$ is the time to correct $N$ tensors using the polar formulation.

For the first step, the Numpy default_rng and random functions were used to generate random 3 × 3 matrices. These matrices were made symmetric by adding the transpose and dividing by two. They were then made deviatoric by subtracting $1/3$ of the trace from each element on the diagonal. This is depicted in Algorithm 2 in Appendix B.

To implement the barycentric corrections, the Numpy eigh function was used. The eigh function uses the LAPACK _syevd routine to calculate the eigenvalues and eigenvectors of a real, symmetric matrix. The eigenvalues were converted to barycentric coordinates according to (2.9). The barycentric coordinates were converted to Cartesian form via (4.6). Once the Cartesian coordinates on the barycentric map were obtained, the correction was calculated with (5.3), converted to eigenvalues through (5.4) and transformed back to standard coordinates via (5.5). All of these steps were combined into a single function. The pseudo-code for the barycentric correction is shown in Algorithm 3 in Appendix B.

The polar invariant corrections were implemented using (3.6b ) and (3.12) to calculate $\alpha _0$ . The norm of $b_{\textit{ij}}$ was calculated using the Numpy norm function. Once $\alpha _0$ was obtained, the correction was applied with (5.8). These steps were combined into a single function. The pseudo-code for the polar correction is shown in Algorithm 4 in Appendix B.

For the test, $N=10\,000$ random matrices were generated with a random seed of 0. The corrections were applied sequentially using separate for loops. The time required for each function evaluation was recorded using the Python time function in the time module. The percent difference in the total time relative to the barycentric formulation was calculated according to (5.17). To reduce the effect of random variations in the performance of the processor, an outer for loop was used to repeat these calculations $M=1000$ times. The percent differences for each run were recorded and averaged to produce the final time estimate. All calculations were performed on an Intel i9–13980HX processor.

Figure 10 shows the distribution of percent reduction in calculation time over the 1000 runs, and compares the initial and corrected points in the polar invariant map for a sample run. The histogram in figure 10(a) shows that the reduction in calculation time across the runs forms an approximately normal distribution. This is the expected result if there are no systemic biases in the variation of processor performance over time. Figure 10(b) confirms that the polar formulation correctly shifts all points to the 1C–2C line. Overall, the results of the test show the polar formulation is equivalent to the barycentric formulation, but with an average reduction in calculation time of 62 %.

Figure 10.

(a) A histogram of the percent reduction in computation time of the polar invariant formulation compared with the barycentric map formulation. The mean is 62 %. (b) A comparison of the initial and corrected points in the polar invariant map for one set of 10 000 random matrices.

5.4. Coefficient constraints: sample calculations

To illustrate how $\alpha$ – $\theta$ coordinates can be used to find the limiting magnitude of turbulence model coefficients, consider the simple case of the quadratic constitutive relation (QCR2000) proposed by Spalart (Reference Spalart2000),

(5.18a) \begin{align}& \tau _{ij,\textit{QCR}} = 2\nu _t S_{\textit{ij}} - C_{cr1}(O_{ik}\tau _{kj} - \tau _{ik}O_{kj}), \end{align}

(5.18b) \begin{align}&\qquad\qquad\quad O_{ik} = \frac {2W_{ik}}{\sqrt {\dfrac {\partial {u_m}}{\partial {x_n}}\dfrac {\partial {u_m}}{\partial {x_n}}}}. \end{align}

Here, $W_{ik}$ is the rate of rotation tensor, $\nu _t$ is the linear eddy viscosity and $C_{cr1}$ is a constant equal to 0.3. The term $\sqrt {({\partial {u_m}}/{\partial {x_n}})( {\partial {u_m}}/{\partial {x_n}}})$ is the velocity gradient magnitude, where $u_m$ is the mean velocity. For brevity, it will be represented as $||\boldsymbol{\nabla }U||$ from this point forward. As an example of how the polar invariant map can be applied to tensor basis models, we will demonstrate how it may be used to derive a limit on the eddy viscosity and $C_{\mu }$ in the $k{-}\varepsilon$ turbulence model.

In the case of QCR2000, the definition of $O_{ik}$ can be substituted back into (5.18a ) to convert it to the form

(5.19) \begin{equation} a_{\textit{ij}} = 2\nu _t S_{\textit{ij}} + \frac {4\nu _t C_{cr1}}{||\boldsymbol{\nabla }U||}(S_{ik} W_{kj} - W_{ik} S_{kj}). \end{equation}

Note that the contribution to $\tau _{\textit{ij}}$ from the turbulent kinetic energy can be added afterwards due to the anti-symmetric tensor in the equation (Spalart Reference Spalart2000). Normalising by $2k$ and factoring out $\nu _t/k$ yields

(5.20) \begin{equation} b_{\textit{ij}} = \frac {\nu _t}{k} \left [S_{\textit{ij}} + \frac {2C_{cr1}}{||\boldsymbol{\nabla }U||}\left (S_{ik} W_{kj} - W_{ik} S_{kj}\right )\right ]. \end{equation}

The terms in parentheses correspond to $T^{(n)}$ from (5.1). Specifically, $S_{\textit{ij}}=T^{(1)}$ and $ (S_{ik} W_{kj} - W_{ik} S_{kj} )=T^{(2)}$ are the first two basis tensors of Pope’s general eddy viscosity formulation (Pope Reference Pope1975).

Comparing (5.20) to (5.11), we can identify $\nu _t/k$ with $\alpha _1/\alpha _0$ and the term in brackets with $\tilde {\mathcal{M}}(\boldsymbol{x})$ . Setting these equal yields

(5.21) \begin{equation} \frac {\nu _t}{k}=\frac {\alpha _1}{\sqrt {6}||\tilde {\mathcal{M}}(\boldsymbol{x})||\cos \left(\theta -\dfrac {\pi }{3}\right)}. \end{equation}

The maximum realisable eddy viscosity occurs when $\alpha _1=1$ , so the limit on the eddy viscosity is

(5.22) \begin{equation} \left |\frac {\nu _t}{k}\right | \leqslant \frac {1}{\sqrt {6}||T^{(1)} + \dfrac {2C_{cr1}}{||\boldsymbol{\nabla }U||}T^{(2)}||\cos {\left(\theta - \dfrac {\pi }{3}\right)}}. \end{equation}

This is an equation that can be solved inexpensively and analytically for any $T^{(1)}$ and $T^{(2)}$ . With the local maximum $\nu _t$ available, it becomes possible to use a model that predicts $\nu _t$ relative to the maximum realisable magnitude.

Next, consider an eddy viscosity as prescribed by the $k{-}\varepsilon$ model (White & Majdalani Reference White and Majdalani2021),

(5.23) \begin{equation} \nu _t = C_{\mu } \frac {k^2}{\varepsilon }. \end{equation}

Here, $C_\mu$ is a proportionality constant. Substituting (5.23) into (5.22) allows a limit to be placed on $C_\mu$ according to

(5.24) \begin{equation} |C_\mu | \leqslant \frac {\sqrt {6}}{6(k/\varepsilon )||T^{(1)} + \dfrac {2C_{cr1}}{||\boldsymbol{\nabla }U||}T^{(2)}||\cos {\left(\theta - \dfrac {\pi }{3}\right)}}. \end{equation}

Equation (5.24) relates the local maximum realisable magnitude to the proportionality coefficient in the $k$ – $\varepsilon$ model. This form of the limiting magnitude would be useful to a turbulence model with adaptive coefficients. Additional research would be necessary to determine how other coefficients could be related to the limiting magnitude.

Consider the case of a high-Reynolds-number turbulent round jet as an example of the limitation on $C_\mu$ . The flow is axisymmetric and non-swirling, so the only non-zero mean velocity components are the axial velocity, $u$ , and the radial velocity, $v$ . By order of magnitude arguments, $v\lt \lt u$ and $\partial /\partial x \lt \lt \partial /\partial r$ (Pope Reference Pope2000). Then, for (5.19), $S_{\textit{ij}}$ and $W_{\textit{ij}}$ can be approximated as having only the components $S_{12} = S_{21} = u_r/2$ and $W_{12} = -W_{21} = u_r/2$ . In this case, $T^{(1)}$ and $T^{(2)}$ are

(5.25) \begin{equation} T^{(1)}=\begin{bmatrix} 0 & u_r/2 & 0 \\[4pt] u_r/2 & 0 & 0 \\[4pt] 0 & 0 & 0 \end{bmatrix}, \quad T^{(2)}=\begin{bmatrix} -u_r^2/2 & 0 & 0 \\[4pt] 0 & u_r^2/2 & 0 \\[4pt] 0 & 0 & 0 \end{bmatrix}. \end{equation}

The velocity gradient magnitude is $|u_r|$ , so $T_{\textit{ij}}$ becomes

(5.26) \begin{equation} T_{\textit{ij}} = \left (T^{(1)} + \frac {2C_{cr1}}{|u_r|}T^{(2)}\right ) = \begin{bmatrix} -\left(u_r^2/|u_r|\right)C_{cr1} & u_r/2 & 0 \\[4pt] u_r/2 & \left(u_r^2/|u_r|\right)C_{cr1} & 0 \\[4pt] 0 & 0 & 0 \end{bmatrix}. \end{equation}

Two simplifications can be made to (5.26). First, $u_r^2/|u_r|$ is strictly non-negative and can simply be replaced with $|u_r|$ . Additionally, it is well known that $u_r$ is strictly non-positive away from the edge of a round jet, so in this specific case, $-u_r$ is strictly positive, so $|u_r|=-u_r$ . Making these simplifications yields

(5.27) \begin{equation} T_{\textit{ij}} = \left (T^{(1)} + \frac {2C_{cr1}}{|u_r|}T^{(2)}\right ) = \begin{bmatrix} u_rC_{cr1} & u_r/2 & 0 \\[4pt] u_r/2 & -u_rC_{cr1} & 0 \\[4pt] 0 & 0 & 0 \end{bmatrix}= u_r \begin{bmatrix} 3/10 & 1/2 & 0 \\[4pt] 1/2 & -3/10 & 0 \\[4pt] 0 & 0 & 0 \end{bmatrix}. \end{equation}

The single zero-column indicates this is a case of plane-strain, so $\theta = \pi /6$ and $\cos {(\theta - \pi /3)} = \sqrt {3}/2$ . The magnitude of $T_{\textit{ij}}$ is $\sqrt {17}/5$ . Substituting these values into (5.22) gives a limit on the eddy viscosity of

(5.28) \begin{equation} \left |\frac {\nu _t}{k}\right | \leqslant \frac {5\sqrt {34}}{51 |u_r|} \approx \frac {0.572}{|u_r|}. \end{equation}

For the $k{-}\varepsilon$ model, this implies the limit on $C_{\mu }$ in a round jet using the quadratic constitutive relation is

(5.29) \begin{equation} |C_{\mu }| \leqslant \frac {5\sqrt {34}}{51 |u_r| (k /\varepsilon )} \approx \frac {0.572}{|u_r| (k /\varepsilon )}. \end{equation}

In general, (5.21) must be solved at each point in the flow field due to differences in the relative weights, $G^{(n)}$ , of (5.1). The above-mentioned example demonstrates the type of analysis involved in applying the polar invariant map to non-linear eddy viscosity models.