The constant that makes a variance negative

What averaging costs is the statement of the problem: averaging the momentum equation leaves a term nobody has an equation for, the Reynolds stress, and the ladder that never closes is what happens to anybody who tries to chase it — each new equation contains a new unknown, for ever.

So the stress has to be modelled, and the model every engineering calculation in the world uses is Boussinesq’s:

$$\langle u_iu_j\rangle = \tfrac23 k\,\delta_{ij} - 2\nu_t S_{ij}, \qquad \nu_t = C_\mu\frac{k^2}{\varepsilon},$$

with $C_\mu = 0.09$.

One number, standing for the entire second-order statistics of a turbulent flow. This essay is about what that number is actually saying.

What the calibration is #

$C_\mu$ was not measured directly. It was inferred from a shear layer in local equilibrium, where production equals dissipation.

With the closure, $P = \nu_t S^2$, so $P = \varepsilon$ gives

$$C_\mu = \left(\frac{\varepsilon}{Sk}\right)^2 \quad\Longleftrightarrow\quad \frac{Sk}{\varepsilon} = \frac{1}{\sqrt{C_\mu}} = \frac{10}{3}.$$

So $C_\mu = 0.09$ is the assertion that $Sk/\varepsilon$ is ten thirds. That is a measurement of one flow, written as a universal constant, and the computed check confirms that production equals dissipation exactly at that strain.

Everything about the model’s behaviour away from that flow follows from what $Sk/\varepsilon$ does there.

There is a second thing hiding in that derivation and it is worth extracting.

$Sk/\varepsilon$ is the strain rate measured on the turbulence’s own time scale $k/\varepsilon$, which is the large-eddy turnover time. So the ratio is asking: does the mean flow deform an eddy appreciably in the time the eddy takes to turn over?

At $10/3$ the answer is “somewhat” — the eddies are being distorted, and a steady state is reached in which the distortion is balanced by the turbulence’s own decorrelation. That balance is what “local equilibrium” means, and it is a special state rather than a general one.

Below the ratio the turbulence has time to forget the strain between distortions and the stress is nearly isotropic. Far above it the turbulence has no time at all: the strain wins, the eddies are stretched into whatever shape the mean flow demands, and the whole notion of a stress proportional to the instantaneous strain is describing a process that has not had time to happen.

The last case is called rapid distortion, it has an exact linear theory of its own, and its answer is not proportional to the strain rate. It is proportional to the accumulated strain, which is a different quantity and does not appear anywhere in the closure.

Where it stops being able to describe a fluid #

Realizability is the requirement that a variance be non-negative — the weakest possible demand on a second-order statistic, and one no measurement can violate.

In a plane strain the closure gives

$$\langle u'^2\rangle = \tfrac23 k - 2C_\mu\frac{k^2}{\varepsilon}S,$$

which is negative as soon as

$$\frac{Sk}{\varepsilon} > \frac{1}{3C_\mu} = 3.704.$$

Three point seven zero four against a calibration at three point three three three: the bound is 11.1 per cent above the point the constant was fixed at, and the computation confirms the crossing is exactly there, positive one per cent below and negative one per cent above.

That is not an edge case at the end of a long extrapolation. It is a tenth of the way past the calibration.

And the strains a real calculation meets #

The strains that matter are far past it.

Ahead of any blunt body there is a stagnation region where the flow decelerates and turns, and $Sk/\varepsilon$ there reaches ten to fifty. A contraction, an impinging jet, the leading edge of an aerofoil at incidence — all of them are in the same range.

At $Sk/\varepsilon = 15$ the closure gives $\langle u'^2\rangle = -2.03$ against a total kinetic energy of 1. That is not a large error in a quantity; it is a quantity of the wrong sign, in a model that is being asked to predict a heat transfer coefficient.

It is worth noting what does not go wrong at the same time, because the failure’s invisibility depends on it.

The turbulent kinetic energy $k$ stays positive: it is a solved variable with its own transport equation, and nothing in that equation can drive it negative. The trace of the modelled stress tensor is $2k$ by construction, so it is positive too. The eddy viscosity is positive. The momentum equation is well posed and its solution is smooth.

Only the individual normal stresses go negative, and only one of the three, and nobody plots them. The model reports a perfectly plausible flow field with a physically impossible statistic underneath it, and every diagnostic a user is likely to look at is fine.

That combination — an impossible internal quantity and an entirely healthy output — is the reason the anomaly was diagnosed from its consequences rather than from its cause.

What actually goes wrong, which is the production #

The symptom everybody knows is not the negative stress — nobody looks at the individual stresses — but its consequence.

$$\frac{P}{\varepsilon} = 2C_\mu\left(\frac{Sk}{\varepsilon}\right)^2,$$

