In the context of flow visualization, a triple decomposition of the velocity gradient into irrotational straining flow, shear flow, and rigid body rotational flow was proposed by Kolář in 2007 [V. Kolář, “Vortex identification: New requirements and limitations,” Int. J. Heat Fluid Flow, **28**, 638–652 (2007)], which has recently received renewed interest. The triple decomposition opens for a refined energy stability analysis of the Navier–Stokes equations, with implications for the mathematical analysis of the structure, computability, and regularity of turbulent flow. We here perform an energy stability analysis of turbulent incompressible flow, which suggests a scenario where at macroscopic scales, any exponentially unstable irrotational straining flow structures rapidly evolve toward linearly unstable shear flow and stable rigid body rotational flow. This scenario does not rule out irrotational straining flow close to the Kolmogorov microscales, since there viscous dissipation stabilizes the unstable flow structures. In contrast to worst case energy stability estimates, this refined stability analysis reflects the existence of stable flow structures in turbulence over extended time.

By the Navier–Stokes equations, turbulent incompressible flow can be described in terms of its velocity and pressure fields. The velocity gradient tensor expresses the local spatial change of the velocity field and can be used, for example, to identify and visualize vortex structures. Yet, the velocity gradient tensor also plays a key role in a flow stability analysis.

In a domain $\Omega \u2282\mathbb{R}3$ over a time interval $[0,T]$, incompressible flow with constant density $\rho >0$ is modeled by the Navier–Stokes equations,

for $(x,t)\u2208\Omega \xd7[0,T]$, together with boundary conditions. Here, $u=u(x,t)=(u1(x,t),u2(x,t),u3(x,t))$ is the velocity vector, and $\rho f=\rho f(x,t)$ is the external force vector. For a Newtonian fluid, the Cauchy stress tensor *σ* is defined by $\sigma =2\mu S(u)\u2212p\xafI,$ with $p\xaf=p\xaf(x,t)$ the mechanical pressure, *μ* the dynamic viscosity, and

is the strain rate tensor. Since $\u2207\xb7u=0$, it follows that:

With the kinematic viscosity $\nu =\mu /\rho $, and the scaled pressure $p=p\xaf/\rho $, we write the Navier–Stokes equations as

If (*u*, *p*) is sufficiently regular to satisfy Eqs. (2)–(3) for each $(x,t)\u2208\Omega \xd7[0,T]$, then it is referred to as a classical solution, or strong solution. Hence, with (*u*, *p*) a classical solution we can take the scalar product of Eq. (2) with *u* and multiply Eq. (3) by *p*, then add them together and integrate over the domain Ω. By the divergence theorem, we obtain the well-known global energy balance

We here neglect the contribution from the boundary $\u2202\Omega $, using a standard assumption that the solution decays toward the boundary. If we include boundary effects in the analysis, additional boundary integrals appear. The energy balance (5) expresses the injection of kinetic energy by external forces and loss of kinetic energy by viscous dissipation. The energy dissipation can equivalently be expressed as

where $\omega =\u2207\xd7u$ is the vorticity vector. The spin tensor,

together with the strain rate tensor (1) constitute the standard double decomposition of the velocity gradient,

The evolution of vorticity is governed by the equation

obtained by taking the curl of the Navier–Stokes equations, and analogous to the balance of kinetic energy (5) we can take the scalar product of the equation with the vorticity vector and integrate, to get

The dissipation term is analogous in the energy balance (5), but there is also a new term with properties determined by the velocity gradient. For example, a velocity gradient, which reflects a straining flow, leads to an exponential growth of vorticity referred to as vortex stretching, which in turbulence is a mechanism that is responsible for the transfer of energy to smaller scales.

Vorticity growth is also used as a criterion to detect blow-up of classical solutions to the Navier–Stokes equations.^{1} Whether or not it always exists a classical solution is one of the Clay Institute $1 million Prize problems.^{2} If blow-up can occur, that would support a mathematical model of turbulence in the form of weak solutions that satisfy the equations only in an average sense,^{3} and which dissipate energy independent of viscosity through the appearances of singularities,^{4} in line with the Onsager conjecture.^{5}

Very similar in structure to the vorticity, Eq. (9) is the system of pertubation equations, which describe the evolution of a velocity perturbation $\phi =\phi (x,t)$ and a pressure perturbation $\theta =\theta (x,t)$ in a flow field (*u*, *p*) governed by the Navier–Stokes equations,

In the vorticity equation, the pressure term from the Navier–Stokes equations disappears because $\u2207\xd7\u2207p=0$, and the vorticity field is divergence free since $\u2207\xb7\omega =\u2207\xb7\u2207\xd7u=0$. The perturbation equations are studied in the fields of hydrodynamic stability and flow control, together with the associated adjoint Navier–Stokes equations.^{6}

