Home / Blogs / Numerical Methods and Stability

Structural Limits of Linear Eddy-Viscosity Models in the Planar Jet

This article examines the validation of planar turbulent jet simulations, highlighting how key metrics like jet spreading and decay laws align with direct numerical simulation data. By achieving low residuals and mass balance stability, the study prompts a discussion on identifying structural errors, model limitations, and the essential role of invariant analysis in standard CFD validation practices.

Share 𝕏 in
Structural Limits of Linear Eddy-Viscosity Models in the Planar Jet

The jet spreads correctly. The decay law matches DNS. The residuals drop six orders. Mass imbalance falls below \(10^{-5}\). And yet the turbulence state at the centerline is mathematically impossible under the closure I’m using.

This isn’t a calibration problem. It’s not a grid issue. It’s not even a turbulence model problem in the way we normally think about them. It’s structural. And it took me five months to understand the difference.

The Paradox

For the past five months, I’ve been running RANS simulations of planar turbulent jets, trying to validate against the Bradbury (1965) experiment at Reynolds number 30,000. The simulations converge cleanly, everything looks textbook.

The spreading rate, that canonical integral metric that every turbulence modeler checks first, comes out to:

\(\frac{d\delta}{dx} = 0.104\)

where \(\delta\) is the jet half-width measured from the centerline. This is within 1% of the DNS value reported by Stanley, Sarkar & Mellado (2002) and within Bradbury’s experimental scatter. By any conventional validation standard, this simulation is a success.

But here’s the problem. When I extract the Reynolds stress tensor at the centerline, where the mean velocity gradient \(\partial U/\partial y = 0\), the model predicts isotropic turbulence:

\(\overline{u’^2} = \overline{v’^2} = \overline{w’^2}\)

The DNS shows something completely different:

\(\frac{\overline{u’^2}}{\overline{v’^2}} \approx 1.47\)

That’s a rod-like anisotropy state: one eigenvalue positive, two equal negatives. The turbulence at the centerline isn’t isotropic. It can’t be isotropic, physically. But my model forces it to be.

And here’s the scary part: this discrepancy persists across \(k\)–\(\epsilon\), across \(k\)–\(\omega\) SST, across grid refinements, across inlet condition variations. It’s not going away.

For a planar jet:

The Mathematical Trap

The Anisotropy Tensor

The Reynolds stress anisotropy tensor is defined as:

\(b_{ij} = \frac{\overline{u_i’ u_j’}}{2k} – \frac{1}{3}\delta_{ij}\)

where \(k = \overline{u_i’u_i’}/2\) is the turbulent kinetic energy. This tensor is symmetric and traceless—it captures exactly how the turbulence deviates from isotropy.

All realizable turbulence states live within a bounded region in invariant space (Lumley & Newman, 1977). The barycentric map (Banerjee et al., 2007) gives us a clean visualization:

\(C_1 = \lambda_1 – \lambda_2,\quad C_2 = 2(\lambda_2 – \lambda_3),\quad C_3 = 3\lambda_3 + 1\)

where \(\lambda_1 \ge \lambda_2 \ge \lambda_3\) are the eigenvalues of \(b_{ij}\). The vertices correspond to:

screenshot 2026 03 04 132139
Barycentric Map introduced by Banerjee et al. (2007)

The Linear Eddy-Viscosity Constraint

For an incompressible planar jet, the mean velocity field is \(U(x, y)\) with \(W=0\). The continuity equation gives:

\(\frac{\partial U}{\partial x} + \frac{\partial V}{\partial y} = 0\)

The strain rate tensor is:

\[
S_{ij}=
\frac{1}{2}
\begin{pmatrix}
2\frac{\partial U}{\partial x} &
\frac{\partial U}{\partial y} + \frac{\partial V}{\partial x} &
0 \\
\frac{\partial U}{\partial y} + \frac{\partial V}{\partial x} &
2\frac{\partial V}{\partial y} &
0 \\
0 & 0 & 0
\end{pmatrix}
\]

This tensor has eigenvalues \((+\lambda, 0, -\lambda)\)—one positive, one negative, one zero.

The linear eddy-viscosity constitutive relation (Boussinesq, 1877) states:

\(b_{ij} = -\frac{\nu_t}{k} S_{ij}\)

Here’s the trap: Because \(\nu_t\) is a scalar, \(b_{ij}\) is forced to have exactly the same eigenstructure as \(S_{ij}\): one positive eigenvalue, one negative, one zero.

And at the centerline, where \(\partial U/\partial z = 0\) and the mean strain vanishes, \(S_{ij} = 0\) gives:

\(b_{ij} = 0 \quad \text{identically}\)

This isn’t an approximation. It’s a mathematical constraint baked into the constitutive assumption itself.

What the DNS Actually Shows

DNS of the planar jet at Re = 30,000 (Stanley et al., 2002; Rosti et al., 2025) gives these anisotropy eigenvalues at the centerline:

