Flow Velocity Gradient-Bounded Constitutive Theory of Incompressible Viscous Fluids

1. Introduction

1.1. The Navier–Stokes Equations and Their Limits

The equations that bear the names of Navier [15] and Stokes [24] have described the motion of incompressible viscous fluids for two centuries. Written in the tensor notation adopted throughout this paper, they read

$$\nabla ·\mathbf{u}=0,\rho \left[{\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·(\nabla \mathbf{u})\right]=-\nabla p+\mu {\nabla }^2\mathbf{u}+\rho \mathbf{b},$$

(1)

where $\mathbf{u}$ is the velocity, p the pressure, $\rho$ the density, $\mu$ the viscosity, $\mathbf{b}$ the body force per unit mass, $\nabla ={\mathbf{e}}_i{\partial }_i$ the gradient operator and $\nabla \mathbf{u}=\nabla \otimes \mathbf{u}=u_{j,i}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ the velocity gradient, so that $\mathbf{u}·(\nabla \mathbf{u})=u_j{u}_{i,j}{\mathbf{e}}_i$ is the convective acceleration. Their success is not in doubt. Direct numerical simulations of (1) reproduce turbulent statistics to within experimental accuracy at every Reynolds number that has so far been computed [14,18], and the equations underpin the whole of engineering fluid mechanics.

Whether (1) is complete, however, has been an open question since Leray [10] constructed weak solutions whose regularity he could not establish. The question was formalised by the Clay Mathematics Institute [7] and remained open until 8 September 2026, when a proof was announced, together with a formalisation in the Lean proof assistant, that a smooth incompressible fluid initially at rest, acted upon by a smooth bounded body force and possessing finite kinetic energy throughout, develops a singularity in finite time: the velocity becomes unbounded while the energy remains finite [17]. Closely related results for the Euler, Boussinesq and porous-medium equations were obtained independently at the same time [3]. The proof establishes statements C and D of the Clay formulation, i.e., blow-up under smooth forcing; the unforced problem remains open. The construction is a collapsing vortex whose core shrinks and accelerates without bound, its strain rate diverging as the singular time is approached.

For the physicist the mathematical statement has an immediate corollary. Infinite velocity and infinite velocity gradient are not attainable by any real liquid or gas: the incompressible equations are valid only for speeds well below the sound speed, and the Newtonian constitutive law from which the viscous term of (1) descends was derived, by Navier and by Stokes alike, under the assumption of a small velocity gradient. A model that can be driven by a bounded, smooth force into a state its own derivation excludes has, at that point, ceased to describe the fluid. The conclusion is not that the Navier–Stokes equations are wrong—they are demonstrably an excellent approximation over an enormous range of conditions—but that they are incomplete, in precisely the sense that they lack the physics that intervenes when the velocity gradient becomes large. The purpose of this paper is to identify that physics from first principles, to formulate the improved governing equations that result, to prove that they have the regularity the Navier–Stokes equations lack, and to work out their consequences in a form that can be tested.

1.2. Where the Equations Can Be Changed

The derivation of (1) contains one, and only one, hypothesis that is not a conservation law. Mass conservation for an incompressible fluid is $\nabla ·\mathbf{u}=0$; momentum conservation for any continuum is Cauchy’s equation

$$\rho \left[{\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·(\nabla \mathbf{u})\right]=\nabla ·\mathbf{\sigma }+\rho \mathbf{b},$$

(2)

with $\mathbf{\sigma }$ the Cauchy stress. Neither admits adjustment. The physical content of a particular fluid enters solely through the constitutive relation

$$\mathbf{\sigma }=-p\mathbf{I}+\mathbf{F}(\nabla \mathbf{u}),$$

(3)

and the choice $\mathbf{F}=\mu (\nabla \mathbf{u}+\mathbf{u}\nabla )$, linear in the velocity gradient, is what turns (2) into (1). This observation was made forcefully by Sun [26], who proposed retaining the next term permitted by the Cayley–Hamilton theorem and obtained $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{D}+4\lambda \mathbf{D}·\mathbf{D}$ with $\mathbf{D}=\frac{1}{2}(\nabla \mathbf{u}+\mathbf{u}\nabla )$ the rate-of-strain tensor and $\lambda$ a new “second-order viscosity”. The model is a particular member of the class of isotropic tensor functions studied by Reiner [20] and Rivlin [21], in which $\lambda$ is the cross-viscosity, and it raises exactly the right question: given that only $\mathbf{F}$ can be changed, what form must $\mathbf{F}$ take?

1.3. Programme and Results

We answer this question in three steps, summarised in Figure 1.

First (§Section 3) we derive the most general $\mathbf{F}$ compatible with material frame indifference, isotropy and the second law. Frame indifference removes the dependence on the spin, so that $\mathbf{F}=\mathbf{F}\left(\mathbf{D}\right)$; the representation theorem of Rivlin Ericksen [22] and the Cayley–Hamilton theorem reduce $\mathbf{F}$ to ${\alpha }_1\mathbf{D}+{\alpha }_2\mathbf{D}·\mathbf{D}$ with coefficients that are functions of the two nontrivial invariants of $\mathbf{D}$; and the second law imposes a dissipation inequality on those coefficients.

Second (§Section 4) we show that the constant-coefficient second-order model cannot be the improvement sought. Three theorems establish that with constant $\lambda \ne 0$ the model violates the second law in three dimensions, that in two dimensions the term $\nabla ·(\mathbf{D}·\mathbf{D})$ is an exact gradient so that the model is dynamically identical to Navier–Stokes, and that a quadratic stress produces a cubic, non-coercive dissipation that cannot bound the velocity gradient. The second of these has a consequence that must be stated plainly: every two-dimensional result obtained from the constant-coefficient model, including the wedge-flow boundary layers of Sun [26], is an artefact of an inconsistent reduction in which the modified pressure was not carried through the transverse momentum balance. The improvement must therefore enter the coefficient ${\alpha }_1$ of the linear term.

Third (§§Section 5Section 6) we formulate the lowest-order admissible law,

$$\mathbf{\sigma }=-p\mathbf{I}+2\mu \left(1+{\tau }^2\mathbf{D}:\mathbf{D}\right)\mathbf{D}+4\lambda \mathbf{D}·\mathbf{D},$$

(4)

in which $\tau$ is a material time—the reciprocal of the strain rate beyond which the fluid stiffens against further deformation—and derive its governing equations. We prove that the second law requires $\tau \ge \left|\lambda \right|/\left(\sqrt{6}\mu \right)$ (theorem 3) and, by an exact decomposition of the tangent stress operator into decoupled blocks, that monotonicity of the stress—the condition for Hadamard well-posedness—requires $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$ (theorem 4), exactly twice the thermodynamic bound. Under the latter condition the momentum operator is strongly monotone with quartic growth and the equations belong to the class introduced by Ladyzhenskaya [9] and Lions [11], for which global existence and uniqueness are known for growth exponent $p\ge 5/2$; here $p=4$. The forced finite-time blow-up of the Navier–Stokes equations therefore cannot occur (theorem 5). The mechanism is transparent: the effective viscosity $\mu (1+{\tau }^2|\mathbf{D}{|}^2)$ grows with the strain rate, so that the viscous stress in a collapsing core eventually grows as ${\left|\mathbf{D}\right|}^3$ and overwhelms the inertial amplification, which grows as ${\left|\mathbf{D}\right|}^2$. The collapse is arrested at strain rate $1/\tau$ and core size $\sqrt{\nu \tau }$; the Navier–Stokes equations, which contain no such scale, can be driven past it.

The remainder of the paper develops the theory to the point of experimental test. Section 7 gives exact solutions in simple shear and in gravity-driven pipe flow, from which a two-measurement protocol determines $\tau$ and $\lambda$ separately and the inequality $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$ provides an internal consistency check; it also presents the stiffened Burgers vortex as a numerical example of arrested collapse. Section 8 develops the boundary-layer theory of the new equations consistently to second order in the layer thickness. Its principal results are that $\lambda$ drops out of the boundary-layer equations identically, that the flat-plate layer has a leading-edge regime with $\delta \propto x^{1/4}$ and ${\tau }_w\propto x^{-3/4}$ that crosses over to Blasius scaling at the material length ${\ell }_{\tau }={\tau }^2{U}^3/\left(2\nu \right)$, and that among the Falkner–Skan wedge flows $U=C{x}^m$ only $m=1/3$ remains exactly self-similar. The non-similar flat-plate equations are solved by marching, the self-similar wedge is tabulated, and both are shown to agree with their analytic limits. Section 9 discusses the relation of the theory to turbulence and to direct numerical simulation, its regime of validity, its extensions, and what has and has not been proved.

2. Kinematics, Balance Laws and Notation

2.1. Tensor Notation

We follow the conventions of Sun [26]. Vectors and second-order tensors are set in boldface; ${\mathbf{e}}_i$ ($i=1,2,3$) are the Cartesian base vectors and summation over repeated indices is implied. The gradient operator is $\nabla ={\mathbf{e}}_i{\partial }_i$ and acts as a vector: the velocity gradient is the outer (dyadic) product

$$\nabla \mathbf{u}=\nabla \otimes \mathbf{u}=u_{j,i}{\mathbf{e}}_i\otimes {\mathbf{e}}_j,\mathbf{u}\nabla =\mathbf{u}\otimes \nabla =u_{i,j}{\mathbf{e}}_i\otimes {\mathbf{e}}_j={(\nabla \mathbf{u})}^{\mathsf{T}}.$$

(5)

The symmetric and antisymmetric parts of the velocity gradient are the rate-of-strain tensor and the spin tensor,

$$\mathbf{D}=\frac{1}{2}\left(\nabla \mathbf{u}+\mathbf{u}\nabla \right)=D_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j,\mathbf{W}=\frac{1}{2}\left(\mathbf{u}\nabla -\nabla \mathbf{u}\right),D_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right).$$

(6)

For second-order tensors $\mathbf{A}=A_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$, $\mathbf{B}=B_{ij}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ the single and double contractions are $\mathbf{A}·\mathbf{B}=A_{ik}{B}_{kj}{\mathbf{e}}_i\otimes {\mathbf{e}}_j$ and $\mathbf{A}:\mathbf{B}=A_{ij}{B}_{ij}$; the divergence of a tensor is the contraction with ∇ on the right, $\nabla ·\mathbf{\sigma }={\sigma }_{ij,j}{\mathbf{e}}_i$, which for a symmetric tensor coincides with the left divergence $\mathbf{\sigma }·\nabla$. The convective acceleration is $\mathbf{u}·(\nabla \mathbf{u})=u_j{u}_{i,j}{\mathbf{e}}_i$; the frequently written $(\mathbf{u}·\nabla )\mathbf{u}$ is avoided because it obscures the fact that the operand is the tensor $\nabla \mathbf{u}$. The Frobenius norm is ${\left|\mathbf{D}\right|}^2=\mathbf{D}:\mathbf{D}$, the Laplacian is ${\nabla }^2=\nabla ·\nabla$, and powers are ${\mathbf{D}}^2=\mathbf{D}·\mathbf{D}$, ${\mathbf{D}}^3=\mathbf{D}·\mathbf{D}·\mathbf{D}$.

The principal invariants of $\mathbf{D}$ are

$$I_1=tr\mathbf{D},I_2=\frac{1}{2}\left[{(tr\mathbf{D})}^2-tr{\mathbf{D}}^2\right],I_3=det\mathbf{D}.$$

(7)

For an incompressible fluid $I_1=\nabla ·\mathbf{u}=0$, so that $I_2=-\frac{1}{2}tr{\mathbf{D}}^2=-\frac{1}{2}{\left|\mathbf{D}\right|}^2$ and, by the Cayley–Hamilton identity ${\mathbf{D}}^3-I_1{\mathbf{D}}^2+I_2\mathbf{D}-I_3\mathbf{I}=\mathbf{0}$ traced, $I_3=\frac{1}{3}tr{\mathbf{D}}^3$. Throughout we use $\left|\mathbf{D}\right|$ and $tr{\mathbf{D}}^3$ as the two independent invariants.

2.2. Balance Laws

For a continuum of constant density $\rho$ the balances of mass and linear momentum are

$$\begin{array}{cc}\hfill \nabla ·\mathbf{u}& =0,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \rho \left[{\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·(\nabla \mathbf{u})\right]& =\nabla ·\mathbf{\sigma }+\rho \mathbf{b},\hfill \end{array}$$

(9)

and the balance of angular momentum, in the absence of body couples, is $\mathbf{\sigma }={\mathbf{\sigma }}^{\mathsf{T}}$. These hold for every continuum and admit no adjustment. Under isothermal conditions the second law, in the form of the Clausius–Duhem inequality with Helmholtz free energy $\psi$ [27], reduces to

$$\mathbf{\sigma }:\mathbf{D}-\rho \dot{\psi }\ge 0,$$

(10)

and for a fluid without elastic energy storage ($\dot{\psi }=0$ in isothermal motion) to the dissipation inequality $\mathbf{\sigma }:\mathbf{D}\ge 0$.

2.3. Constitutive Relation

To close (8)–() a constitutive relation must be supplied. For a fluid without memory the stress is a function of the present velocity gradient,

$$\mathbf{\sigma }=-p\mathbf{I}+\mathbf{F}(\nabla \mathbf{u}),$$

(11)

in which p is the constraint pressure enforcing (8) and $\mathbf{F}$ is the extra stress. The Newtonian choice $\mathbf{F}=2\mu \mathbf{D}$ gives (1), because $\nabla ·\mathbf{D}=\frac{1}{2}{\nabla }^2\mathbf{u}+\frac{1}{2}\nabla (\nabla ·\mathbf{u})=\frac{1}{2}{\nabla }^2\mathbf{u}$ for solenoidal $\mathbf{u}$. Every other property of a particular fluid enters through $\mathbf{F}$, and it is to the determination of $\mathbf{F}$ that we now turn.

3. The Admissible Form of the Constitutive Relation

Three principles constrain $\mathbf{F}$: material frame indifference, material symmetry (isotropy), and the second law. We take them in turn. None of the results of this section is new in continuum mechanics [16,22,27], but assembling them in the present context is what determines the structure of the improved equations, and the lemma on $tr{\mathbf{D}}^3$ that closes the section is used repeatedly.

3.1. Material Frame Indifference

Under a superposed rigid motion ${\mathbf{x}}^{*}=\mathbf{Q}\left(t\right)·\mathbf{x}+\mathbf{c}\left(t\right)$, with $\mathbf{Q}$ orthogonal, the velocity gradient transforms as

$${(\nabla \mathbf{u})}^{*}=\mathbf{Q}·(\nabla \mathbf{u})·{\mathbf{Q}}^{\mathsf{T}}+{\Omega }^{\mathsf{T}},\Omega =\dot{\mathbf{Q}}·{\mathbf{Q}}^{\mathsf{T}},$$

(12)

where $\Omega$ is antisymmetric. The symmetric part transforms objectively, ${\mathbf{D}}^{*}=\mathbf{Q}·\mathbf{D}·{\mathbf{Q}}^{\mathsf{T}}$, while the spin acquires the frame rotation, ${\mathbf{W}}^{*}=\mathbf{Q}·\mathbf{W}·{\mathbf{Q}}^{\mathsf{T}}+\Omega$. Since $\Omega$ is arbitrary, the requirement ${\mathbf{\sigma }}^{*}=\mathbf{Q}·\mathbf{\sigma }·{\mathbf{Q}}^{\mathsf{T}}$ can hold for all observers only if $\mathbf{F}$ does not depend on $\mathbf{W}$ [16]:

$$\mathbf{F}(\nabla \mathbf{u})=\mathbf{F}\left(\mathbf{D}\right),\mathbf{F}(\mathbf{Q}·\mathbf{D}·{\mathbf{Q}}^{\mathsf{T}})=\mathbf{Q}·\mathbf{F}\left(\mathbf{D}\right)·{\mathbf{Q}}^{\mathsf{T}}.$$

(13)

This is the first structural fact: whatever nonlinearity is added to the Newtonian law must be built from $\mathbf{D}$ alone. Models in which the stress depends on $\nabla \mathbf{u}$ through combinations other than $\mathbf{D}$ are not constitutive relations for a fluid but artefacts of a particular observer.

3.2. Isotropy and the Cayley–Hamilton Theorem

For an isotropic fluid the second condition in (13) holds for all orthogonal $\mathbf{Q}$, so that $\mathbf{F}$ is an isotropic tensor function of the single symmetric tensor $\mathbf{D}$. The representation theorem of Rivlin Ericksen [22] states that every such function is a polynomial in $\mathbf{D}$ with coefficients that are functions of the invariants, and the Cayley–Hamilton theorem [5],

$${\mathbf{D}}^3=I_1{\mathbf{D}}^2-I_2\mathbf{D}+I_3\mathbf{I},$$

(14)

reduces every power ${\mathbf{D}}^n$, $n\ge 3$, to a combination of $\mathbf{I}$, $\mathbf{D}$ and ${\mathbf{D}}^2$. Hence

$$\mathbf{F}\left(\mathbf{D}\right)={\alpha }_0(I_2,I_3)\mathbf{I}+{\alpha }_1(I_2,I_3)\mathbf{D}+{\alpha }_2(I_2,I_3)\mathbf{D}·\mathbf{D},$$

(15)

with $I_1=0$ for an incompressible fluid. The term ${\alpha }_0\mathbf{I}$ is an isotropic stress and is absorbed into the constraint pressure p. Two remarks are essential for what follows.

First, (15) is exhaustive: there is no further term to add. The Navier–Stokes equations correspond to ${\alpha }_1=2\mu$, ${\alpha }_2=0$; the model of Sun [26] to ${\alpha }_1=2\mu$, ${\alpha }_2=4\lambda$ with both constant; the fluid of Reiner [20] and Rivlin [21], of which the latter is a special case, to arbitrary ${\alpha }_1,{\alpha }_2$.

Second, the coefficients are functions of the invariants, not constants. The historical reason for treating them as constants is that in the small-gradient limit ${\alpha }_1\to 2\mu$ and the invariant dependence is a higher-order correction. Precisely because the problem before us is the large-gradient limit, the invariant dependence of ${\alpha }_1$ is not a refinement but the heart of the matter.

3.3. The Second Law

With (15) and $I_1=0$ the stress power is

$$\mathbf{\sigma }:\mathbf{D}=-ptr\mathbf{D}+{\alpha }_1\mathbf{D}:\mathbf{D}+{\alpha }_2(\mathbf{D}·\mathbf{D}):\mathbf{D}={\alpha }_1{\left|\mathbf{D}\right|}^2+{\alpha }_{2}tr{\mathbf{D}}^3,$$

(16)

since $(\mathbf{D}·\mathbf{D}):\mathbf{D}=D_{ik}{D}_{kj}{D}_{ij}=tr{\mathbf{D}}^3$ by symmetry. The dissipation inequality therefore reads

$$\Phi \left(\mathbf{D}\right)\equiv {\alpha }_1(I_2,I_3){\left|\mathbf{D}\right|}^2+{\alpha }_2(I_2,I_3)tr{\mathbf{D}}^3\ge 0\mathrm{for}\mathrm{all}\mathrm{traceless}\mathrm{symmetric}\mathbf{D}.$$

(17)

Whatever ${\alpha }_1,{\alpha }_2$ are chosen, they must satisfy (17) identically. The following lemma converts (17) into an explicit condition.

Lemma 1 

(bound on the cubic invariant). For every traceless symmetric $\mathbf{D}$ in three dimensions,

$$\left|tr{\mathbf{D}}^3\right|\le {\displaystyle \frac{{\left|\mathbf{D}\right|}^3}{\sqrt{6}}},$$

(18)

with equality if and only if $\mathbf{D}$ has eigenvalues proportional to $(2,-1,-1)$, i.e., is an axisymmetric extension or compression.

Proof. 

In the eigenbasis of $\mathbf{D}$, with eigenvalues $d_1,d_2,d_3$ and $d_1+d_2+d_3=0$, one has $tr{\mathbf{D}}^3=d_1^3+d_2^3+d_3^3=3d_1{d}_2{d}_3$, by the identity $a^3+b^3+c^3-3abc=(a+b+c)(a^2+b^2+c^2-ab-bc-ca)$. It suffices to maximise $d_1{d}_2{d}_3$ on the circle $d_1^2+d_2^2+d_3^2=1$, $d_1+d_2+d_3=0$. The Lagrange conditions $d_j{d}_k=2\alpha d_i+\beta$ for cyclic $(i,j,k)$ imply, on subtracting pairs, $(d_i-d_j)(d_k+2\alpha )=0$, so that at an extremum at least two eigenvalues coincide. With $d_2=d_3=-d_1/2$ and $d_1^2+d_1^2/2=1$ one finds $d_1=\pm \sqrt{2/3}$ and $3d_1{d}_2{d}_3=\pm \frac{3}{4}{d}_1^3=\pm \frac{3}{4}{(2/3)}^{3/2}=\pm 1/\sqrt{6}$. Homogeneity of degree three in $\mathbf{D}$ gives the general statement. □

Corollary 1. 

For traceless symmetric $\mathbf{D}$ the largest eigenvalue satisfies $d_{max}\le \sqrt{2/3}\left|\mathbf{D}\right|$, with equality in the same configuration.

Proof. 

Maximise $d_1$ subject to $\sum d_i=0$, $\sum d_i^2={\left|\mathbf{D}\right|}^2$: by Cauchy–Schwarz $d_2^2+d_3^2\ge \frac{1}{2}{(d_2+d_3)}^2=\frac{1}{2}{d}_1^2$, so ${\left|\mathbf{D}\right|}^2\ge \frac{3}{2}{d}_1^2$. □

3.4. Consequence: the Linear Coefficient Must Depend on the Strain Rate

Inequality (17) with lemma 1 gives the necessary and sufficient condition

$${\alpha }_1\ge {\displaystyle \frac{|{\alpha }_2|}{\sqrt{6}}}\left|\mathbf{D}\right|\mathrm{at}\mathrm{every}\mathbf{D},$$

(19)

because the cubic invariant can be made to take either sign at the extremal orientation. Equation (19) is the key structural constraint of the theory. If ${\alpha }_2$ is a nonzero constant, then ${\alpha }_1$ must grow at least linearly with $\left|\mathbf{D}\right|$; a constant ${\alpha }_1=2\mu$ is excluded. In other words, one cannot add a second-order term to the Newtonian law without at the same time modifying the coefficient of the first-order term, and the modification must make the fluid stiffer at large strain rate. The next section makes the failure of the constant-coefficient model explicit; §Section 5 then chooses the simplest ${\alpha }_1$ that satisfies (19).

4. The Constant-Coefficient Second-Order Model

The model $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{D}+4\lambda \mathbf{D}·\mathbf{D}$ with constant $\mu ,\lambda$ [26] is the natural first candidate for an improved law, and it is instructive to see precisely why it fails. The three theorems of this section are independent of one another; each alone would disqualify the model as the sought-for completion of the Navier–Stokes equations.

4.1. Violation of the Second Law

Theorem 1 

(thermodynamic inadmissibility). For constant $\mu >0$ and constant $\lambda \ne 0$, the dissipation $\Phi ={2\mu \left|\mathbf{D}\right|}^2+4\lambda tr{\mathbf{D}}^3$ is negative for some traceless symmetric $\mathbf{D}$. Specifically, with $\mathbf{D}=s\mathrm{sgn}(-\lambda )diag(2,-1,-1)/\sqrt{6}$,

$$\Phi =2\mu s^2-{\displaystyle \frac{4\left|\lambda \right|}{\sqrt{6}}}{s}^3<0\mathrm{whenever}s>{\displaystyle \frac{\sqrt{6}\mu }{2\left|\lambda \right|}}.$$

(20)

Proof. 

Immediate from lemma 1: at the extremal configuration $tr{\mathbf{D}}^3=\mp s^3/\sqrt{6}$, the sign being chosen opposite to that of $\lambda$. □

The failure is not marginal. For any fluid with $\lambda \ne 0$ there exists an axisymmetric straining motion of finite strain rate $\sqrt{6}\mu /\left(2\right|\lambda \left|\right)$ beyond which the model generates mechanical energy from nothing. This is the well-known thermodynamic defect of the Reiner–Rivlin fluid with constant cross-viscosity [6,27]. It is worth noting why it escaped detection in Sun [26]: the entropy inequality was verified there only for two-dimensional flow, Equation (25) of that paper, where $tr{\mathbf{D}}^3=u_{,x}^3+v_{,y}^3$ vanishes identically by continuity ($v_{,y}=-u_{,x}$). The check is vacuous in two dimensions and fails in three.

4.2. Dynamical Invisibility in Two Dimensions

Theorem 2 

(planar flows). For every planar incompressible flow $\mathbf{u}=u(x,y,t){\mathbf{e}}_x+v(x,y,t){\mathbf{e}}_y$ with $u_{,x}+v_{,y}=0$,

$$\nabla ·(\mathbf{D}·\mathbf{D})=\nabla \left(\frac{1}{2}\mathbf{D}:\mathbf{D}\right).$$

(21)

Consequently the equations of motion of the constant-coefficient model coincide, in planar flow, with the Navier–Stokes equations under the substitution $p\to \tilde{p}=p-2\lambda \mathbf{D}:\mathbf{D}$.

Proof. 

The in-plane part ${\mathbf{D}}_2$ of $\mathbf{D}$ is a traceless symmetric $2\times 2$ tensor, for which the Cayley–Hamilton theorem reads ${\mathbf{D}}_2·{\mathbf{D}}_2=-(det{\mathbf{D}}_2){\mathbf{I}}_2$; since $det{\mathbf{D}}_2=-\frac{1}{2}tr{\mathbf{D}}_2^2=-\frac{1}{2}\mathbf{D}:\mathbf{D}$ when $tr{\mathbf{D}}_2=0$, one has ${\mathbf{D}}_2·{\mathbf{D}}_2=\frac{1}{2}(\mathbf{D}:\mathbf{D}){\mathbf{I}}_2$. In components, with $u_{,x}=-v_{,y}$ and $\gamma \equiv u_{,y}+v_{,x}$, ${(\mathbf{D}·\mathbf{D})}_{xx}={(\mathbf{D}·\mathbf{D})}_{yy}=u_{,x}^2+\frac{1}{4}{\gamma }^2$ and ${(\mathbf{D}·\mathbf{D})}_{xy}=\frac{1}{2}(u_{,x}+v_{,y})\gamma =0$. The components involving ${\mathbf{e}}_z$ vanish. Taking the divergence, with ${\partial }_z=0$, gives (21) with $\frac{1}{2}\mathbf{D}:\mathbf{D}=u_{,x}^2+\frac{1}{4}{\gamma }^2$. Substituting $\mathbf{\sigma }=-p\mathbf{I}+2\mu \mathbf{D}+4\lambda \mathbf{D}·\mathbf{D}$ in () then gives $\nabla ·\mathbf{\sigma }=-\nabla (p-2\lambda \mathbf{D}:\mathbf{D})+\mu {\nabla }^2\mathbf{u}$. □

The identity was also confirmed by symbolic computation for an arbitrary stream function (appendix ). Its consequence for the two-dimensional results of Sun [26] is unavoidable. The wedge-flow boundary-layer equation of that paper, its Equation (27), retains a term $2\kappa u_{,y}{u}_{,xy}$ in the streamwise momentum balance while taking the pressure to be the outer pressure $p_e\left(x\right)$. But at the same order in the layer thickness the transverse momentum balance reads

$$0=-{\displaystyle \frac{1}{\rho }}{p}_{,y}+\kappa {\left[{\left(u_{,y}\right)}^2\right]}_{,y}⟹p=p_e\left(x\right)+\rho \kappa {\left(u_{,y}\right)}^2,$$

(22)

where $\kappa =\lambda /\rho$, and substituting this into the streamwise balance gives $-{\rho }^{-1}{p}_{,x}=U{U}_{,x}-2\kappa u_{,y}{u}_{,xy}$, which cancels the retained term identically. The consistent boundary-layer equations are Prandtl’s, unchanged; the wall shear parameter is $f_{,\eta \eta }\left(0\right)=1.233$ for every value of the parameter $\Lambda$ of that paper, and the “turbulent” profiles of its Figure 2 and Table 1 are not predictions of the model. We return to the boundary-layer reduction in §Section 8, where it is carried out consistently for the improved law.

The physical content of theorem 2 is that in planar flow the second-order stress $4\lambda \mathbf{D}·\mathbf{D}$ is isotropic in the plane: the first normal-stress difference of a Reiner–Rivlin fluid vanishes identically [2]. The term is dynamically active only through genuinely three-dimensional straining, in which the three eigenvalues of $\mathbf{D}$ are distinct.

4.3. Failure to Bound the Velocity Gradient

The third defect is not thermodynamic but analytic. Even where $\Phi \ge 0$, the dissipation of the constant-coefficient model is at most cubic in $\left|\mathbf{D}\right|$ and its sign is indefinite in the cubic term; it therefore provides no better a priori control of $\nabla \mathbf{u}$ than the quadratic dissipation ${2\mu \left|\mathbf{D}\right|}^2$ of the Navier–Stokes equations, which is exactly the control that fails to prevent blow-up. A term quadratic in $\mathbf{D}$ in the stress can never be coercive, because $\mathbf{D}\mapsto \mathbf{D}·\mathbf{D}$ is not monotone: the tangent map $\mathbf{E}\mapsto \mathbf{D}·\mathbf{E}+\mathbf{E}·\mathbf{D}$ has $\mathbf{E}:(\mathbf{D}·\mathbf{E}+\mathbf{E}·\mathbf{D})=2tr(\mathbf{D}·\mathbf{E}·\mathbf{E})$, which takes both signs.

Taken together, theorems 1 and 2 and the observation just made show that the Cayley–Hamilton term $\mathbf{D}·\mathbf{D}$ cannot by itself supply the missing physics. Condition (19) locates where the physics must go: in the coefficient of $\mathbf{D}$.

5. The Improved Constitutive Relation and the Governing Equations

5.1. Choice of the Coefficients

Condition (19) requires ${\alpha }_1$ to grow at least linearly in $\left|\mathbf{D}\right|$ if ${\alpha }_2\ne 0$, and the regularity theory of §Section 6 will require, independently of ${\alpha }_2$, that ${\alpha }_1$ grow at least as ${\left|\mathbf{D}\right|}^{1/2}$. The coefficient ${\alpha }_1$ is a scalar isotropic function of $\mathbf{D}$ and may depend on both invariants. Dependence on $tr{\mathbf{D}}^3$, which is odd under $\mathbf{D}\to -\mathbf{D}$, would make the linear coefficient different for extension and compression; we exclude it on the grounds that a fluid with no preferred direction should resist a straining motion and its reverse equally, and retain only dependence on ${\left|\mathbf{D}\right|}^2$. The lowest-order analytic correction is then

$${\alpha }_1=2\mu \left(1+{\tau }^2\mathbf{D}:\mathbf{D}\right),{\alpha }_2=4\lambda ,$$

(23)

which gives the constitutive relation

$$\begin{array}{|c|}\hline \mathbf{\sigma }=-p\mathbf{I}+2\mu \left(1+{\tau }^2\mathbf{D}:\mathbf{D}\right)\mathbf{D}+4\lambda \mathbf{D}·\mathbf{D}\\ \hline\end{array}$$

(24)

Three material constants appear: the Newtonian viscosity $\mu$ (dimension M/LT), the cross-viscosity $\lambda$ (dimension M/L, the “second-order viscosity coefficient” of Sun [26]), and a material time $\tau$ (dimension T). The combination $1/\tau$ is the strain rate at which the effective viscosity ${\mu }_{\mathrm{eff}}=\mu (1+{\tau }^2{\left|\mathbf{D}\right|}^2)$ has doubled; below it the fluid is Newtonian to $O\left({\tau }^2\right|\mathbf{D}{|}^2)$, above it the fluid stiffens. Higher even powers, ${\alpha }_1=2\mu (1+{\tau }^2{\left|\mathbf{D}\right|}^2+{\tau }_4^4{\left|\mathbf{D}\right|}^4+\dots )$, are admissible and would refine (24) at still larger strain rate; nothing below depends on their absence except the numerical values of the constants in the regularity estimates, and the theory of §Section 6 applies to any ${\alpha }_1$ with growth exponent $p-2\ge 1/2$.

5.2. Thermodynamic Admissibility

Theorem 3 

(second law). The dissipation of (24),

$$\Phi \left(\mathbf{D}\right)={2\mu \left|\mathbf{D}\right|}^2+2\mu {\tau }^2{\left|\mathbf{D}\right|}^4+4\lambda tr{\mathbf{D}}^3,$$

(25)

is non-negative for every traceless symmetric $\mathbf{D}$ if and only if

$$\tau \ge {\tau }_{\mathrm{th}}\equiv {\displaystyle \frac{\left|\lambda \right|}{\sqrt{6}\mu }}.$$

(26)

Proof. 

Let $s=\left|\mathbf{D}\right|$. By lemma 1, $\Phi \ge s^2\left[2\mu +2\mu {\tau }^2{s}^2-4\left|\lambda \right|s/\sqrt{6}\right]$, and by the inequality of arithmetic and geometric means $2\mu +2\mu {\tau }^2{s}^2\ge 4\mu \tau s$, so that $\Phi \ge 4s^3(\mu \tau -|\lambda |/\sqrt{6})\ge 0$ under (26). Conversely, at the extremal $\mathbf{D}$ of lemma 1 with the sign of $tr{\mathbf{D}}^3$ opposite to that of $\lambda$, $\Phi =s^2\left[2\mu -4\left|\lambda \right|s/\sqrt{6}+2\mu {\tau }^2{s}^2\right]$, a quadratic in s whose discriminant is $16{\lambda }^2/6-16{\mu }^2{\tau }^2$; it is negative for some s precisely when $\tau <\left|\lambda \right|/\left(\sqrt{6}\mu \right)$. At $\tau ={\tau }_{\mathrm{th}}$ the bracket is the perfect square ${(\lambda s-\sqrt{6}\mu )}^2/\left(3\mu \right)$, vanishing at $s=\sqrt{6}\mu /\left|\lambda \right|$. □

The theorem couples the two new constants. A fluid with cross-viscosity $\lambda$ cannot exist without a gradient-stiffening time $\tau$ of at least $\left|\lambda \right|/\left(\sqrt{6}\mu \right)$; the model of Sun [26] is the limit $\tau \to 0$, excluded unless $\lambda =0$. Figure 2(a) shows the admissible region in the $\left(\right|\lambda |/\mu ,\tau )$ plane; Figure 2(b) shows $\Phi$ along the worst orientation for several $\tau$. A second, stronger boundary also appears in the figure; it is the monotonicity condition of theorem 4, to which we come in §Section 6.

5.3. Governing Equations

Substituting (24) into () and using $\nabla ·\mathbf{D}=\frac{1}{2}{\nabla }^2\mathbf{u}$ for the Newtonian part, with $\nu =\mu /\rho$ and $\kappa =\lambda /\rho$ (dimension L²), we obtain the improved governing equations

$$\begin{array}{|c|}\hline \begin{array}{cc}\hfill \nabla ·\mathbf{u}& =0,\hfill \\ \hfill {\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·(\nabla \mathbf{u})& =\mathbf{b}-{\displaystyle \frac{1}{\rho }}\nabla p+\nu {\nabla }^2\mathbf{u}+2\nu {\tau }^2\nabla ·\left[(\mathbf{D}:\mathbf{D})\mathbf{D}\right]+4\kappa \nabla ·(\mathbf{D}·\mathbf{D}).\hfill \end{array}\\ \hline\end{array}$$

(27)

In Cartesian components,

$$u_{i,t}+u_j{u}_{i,j}=b_i-{\displaystyle \frac{1}{\rho }}{p}_{,i}+\nu u_{i,kk}+2\nu {\tau }^2{\left[\left(D_{mn}{D}_{mn}\right)D_{ij}\right]}_{,j}+4\kappa {\left[D_{ik}{D}_{kj}\right]}_{,j},$$

(28)

where $D_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})$. The last term is identical to that of Sun [26], Equation (16) there, since $4D_{ik}{D}_{kj}=(u_{i,k}+u_{k,i})(u_{k,j}+u_{j,k})$; the penultimate term is new. Expanding the derivatives,

$$2\nu {\tau }^2{\left[\left(D_{mn}{D}_{mn}\right)D_{ij}\right]}_{,j}=\nu {\tau }^2\left[\left(D_{mn}{D}_{mn}\right)u_{i,kk}+2D_{mn}{D}_{mn,j}(u_{i,j}+u_{j,i})\right],$$

(29)

which exhibits the structure of the stiffening term: an isotropic viscosity enhancement $\nu {\tau }^2{\left|\mathbf{D}\right|}^2{\nabla }^2\mathbf{u}$ plus a force directed along the gradient of ${\left|\mathbf{D}\right|}^2$. Explicit Cartesian and cylindrical component forms are given in appendix .

5.4. Dimensionless Form

With a length L and a velocity U, (27) becomes

$${\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\mathbf{u}·(\nabla \mathbf{u})=-\nabla p+{\displaystyle \frac{1}{\mathit{Re}}}\left\{{\nabla }^2\mathbf{u}+2{\mathcal{G}}^2\nabla ·\left[(\mathbf{D}:\mathbf{D})\mathbf{D}\right]\right\}+{\displaystyle \frac{4}{\mathcal{K}}}\nabla ·(\mathbf{D}·\mathbf{D}),$$

(30)

$$\mathit{Re}={\displaystyle \frac{UL}{\nu }},\mathcal{G}={\displaystyle \frac{\tau U}{L}},\mathcal{K}={\displaystyle \frac{L^2}{\kappa }}.$$

(31)

$\mathcal{K}$ is the parameter introduced by Sun [26]; it depends on the length scale alone, as noted there. $\mathcal{G}$ is a gradient number, the ratio of the flow strain rate $U/L$ to the material strain rate $1/\tau$; it plays the role that the Weissenberg number plays for viscoelastic liquids, although the physics is different, since (24) has no memory. The admissibility conditions read $\mathcal{G}\ge \mathit{Re}/\left(\sqrt{6}\mathcal{K}\right)$ (second law) and $\mathcal{G}\ge \sqrt{2/3}\mathit{Re}/\mathcal{K}$ (monotonicity). Relative to the Newtonian viscous term the stiffening term is ${\mathcal{G}}^2{\left|\mathbf{D}\right|}^2$, with $\left|\mathbf{D}\right|$ the dimensionless local strain rate: it is negligible where the local strain rate is small compared with $1/\tau$ and dominant where it is large, irrespective of the Reynolds number. This is the decisive difference from the term $4{\mathcal{K}}^{-1}\nabla ·(\mathbf{D}·\mathbf{D})$, which is fixed by the length scale.

5.5. Energy Identity

Taking the inner product of the momentum equation with $\rho \mathbf{u}$, integrating over a domain $\Omega$ with no-slip walls or periodic boundaries, and using $\nabla ·\mathbf{u}=0$ and the symmetry of $\mathbf{\sigma }$,

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\int }_{\Omega }{\displaystyle \frac{\rho }{2}}{\left|\mathbf{u}\right|}^2{\mathrm{d}}^{3}x+{\int }_{\Omega }\left[{2\mu \left|\mathbf{D}\right|}^2+2\mu {\tau }^2{\left|\mathbf{D}\right|}^4+4\lambda tr{\mathbf{D}}^3\right]{\mathrm{d}}^{3}x={\int }_{\Omega }\rho \mathbf{b}·\mathbf{u}{\mathrm{d}}^{3}x.$$

(32)

By theorem 3, whenever $\tau >{\tau }_{\mathrm{th}}$ the integrand is bounded below by $c_0{\mu \left|\mathbf{D}\right|}^2+c_1\mu {\tau }^2{\left|\mathbf{D}\right|}^4$ with $c_0,c_1>0$ (take $c_0=c_1=1-{\tau }_{\mathrm{th}}/\tau$ and use the arithmetic–geometric bound with the remaining fractions). Hence for any bounded forcing on $[0,T]$,

$$\underset{t\le T}{sup}{\parallel \mathbf{u}\parallel }_{L^2}^2+{\int }_0^T{\int }_{\Omega }{\left|\mathbf{D}\right|}^4{\mathrm{d}}^{3}x\mathrm{d}t\le C(T,{\mathbf{u}}_0,\mathbf{b}).$$

(33)

This is the estimate that the Navier–Stokes equations lack, and it is the source of everything in §Section 6.

5.6. Limiting Cases

1.

$\tau \to 0,\lambda \to 0$: the Navier–Stokes equations. More precisely, (27) reduces to (1) with relative error $O({\tau }^2{\left|\mathbf{D}\right|}^2)+O(\kappa \left|\mathbf{D}\right|/\nu ·{\ell }^{-1}\dots )$ wherever $\tau \left|\mathbf{D}\right|\ll 1$. Every flow for which direct numerical simulation of (1) agrees with experiment lies in this regime; the improved equations do not disturb that agreement.

2.

$\lambda =0$: the fluid of Ladyzhenskaya [9], $\mathbf{\sigma }=-p\mathbf{I}+2\mu (1+{\tau }^2|\mathbf{D}{|}^2)\mathbf{D}$, for which global well-posedness is classical.

3.

$\tau =0$, $\lambda \ne 0$: the constant-coefficient model of §Section 4, excluded by theorem 3.

4.

Planar flow: by theorem 2 the $\lambda$ term is absorbed into the pressure and the dynamics is that of the Ladyzhenskaya fluid.

6. Monotonicity, Well-Posedness and the Removal of Blow-Up

Thermodynamic admissibility guarantees that the fluid dissipates energy, but not that the equations are well posed. The additional requirement is that the extra stress

$$\mathbf{S}\left(\mathbf{D}\right)\equiv 2\mu (1+{\tau }^2|\mathbf{D}{|}^2)\mathbf{D}+4\lambda \mathbf{D}·\mathbf{D}$$

(34)

be a monotone function of $\mathbf{D}$: for all traceless symmetric ${\mathbf{D}}_1,{\mathbf{D}}_2$,

$$\left[\mathbf{S}\left({\mathbf{D}}_1\right)-\mathbf{S}\left({\mathbf{D}}_2\right)\right]:({\mathbf{D}}_1-{\mathbf{D}}_2)\ge 0.$$

(35)

Monotonicity is what makes the viscous operator elliptic in the sense required for the initial-value problem to be well posed in the sense of Hadamard: where it fails, the linearised equations admit modes that grow without bound at arbitrarily small scales, and the constitutive relation is unstable [8]. For the differentiable $\mathbf{S}$ at hand, (35) holds if and only if the tangent (Gâteaux) operator is positive semidefinite,

$$\mathcal{Q}(\mathbf{D};\mathbf{E})\equiv \mathbf{E}:\partial \mathbf{S}\left[\mathbf{D}\right]·\mathbf{E}\ge 0\mathrm{for}\mathrm{all}\mathrm{traceless}\mathrm{symmetric}\mathbf{D},\mathbf{E},$$

(36)

since $[\mathbf{S}\left({\mathbf{D}}_1\right)-\mathbf{S}\left({\mathbf{D}}_2\right)]:({\mathbf{D}}_1-{\mathbf{D}}_2)={\int }_0^1\mathcal{Q}({\mathbf{D}}_2+\theta ({\mathbf{D}}_1-{\mathbf{D}}_2);{\mathbf{D}}_1-{\mathbf{D}}_2)\mathrm{d}\theta$.

6.1. The Tangent Operator

Differentiating,

$$\partial \mathbf{S}\left[\mathbf{D}\right]·\mathbf{E}=2\mu (1+{\tau }^2{\left|\mathbf{D}\right|}^2)\mathbf{E}+4\mu {\tau }^2(\mathbf{D}:\mathbf{E})\mathbf{D}+4\lambda (\mathbf{D}·\mathbf{E}+\mathbf{E}·\mathbf{D}),$$

(37)

so that, using $\mathbf{E}:(\mathbf{D}·\mathbf{E}+\mathbf{E}·\mathbf{D})=2tr(\mathbf{D}·\mathbf{E}·\mathbf{E})$,

$$\mathcal{Q}(\mathbf{D};\mathbf{E})=2\mu (1+{\tau }^2{\left|\mathbf{D}\right|}^2{\left)\right|\mathbf{E}|}^2+4\mu {\tau }^2{(\mathbf{D}:\mathbf{E})}^2+8\lambda tr(\mathbf{D}·{\mathbf{E}}^2).$$

(38)

The first two terms are non-negative; the third is of indefinite sign and is the only route by which monotonicity can fail.

6.2. The Sharp Monotonicity Theorem

Theorem 4 

(monotonicity). The extra stress $\mathbf{S}$ of (24) is monotone on traceless symmetric tensors if and only if

$$\tau \ge {\tau }_{\mathrm{m}}\equiv \sqrt{\frac{2}{3}}{\displaystyle \frac{\left|\lambda \right|}{\mu }}=2{\tau }_{\mathrm{th}}.$$

(39)

For $\tau >{\tau }_{\mathrm{m}}$ it is strongly monotone: for all $\mathbf{D},\mathbf{E}$,

$$\mathcal{Q}(\mathbf{D};\mathbf{E})\ge 2\mu \left(1-{\displaystyle \frac{{\tau }_{\mathrm{m}}^2}{{\tau }^2}}\right){\left|\mathbf{E}\right|}^2.$$

(40)

Proof. 

Fix $\mathbf{D}$ and work in its eigenbasis, $\mathbf{D}=diag(d_1,d_2,d_3)$, $\sum d_i=0$, $s=\left|\mathbf{D}\right|$. Write $\mathbf{E}={\mathbf{E}}_{\mathrm{d}}+{\mathbf{E}}_{\mathrm{o}}$, with ${\mathbf{E}}_{\mathrm{d}}=diag(e_1,e_2,e_3)$, $\sum e_i=0$, and ${\mathbf{E}}_{\mathrm{o}}$ off-diagonal with entries $E_{ij}=E_{ji}$, $i\ne j$. The three terms of (38) decouple exactly between the diagonal and off-diagonal parts: ${\left|\mathbf{E}\right|}^2=|{\mathbf{E}}_{\mathrm{d}}{|}^2+{\left|{\mathbf{E}}_{\mathrm{o}}\right|}^2$; $\mathbf{D}:\mathbf{E}=\sum d_i{e}_i$ involves only ${\mathbf{E}}_{\mathrm{d}}$; and $tr(\mathbf{D}·{\mathbf{E}}^2)={\sum }_i{d}_i{\left({\mathbf{E}}^2\right)}_{ii}={\sum }_i{d}_i{\sum }_j{E}_{ij}^2={\sum }_i{d}_i{e}_i^2+{\sum }_{i<j}(d_i+d_j)E_{ij}^2$ contains no cross term. Hence $\mathcal{Q}={\mathcal{Q}}_{\mathrm{o}}+{\mathcal{Q}}_{\mathrm{d}}$ with

$$\begin{array}{cc}\hfill {\mathcal{Q}}_{\mathrm{o}}& =\sum _{i<j}{E}_{ij}^2\left[4\mu (1+{\tau }^2{s}^2)+8\lambda (d_i+d_j)\right]=\sum _{i<j}{E}_{ij}^2\left[4\mu (1+{\tau }^2{s}^2)-8\lambda d_k\right],\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\mathcal{Q}}_{\mathrm{d}}& =2\mu (1+{\tau }^2{s}^2){\left|\mathbf{e}\right|}^2+4\mu {\tau }^2{(\mathbf{d}·\mathbf{e})}^2+8\lambda \sum _i{d}_i{e}_i^2,\hfill \end{array}$$