which with a constant $C_\mu$ grows as the square of the strain, without bound. The exponent is measured back out at exactly 2.

At $Sk/\varepsilon = 60$ it is 648. The model reports the turbulence generating energy six hundred times faster than it is destroying it, $k$ grows enormously in the stagnation region, and everything downstream — the boundary layer, the separation point, the heat transfer — is computed from a turbulent energy that is an artefact.

That is the stagnation-point anomaly, and it is one of the best-known failures in engineering turbulence modelling. Its cause is one line of algebra.

What the runaway does downstream #

The chain from the algebra to a wrong answer is worth following once, because it explains why this particular failure has a name.

In the stagnation region the production is enormous, so $k$ grows to a value far above anything measured. That inflated $k$ is convected downstream, into the boundary layer that forms on the body’s surface, where it raises the eddy viscosity.

A boundary layer with too much eddy viscosity is too full: its profile is too resistant to an adverse pressure gradient, it separates too late or not at all, and its heat transfer is too high. How much uphill a layer can take is the essay about the quantity being got wrong there, and the error arrives from an entirely different part of the flow.

Measured heat transfer at the stagnation point of a cylinder in a turbulent stream is over-predicted by a factor of two or three by an uncorrected model. That is the number that made the problem famous, and it is downstream — literally — of a normal stress nobody looked at.

The fix, derived rather than borrowed #

The bound says what has to happen: $C_\mu$ must fall as $1/(3\eta)$ at large $\eta = Sk/\varepsilon$, or the stress goes negative.

So take $C_\mu = 1/(A_0 + 3\eta)$ and fix $A_0$ by requiring the equilibrium value:

$$\frac{1}{A_0 + 3 \times \tfrac{10}{3}} = 0.09 \quad\Longrightarrow\quad A_0 = \frac{1}{0.09} - 10 = 1.111.$$

The normal stress is then $\tfrac23 - 2\eta/(A_0 + 3\eta)$, which is positive at every strain and tends to zero — marginally realizable in the limit, which is the most a model of this form can be. The measured exponent of $C_\mu$ against $\eta$ at large strain is $-0.991$.

That is Shih’s realizable $k$–$\varepsilon$ model, arrived at from the bound rather than from its author’s algebra. The point of deriving it that way is that the form is forced: any model that keeps a variance positive at large strain has $C_\mu \propto 1/\eta$ there.

The limit being taken, which is a small one #

The reason this belongs beside a run of essays about limits is that Boussinesq’s closure is a limit, and it is easy to forget.

Writing the stress as a linear function of the strain rate is the leading term of an expansion in the strain measured on the turbulence’s own time scale — that is, in $Sk/\varepsilon$. Every term after it has been dropped: the quadratic ones, which produce normal-stress anisotropy in a simple shear; the history terms, which carry the fact that turbulence takes time to respond; and the whole of the alignment between the stress and the strain, which is assumed rather than computed.

The expansion parameter is $Sk/\varepsilon$, the calibration is at $10/3$, and the model is used at 50.

So the residue of this limit is not subtle. It is every term that was dropped, at a parameter value where the expansion has no business being truncated — and the dropped terms carry a definite sign, because what they do is stop a variance going negative.

What a second-moment closure does instead #

The alternative is to abandon the eddy viscosity and solve transport equations for the stresses themselves, which is what a Reynolds-stress model does.

That does not close the problem — the ladder that never closes still applies, and the stress equations contain triple correlations and a pressure-strain term that have to be modelled in turn — but it moves the modelling to a place where realizability can be imposed directly, because the stresses are the variables rather than being reconstructed from a strain.

The price is six equations instead of one, stiffness, and a pressure-strain model with rather more constants in it than $C_\mu$. Which is why, thirty years after the problem was understood, most industrial computations still use a two-equation model with a corrected $C_\mu$: the correction fixes the specific failure at a cost of nothing.

Why the failure is invisible in the usual tests #

It is worth asking why a model with a hole this large survives, and the answer is the same as the answer to most such questions in this subject.

The flows a model is validated on are the flows the constants were fitted to. A boundary layer, a channel, a free shear layer, a pipe — all of them are near-equilibrium shear flows with $Sk/\varepsilon$ within a factor of two of $10/3$. Over that range the closure is a perfectly good interpolation and every prediction it makes is right.

The failures appear in flows with normal strain: stagnation, impingement, contraction, strong curvature. Those are less often used for validation, they are harder to measure, and their failure mode — an over-predicted turbulent energy — shows up as a wrong heat transfer or a wrong separation point rather than as an obviously absurd number.