\(\lambda_1 = 0.12,\quad \lambda_2 = -0.06,\quad \lambda_3 = -0.06\)

The barycentric coordinates are:

\(C_1 = 0.18,\quad C_2 = 0.12,\quad C_3 = 0.70\)

This state lies inside the realizability triangle, away from the plane-strain line and distinct from the isotropic vertex. The turbulence at the centerline remembers its anisotropic history—it’s been transported from the shear layer, and pressure-strain redistribution hasn’t fully isotropized it.

The Structural Error Quantified

Centerline Comparison

Source\(C_1\)\(C_2\)\(C_3\)State
DNS0.180.120.70Interior of triangle
LEVM0.000.001.00Isotropic vertex

The Frobenius norm of the anisotropy tensor error:

\(|b_{ij}^{\text{LEVM}} – b_{ij}^{\text{DNS}}|_F = 0.17\)

For context, the maximum possible Frobenius norm (at the one-component limit) is \(\sqrt{2/3} \approx 0.82\).

Across the Jet

\(y/\delta\)0.00.51.01.52.0
Error0.170.120.090.080.10

The error is not random. It is maximal precisely where the closure’s constitutive assumption becomes degenerate.

What This Looks Like in Invariant Space

The DNS trajectory traverses the interior of the barycentric triangle, reaching a rod-like state (\(C_1 > C_2\), \(C_3 < 0.8\)) in the shear layer and retaining significant anisotropy at the centerline.

The LEVM trajectory is confined to the plane-strain line in the shear layer, unable to access states with simultaneous \(C_1 > 0\) and \(C_2 > 0\), and collapses to the isotropic vertex at the centerline.

The geometric restriction is absolute: linear eddy-viscosity models cannot represent turbulence states with \(C_2 > 0\) and \(C_1 > 0\) simultaneously away from the plane-strain line, regardless of calibration.

Structural Error Cannot Be Tuned Away

This is the critical distinction I had to learn:

In this case:

Because the closure constrains \(b_{ij} \propto S_{ij}\). When \(S_{ij} = 0\), the model must predict isotropy. The physics doesn’t care.

Why Does DNS Show Anisotropy at the Centerline?

This is where understanding turbulence physics matters more than tensor algebra.

The centerline isn’t strain-free in a Lagrangian sense. Turbulence there has been advected from the shear layer, where it was highly anisotropic. The pressure-strain redistribution term—the mechanism that returns turbulence toward isotropy, acts on a finite time scale. By the time fluid parcels reach the centerline, they haven’t fully isotropized.

The anisotropy persists because:

Linear eddy-viscosity models have no mechanism to represent this. They assume local equilibrium—that turbulence adjusts instantaneously to the local mean strain. At the centerline, where mean strain vanishes, that assumption forces isotropy.

When Does This Error Actually Matter?

Here’s the uncomfortable truth: my spreading rate matched DNS to within 1%. The velocity profiles collapsed perfectly. By every integral metric, the simulation succeeded.

The spreading rate is remarkably insensitive to centerline anisotropy. Both linear and quadratic EVMs predict spreading rates within 2% of DNS (0.0993 vs 0.098–0.100), despite factor-of-two errors in \(b_{11}\) at the centerline.

So why should anyone care?

Because other quantities depend directly on the Reynolds stress tensor:

The Model Selection Framework This Suggests

This experience clarified for me that model adequacy must be evaluated relative to the quantity of interest:

ApplicationLEVM AdequacyWhy
Spreading rate prediction✓ AdequateIntegral quantity, insensitive to anisotropy
Mean velocity profiles✓ AdequateConstrained by momentum integral
Centerline turbulence intensity✗ InadequateDirectly mispredicted by factor 1.5
Scalar mixing rate⚠️ Context-dependentError propagates through turbulent diffusivity
Aeroacoustics✗ InadequateSensitive to two-point correlations
Combustion with strained flames✗ InadequateScalar dissipation couples to anisotropy

The geometric framing—barycentric coordinates, invariant maps, admissible manifolds, provides the language to articulate this distinction precisely.

What I’d Do Differently Now

For validation:

For model selection:

For diagnosing structural error:

What This Teaches Us About CFD:

Model adequacy is application-dependent. The right question isn’t “is this model good?” but “is this model adequate for my quantity of interest?”

Integral metrics can lie. A perfectly matched spreading rate doesn’t mean your turbulence state is physical.

Structural error is different from parametric error. If it persists across models, grids, and schemes, it’s in the constitutive assumption.

Invariant analysis should be standard validation practice. Project your results into barycentric space. Ask where your model lives relative to DNS.

References

Some more papers of interest


Discussion:

  1. Have you encountered structural errors in your own simulations that integral metrics couldn’t reveal?
  2. How do you determine whether a model’s limitations matter for your specific application?
  3. What role should invariant analysis play in standard CFD validation practice?

Share your experiences below!

Leave a Reply

Your email address will not be published. Required fields are marked *