(42)

where in (41) k is the index complementary to $(i,j)$ and $d_i+d_j=-d_k$.

Off-diagonal block.${\mathcal{Q}}_{\mathrm{o}}\ge 0$ for all ${\mathbf{E}}_{\mathrm{o}}$ if and only if $\mu (1+{\tau }^2{s}^2)\ge 2\lambda d_k$ for each k, i.e., $\mu (1+{\tau }^2{s}^2)\ge 2\left|\lambda \right|d_{max}$ over both signs of $\lambda$ (replace $\mathbf{D}$ by $-\mathbf{D}$). By corollary 1, $d_{max}\le \sqrt{2/3}s$ with equality attainable, so the condition for all $\mathbf{D}$ of norm s is

$$A\left(s\right)\equiv 2\mu (1+{\tau }^2{s}^2)\ge k\left(s\right)\equiv {\displaystyle \frac{8\left|\lambda \right|s}{\sqrt{6}}}.$$

(43)

Diagonal block. Parametrise the traceless diagonal tensors by an angle: $\mathbf{d}=s\sqrt{2/3}(cos\varphi ,cos(\varphi -\frac{2\pi }{3}),cos(\varphi +\frac{2\pi }{3}))$, and let $\widehat{\mathbf{d}}=\mathbf{d}/s$ and ${\widehat{\mathbf{e}}}_{\perp }=\sqrt{2/3}(sin\varphi ,sin(\varphi -\frac{2\pi }{3}),sin(\varphi +\frac{2\pi }{3}))$ be the orthonormal basis of the two-dimensional space of traceless diagonal tensors. Using ${\sum }_{i}cos{\varphi }_i={\sum }_{i}sin{\varphi }_i=0$, ${\sum }_{i}cos3{\varphi }_i=3cos3\varphi$, ${\sum }_{i}sin3{\varphi }_i=3sin3\varphi$ and the triple-angle formulae,

