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:
- \(x\) = streamwise direction (direction of the jet exit)
- \(y\) = cross-stream (lateral) direction (perpendicular to the jet centerline)
- \(z\) = spanwise direction (homogeneous direction, infinite extent)
- \(U\) = mean velocity in streamwise (\(x\)) direction
- \(V\) = mean velocity in cross-stream (\(y\)) direction
- \(W\) = mean velocity in spanwise (\(z\)) direction
- \(u’\) = velocity fluctuation in streamwise direction
- \(v’\) = velocity fluctuation in cross-stream direction
- \(w’\) = velocity fluctuation in spanwise direction
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:
- \(C_1 = 1\): one-component turbulence (rod-like)
- \(C_2 = 1\): two-component turbulence (disk-like)
- \(C_3 = 1\): three-component isotropic turbulence
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 |
|---|---|---|---|---|
| DNS | 0.18 | 0.12 | 0.70 | Interior of triangle |
| LEVM | 0.00 | 0.00 | 1.00 | Isotropic 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.0 | 0.5 | 1.0 | 1.5 | 2.0 |
|---|---|---|---|---|---|
| Error | 0.17 | 0.12 | 0.09 | 0.08 | 0.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:
- Parametric error:Â Wrong constants, wrong inlet conditions, wrong grid resolution. These can be fixed by calibration.
- Structural error:Â The constitutive relation itself cannot represent the required physics. No amount of tuning helps.
In this case:
- Changing the grid does nothing
- Changing inlet turbulence does nothing
- Switching from \(k\)–\(\epsilon\) to \(k\)–\(\omega\) SST does nothing
- Adjusting model constants does nothing
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:
- Transport:Â Turbulent kinetic energy and stresses are convected from production regions
- Finite redistribution rate:Â Pressure-strain correlation operates on a turbulent time scale comparable to the mean flow time scale
- History effects:Â Turbulence retains memory of its anisotropic generation
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:
- Turbulent kinetic energy distribution:Â Affects mixing, entrainment, and shear layer growth
- Scalar transport:Â For a passive scalar with source at the centerline, the scalar variance profile would differ by approximately 20% between DNS and LEVM predictions
- Aeroacoustics:Â Sound generation depends on two-point correlations of velocity fluctuations
- Conjugate heat transfer:Â Wall heat flux depends on turbulent Prandtl number assumptions that cascade from stress anisotropy
- Reynolds-stress-sensitive combustion models:Â Chemical source terms depend on scalar dissipation rates that couple to anisotropy
For my current work on film cooling—where scalar transport depends on accurate stress distributions—LEVMs are structurally insufficient. For applications requiring only integral quantities, they may be perfectly adequate.
The Model Selection Framework This Suggests
This experience clarified for me that model adequacy must be evaluated relative to the quantity of interest:
| Application | LEVM Adequacy | Why |
|---|---|---|
| Spreading rate prediction | âś“ Adequate | Integral quantity, insensitive to anisotropy |
| Mean velocity profiles | âś“ Adequate | Constrained by momentum integral |
| Centerline turbulence intensity | âś— Inadequate | Directly mispredicted by factor 1.5 |
| Scalar mixing rate | ⚠️ Context-dependent | Error propagates through turbulent diffusivity |
| Aeroacoustics | âś— Inadequate | Sensitive to two-point correlations |
| Combustion with strained flames | âś— Inadequate | Scalar 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:
- Don’t stop at integral quantities
- Project results into invariant space
- Ask: “Is my predicted turbulence state even realizable under my closure?”
For model selection:
- Map the quantity of interest’s sensitivity to anisotropy
- If it’s insensitive, LEVM might be fine
- If it’s sensitive, consider EARSM or full RSM
For diagnosing structural error:
- If the error persists across grids, schemes, and model variants—it’s structural
- If the error occurs where constitutive assumptions become degenerate—it’s structural
- If no amount of calibration fixes it—it’s structural
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
Banerjee, S., Krahl, R., Durst, F., & Zenger, C. (2007). Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches. Journal of Turbulence, 8, N32.
Boussinesq, J. (1877). Essai sur la thĂ©orie des eaux courantes. MĂ©moires prĂ©sentĂ©s par divers savants Ă l’AcadĂ©mie des Sciences, 23(1), 1–680.
Bradbury, L. J. S. (1965). The structure of a self-preserving turbulent plane jet. Journal of Fluid Mechanics, 23(1), 31–64.
Lumley, J. L., & Newman, G. R. (1977). The return to isotropy of homogeneous turbulence. Journal of Fluid Mechanics, 82(1), 161–178.
Rosti, M. E., Soligo, G., & Chiarini, A. (2025). Reynolds number effect on the flow statistics and turbulent-non-turbulent interface of a planar jet. Journal of Fluid Mechanics, in press.
Stanley, S. A., Sarkar, S., & Mellado, J. P. (2002). A study of the flow-field evolution and mixing in a planar turbulent jet using direct numerical simulation. Journal of Fluid Mechanics, 450, 377–407.
A natural extension of this work is to quantify how anisotropy misprediction propagates into scalar flux modeling (\(\overline{u_i’\theta’}\)), particularly in configurations where Reynolds stress–scalar gradient alignment is critical for film cooling effectiveness. This remains an open question in the turbulence modeling community.
Discussion:
- Have you encountered structural errors in your own simulations that integral metrics couldn’t reveal?
- How do you determine whether a model’s limitations matter for your specific application?
- What role should invariant analysis play in standard CFD validation practice?
Share your experiences below!