In the context of adaptive numerical methods, the adjoint Navier–Stokes equations are used for goal-oriented *a posteriori* error estimation^{7} and take the following form:

with $\lambda =\lambda (x,t)$ the adjoint velocity vector, and $\pi =\pi (x,t)$ the scalar adjoint pressure. Here, $U=U(x,t)$ is a computed approximation of the velocity, and $\u2207UT$ is the transpose of $\u2207U$. The goal of the computation is an output of interest, which is defined in terms of sources or boundary conditions in the adjoint equations.^{8} Note that the adjoint equations evolve backward in time from a final condition at *t *=* T*, but with the change of variables $s=T\u2212t$ we instead obtain a problem that evolves from *s *=* *0 to *s *=* T*,

The similarities between the perturbation equations, the adjoint equations, and the vorticity equation are clear. Specifically, in an energy analysis they all share the same stability properties. From the perturbation equations, we obtain an energy balance for the perturbations, the Reynolds–Orr energy equation,^{9}

and similar for the adjoint equations (12),

We recall that $\u2207U\u2248\u2207u$, and by the definition of the scalar product, we have that

Therefore, the key to any stability analysis is the term

To analyze *P*(*u*), a natural starting point is the classical decomposition (8) into the symmetric strain rate tensor *S*(*u*) and the skew-symmetric spin tensor $\Omega (u)$. By Eqs. (8) and (15), it follows that:

That is, only *S*(*u*) contributes to the critical term (16).

To better understand the stability of different flow structures, we recall three simple archetypical solutions to the Navier–Stokes equations. Let $(urr,prr)$ denote rigid body rotational flow about the *x*_{1} coordinate axis,

an example for which $S(urr)=0$ but $\Omega (urr)\u22600$. In contrast, irrotational straining flow (elongation and compression) in the $x1x2$ plane is an example for which $\Omega (uel)=0$ but $S(uel)\u22600$,

This follows from inspection of the velocity gradients:

One might then be lead to believe that Eq. (8) decomposes the flow in every point of the domain into a sum of rigid body rotational flow and irrotational straining flow. This, however, does not take shear flow into consideration. An example of shear flow in the $x2x3$ plane is

with gradient

which contributes to both the strain rate tensor and the spin tensor,

Hence, the decomposition (8) is not able to distinguish between the three archetypes, each with its own unique properties. From Eq. (17), we conclude that *u _{rr}* does not contribute to the critical term (16); therefore,

*u*is stable in the sense that the only active terms in the energy balances (10), (13), and (14) are the dissipative terms. The critical term (16) for the three examples takes the forms

_{rr}To determine the stability of *u _{el}* and

*u*, it is instructive to analyze the energy in each vector component individually. In the case of the vorticity, we take the scalar product of Eq. (9) with each of the vectors

_{sh}and then integrate over the domain Ω, to get the energy balance for each component of the vorticity vector. For the irrotational straining flow (19), this leads to

from which we detect exponential growth of the energy in *ω*_{1}, if it is not damped by the viscous dissipation. The growth is aligned with the direction of elongation in *u _{el}*, reflecting the well-known vortex stretching mechanism. Note that

*ω*

_{1}may represent either rigid body rotational flow or shear flow, or a combination thereof.

A similar componentwise analysis of the perturbations gives

with two notable differences from the analysis of the vorticity. First, the critical term (16) has the opposite sign; hence, exponential growth here appears in the second component $\phi 2$. Second, the individual components of the pressure perturbation gradient appear, even though they cancel when summed up in the Reynolds-Orr equation due to the fact that $\phi $ is divergence free.

In the case of the shear flow (20), the analogous stability analysis gives

and

from which we observe that no exponential growth takes place, but instead potential linear growth of the energy in *ω*_{2} and $\phi 2$, with *x*_{2} the main flow direction in Eq. (20).

We conclude that the irrotational straining flow (19) is exponentially unstable, that the shear flow (20) is linearly unstable, and that the rigid body rotational flow (18) is stable. These archetypical solutions to the Navier–Stokes equations are clearly very simple, but we will find that they in their simplicity hold the essence of the stability properties of a real turbulent flow.

The discrepancy between the double decomposition (8) and the three archetypical flow structures illustrates a problem of both practical and pedagogical nature. The spin tensor does not isolate rigid body rotational flow but may also include shear flow, which is a well-known complication in the context of flow visualization.^{10–12} For example, in a shear layer the spin tensor is large, while there may be no rigid body rotational flow what so ever. A pedagogical problem arises in the understanding of vortex lines, tubes, filaments, and sheets, all defined in terms of vorticity that can mean either shear flow or rigid body rotational flow, or a combination of the two.