$$\sum _i{d}_i{\widehat{d}}_i^2={\displaystyle \frac{scos3\varphi }{\sqrt{6}}},\sum _i{d}_i{\widehat{e}}_{\perp i}^2=-{\displaystyle \frac{scos3\varphi }{\sqrt{6}}},\sum _i{d}_i{\widehat{d}}_i{\widehat{e}}_{\perp i}={\displaystyle \frac{ssin3\varphi }{\sqrt{6}}}.$$

(44)

With $\mathbf{e}=a\widehat{\mathbf{d}}+b{\widehat{\mathbf{e}}}_{\perp }$, () becomes the quadratic form

$${\mathcal{Q}}_{\mathrm{d}}=\left(\begin{array}{cc}a& b\end{array}\right)\left(\begin{array}{cc}A+4\mu {\tau }^2{s}^2+k^{\prime }c& k^{\prime }\varsigma \\ k^{\prime }\varsigma & A-k^{\prime }c\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right),$$

(45)

where $c=cos3\varphi$, $\varsigma =sin3\varphi$ and $k^{\prime }=8\lambda s/\sqrt{6}$. Its trace $2A+4\mu {\tau }^2{s}^2$ is positive and its determinant is

$$\begin{array}{cc}\hfill (A+4\mu {\tau }^2{s}^2+k^{\prime }c)(A-k^{\prime }c)-k^{\prime 2}{\varsigma }^2& =A^2-k^{\prime 2}+4\mu {\tau }^2{s}^2(A-k^{\prime }c)\hfill \\ \hfill & =(A-|k^{\prime }\left|\right)(A+|k^{\prime }\left|\right)+4\mu {\tau }^2{s}^2(A-k^{\prime }c).\hfill \end{array}$$

