Energy stability analysis of turbulent incompressible flow based on the triple decomposition of the velocity gradient tensor

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 Kolar in 2007 [V. Kolar, International journal of heat and fluid flow, 28, 638, (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 towards 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. But the velocity gradient tensor also plays a key role in a flow stability analysis.
If (u, p) is sufficiently regular to satisfy equations (2)-(3) for each (x,t) ∈ Ω × [0, T ], then it is referred to as a classical solution, or strong solution. Hence, with (u, p) a classical a) http://www.kth.se/profile/~jhoffman solution we can take the scalar product of equation (2) with u and multiply equation (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 ∂ Ω, using a standard assumption that the solution decays towards 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 ω = ∇ × u 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 equatioṅ ω(x, 0) = ω 0 (x), 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 here 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 blowup 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 equation (9) is the system of pertubation equations, which describe the evolution of a velocity perturbation ϕ = ϕ(x,t) and a pressure perturbation θ = θ (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 ∇ × ∇p = 0, and the vorticity field is divergence free since ∇ · ω = ∇ · ∇ × u = 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 λ = λ (x,t) the adjoint velocity vector, and π = π(x,t) the scalar adjoint pressure. Here U = U(x,t) is a computed approximation of the velocity, and ∇U T is the transpose of ∇U. 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 backwards in time from a final condition at t = T , but with the change of variables s = T −t we instead obtain a problem that evolves from s = 0 to s = T , λ (x, 0) = λ T (x).
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 , (13) and similar for the adjoint equations (12), Recall that ∇U ≈ ∇u, 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 Ω(u). By equation (8) and equation (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 (u rr , p rr ) denote rigid body rotational flow about the x 1 coordinate axis, an example for which S(u rr ) = 0 but Ω(u rr ) = 0. In contrast, irrotational straining flow (elongation and compression) in the x 1 x 2 plane is an example for which Ω(u el ) = 0 but S(u el ) = 0, This follows from inspection of the velocity gradients, One might then be lead to believe that equation (8) decomposes the flow in every point of the domain into a sum of rigid body rotational motion and irrotational straining motion. This, however, does not take shear flow into consideration. An example of shear flow in the x 2 x 3 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 equation (17) we conclude that u rr does not contribute to the critical term (16), therefore, u rr 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 To determine the stability of u el and u sh , it is instructive to analyze the energy in each vector component individually. In the case of the vorticity, we take the scalar product of equation (9) with each of the vectors 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 perturbation equations (11) 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 ϕ 2 .
Second, the individual components of the pressure perturbation gradient appear, even though they cancel when summed up in the Reynolds-Orr equation (13) due to the fact that ϕ is divergence free.
In the case of the shear flow (20), the analogous stability analysis gives from which we observe that no exponential growth takes place, but instead potential linear growth of the energy in ω 2 and ϕ 2 , with x 2 the main flow direction in equation (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][11][12] . For example, in a shear layer the spin tensor is large while there may be no rigid rotation 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 Q T = Q −1 to transform the velocity gradient tensor into a frame of reference where a shear flow tensor (∇u) SH can be isolated, after which the remaining residual tensor (∇u) RES is split into a symmetric tensor (∇u) EL and a skew-symmetric tensor (∇u) RR , corresponding to irrotational strain and rigid body rotation, respectively, This amounts to the triple decomposition which states that the velocity gradient tensor can be decomposed into a sum of irrotational strain, rigid body rotation and shear flow. Equation (25) follows the standard double decomposition (8), whereas the decomposition in equation (24) needs to be defined. Kolář gave the following definition of the residual tensor, with [·] i j each tensor component and ∇u = Q ∇u Q T , where Q should be determined such that the residual tensor is minimized in the Frobenius norm · 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 . 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 equation (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 −1 , T an upper triangular matrix with the eigenvalues of A on its diagonal, D the diagonal part of T , and R = T − D its strictly upper triangular part. The eigenvalues λ 1 , λ 2 , λ 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 equation (28), with σ 1 , σ 2 , σ 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 (∇u) SH = R, the Schur decomposition can be used to construct the triple decomposition 17 .
(∇u) SH = R also satisfies equation (27), possibly using complex arithmetic if the velocity gradient tensor has complex eigenvalues, corresponding to its skew-symmetric part being non-zero (∇u) RR = 0. Complex eigenvalues of a real matrix always exist as pairs of complex conjugates (λ ,λ ) = (a + ib, a − ib). Complex arithmetic can be avoided by using the the real Schur decomposition 18 of the velocity gradient tensor, where each pair of complex eigenvalues are represented by the real 2 × 2 block on the diagonal of the matrix D. Equation (27) is satisfied for the triple decomposition based on the real Schur decomposition, since the block (30) is symmetric in absolute value. If (∇u) RR = 0 the complex and real Schur decompositions are equivalent. The Schur decomposition, real or complex, can be constructed e.g. 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 (∇u) RR = 0 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 (∇u) RR = 0 and the velocity gradient tensor is normal, the frame of reference is given by eigenvectors of the three real eigenvalues corresponding to (∇u) EL .
In light of the triple decomposition (26), the critical term (16) can be decomposed as follows, Both (∇u) RR and (∇u) SH have zero diagonal components, and (∇u) 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 (∇u) 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, (∇u) SH is a non-normal strictly triangular matrix, which generates linear growth. The tensor (∇u) RR is skew-symmetric with purely imaginary eigenvalues, which corresponds to a stable vortex.
By Taylor's formula and equation (26), near each interior point x 0 ∈ Ω 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, towards 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 scale, 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 supports 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][24][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 scale, 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][27][28][29] , and help to solve the outstanding mathe-matical 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 know 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 scale. 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][33][34][35][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 strain, rigid rotation, shear flow and viscous dissipation 38 .

I. ACKNOWLEDGEMENTS
The following article has been accepted for publication in Physics of Fluids. This work is supported by the Swedish Research Council (contract number 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.