In 2007, Kolář^{13} proposed a triple decomposition of the velocity gradient tensor, which corresponds to the three archetypes of rigid body rotational flow, irrotational straining flow, and shear flow. The key idea is to first use an orthogonal matrix *Q* and its transpose $QT=Q\u22121$ to transform the velocity gradient tensor into a frame of reference where a shear flow tensor $(\u2207u)SH$ can be isolated

after which the remaining residual tensor $(\u2207u)RES$ is split into a symmetric tensor $(\u2207u)EL$ and a skew-symmetric tensor $(\u2207u)RR$, corresponding to irrotational straining flow and rigid body rotational flow, respectively,

This amounts to the triple decomposition

which states that the velocity gradient tensor can be decomposed into a sum of irrotational straining flow, rigid body rotational flow, and shear flow.

Equation (25) follows the standard double decomposition (8), whereas the decomposition in Eq. (24) needs to be defined. Kolář gave the following definition of the residual tensor,

with $[\xb7]ij$ each tensor component and $\u2207u\xaf=Q\u2009\u2207u\u2009QT$, where *Q* should be determined such that the residual tensor is minimized in the Frobenius norm $||\xb7||F$, referred to as the *basic reference frame*. By the properties of the norm and the orthogonal matrix *Q*, minimization of the residual tensor is equivalent to maximization of the quantity

since $\u2207u$ is independent of *Q* and

Kolář first solved this optimization problem for two-dimensional flow^{13} and later also for three-dimensional flow.^{14} More recently, new algorithms have been proposed to solve the same optimization problem.^{15,16}

Another strategy is to identify the shear flow tensor with the non-normal part of the velocity gradient tensor in which case Eq. (24) follows from the Schur decomposition of the velocity gradient tensor. The Schur decomposition of a general (possibly complex) 3 × 3 matrix *A* takes the form

with *U* a unitary matrix with adjoint $U*=U\u22121$, *T* an upper triangular matrix with the eigenvalues of *A* on its diagonal, *D* the diagonal part of *T*, and $R=T\u2212D$ its strictly upper triangular part. The eigenvalues $\lambda 1,\lambda 2,\lambda 3$ are unique, but how they are ordered on the diagonal of *T* is not. *U* is not unique, and neither is *R* although the Frobenius norm of *R* is, since by Eq. (28),

with $\sigma 1,\sigma 2,\sigma 3$ the singular values of *A*. If *A* is a normal matrix then *R *=* *0. Hence, the matrix $U*DU$ represents the normal part of *A*, and $U*RU$ the non-normal part. By identifying the non-normal part of the velocity gradient tensor with the shear flow tensor $(\u2207u)SH=R$, the Schur decomposition can be used to construct the triple decomposition.^{17}

$(\u2207u)SH=R$ also satisfies Eq. (27), possibly using complex arithmetic if the velocity gradient tensor has complex eigenvalues, corresponding to its skew-symmetric part being non-zero $(\u2207u)RR\u22600$. Complex eigenvalues of a real matrix always exist as pairs of complex conjugates $(\lambda ,\lambda \xaf)=(a+ib,a\u2212ib)$. Complex arithmetic can be avoided by using the real Schur decomposition^{18} of the velocity gradient tensor, where each pair of complex eigenvalues is represented by the real 2 × 2 block

on the diagonal of the matrix *D*. Eq. (27) is satisfied for the triple decomposition based on the real Schur decomposition, since the block (30) is symmetric in absolute value. If $(\u2207u)RR=0$, the complex and real Schur decompositions are equivalent. The Schur decomposition, real or complex, can be constructed, for example, by the QR algorithm, which, therefore, offers an alternative to methods based on the optimization problem.^{17,19} But the relation between the Schur transformations and the basic reference frame is still an open question.

The non-uniqueness of the Schur decomposition is not as serious as one may first think. If $(\u2207u)RR\u22600$, the flow exhibits rigid body rotation and the velocity gradient tensor has a pair of complex conjugate eigenvalues, and an eigenvector of the remaining real eigenvalue defines the axis of rotation. Hence, the non-uniqueness of the frame of reference is merely a matter of choosing a coordinate system aligned with the axis of rotation.^{20} In the case $(\u2207u)RR=0$ and the velocity gradient tensor is normal, the frame of reference is given by eigenvectors of the three real eigenvalues corresponding to $(\u2207u)EL$.

Both $(\u2207u)RR$ and $(\u2207u)SH$ have zero diagonal components, and $(\u2207u)EL$ is a symmetric tensor with zero trace, which follows from the divergence free condition of incompressible flow, the cyclic property of the trace operation, and the definition of an orthogonal matrix,

Hence, the eigenvalues of the matrix $(\u2207u)EL$ are real and sum to zero, which means that either all eigenvalues are zero, or there will always be at least one eigenvalue that will generate exponential growth. With the Schur decomposition, $(\u2207u)SH$ is a non-normal strictly triangular matrix, which generates linear growth. The tensor $(\u2207u)RR$ is skew-symmetric with purely imaginary eigenvalues, which corresponds to a stable vortex.