(46)

Since $|k^{\prime }c|\le |k^{\prime }|=k$, the determinant is non-negative for all $\varphi$ whenever $A\ge k$, and it is negative at $c=\mathrm{sgn}\left(\lambda \right)$ when $A<k$. Thus the diagonal block is positive semidefinite for all $\mathbf{D}$ of norm s under exactly the same condition (43) as the off-diagonal block.

Minimisation over s. Condition (43) for all $s\ge 0$ reads $2\mu +2\mu {\tau }^2{s}^2-8\left|\lambda \right|s/\sqrt{6}\ge 0$; the minimum of the left side over s is $2\mu -4{\lambda }^2/\left(3\mu {\tau }^2\right)$, which is non-negative if and only if ${\tau }^2\ge 2{\lambda }^2/\left(3{\mu }^2\right)$, i.e., (39). Finally, both blocks satisfy $\mathcal{Q}\ge {(A-k)\left|\mathbf{E}\right|}^2$ (for the diagonal block, the smaller eigenvalue of the matrix is at least $A-k$), and $A-k\ge 2\mu -4{\lambda }^2/\left(3\mu {\tau }^2\right)=2\mu (1-{\tau }_{\mathrm{m}}^2/{\tau }^2)$, which is (40). □

The theorem is sharp: at $\tau ={\tau }_{\mathrm{m}}$ the form $\mathcal{Q}$ vanishes when $\mathbf{D}$ is the axisymmetric tensor $\pm s^{*}diag(2,-1,-1)/\sqrt{6}$ of norm $s^{*}=\sqrt{3/2}\mu /\left|\lambda \right|$ and $\mathbf{E}$ is any shear in a plane containing the axis. In physical terms, for an axisymmetric strain $\mathbf{D}=\gamma diag(-\frac{1}{2},-\frac{1}{2},1)$ and a perturbing shear $\mathbf{E}$ in a plane containing the axis, $(\mathbf{D}·\mathbf{D})$ contributes $-\gamma \mathbf{E}$ to the shear stress, an effective viscosity reduction$-2\lambda \gamma$; the condition $\mu (1+\frac{3}{2}{\tau }^2{\gamma }^2)-2\lambda \gamma \ge 0$ for all $\gamma$ is precisely $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$. This configuration—a stretched vortex—is exactly the geometry of the blow-up construction, and we return to it in §Section 7.4. An independent numerical minimisation of $\mathcal{Q}/{\left|\mathbf{E}\right|}^2$ over all $(\mathbf{D},\mathbf{E})$ (appendix ) located the threshold at $\tau /\left(\right|\lambda |/\mu )=0.82\pm 0.02$, in agreement with $\sqrt{2/3}=0.8165$, and reproduced the bound (40) to four significant figures at $\tau /{\tau }_{\mathrm{m}}=1.04$ and $1.22$.

Corollary 2 

(coercivity and growth). For $\tau >{\tau }_{\mathrm{m}}$ there are constants $c_1,c_2,C>0$ depending only on $\mu ,\lambda ,\tau$ such that for all traceless symmetric $\mathbf{D},{\mathbf{D}}_1,{\mathbf{D}}_2$,

$$\begin{array}{cc}\hfill \mathbf{S}\left(\mathbf{D}\right):\mathbf{D}& \ge c_1{\left|\mathbf{D}\right|}^2+c_2{\tau }^2{\left|\mathbf{D}\right|}^4,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill \left|\mathbf{S}\right(\mathbf{D}\left)\right|& \le C(1+{\tau }^2|\mathbf{D}{|}^2\left)\right|\mathbf{D}|,\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill \mathbf{S}\left({\mathbf{D}}_1\right)-\mathbf{S}\left({\mathbf{D}}_2\right)]:({\mathbf{D}}_1-{\mathbf{D}}_2)& \ge \left[c_1+c_2{\tau }^2\left(|{\mathbf{D}}_1{|}^2+{\left|{\mathbf{D}}_2\right|}^2\right)\right]{|{\mathbf{D}}_1-{\mathbf{D}}_2|}^2.\hfill \end{array}$$

(49)