What a code says to a wall is the same shape of finding one level down: a wall function is exact between $y^+ = 30$ and $10{,}000$, sixty per cent wrong at $y^+ = 1$, and the sixty per cent is not visible in the output either.

What a constant is, in this subject #

There is a general point here that this collection keeps arriving at, and it is worth stating plainly.

$\kappa = 0.41$ is a constant of a dimensional argument: the overlap region has no length in it, so the gradient is $u_\tau/\kappa y$, and $\kappa$ is the coefficient of a relation that follows from what is forbidden. The layer with no length in it is careful to say that it has never been derived from anything.

$C_\mu = 0.09$ is not that. It is the value of a ratio in one flow, promoted to a constant by being written as a coefficient. There is no argument that says $Sk/\varepsilon$ should be the same in every flow — it is not, and the whole of this essay is about where.

Telling the two kinds apart is most of the skill in reading a turbulence model. A constant whose dimensional argument forbids the alternative is on solid ground; a constant that is a measured ratio under a different name is a calibration with a domain of validity.

The arithmetic that makes it memorable #

Three numbers, and the relation between them is the whole essay.

$1/\sqrt{C_\mu} = 3.333$ is where the model was calibrated. $1/(3C_\mu) = 3.704$ is where it stops describing a fluid. Their ratio is $\sqrt{C_\mu}/(3C_\mu)\cdot\ldots$ — more simply, $\tfrac{1}{3C_\mu} \big/ \tfrac{1}{\sqrt{C_\mu}} = \tfrac{1}{3\sqrt{C_\mu}} = 1.111$ — exactly ten ninths, independent of the value of $C_\mu$.

The margin is always eleven per cent, whatever the constant is. Choosing a different $C_\mu$ moves both the calibration and the bound and does not widen the gap between them. That is worth knowing, because it says the problem cannot be tuned away: any model of this form fails at a tenth past its own calibration point, and only changing the form helps.

What the same failure looks like in a wall model #

The pattern here — a relation calibrated in one state and applied in another — recurs one level down, and the collection has the case already.

What a code says to a wall is about a wall function: the law of the wall applied as a boundary condition at one grid point, on the assumption that the point lies in the overlap region. Where it does, the answer is exact. Where it does not — at $y^+ = 1$, which is where a fine grid puts the first point — the friction is sixty per cent wrong, and refining the grid makes it worse rather than better.

Both failures have the same shape: an asymptotic relation with a domain of validity, used outside it, with no diagnostic firing. And both have the same fix, which is to blend the relation with something valid in the other regime and to accept a model that is approximate everywhere rather than exact somewhere.

That trade — exact in one place and wrong elsewhere, against approximate everywhere — is the central choice in engineering turbulence modelling, and it is worth recognising as a choice rather than as a technical detail.

Limits recorded rather than smoothed over #

Only the plane strain is examined. In a simple shear the closure’s normal stresses are all equal to $\tfrac23 k$, which is wrong — measured anisotropies are large — but not negative, so realizability is not violated there. The bound computed here is the plane-strain one, which is the relevant case for the stagnation-point anomaly.

The realizable form is one of several. Shih’s is the one derived here; there are others, including simply limiting the production, and they differ in the transitional range. All of them agree that $C_\mu$ must fall as $1/\eta$ at large strain, because the bound says so.

Nothing here is a computation of a flow. The strains quoted for stagnation regions are from the literature, and the whole essay is algebra on a closure rather than a solution of anything.

And realizability is necessary, not sufficient. A model that keeps every variance positive can still be wrong in every other way. What has been established is that the standard model fails a test no model can be allowed to fail, not that a model passing it is right.

What the collection has said about closures #

Four essays, and together they are the whole account.

What averaging costs states the problem: averaging leaves a term with no equation for it.

The ladder that never closes shows that chasing it does not terminate — each new equation brings a new unknown, for ever, so a model is not optional.

A guess with a constant in it supplies the simplest model there is, Prandtl’s mixing length, and is careful about what its agreement with a wall profile does and does not establish.

And this one asks what the standard two-equation model’s constant is a statement about, and finds that it is a statement about one flow with a domain of validity eleven per cent wide.

The sequence is not an argument that modelling is hopeless. It is an argument that a model’s constants should be read as measurements of particular flows, and that knowing which flow is most of knowing when the model applies. Ideal against real makes the same point about the exact theory, from the other end of the subject.

The residue #

The limit is the expansion in $Sk/\varepsilon$, truncated at the first term.

What survives it, at the strains real calculations meet, is everything that was dropped — and the dropped terms are not a small correction to a value. They are what makes a variance positive, and without them the model reports a mean square below zero and a production six hundred times the dissipation while every equation it is solving remains perfectly well behaved.

The limit was taken at a ratio of ten thirds. It is used at fifty.