By Taylor's formula and Eq. (26), near each interior point $x0\u2208\Omega $ we can construct a linear approximation of the velocity field,

a decomposition of the local velocity field defined by

where each part represents constant flow, irrotational straining flow, rigid body rotational flow, and shear flow, respectively. Each part has its own specific stability properties, analogous to the archetypes (19), (18), and (20).

In the energy balances (10), (13), and (14), the growth generated by strain and shear may be damped by viscous dissipation, but at macroscopic turbulent scales the critical term (31) will dominate. Hence, the stability analysis suggests that at macroscopic turbulent scales, rigid body rotational flow can exist for a long time, shear flow for a short time, and irrotational straining flow can only appear as intermittent very short-lived events. The decomposition (32) always exists, assuming the velocity field is regular enough, and therefore implies that the flow evolves from exponentially unstable irrotational straining flow and linearly unstable shear flow, toward stable rigid body rotational flow. Mechanisms for this evolution include vortex stretching by strain, and shear layer roll up by the Kelvin–Helmholtz instability.^{15,21,22}

At microscopic turbulent scales, close to the Kolmogorov microscales, viscous dissipation is significant. A detailed understanding of how energy is dissipated in turbulent flow is still lacking, but the stability analysis suggests that irrotational straining flow can be sustained longer at the microscopic scales, since it is stabilized by viscous dissipation. Experimental data also support the prediction that flow structures at microscopic turbulent scales involve strain, for example, elongation of vortices (vortex stretching) and compression of shear layers (strain self-amplification).^{15,21,23–25}

The triple decomposition of the velocity gradient tensor (26) together with the linearization (32) gives a characterization of the local velocity field as a sum of a constant part, a strain part, a rigid body rotational part, and a shear part. In the energy stability analysis, the critical term (31) exhibits the stability property of each part, where strain corresponds to exponential growth and shear to linear growth. At macroscopic turbulent scales viscous dissipation can be neglected, which suggests that straining flow cannot exist other than for brief instants and shear flow only for a short time. Rigid body rotation, on the other hand, is a stable flow structure. In contrast, at microscopic turbulent scales viscous dissipation permits flow structures that include strain and shear to exist over extended time.

Experimental observations appear to support the predictions of the stability analysis that strain and shear flow are short-lived at macroscopic turbulent scales, to instead evolve into vortex structures through vortex stretching and shear layer roll up.^{15,21,22} Close to the Kolmogorov microscales, experimental data recently reported identifies intense vortex stretching.^{23,24} A detailed characterization of the finest turbulent flow structures could give clues to how energy is dissipated^{26–29} and help to solve the outstanding mathematical questions of Onsager's conjecture^{5} and the blow-up problem.^{2}

The focus of this note is turbulent structures away from boundaries, but the decomposition (32) is valid also near boundaries. It is well known that as the Reynolds number increases, shear flow in laminar boundary layers transition into turbulent boundary layers populated by vortex structures.

In a large eddy simulation, the turbulent scales are truncated by the resolution of the computational mesh, typically in the inertial range far from the Kolmogorov microscales. Hence, the stability analysis predicts that the resolved turbulent flow will be dominated by a collection of stable rigid body rotational flow structures, vortex tubes; and linearly unstable shear flow structures, shear pancakes, which either dissipate or roll up to form vortex tubes. Strain is short-lived and only appears intermittently in scenarios where the other two types of flow structures are intensified. Computability of mean quantities in turbulent flow can be formulated in terms of *a posteriori* error analysis and the adjoint equations.^{30,31} In the context of adaptive finite element methods, stable solutions to the adjoint equations have been computed for a wide variety of examples of fully developed turbulent flow.^{8,32–36}

Based on the stability analysis in this note, we conjecture that turbulence is a less chaotic phenomenon than is commonly believed, in the sense that its elements are robust and predictable, vortex tubes and shear pancakes at macroscopic scales, which are intensified by strain to be dissipated at microscopic scales. Therefore, it is not surprising that mean quantities that do not depend on the detailed dynamics of these turbulent flow structures, such as the drag of a car, can be predicted in computer simulations, which is also reflected in the adjoint equations.

Apart from this general characterization of turbulent flow, the triple decomposition can be used in an energy stability analysis of model solutions to the Navier–Stokes equations, such as a Burgers vortex,^{37} which exhibits the stability properties of combined straining flow, rigid body rotational flow, shear flow and viscous dissipation.^{38}

## ACKNOWLEDGMENTS

This work is supported by the Swedish Research Council (Contract No. 2018-04854). The author would like to thank Dr. Václav Kolář for a helpful discussion regarding the optimization problem definition of the triple decomposition.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.