Proof. (47) is theorem 3 with a fraction $1-{\tau }_{\mathrm{m}}/\left(2\tau \right)$ of the quartic term retained after the arithmetic–geometric step; () follows from $|\mathbf{D}·\mathbf{D}|\le {\left|\mathbf{D}\right|}^2$; () follows by integrating (40) along the segment from ${\mathbf{D}}_2$ to ${\mathbf{D}}_1$ after retaining, in the proof of theorem 4, the fraction of $2\mu {\tau }^2{s}^2{\left|\mathbf{E}\right|}^2$ not consumed by the arithmetic–geometric bound (for $\tau >{\tau }_{\mathrm{m}}$ a positive fraction is available). □

Inequalities (47)–() are the hypotheses of the theory of generalised Newtonian fluids with p-structure, with growth exponent $p=4$ [9,11,13].

6.3. Global Existence, Uniqueness and Regularity

Theorem 5 

(global well-posedness). Let $\tau >{\tau }_{\mathrm{m}}$, let $\Omega ={\mathbb{R}}^3/{\mathbb{Z}}^3$, and let ${\mathbf{u}}_0$ be a smooth divergence-free initial velocity and $\mathbf{b}$ a smooth bounded body force. Then for every $T>0$ the system (27) possesses a unique weak solution $\mathbf{u}$ on $[0,T]$ with

$$\mathbf{u}\in L^{\infty }(0,T;L^2)\cap L^4(0,T;W^{1,4}),{\partial }_t\mathbf{u}\in L^{4/3}(0,T;{\left(W_{\mathrm{div}}^{1,4}\right)}^{*}),$$

(50)

satisfying the energy identity (32). Moreover the solution is strong,

$$\nabla \mathbf{u}\in L^{\infty }(0,T;L^2),{\nabla }^2\mathbf{u}\in L^2(0,T;L^2),{\int }_0^T{\int }_{\Omega }(1+{\tau }^2{\left|\mathbf{D}\right|}^2{\left)\right|\nabla \mathbf{D}|}^2{\mathrm{d}}^{3}x\mathrm{d}t<\infty ,$$

(51)

so that the enstrophy $\frac{1}{2}\int {|\nabla \times \mathbf{u}|}^2{\mathrm{d}}^{3}x$ is bounded uniformly on $[0,T]$ for every T.

Proof. 

The proof follows the monotone-operator method of Ladyzhenskaya [9] and Lions [11], chapter 2, §5, and we indicate only where the hypotheses enter. Existence. Galerkin approximations ${\mathbf{u}}_n$ satisfy the energy identity, and by (47) they are bounded in $L^{\infty }\left(L^2\right)\cap L^4\left(W^{1,4}\right)$ (Korn’s inequality on the torus gives ${\parallel \nabla \mathbf{u}\parallel }_{L^4}\le C{\parallel \mathbf{D}\parallel }_{L^4}$). By (), $\mathbf{S}\left({\mathbf{D}}_n\right)$ is bounded in $L^{4/3}$, and the convective term ${\mathbf{u}}_n·(\nabla {\mathbf{u}}_n)$ is bounded in $L^{4/3}\left(L^{4/3}\right)$ because ${\mathbf{u}}_n\in L^4\left(L^{\infty }\right)$ by the Sobolev embedding $W^{1,4}({\mathbb{R}}^3/{\mathbb{Z}}^3)\subset L^{\infty }$. Hence ${\partial }_t{\mathbf{u}}_n$ is bounded in $L^{4/3}\left({\left(W_{\mathrm{div}}^{1,4}\right)}^{*}\right)$, the Aubin–Lions lemma gives strong convergence in $L^2\left(L^2\right)$, the convective term passes to the limit, and monotonicity (35) with Minty’s trick identifies the weak limit of $\mathbf{S}\left({\mathbf{D}}_n\right)$ as $\mathbf{S}\left(\mathbf{D}\right)$. Uniqueness. For two solutions let $\mathbf{w}={\mathbf{u}}_1-{\mathbf{u}}_2$. Subtracting the equations, testing with $\mathbf{w}$ and using (),

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\parallel \mathbf{w}\parallel }_2^2+c_1\parallel {\mathbf{D}}_{\mathbf{w}}{\parallel }_2^2+c_2{\tau }^2\parallel \left(\right|{\mathbf{D}}_1|+|{\mathbf{D}}_2\left|\right){\mathbf{D}}_{\mathbf{w}}{\parallel }_2^2\le {\int }_{\Omega }\left|\mathbf{w}\right||\nabla {\mathbf{u}}_2\left|\right|\mathbf{w}|{\mathrm{d}}^{3}x.$$

(52)

The right side is bounded by $\parallel \nabla {\mathbf{u}}_2{\parallel }_4{\parallel \mathbf{w}\parallel }_{8/3}^2\le \parallel \nabla {\mathbf{u}}_2{\parallel }_4{\parallel \mathbf{w}\parallel }_2^{5/4}{\parallel \nabla \mathbf{w}\parallel }_2^{3/4}\le {ϵ\parallel \nabla \mathbf{w}\parallel }_2^2+C_{ϵ}\parallel \nabla {\mathbf{u}}_2{\parallel }_4^{8/5}{\parallel \mathbf{w}\parallel }_2^2$ by interpolation between $L^2$ and $L^6$ and Young’s inequality; since ${\int }_0^T{\parallel \nabla {\mathbf{u}}_2\parallel }_4^{8/5}\mathrm{d}t\le T^{3/5}{\left({\int }_0^T{\parallel \nabla {\mathbf{u}}_2\parallel }_4^4\mathrm{d}t\right)}^{2/5}<\infty$, Gronwall’s lemma gives $\mathbf{w}=\mathbf{0}$. (For general growth exponent p the same argument closes for $p\ge 5/2$; here $p=4$.) Regularity. On the torus one may test the equation with $-{\nabla }^2\mathbf{u}$; the viscous term yields, by () applied pointwise to difference quotients, the positive quantity $\int (c_1+c_2{\tau }^2{\left|\mathbf{D}\right|}^2{\left)\right|\nabla \mathbf{D}|}^2$, and the convective term is bounded by ${\parallel \mathbf{u}\parallel }_{\infty }{\parallel \nabla \mathbf{u}\parallel }_2\parallel {\nabla }^2{\mathbf{u}\parallel }_2\le {C\parallel \nabla \mathbf{u}\parallel }_4{\parallel \nabla \mathbf{u}\parallel }_2{\parallel {\nabla }^2\mathbf{u}\parallel }_2$, which is absorbed for $p\ge 5/2$ exactly as in Málek et al. [12]; see also Málek et al. [13], chapter 5. The estimate closes with $\nabla \mathbf{u}\in L^{\infty }\left(L^2\right)$ and ${\nabla }^2\mathbf{u}\in L^2\left(L^2\right)$. □

Corollary 3 

(no finite-time blow-up). Under the hypotheses of theorem 5, the velocity gradient of the improved equations satisfies ${\int }_0^T{\int }_{\Omega }{|\nabla \mathbf{u}|}^4{\mathrm{d}}^{3}x\mathrm{d}t<\infty$ and ${sup}_{t\le T}{\int }_{\Omega }{|\nabla \mathbf{u}|}^2{\mathrm{d}}^{3}x<\infty$ for every $T>0$. The finite-time singularity of the Navier–Stokes equations under smooth bounded forcing [17], in which the velocity becomes unbounded, cannot occur: for the Navier–Stokes equations bounded enstrophy implies regularity by the Ladyzhenskaya–Prodi–Serrin criterion, so any Navier–Stokes blow-up is accompanied by unbounded enstrophy, whereas the enstrophy of the improved equations is bounded on every finite interval.

6.4. Physical Mechanism and the Regularisation Scale

The mathematics has a simple physical reading, sketched in Figure 3. In the blow-up solution a vortex core of size ℓ collapses under an axial strain, and its azimuthal velocity v grows so that $\Gamma \sim v\ell$ stays finite. The strain rate in the core is $\left|\mathbf{D}\right|\sim v/\ell \sim \Gamma /{\ell }^2$, the inertial (stretching) amplification of vorticity is proportional to ${\left|\mathbf{D}\right|}^2$, and the Newtonian viscous stress that opposes it is proportional to $\left|\mathbf{D}\right|$: viscosity loses at small scale, and nothing in (1) stops the collapse. In (27) the viscous stress is $\mu (1+{\tau }^2|\mathbf{D}{|}^2\left)\right|\mathbf{D}|$, which grows as ${\left|\mathbf{D}\right|}^3$ and overtakes the amplification once $\tau \left|\mathbf{D}\right|\sim 1$. The collapse is arrested at strain rate $\left|\mathbf{D}\right|\sim 1/\tau$ and the core reaches a minimum size

$${\ell }_{min}\sim \sqrt{\nu \tau },$$

(53)

the length at which the Newtonian viscous time ${\ell }^2/\nu$ equals the material time $\tau$. This is the regularisation scale of the theory. It is a property of the fluid, not of the flow, and it plays for the improved equations the role that the Kolmogorov scale $\eta ={({\nu }^3/\epsilon )}^{1/4}$ plays for a turbulent flow: whereas $\eta$ can be made arbitrarily small by increasing the dissipation rate $\epsilon$, ${\ell }_{min}$ cannot. In a turbulent flow with $\eta \gg {\ell }_{min}$ the stiffening term is negligible everywhere and (27) is indistinguishable from (1); it acts only where the cascade, or a forced collapse, would otherwise drive the strain rate past $1/\tau$.

7. Exact Solutions and the Determination of the Material Constants

The theory contains two constants beyond the Newtonian viscosity. Both must be measured, and the measurements must be consistent with theorem 4. This section gives the exact solutions on which the measurements can be based, and a numerical example—the stiffened Burgers vortex—that exhibits the arrest of collapse in a genuinely three-dimensional flow.

7.1. Simple Shear and Viscometry

For the steady simple shear $\mathbf{u}=\dot{\gamma }y{\mathbf{e}}_x$ the velocity gradient is $\nabla \mathbf{u}=\dot{\gamma }{\mathbf{e}}_y\otimes {\mathbf{e}}_x$, so

$$\mathbf{D}={\displaystyle \frac{\dot{\gamma }}{2}}\left({\mathbf{e}}_x\otimes {\mathbf{e}}_y+{\mathbf{e}}_y\otimes {\mathbf{e}}_x\right),\mathbf{D}:\mathbf{D}={\displaystyle \frac{{\dot{\gamma }}^2}{2}},\mathbf{D}·\mathbf{D}={\displaystyle \frac{{\dot{\gamma }}^2}{4}}\left({\mathbf{e}}_x\otimes {\mathbf{e}}_x+{\mathbf{e}}_y\otimes {\mathbf{e}}_y\right).$$

(54)

The stress (24) is therefore

$${\sigma }_{xy}=\mu \left(1+\frac{1}{2}{\tau }^2{\dot{\gamma }}^2\right)\dot{\gamma },{\sigma }_{xx}={\sigma }_{yy}=-p+\lambda {\dot{\gamma }}^2,{\sigma }_{zz}=-p,$$

(55)

and the viscometric functions are

$$\eta \left(\dot{\gamma }\right)=\mu \left(1+\frac{1}{2}{\tau }^2{\dot{\gamma }}^2\right),N_1={\sigma }_{xx}-{\sigma }_{yy}=0,N_2={\sigma }_{yy}-{\sigma }_{zz}=\lambda {\dot{\gamma }}^2.$$

(56)

Figure 4 displays them. Three remarks follow. (i) The shear viscosity increases quadratically with the shear rate; the material time is read directly from the rate ${\dot{\gamma }}_{*}=\sqrt{2}/\tau$ at which it has doubled, or from the curvature of ${\sigma }_{xy}\left(\dot{\gamma }\right)$ at small rate. (ii) The first normal-stress difference vanishes identically—this is theorem 2 again—and the second is quadratic in the rate with coefficient $\lambda$. The measurement of $N_2$ therefore gives $\lambda$ including its sign. (iii) The consistency condition $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$ becomes, in terms of measurable quantities,

$${\displaystyle \frac{\eta \left(\dot{\gamma }\right)-\mu }{\mu }}\ge {\displaystyle \frac{1}{3}}{\displaystyle \frac{N_2{\left(\dot{\gamma }\right)}^2}{{\mu }^2{\dot{\gamma }}^2}}\mathrm{for}\mathrm{all}\dot{\gamma },$$

(57)

which is a falsifiable statement about a single rheogram. A fluid that showed a measurable $N_2$ without shear thickening at the rate $\sqrt{3/2}\mu /\left|\lambda \right|$ would violate the theory.

7.2. Gravity-Driven Flow in an Inclined Pipe

Consider the steady flow of the fluid down an infinitely long circular pipe of radius R inclined at an angle $\alpha$ to the horizontal, with no imposed pressure gradient—the configuration proposed by Sun [26] for the determination of $\lambda$. We use cylindrical coordinates $(r,\theta ,z)$ with ${\mathbf{e}}_z$ along the pipe axis in the direction of the flow, so that the gravitational body force is $\mathbf{b}=gsin\alpha {\mathbf{e}}_z-gcos\alpha (cos\theta {\mathbf{e}}_r-sin\theta {\mathbf{e}}_{\theta })$, and seek by the semi-inverse method a solution $\mathbf{u}=w\left(r\right){\mathbf{e}}_z$. Then

$$\mathbf{D}={\displaystyle \frac{w^{\prime }}{2}}\left({\mathbf{e}}_r\otimes {\mathbf{e}}_z+{\mathbf{e}}_z\otimes {\mathbf{e}}_r\right),\mathbf{D}:\mathbf{D}={\displaystyle \frac{w^{\prime 2}}{2}},\mathbf{D}·\mathbf{D}={\displaystyle \frac{w^{\prime 2}}{4}}\left({\mathbf{e}}_r\otimes {\mathbf{e}}_r+{\mathbf{e}}_z\otimes {\mathbf{e}}_z\right),$$

(58)

and the nonzero stress components are

$${\sigma }_{rz}=\mu \left(1+\frac{1}{2}{\tau }^2{w}^{\prime 2}\right)w^{\prime },{\sigma }_{rr}={\sigma }_{zz}=-p+\lambda w^{\prime 2},{\sigma }_{\theta \theta }=-p.$$

(59)

Continuity is satisfied identically. The three components of () are

$$\begin{array}{cc}\hfill & r:-p_{,r}+{\displaystyle \frac{\lambda }{r}}{\left(r{w}^{\prime 2}\right)}^{\prime }-\rho gcos\alpha cos\theta =0,\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill & \theta :-{\displaystyle \frac{1}{r}}{p}_{,\theta }+\rho gcos\alpha sin\theta =0,\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill & z:-p_{,z}+{\displaystyle \frac{1}{r}}{\left[r\mu \left(1+\frac{1}{2}{\tau }^2{w}^{\prime 2}\right)w^{\prime }\right]}^{\prime }+\rho gsin\alpha =0.\hfill \end{array}$$

(62)

Since the pipe is infinite and the ends are at equal pressure, $p_{,z}=0$; the $\theta$ and r equations then give

$$p=p_0-\rho gcos\alpha rcos\theta +\lambda {\int }_0^r{\displaystyle \frac{1}{\varrho }}{\left(\varrho w^{\prime 2}\right)}^{\prime }\mathrm{d}\varrho ,$$

(63)

where $p_0$ is the pressure on the axis; the hydrostatic term vanishes in a vertical pipe ($\alpha =\pi /2$), and we write $p_R$ for the wall pressure in that case. The z equation integrates once, with regularity at the axis, to

$$\mu \left(1+\frac{1}{2}{\tau }^2{w}^{\prime 2}\right)w^{\prime }=-\frac{1}{2}\rho gsin\alpha r.$$

(64)

Equation (64) is a cubic for $w^{\prime }\left(r\right)$ whose left side is strictly increasing in $w^{\prime }$; it has exactly one real root for every r, and the profile is obtained by quadrature from $w\left(R\right)=0$. Note the division of labour between the two constants: $\tau$ controls the velocity profile through (64) and $\lambda$ does not enter it, because ${(\mathbf{D}·\mathbf{D})}_{rz}=0$; $\lambda$ controls the transverse pressure distribution through (60) and $\tau$ enters that only through $w^{\prime }$.

With the wall shear rate of Poiseuille flow $s_0=\rho gRsin\alpha /\left(2\mu \right)$, the dimensionless variables $\zeta =r/R$, $W=w^{\prime }/s_0$ and the parameter

$$\Theta =\frac{1}{2}{\tau }^2{s}_0^2,$$

(65)

(64) reads $(1+\Theta W^2)W=-\zeta$, and the flow rate and wall–axis pressure difference in the vertical pipe ($\alpha =\pi /2$) are

$${\displaystyle \frac{Q}{Q_{\mathrm{P}}}}=8{\int }_0^1\zeta {\int }_{\zeta }^1(-W)\mathrm{d}{\zeta }^{\prime }\mathrm{d}\zeta ,p_R-p_0=\lambda s_0^2\mathcal{P}\left(\Theta \right),\mathcal{P}\left(\Theta \right)={\int }_0^1{\displaystyle \frac{1}{\zeta }}{\left(\zeta W^2\right)}^{\prime }\mathrm{d}\zeta ,$$

(66)

with $Q_{\mathrm{P}}=\pi \rho g{R}^{4}sin\alpha /\left(8\mu \right)$ the Poiseuille flow rate. For $\Theta =0$, $W=-\zeta$, $Q=Q_{\mathrm{P}}$ and $\mathcal{P}=3/2$, which recovers Equation (47) of Sun [26], $\lambda =\frac{8}{3}{(\mu /\rho g)}^2(p_R-p_0)/R^2$. For $\Theta \to \infty$, $W\to -{(\zeta /\Theta )}^{1/3}$ and $Q/Q_{\mathrm{P}}\to \frac{6}{5}{\Theta }^{-1/3}$: the profile approaches that of a power-law fluid of index $n=3$ and the wall shear rate falls as ${\Theta }^{-1/3}$. Table 2 and Figure 5 give the intermediate values.

7.3. Protocol for the Determination of $\tau$ and $\lambda$

The two exact solutions lead to the following procedure (Figure 6).

1.

Measure $\mu$ in the Newtonian range, by any standard method.

2.

In a vertical pipe of radius R, measure the flow rate Q under gravity alone, over a range of R (or of $gsin\alpha$) so that $s_0=\rho gRsin\alpha /\left(2\mu \right)$ spans the range where $Q/Q_{\mathrm{P}}$ departs from unity. Invert the first of (66), or read Table 2, to obtain $\Theta$ and hence $\tau =\sqrt{2\Theta }/s_0$. Every pipe must return the same $\tau$; the collapse of $Q/Q_{\mathrm{P}}$ against ${\tau }^2{s}_0^2$ is the first test of the theory.

3.

In the same pipe, measure the pressure difference between the wall and the axis in a cross-section, $p_R-p_0$. With $\Theta$ known, the second of (66) gives $\lambda =(p_R-p_0)/\left[s_0^2\mathcal{P}\left(\Theta \right)\right]$, including its sign.

4.

Verify $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$ (theorem 4). Failure would indicate that the fluid is outside the class described by (24), for instance because it has memory.

5.

Independently, in a rotational rheometer, measure $\eta \left(\dot{\gamma }\right)$ and $N_2\left(\dot{\gamma }\right)$ and compare with (56). The agreement of the two determinations of $\tau$ and of $\lambda$ is the second test.

For low-viscosity liquids such as water no departure from Newtonian behaviour has been resolved by rheometry up to shear rates of order ${10}^5$–${10}^6{\mathrm{s}}^{-1}$, which places $\tau$ below roughly ${10}^{-6}$ s and ${\ell }_{min}=\sqrt{\nu \tau }$ below about a micrometre; the theory is therefore not in conflict with any existing measurement on simple liquids, and its regularisation scale, while small, remains three orders of magnitude above the molecular scale at which the continuum description itself would fail. For suspensions and other complex fluids that shear-thicken at accessible rates, all of the quantities in steps 2–5 are measurable with existing apparatus.

7.4. A Strained Vortex: Arrest of Collapse

The Burgers vortex [4] is the exact steady Navier–Stokes solution in which an axisymmetric vortex $\mathbf{u}=v\left(r\right){\mathbf{e}}_{\theta }$ is held against viscous diffusion by the irrotational strain $-\frac{1}{2}\gamma r{\mathbf{e}}_r+\gamma z{\mathbf{e}}_z$. It is the simplest model of the stretched vortex tubes that populate turbulent flows and, in the geometry of OpenAI [17], of the collapsing core of the blow-up solution, in which $\gamma$ grows without bound. We use it here as a three-dimensional test of the improved law with all terms active.

With $\mathbf{u}=-\frac{1}{2}\gamma r{\mathbf{e}}_r+v\left(r\right){\mathbf{e}}_{\theta }+\gamma z{\mathbf{e}}_z$ the rate of strain is

$$\mathbf{D}=\gamma \left(-\frac{1}{2}{\mathbf{e}}_r\otimes {\mathbf{e}}_r-\frac{1}{2}{\mathbf{e}}_{\theta }\otimes {\mathbf{e}}_{\theta }+{\mathbf{e}}_z\otimes {\mathbf{e}}_z\right)+{\displaystyle \frac{\varsigma }{2}}\left({\mathbf{e}}_r\otimes {\mathbf{e}}_{\theta }+{\mathbf{e}}_{\theta }\otimes {\mathbf{e}}_r\right),\varsigma \equiv v^{\prime }-{\displaystyle \frac{v}{r}},$$

(67)

so that ${\left|\mathbf{D}\right|}^2=\frac{3}{2}{\gamma }^2+\frac{1}{2}{\varsigma }^2$ and, because the strain part has $D_{rr}+D_{\theta \theta }=-\gamma$, ${(\mathbf{D}·\mathbf{D})}_{r\theta }=-\frac{1}{2}\gamma \varsigma$. The azimuthal shear stress is therefore

$${\sigma }_{r\theta }=\left[\mu \left(1+{\tau }^2{\left|\mathbf{D}\right|}^2\right)-2\lambda \gamma \right]\varsigma ,$$

(68)

which makes explicit the observation following theorem 4: in a stretched vortex the cross-viscosity reduces the effective viscosity by $2\lambda \gamma$, and it is the term $\frac{3}{2}\mu {\tau }^2{\gamma }^2$ that keeps the effective viscosity positive for all $\gamma$ precisely when $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$. The azimuthal momentum equation, with the radial and axial equations satisfied by a pressure that is a function of r and z alone, is

$$-{\displaystyle \frac{\gamma r}{2}}\left(v^{\prime }+{\displaystyle \frac{v}{r}}\right)={\displaystyle \frac{1}{\rho r^2}}{\left(r^2{\sigma }_{r\theta }\right)}^{\prime },$$

(69)

with $v\left(0\right)=0$ and $v\to \Gamma /\left(2\pi r\right)$ as $r\to \infty$. In terms of the Burgers radius $r_B=\sqrt{4\nu /\gamma }$, the scaled variables $R=r/r_B$, $V=2\pi r_{B}v/\Gamma$, and the parameters $\tau \gamma$, $\lambda \gamma /\mu$ and ${\mathit{Re}}_{\Gamma }=\Gamma /\left(2\pi \nu \right)$, this becomes

$$-2R\left(V^{\prime }+{\displaystyle \frac{V}{R}}\right)={\displaystyle \frac{1}{R^2}}{\left[R^{2}M\left(V^{\prime }-{\displaystyle \frac{V}{R}}\right)\right]}^{\prime },M=1+{\tau }^2{\gamma }^2\left[\frac{3}{2}+\frac{1}{32}{\mathit{Re}}_{\Gamma }^2{\left(V^{\prime }-{\displaystyle \frac{V}{R}}\right)}^2\right]-2{\displaystyle \frac{\lambda \gamma }{\mu }},$$

(70)

whose Newtonian solution is $V=(1-{\mathrm{e}}^{-R^2})/R$ with peak vorticity $2\pi r_B^2\omega \left(0\right)/\Gamma =2$. Equation (70) was solved by a finite-difference Newton iteration in the regular variable $V/R$ (appendix ); Figure 7 shows the result for ${\mathit{Re}}_{\Gamma }=100$.

The peak vorticity of the Burgers vortex is $\gamma \Gamma /\left(4\pi \nu \right)$: for a given circulation it grows linearly with the strain rate, without limit, which is the elementary form of the collapse mechanism. In the improved law it is reduced by 2% at $\tau \gamma =0.02$, 7% at $0.04$, 18% at $0.08$ and 28% at $0.12$, and the core radius (the position of the velocity maximum) grows from $1.12r_B$ to $1.32r_B$ over the same range. The reduction sets in when $\tau \left|\mathbf{D}\right|$ in the core reaches order one; since $\left|\mathbf{D}\right|$ in the core is dominated by the azimuthal shear, of order $\Gamma \gamma /\left(8\pi \nu \right)=\frac{1}{4}{\mathit{Re}}_{\Gamma }\gamma$, the onset is at $\tau \gamma \sim 4/{\mathit{Re}}_{\Gamma }$, i.e., $0.04$ for the case shown, which is what the figure displays. The cross-viscosity, with the sign that reduces the effective viscosity, slightly intensifies the vortex at small $\tau \gamma$ but is overpowered by the stiffening as $\tau \gamma$ grows. Thus a vortex subjected to ever-increasing strain does not collapse: its core saturates at the size $\sqrt{\nu \tau }$ of (53) and its peak vorticity at the value $\Gamma /\left(2\pi \nu \tau \right)$ set by the material time, rather than growing without bound with $\gamma$.

8. Boundary-Layer Theory

The boundary layer is where velocity gradients are largest in a high-Reynolds-number flow, and it is where Prandtl [19] showed that a viscous term, however small its coefficient, cannot be neglected. It is therefore the natural setting in which the improved law acts, and it is also the setting most amenable to experimental verification, since wall shear stress and boundary-layer profiles are routinely measured. This section derives the boundary-layer equations of (27) consistently, solves them for the flat plate and for wedge flows, and identifies the measurements that would confirm or refute the theory.

8.1. Order-of-Magnitude Analysis

Consider steady two-dimensional flow past a wall $y=0$ with outer velocity $U\left(x\right)$, and let $\delta$ be the layer thickness, L the streamwise length, $ϵ=\delta /L\ll 1$. Prandtl’s estimates, listed in table I of Sun [26], are $u,u_{,x},v_{,y}=O\left(1\right)$; $v,v_{,x}=O\left(ϵ\right)$; $u_{,y}=O(1/ϵ)$; $u_{,yy}=O(1/{ϵ}^2)$, in units of U and L. In planar flow

$$\mathbf{D}:\mathbf{D}=2u_{,x}^2+\frac{1}{2}{\left(u_{,y}+v_{,x}\right)}^2=\frac{1}{2}{u}_{,y}^2\left[1+O\left({ϵ}^2\right)\right],$$

(71)

and the components of the two new stress operators are, to leading order,

$$\begin{array}{cc}\hfill 2\nu {\tau }^2\nabla ·\left[(\mathbf{D}:\mathbf{D})\mathbf{D}\right]·{\mathbf{e}}_x& =\frac{1}{2}\nu {\tau }^2{\left(u_{,y}^3\right)}_{,y}+\dots =O\left({\displaystyle \frac{\nu {\tau }^2{U}^3}{{\delta }^4}}\right),\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill 2\nu {\tau }^2\nabla ·\left[(\mathbf{D}:\mathbf{D})\mathbf{D}\right]·{\mathbf{e}}_y& =\nu {\tau }^2{\left(u_{,y}^2{v}_{,y}\right)}_{,y}+\dots =O\left({\displaystyle \frac{\nu {\tau }^2{U}^3}{{\delta }^{2}L}}\right),\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill 4\kappa \nabla ·(\mathbf{D}·\mathbf{D})·{\mathbf{e}}_x& =\kappa {\left[u_{,y}^2\right]}_{,x}+\dots =O\left({\displaystyle \frac{\kappa U^2}{{\delta }^{2}L}}\right),\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill 4\kappa \nabla ·(\mathbf{D}·\mathbf{D})·{\mathbf{e}}_y& =\kappa {\left[u_{,y}^2\right]}_{,y}+\dots =O\left({\displaystyle \frac{\kappa U^2}{{\delta }^3}}\right).\hfill \end{array}$$

(75)

The y-momentum equation determines the pressure variation across the layer. The term () produces $p=p_e\left(x\right)+\rho \kappa u_{,y}^2+\dots$, whose streamwise gradient is of the same order as () and, as shown in §Section 4, cancels it identically—this is theorem 2 in boundary-layer form. The term () produces a pressure variation $\Delta p/\rho =O(\nu {\tau }^2{U}^3/\left(\delta L\right))$ whose streamwise gradient, $O(\nu {\tau }^2{U}^3/\left(\delta L^2\right))$, is smaller than (72) by a factor ${ϵ}^2$. We conclude:

Proposition 1 

(boundary-layer equations). To leading order in ϵ, the boundary-layer equations of the improved system (27) are

$$u_{,t}+u{u}_{,x}+v{u}_{,y}=U{U}_{,x}+\nu {\left[\left(1+\frac{1}{2}{\tau }^2{u}_{,y}^2\right)u_{,y}\right]}_{,y},u_{,x}+v_{,y}=0,$$

(76)

with $u=v=0$ at $y=0$ and $u\to U\left(x\right)$ as $y/\delta \to \infty$. The cross-viscosity λ does not appear; the pressure across the layer is $p=p_e\left(x\right)+\rho \kappa u_{,y}^2+O\left({ϵ}^2\right)$. The wall shear stress is

$${\tau }_w=\mu \left(1+\frac{1}{2}{\tau }^2{u}_{,y}^2\right)u_{,y}{|}_{y=0}.$$

(77)

Equation (76) is Prandtl’s equation with the shear viscosity replaced by the viscometric function $\eta \left(u_{,y}\right)$ of (56). This is as it should be: within a thin layer the flow is locally a simple shear, and the only material function a simple shear can probe is $\eta \left(\dot{\gamma }\right)$. The corresponding equation for a power-law fluid was studied by Acrivos et al. [1]; the present case differs in possessing a Newtonian core at small rate, which produces the two-regime structure described next. The retained term is the largest in the equation wherever $\tau u_{,y}=O\left(1\right)$, and the parameter that measures its importance is

$$\tau u_{,y}\sim \tau {\displaystyle \frac{U}{\delta }}={\displaystyle \frac{\tau U}{L}}\sqrt{\mathit{Re}}=\mathcal{G}\sqrt{\mathit{Re}}:$$

(78)

the correction is governed by $\mathcal{G}\sqrt{\mathit{Re}}$, not by $\mathit{Re}$ alone, and increases with Reynolds number at fixed gradient number because the layer thins.

8.2. The Flat Plate: Non-Similar Equations and the Material Length

For $U=\mathrm{const}$ introduce the Blasius variables $\psi =\sqrt{\nu Ux}f(x,\eta )$, $\eta =y\sqrt{U/\left(\nu x\right)}$, so that $u=U{f}_{,\eta }$ and $u_{,y}=U\sqrt{U/\left(\nu x\right)}{f}_{,\eta \eta }$. Then $\frac{1}{2}{\tau }^2{u}_{,y}^2=({\tau }^2{U}^3/2\nu x)f_{,\eta \eta }^2$, and the natural streamwise scale is the material length

$${\ell }_{\tau }\equiv {\displaystyle \frac{{\tau }^2{U}^3}{2\nu }},X\equiv {\displaystyle \frac{x}{{\ell }_{\tau }}}.$$

(79)

Substituting into (76), the standard non-similar form of the convective terms [23] gives

$${\left[\left(1+{\displaystyle \frac{f_{,\eta \eta }^2}{X}}\right)f_{,\eta \eta }\right]}_{,\eta }+\frac{1}{2}f{f}_{,\eta \eta }=X\left(f_{,\eta }{f}_{,\eta X}-f_{,\eta \eta }{f}_{,X}\right),$$

(80)

with $f=f_{,\eta }=0$ at $\eta =0$ and $f_{,\eta }\to 1$ as $\eta \to \infty$. For $X\to \infty$ the stiffening term vanishes and (80) reduces to the Blasius equation; the flat plate “forgets” the material time downstream, because the wall shear rate $u_{,y}\propto x^{-1/2}$ falls below $1/\tau$ at $x\sim {\ell }_{\tau }$. Upstream of ${\ell }_{\tau }$ the stiffening term dominates and a different similarity structure emerges.

8.2.1. Leading-Edge Similarity Solution

Seek $f=X^{-1/4}F\left(\zeta \right)$ with $\zeta =\eta X^{1/4}$. Then $f_{,\eta }=F^{\prime }$, $f_{,\eta \eta }=X^{1/4}{F}^{\prime \prime }$, $f_{,\eta \eta \eta }=X^{1/2}{F}^{\prime \prime \prime }$, and one finds $X{f}_{,\eta }{f}_{,\eta X}=\frac{1}{4}\zeta F^{\prime }{F}^{\prime \prime }$ and $X{f}_{,\eta \eta }{f}_{,X}=\frac{1}{4}{F}^{\prime \prime }(\zeta F^{\prime }-F)$, so that the right side of (80) is $\frac{1}{4}F{F}^{\prime \prime }$ exactly. The term ${[f_{,\eta \eta }^3/X]}_{,\eta }={\left(F^{\prime \prime 3}\right)}^{\prime }$ is $O\left(1\right)$ while the Newtonian term $f_{,\eta \eta \eta }=X^{1/2}{F}^{\prime \prime \prime }$ vanishes as $X\to 0$. The leading-edge problem is therefore

$${\left(F^{\prime \prime 3}\right)}^{\prime }+\frac{1}{4}F{F}^{\prime \prime }=0,F\left(0\right)=F^{\prime }\left(0\right)=0,F^{\prime }(\infty )=1,$$

(81)

which is the boundary-layer equation of a power-law fluid of index $n=3$ in the form of Acrivos et al. [1]. Its numerical solution (appendix ) gives

$$F^{\prime \prime }\left(0\right)=0.46239,F^{\prime \prime }{\left(0\right)}^3=0.09886,{\int }_0^{\infty }(1-F^{\prime })\mathrm{d}\zeta =1.1300.$$

(82)

The corresponding physical scalings are

$$\delta \propto {\ell }_{\tau }^{1/4}{\left({\displaystyle \frac{\nu x}{U}}\right)}^{1/4}\propto x^{1/4},{\displaystyle \frac{{\tau }_w}{\rho U^2}}\sqrt{{\mathit{Re}}_x}=F^{\prime \prime }\left(0\right)X^{1/4}+F^{\prime \prime }{\left(0\right)}^3{X}^{-1/4}(X\ll 1),$$

(83)

where ${\mathit{Re}}_x=Ux/\nu$. The wall shear stress therefore varies as $x^{-3/4}$ near the leading edge, more steeply than the Blasius $x^{-1/2}$, and the layer grows as $x^{1/4}$, more slowly than $x^{1/2}$: near the leading edge the shear-stiffened fluid behaves, in the boundary layer, like a fluid of higher viscosity that thins the layer and raises the friction. Figure 8 sketches the two regimes.

8.2.2. Numerical Solution and Crossover

Equation (80) was integrated by marching in X from $X={10}^{-3}$, where the leading-edge solution (81) supplies the initial profile, to $X=200$, with the streamwise derivatives discretised backward and the two-point boundary-value problem in $\eta$ at each station solved by collocation (appendix ). Table 1 and Figure 9 summarise the results. The skin-friction function $C\left(X\right)\equiv \frac{1}{2}{c}_f\sqrt{{\mathit{Re}}_x}=(1+f_{,\eta \eta }^2\left(0\right)/X)f_{,\eta \eta }\left(0\right)$ follows the leading-edge law (83) to within $1\%$ for $X<0.01$ and approaches the Blasius value $0.3321$ from below for $X>1$, reaching $0.3313$ at $X=100$. Between the two it passes through a shallow minimum, $C\simeq 0.313$ at $X\simeq 0.3$, some $6\%$ below Blasius: in the crossover region the wall shear rate has fallen below $1/\tau$, so that the stress is nearly Newtonian, while the profile is still the fuller profile inherited from the stiffened regime upstream, with the smaller wall gradient that a fuller profile of given thickness implies. The displacement thickness is monotone and reaches its Blasius value $1.7208$ (in units of $x/\sqrt{{\mathit{Re}}_x}$) from above.

The total drag of a plate of length L, $C_D=(1/L){\int }_0^L{c}_f\mathrm{d}x$, is obtained by quadrature; in the last column of Table 1 it is compared with the Blasius value $1.328{\mathit{Re}}_L^{-1/2}$. A plate shorter than ${\ell }_{\tau }$ experiences markedly more drag than Blasius—$36\%$ more at $L=0.1{\ell }_{\tau }$ and a factor $2.1$ at $L=0.01{\ell }_{\tau }$—while a plate longer than $10{\ell }_{\tau }$ is within $1\%$ of Blasius. Since drag on a plate is among the simplest fluid-mechanical measurements, this provides a direct test: plotted as $C_D\sqrt{{\mathit{Re}}_L}$ against $L/{\ell }_{\tau }=2\nu L/\left({\tau }^2{U}^3\right)$, data from plates of different lengths at different speeds must collapse onto the curve of Table 1. The collapse requires varying L and U independently; a set of measurements at fixed L and varying U would show $C_D\sqrt{{\mathit{Re}}_L}$ rising above $1.328$ as $U^3$ exceeds $2\nu L/{\tau }^2$.

8.3. Wedge Flows

For the outer flow $U=C{x}^m$ the Falkner–Skan variables are $\eta =y\sqrt{(m+1)U/\left(2\nu x\right)}$ and $\psi =\sqrt{2\nu Ux/(m+1)}f$, so that the wall-normal gradient is $u_{,y}=U\sqrt{(m+1)U/\left(2\nu x\right)}{f}_{,\eta \eta }$, and hence

$$\frac{1}{2}{\tau }^2{u}_{,y}^2=\Lambda \left(x\right)f_{,\eta \eta }^2,\Lambda \left(x\right)={\displaystyle \frac{(m+1){\tau }^2{C}^3}{4\nu }}{x}^{3m-1}.$$

(84)

Three cases arise.

Decelerating and weakly accelerating flows, $m<1/3$.$\Lambda \to \infty$ at the leading edge and $\Lambda \to 0$ downstream: the flat-plate structure of two regimes and a crossover at $x_{\tau }={[4\nu /\left((m+1){\tau }^2{C}^3\right)]}^{1/(1-3m)}$ is generic.

Strongly accelerating flows, $m>1/3$.$\Lambda$ grows downstream as $x^{3m-1}$: the flow is Newtonian near the origin and stiffens downstream, the reverse of the flat-plate ordering. The case $m=1$ (plane stagnation point, wedge half-angle $\pi /2$) is the one studied by Sun [26]; there $\Lambda \propto x^2$, and the boundary layer is not self-similar under the improved law. The self-similarity claimed for that case in Sun [26] was a consequence of the $\lambda$ term, which we have shown drops out; with the physically active $\tau$ term the Hiemenz flow at the stagnation point is unchanged and departures from it grow quadratically with distance from the stagnation line. The correct equation is

$${\left[\left(1+{\Lambda }_1{x}^2{f}_{,\eta \eta }^2\right)f_{,\eta \eta }\right]}_{,\eta }+f{f}_{,\eta \eta }+1-f_{,\eta }^2=x\left(f_{,\eta }{f}_{,\eta x}-f_{,\eta \eta }{f}_{,x}\right),{\Lambda }_1={\displaystyle \frac{{\tau }^2{C}^3}{2\nu }},$$

(85)

which replaces Equation (36) of that paper.

The exactly self-similar wedge, $m=1/3$. Equation (84) shows that $\Lambda$ is independent of x if and only if $m=1/3$, which corresponds to the Falkner–Skan pressure-gradient parameter $\beta =2m/(m+1)=1/2$, a wedge of half-angle $\pi \beta /2=\pi /4$. For this flow, and for no other wedge, the improved boundary layer is exactly self-similar:

$${\left[\left(1+\Lambda f_{,\eta \eta }^2\right)f_{,\eta \eta }\right]}_{,\eta }+f{f}_{,\eta \eta }+\frac{1}{2}\left(1-f_{,\eta }^2\right)=0,\Lambda ={\displaystyle \frac{{\tau }^2{C}^3}{3\nu }},$$

(86)

with $f\left(0\right)=f_{,\eta }\left(0\right)=0$, $f_{,\eta }(\infty )=1$. The physical reason is that in this flow the wall shear rate $u_{,y}{|}_w\propto x^{(3m-1)/2}$ is constant along the wall, so that the local ratio of the strain rate to $1/\tau$ does not change with x. Table 3 and Figure 10 give the solution. At $\Lambda =0$ the Falkner–Skan value $f_{,\eta \eta }\left(0\right)=0.9277$ is recovered. As $\Lambda$ increases the wall value of $f_{,\eta \eta }$ decreases while the wall stress $(1+\Lambda f_{,\eta \eta }^2)f_{,\eta \eta }$ increases, the profiles become fuller and the displacement thickness grows. Table 3 is the correct counterpart, for the improved theory, of table II in Sun [26]; unlike that table it is derived from equations that are consistent in both momentum components, and it applies to the $\pi /4$ wedge rather than to the stagnation-point flow.

8.4. Unsteady Boundary Layers

The similarity transformation of Sun [25] for unsteady laminar layers, $\psi =U\delta f(\eta ,\varsigma )$, $\eta =y/\delta \left(x\right)$, $\varsigma =\nu t/{\delta }^2$, applies to (76) without modification other than the replacement of $f_{,\eta \eta \eta }$ by ${\left[(1+\Lambda \left(x\right)f_{,\eta \eta }^2)f_{,\eta \eta }\right]}_{,\eta }$ with $\Lambda \left(x\right)=\frac{1}{2}{\tau }^2{U}^2/{\delta }^2$. For the $m=1/3$ wedge, $\delta \propto x^{1/3}$ and $\Lambda$ is constant, so that the unsteady problem ${\left[(1+\Lambda f_{,\eta \eta }^2)f_{,\eta \eta }\right]}_{,\eta }+f{f}_{,\eta \eta }+\frac{1}{2}(1-f_{,\eta }^2)=f_{,\varsigma \eta }$, with $f(\varsigma ,0)=0$, remains a problem in two variables. The impulsively started $\pi /4$ wedge is thus the simplest unsteady problem in which the material time enters, and its early-time behaviour—a Rayleigh layer $u_{,y}\sim U/\sqrt{\nu t}$ that is stiffened for $t<{\tau }^2{U}^2/\nu$—can be worked out from the same equation. We do not pursue it here.

8.5. Experimental Tests

The boundary-layer results reduce to a small number of measurable statements.

1.

In a shear-thickening liquid of known $\mu$ and $\tau$, the skin-friction distribution on a flat plate, plotted as $\frac{1}{2}{c}_f\sqrt{{\mathit{Re}}_x}$ against $x/{\ell }_{\tau }$, must follow the curve of Figure 9(b): $x^{-1/4}$ near the leading edge, a $6\%$ dip near $x=0.3{\ell }_{\tau }$, and Blasius beyond. Plates of different lengths at different speeds must collapse.

2.

The total drag of a plate, $C_D\sqrt{{\mathit{Re}}_L}$, must follow the last column of Table 1 as a function of $L/{\ell }_{\tau }$; this requires only a force balance.

3.

Velocity profiles measured in the stiffened regime must be fuller than Blasius when scaled by ${\delta }_{99}$, as in Figure 9(a).

4.

On a $\pi /4$ wedge the wall shear stress must be independent of x and equal to $\rho \sqrt{2\nu C^3/3}(1+\Lambda f_{,\eta \eta }^2)f_{,\eta \eta }{|}_0$ with the value taken from Table 3; on a plane stagnation point it must follow Hiemenz near the stagnation line and depart quadratically with distance from it.

5.

No dependence on $\lambda$ may be detected in any boundary-layer quantity; $\lambda$ is measurable only through the transverse pressure or the normal-stress difference $N_2$.

The last statement is a strong one. It says that the boundary layer is blind to half of the new physics, and that the search for the cross-viscosity must be conducted in pipe flow or in a rheometer, not on a flat plate.

9. Discussion

9.1. Relation to Turbulence and to Direct Numerical Simulation

The improved equations are not a theory of turbulence, and it is important to say precisely what they are and are not. Turbulence at the Reynolds numbers of laboratory and engineering interest is reproduced by direct numerical simulation of the Navier–Stokes equations to within experimental accuracy [14]. In such flows the smallest strain rates are of order ${(\epsilon /\nu )}^{1/2}$, with $\epsilon$ the dissipation rate, and the largest gradients, at the Kolmogorov scale $\eta$, are still far from any strain rate at which a simple liquid departs from Newtonian behaviour. The stiffening term of (27) is then negligible everywhere, the cross-viscosity term is a small correction to the pressure field, and the improved equations reproduce every result of Navier–Stokes DNS. The theory does not predict a change in the turbulent statistics of water or air at accessible Reynolds numbers, and it should not be read as doing so.

What the theory changes is the small-scale end of the cascade in the limit that the Navier–Stokes equations cannot handle. Whether by forcing, as in OpenAI [17], or possibly by the internal dynamics of the cascade at Reynolds numbers far beyond present reach, the Navier–Stokes equations can be driven to strain rates beyond any bound. The improved equations cannot: they possess a fixed regularisation scale ${\ell }_{min}=\sqrt{\nu \tau }$ at which the effective viscosity has grown enough to arrest the cascade, and they provide the a priori bound (33) that the Navier–Stokes equations lack. In this sense the improved theory stands to the Navier–Stokes equations as the Navier–Stokes equations stand to the Euler equations: the added term is negligible almost everywhere and decisive in the thin regions where the lower-order theory becomes singular.

The Reynolds-averaging problem is unaffected. The closure problem of turbulence is the problem of predicting means without resolving the scales that produce them, and no change in the constitutive relation removes it. If anything, the improved equations sharpen it, since they show that the small-scale dynamics is well defined and bounded, so that the difficulty is one of computation and statistics, not of the existence of the object being computed. We differ here from Sun [26], who proposed that turbulence arises from the second-order stress at a critical $\mathcal{K}$; theorem 2 shows that in two dimensions that stress does nothing, and in three dimensions the transition to turbulence in a Newtonian liquid is, on all the evidence, a property of the Navier–Stokes equations themselves.

9.2. Regime of Validity and the Physical Origin of $\tau$

The improved law is a continuum constitutive relation and shares the assumptions of any such relation: local equilibrium, scale separation between the molecular and the hydrodynamic, and absence of memory. The material time $\tau$ is the strain rate scale at which the fluid’s microstructure—for a simple liquid, its molecular arrangement; for a suspension, the particle configuration—can no longer relax within a strain period, so that the stress rises faster than linearly. In dense suspensions this is the well-known shear-thickening transition, at rates of 10–${10}^3$ s⁻¹, and there all of the predictions of §§Section 7Section 8 are directly testable. In simple liquids $\tau$ is bounded above by roughly ${10}^{-6}$ s by existing rheometry, and the corresponding regularisation length $\sqrt{\nu \tau }\lesssim 1\mu$m is well above the molecular scale; the continuum description therefore remains valid at ${\ell }_{min}$, and the cut-off it provides intervenes before compressibility or molecular-scale effects, which set in at scales of ${10}^{-9}$ m and speeds of ${10}^3$ m s⁻¹. Whether the actual cut-off in a given liquid is set by the constitutive law, as modelled here, or by compressibility or molecular effects intervening first, is an experimental question that the present theory cannot settle; what the theory establishes is that a cut-off of the constitutive kind, if present, is sufficient to remove the singularity, and it gives the form that such a cut-off must take.

Memory effects have been excluded. A fluid whose stress depends on the history of $\mathbf{D}$ (a viscoelastic liquid) is described by the Rivlin–Ericksen or Oldroyd families, in which the second-order term involves the convected derivative of $\mathbf{D}$ and $N_1\ne 0$. For such fluids the present law is the zero-relaxation-time limit, and the inequality $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$ should be regarded as a test of the applicability of the memoryless description.

9.3. Extensions

The developments of Sun [26], appendices B and C, for compressible flow and for heat transfer carry over. For a compressible fluid the incompressibility constraint is removed, $I_1=\nabla ·\mathbf{u}\ne 0$, and (15) acquires the bulk term ${\alpha }_0(I_1,I_2,I_3)\mathbf{I}$; the natural generalisation of (24) is

$$\mathbf{\sigma }=-p\mathbf{I}+\left(\zeta -\frac{2}{3}\mu \right)(tr\mathbf{D})\mathbf{I}+2\mu \left(1+{\tau }^2{\mathbf{D}}^{\prime }:{\mathbf{D}}^{\prime }\right)\mathbf{D}+4\lambda {\mathbf{D}}^{\prime }·{\mathbf{D}}^{\prime },$$

(87)

with ${\mathbf{D}}^{\prime }$ the deviator of $\mathbf{D}$ and $\zeta$ the bulk viscosity, and theorems 3 and 4 apply to the deviatoric part unchanged. In the energy equation the dissipation function is $\Phi$ of (25), positive by theorem 3, so that the temperature equation is

$$\rho c_p\left[{\displaystyle \frac{\partial T}{\partial t}}+\mathbf{u}·\nabla T\right]=\nabla ·(k\nabla T)+2\mu \left(1+{\tau }^2{\left|\mathbf{D}\right|}^2\right){\left|\mathbf{D}\right|}^2+4\lambda tr{\mathbf{D}}^3;$$

(88)

in a boundary layer the viscous heating is $\mu (1+\frac{1}{2}{\tau }^2{u}_{,y}^2)u_{,y}^2$, enhanced in the stiffened regime. The Reynolds-averaged form of (27) contains, in addition to the Reynolds stress, the mean of the nonlinear viscous stress, $\overline{{\left|\mathbf{D}\right|}^2\mathbf{D}}\ne {\left|\overline{\mathbf{D}}\right|}^2\overline{\mathbf{D}}$; in flows where $\tau \left|\mathbf{D}\right|\ll 1$ this is negligible and the standard RANS closures apply unchanged.

The mathematical theory can be extended in several directions that we have not pursued: to bounded domains with no-slip walls, for which the results of Málek et al. [13] and subsequent work on p-structure fluids apply with additional technicalities in the regularity argument; to the whole space ${\mathbb{R}}^3$, the setting of statement C of the Clay problem; and to general growth ${\alpha }_1\sim {\left|\mathbf{D}\right|}^{p-2}$ with $5/2\le p<4$, for which theorem 4 must be re-proved with the appropriate constants.

9.4. What Has and Has Not Been Proved

We have been careful to distinguish the status of the various results. Lemma 1, corollary 1 and theorems 1, 2, 3 and 4 are proved completely in the text, and the last has been verified independently by numerical minimisation. The exact solutions of §Section 7 and the reductions of §Section 8 are elementary. Theorem 5 is proved by verifying that the hypotheses of the classical monotone-operator theory hold and by reference to that theory for the passage to the limit and for the strong-solution estimate; the uniqueness argument is given in full. A reader who wishes to see the theorem proved without reference to the literature will find every step in Lions [11] and Málek et al. [13]; a formalisation in a proof assistant, in the manner of OpenAI [17], would be a natural next step and we would welcome it. The numerical results of §§Section 7Section 8 have been checked against their analytic limits (Poiseuille, Burgers, Blasius, Falkner–Skan, the leading-edge law) and are believed accurate to the figures quoted.

We have not proved that any real fluid obeys (24); that is a matter for experiment, and §§Section 7.3 and Section 8 specify the experiments. We have not proved that the unforced Navier–Stokes equations blow up, which remains open, nor that the improved equations are the only admissible completion—higher even powers of $\left|\mathbf{D}\right|$ in ${\alpha }_1$ are equally admissible and differ only at strain rates beyond $1/\tau$. What we have proved is that the constant-coefficient second-order model is not admissible, that any admissible completion must stiffen the linear coefficient, that the lowest-order such completion is globally well posed, and that its consequences for pipe flow, rheometry, vortex dynamics and boundary layers are definite and measurable.

10. Conclusions

Starting from the Cauchy momentum equation and asking what constitutive relation a fluid must obey if its velocity gradients are to remain bounded, we have arrived at

$$\mathbf{\sigma }=-p\mathbf{I}+2\mu \left(1+{\tau }^2\mathbf{D}:\mathbf{D}\right)\mathbf{D}+4\lambda \mathbf{D}·\mathbf{D},\tau \ge \sqrt{\frac{2}{3}}{\displaystyle \frac{\left|\lambda \right|}{\mu }},$$

(89)

and the governing equations (27). The principal results are the following.

1.

Frame indifference, isotropy and the Cayley–Hamilton theorem restrict the extra stress to ${\alpha }_1\mathbf{D}+{\alpha }_2\mathbf{D}·\mathbf{D}$ with invariant-dependent coefficients, and the second law requires ${\alpha }_1\ge |{\alpha }_2\left|\right|\mathbf{D}|/\sqrt{6}$: a second-order term cannot be added without stiffening the first-order term.

2.

The constant-coefficient second-order model violates the second law in three dimensions, is an exact gradient in two dimensions, and cannot bound the velocity gradient. Its two-dimensional boundary-layer predictions are artefacts of an inconsistent reduction.

3.

For the improved law, thermodynamic admissibility requires $\tau \ge \left|\lambda \right|/\left(\sqrt{6}\mu \right)$ and monotonicity requires $\tau \ge \sqrt{2/3}\left|\lambda \right|/\mu$; the second bound is sharp, is exactly twice the first, and is attained in a stretched vortex.

4.

Under the monotonicity condition the momentum operator is strongly monotone with growth exponent $p=4$, the equations have unique global strong solutions on the torus, the enstrophy is bounded on every finite interval, and the forced finite-time blow-up of the Navier–Stokes equations cannot occur. The regularisation scale is $\sqrt{\nu \tau }$.

5.

In simple shear $\eta =\mu (1+\frac{1}{2}{\tau }^2{\dot{\gamma }}^2)$, $N_1=0$, $N_2=\lambda {\dot{\gamma }}^2$. In gravity-driven pipe flow the velocity profile depends on $\tau$ alone and the transverse pressure on $\lambda$; the two are determined by two measurements in one apparatus and must satisfy the monotonicity inequality.

6.

In a stretched vortex the cross-viscosity reduces and the stiffening increases the effective viscosity; the peak vorticity saturates instead of growing linearly with the strain.

7.

The boundary-layer equations are Prandtl’s with the viscometric viscosity; $\lambda$ drops out identically. The flat-plate layer has a leading-edge regime, $\delta \propto x^{1/4}$, ${\tau }_w\propto x^{-3/4}$, with $F^{\prime \prime }\left(0\right)=0.4624$, crossing over at ${\ell }_{\tau }={\tau }^2{U}^3/\left(2\nu \right)$ to Blasius. Among wedge flows only $U\propto x^{1/3}$ is self-similar, with the wall stress function tabulated. The plane stagnation point is not self-similar.

The theory is complete in the sense that it contains no free functions, only two material constants beyond the viscosity, both measurable by the protocols given, and it makes predictions—for plate drag, for skin friction near a leading edge, for the flow rate and transverse pressure in a pipe, for the wall stress on a $\pi /4$ wedge—that can confirm or refute it with existing apparatus.