The triple decomposition of a velocity gradient tensor provides an analysis tool in fluid mechanics by which the flow can be split into a sum of irrotational straining flow, shear flow, and rigid body rotational flow. In 2007, Kolář formulated an optimization problem to compute the triple decomposition [V. Kolář, “Vortex identification: New requirements and limitations,” Int. J. Heat Fluid Flow **28**, 638–652 (2007)], and more recently, the triple decomposition has been connected to the Schur form of the associated matrix. We show that the standardized real Schur form, which can be computed by state of the art linear algebra routines, is a solution to the optimization problem posed by Kolář. We also demonstrate why using the standardized variant of the real Schur form makes computation of the triple decomposition more efficient. Furthermore, we illustrate why different structures of the real Schur form correspond to different alignments of the coordinate system with the fluid flow and may, therefore, lead to differences in the resulting triple decomposition. Based on these results, we propose a new, simplified algorithm for computing the triple decomposition, which guarantees consistent results.

In recent years, the triple decomposition of a velocity gradient tensor has emerged as a new analysis tool in fluid mechanics, originally proposed in the context of flow visualization^{1} and later in the analysis of the Navier–Stokes equations and turbulence.^{2} The triple decomposition is used to split the velocity gradient tensor into a straining flow part, a rigid body rotation part, and a shear flow part, thereby revealing useful information about the fluid flow. Applications include vortex identification, where it has been demonstrated to prevent shear contamination that other methods suffer from;^{3} statistical analysis of small-scale properties of the velocity gradient tensor in turbulence;^{4} and analysis of shear in blood flow simulations, to improve assessment of risk for thrombosis formation.^{5}

While the usefulness of the triple decomposition has been clearly demonstrated, different interpretations of the triple decomposition have been suggested. Specifically, we can distinguish between the triple decomposition, on the one hand, being defined in terms of an optimization problem^{1,6,7} and, on the other hand, as a Schur form of the velocity gradient tensor.^{8,9} As the first main result of this paper, we show that, in fact, the two interpretations are equivalent in the sense that a standardized real Schur form of the velocity gradient tensor is a solution to the optimization problem. Furthermore, we show that the standardized real Schur form can be directly decomposed into the triple components, without requiring additional rotation, used, e.g., in the Rortex algorithm by Liu *et al.*^{8}

In the second part of this paper, we demonstrate how different structures of the real Schur form have different physical interpretations. The real Schur form always aligns the coordinate system with the fluid flow, but the alignment can be either with the flow rotation axis or with the flow rotation plane. The magnitude of the rigid body rotation measure in these two cases can differ. We argue that the real Schur form corresponding to alignment with the flow rotation axis should always be used in order to guarantee consistent results. Finally, we conclude by suggesting a new simplified algorithm for computing the triple decomposition.

To explain why the real Schur form also provides a solution to the optimization problem, consider the velocity gradient tensor $\u2207u$, which expresses the local spatial change of the velocity field, in the form of a 3 × 3 matrix defined by its components

The classical double decomposition of the velocity gradient tensor takes the following form:

defined in terms of the symmetric strain rate tensor,

and the skew-symmetric spin tensor,

The operation is purely algebraic, valid for any square matrix.

Although the double decomposition is used extensively in fluid mechanics to distinguish between straining and rotating flows, it has a significant flaw: it does not separate shear flow from strain and rotation. In fact, shear flow may contaminate both the strain and rotation components.^{3} This has important consequences, both for the analysis of the Navier–Stokes equations^{2} and for the processing of fluid dynamics data. To address this flaw, in the context of flow visualization, Kolář proposed a triple decomposition of the velocity gradient tensor^{1} into the components of irrotational straining flow, rigid body rotational flow, and shear flow. This triple decomposition was shown to improve vortex identification, which otherwise can be a challenge in parts of the flow with high shear.

The key idea of Kolář was to first construct an orthogonal matrix *Q* to transform the velocity gradient tensor into a frame of reference where a shear flow tensor $(\u2207u)SH$ can be isolated,

Note that this transformation is not unique. In a second step, the remaining residual tensor $(\u2207u)RES$ is split into a symmetric tensor $(\u2207u)EL$ (_{EL} for elongation) and a skew-symmetric tensor $(\u2207u)RR$, which correspond to the irrotational straining flow tensor and the rigid body rotational flow tensor, respectively,

This leads to the triple decomposition

Equation (6) follows the construction of the double decomposition (2), but the determination of the orthogonal matrix *Q* in Eq. (5) is less obvious. The frame of reference defined by *Q* is denoted the basic reference frame (BRF), and Kolář formulated the following optimization problem to determine *Q*: find *Q* such that $(\u2207u)RES$ is minimized in the Frobenius norm $||\xb7||F$, where the residual tensor is defined by

with $\u2207u\xaf=QT\u2009\u2207u\u2009Q$. It follows that minimization of $(\u2207u)RES$ is equivalent to maximization of the quantity

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

This optimization problem was solved by Kolář, first for two-dimensional flow^{1} and later for three-dimensional flow.^{10} New solution algorithms have also been proposed more recently.^{6,7}

Another strategy, independently proposed by Keylock^{9} and Liu *et al.*,^{8} connects the triple decomposition to the Schur form of the velocity gradient tensor. The Schur form of a general (possibly complex) 3 × 3 matrix *A* is given by

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

where $\sigma 1,\sigma 2,\sigma 3$ are the singular values of *A*. If *A* is a normal matrix then *R* = 0. Therefore, we interpret $UDU*$ as the normal part of the matrix *A* and $URU*$ as 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 form (11) can be used to construct a triple decomposition, analogous to Eq. (5). Complex eigenvalues of the velocity gradient tensor correspond to a rotating flow. Complex eigenvalues of a real matrix always exist as a pair of complex conjugates $(\lambda ,\lambda *)=(a+ib,a\u2212ib)$. To avoid complex arithmetic, we can use the real Schur form,^{11} in which each pair of complex eigenvalues is represented by the real 2 × 2 block

on the diagonal of the matrix *D*. Consequently, *D* is then quasi-diagonal rather than diagonal. The real Schur form is said to be standardized^{12} if

in which case

Note that since we have complex eigenvalues, we know that $\beta \gamma <0$. If we let *λ*_{1} be the real eigenvalue of the velocity gradient tensor, positioned in the top left of the resulting matrix, and we number the complex eigenvalues $\lambda 2=\lambda $ and $\lambda 3=\lambda *$, we have the standardized real Schur form of the velocity gradient tensor

with $(x\xaf1,x\xaf2,x\xaf3)$ being the coordinates in the frame of reference defined by the orthogonal (real) Schur transformation matrix *U* and $(u\xaf1,u\xaf2,u\xaf3)$ being the corresponding velocities.

Liu *et al.*^{8} defined the existence of rotation around a coordinate axis $x\xafi$ by the condition

with $i\u2260j\u2260k$. The fact that

implies that rotation in this sense is restricted to the coordinate axis $x\xaf1$.^{8} This standardized real Schur form with *λ*_{1} top left defines a frame of reference in which the single axis of rotation is the $x\xaf1$ axis, about which rotation is characterized by the complex conjugate eigenvalues $a\xb1\beta \gamma $. If $\beta =\u2212\gamma $, the rotation is circular, otherwise it is elliptic due to shear. Furthermore, Eq. (18) implies that the other possible shear directions are across the $x\xaf1x\xaf2$ and $x\xaf1x\xaf3$ planes.

To extract the shear tensor $(\u2207u)SH$ from the real Schur form, it is natural to use the same approach as in Eq. (8). That is, to subtract the non-skew-symmetric part of the off diagonal components of the standardized real Schur form,

The corresponding matrices for strain and rotation are, respectively, the symmetric part,

and the skew-symmetric part,

The main purpose of this note is to clarify the relation between the standardized Schur transformation *U* and the basic reference frame defined by *Q*. We will show that the Schur transformation *U* is a solution to the optimization problem defined by Kolář and, therefore, defines a basic reference frame. A basic reference frame *Q* maximizes the quantity $\Delta (Q)$, defined in Eq. (9). If there is no rotation, corresponding to three real eigenvalues of the velocity gradient tensor, the Schur transformation leads to a triangular transformed velocity gradient tensor $\u2207u\xaf=QT\u2207uQ$, which maximizes $\Delta (Q)$ since all subdiagonal components are zero.

In flow with non-zero rotation, a real Schur form with this structure defines a single axis of rotation $x\xaf1$ that maximizes two out of three of the terms in $\Delta (Q)$, and which reduces the maximization problem to the 2D plane $x\xaf2x\xaf3$. However, the solution to this maximization problem is trivial, which leads us to the first main conclusion of this paper:

**Theorem 1.** *A standardized real Schur form is a solution to the optimization problem posed by Kolář*.

*Proof.* We have shown that two out of three terms in Eq. (9) are optimal, and it remains to prove that the remaining term is optimal. In the plane $x\xaf2x\xaf3$, we can define the maximization problem as finding the critical point of the subdiagonal component in the matrix (16) under a Givens rotation

where $\theta \u2208[0,2\pi )$ is the counterclockwise rotation angle in the $x\xaf2x\xaf3$ plane. The Givens transformation is expressed as

with $[\u2207u\xaf]x\xaf2x\xaf3$ being the 2 × 2 block corresponding to the $x\xaf2x\xaf3$ plane.

The optimization problem can now be stated as to find the critical point of the sub-diagonal component of the Givens transformation as a function of the rotation angle *θ*,

with derivative

It follows that for a critical point $f\u2032(\theta )=0$, the rotation angle has four possible solutions:

This corresponds to the optimal value of the subdiagonal element being either $f(\theta )=\gamma $ or $f(\theta )=\u2212\beta $. Thus, the standardized Schur form (16) is already optimal up to a right angle rotation, which in the block $[\u2207u\xaf]x\xaf2x\xaf3$ only results in switching the positions and signs of *γ* and *β*. Such a matrix rotation does not change the value of $\Delta (Q)$ given in Eq. (9). Hence, it is clear that the standardized Schur transformation *U* is also a basic reference frame *Q* as defined by the optimization problem stated by Kolář.□

The link between the Schur form and the triple decomposition of a velocity gradient tensor has been established before.^{8,9} Proving that the standardized real Schur form is indeed a solution to the original optimization problem posed by Kolář solidifies the theory behind the triple decomposition of the velocity gradient tensor, while also offering improvements to algorithms for computing the triple decomposition; a real Schur form on non-standardized form requires the Givens rotation in Eq. (23) around the single axis of rotation of the flow to maximize the non-zero term of $\Delta (Q)$ in Eq. (9), see, e.g., the algorithm for computing the Rortex by Liu *et al.,*^{8} whereas the standardized real Schur form is already optimal. The standardized real Schur form can be computed by using state of the art linear algebra routines,^{12} for example, the function *schur* in MATLAB or *DGEES* in LAPACK.

While standardizing the real Schur form makes the computation of the triple decomposition more straightforward, the sorting of the eigenvalues on the quasi-diagonal of the real Schur form influences the outcome of the triple decomposition. In the following part of the paper, we will show how different sorting of eigenvalues corresponds to different alignments of the coordinate system with the flow. Understanding this difference is important for explaining why different triple decomposition algorithms may end up with different results. Together with the standardization of the real Schur form, the following discussion enables us to arrive at a simple and reliable algorithm for computing the triple decomposition.

Above we assume that the real Schur form results in a quasi-triangular matrix *T* with the real eigenvalue as the top left element, and that the 2 × 2 block is placed in the bottom right corner. Furthermore, we state that this aligns the $x\xaf1$ axis of the coordinate system with the rotation axis of the flow. However, this is not obvious from the conditions of Liu *et al.*^{8} Since $\u2202u\xaf2/\u2202x\xaf1=\u2202u\xaf3/\u2202x\xaf1=0$, we have no rotation around the $x\xaf2$ or $x\xaf3$ axes, according to the condition (17). This is not the same as the $x\xaf1$ axis being aligned with the rotation axis.

In 1776, Euler stated his rotation theorem as follows: “When a sphere is moved around its center it is always possible to find a diameter whose direction in the displaced position is the same as in the initial position.”^{13} In other words, we may define the rotation axis of a rigid body rotation transformation as a line on which all points are mapped onto the same line under the transformation. Specifically, for a rigid body rotation transformation *R* and a (real) rotation axis vector *v*,

i.e., *v* is a real eigenvector of *R* and *λ* is the corresponding real eigenvalue.

In the real Schur form above, if the $x\xaf1$ axis is aligned with the rotation axis, it thereby has to be an eigenvector of $\u2207u\xaf$. To check this, let again $\u2207u$ be a velocity gradient tensor, and *Q* and $\u2207u\xaf$ be given by the real Schur form, i.e.,

where $\u2207u\xaf$ has a structure defined in Eq. (16), i.e., with $\u2202u\xaf2/\u2202x\xaf1=\u2202u\xaf3/\u2202x\xaf1=0$. Let *e*_{1} be the standard basis vector $e1=(1,0,0)T$, and

Multiplication by *Q ^{T}* from the left gives

and thereby,

Thus, *v _{R}* is an eigenvector corresponding to the real eigenvalue

*λ*

_{1}for $\u2207u$, and $e1=QTvR$ is an eigenvector of the real eigenvalue for $\u2207u\xaf$. Hence, in accordance with Euler's rotation theorem, the coordinate axis $x\xaf1$ is aligned with the rotation axis.

If instead the order of the eigenvalues on the quasi-diagonal of *T* is different, i.e., the 2 × 2 block is placed top left, and the real eigenvalue is located in the bottom right of *T*, a different structure is obtained. Let

with $\u2207u\u0302$ also given by a real Schur form

In this case, in accordance with the condition in Eq. (17), there is neither rotation around the $x\u03021$ axis nor the $x\u03022$ axis, since

It remains to check whether the rotation axis is aligned with the $x\u03023$ axis. Let

with $e3=(0,0,1)T$. Applying the same procedure as above gives

Thereby, we conclude that an eigenvector corresponding to *λ*_{1} for $\u2207u$ is

and a real eigenvector of $\u2207u\u0302$ cannot coincide with *e*_{3} unless $\zeta \u0302=\epsilon \u0302=0$, i.e., when there is no shear outside the rotation plane $x\u03021x\u03022$. Consequently, a real Schur form with the eigenvalues ordered this way does not align any of the axes of the coordinate system with the rotation axis. Instead, $\u2202u\u03023/\u2202x\u03021=\u2202u\u03023/\u2202x\u03022=0$ implies that there is no shear in the $x\u03023$ direction, so the only possible shear direction is across the $x\u03021x\u03022$ plane. The rotation around the $x\u03021$ and $x\u03022$ axes is also zero, meaning that this Schur form aligns the $x\u03021x\u03022$ plane with the rotation plane of the flow, rather than aligning the $x\u03023$ axis with the rotation axis. In the original Schur form $\u2207u\xaf$ with the real eigenvalue top left, there is no such alignment of any coordinate plane with the rotation plane.

We define the rotation plane of a matrix analogously to the definition of the rotation axis: all points in the rotation plane are mapped back onto the same plane under the transformation. For $\u2207u\u0302$ given in Eq. (32), and a point $(\xi \u0302,\nu \u0302,0)T$ in the $x\u03021x\u03022$ plane, we can easily confirm that this condition is fulfilled

In the case of $\u2207u\xaf$, for which the $x\xaf1$ axis is aligned with the rotation axis, there is no such alignment of the $x\xaf2x\xaf3$ plane with the rotation plane

Going back to the definition of the BRF, it is clear that both *Q* and $Q\u0302$ fulfill the requirement to maximize Eq. (9). By this requirement alone, they could both be applicable for calculating a triple decomposition, but since they correspond to different alignments of the coordinate system with the flow, they may lead to different outcomes. Specifically, since the rigid body rotation has different interpretations in the two cases, its magnitude in the Frobenius norm may differ depending on which Schur form is used.

This is illustrated in Figs. 1(a) and 1(b), which show snapshots of the magnitude of rigid body rotation computed from the two different forms $\u2207u\xaf$ and $\u2207u\u0302$ in a simulation of the flow of blood into the left ventricle of a human heart. The data are from a previous study conducted by our research group.^{5} Figure 1(c) shows the difference between the two rotation fields. The two fields appear quite similar but with significant differences in some areas. The difference is quantified in Table I, where the mean values of the magnitudes are shown, as well as the mean of the absolute value of the difference. The average difference amounts to $17.71%$ of the average rigid body rotation magnitude for the Schur form $\u2207u\xaf$ with the real eigenvalue top left, indicating a clear statistical difference between the two cases.

Mean $(||(\u2207u\xaf)RR||F)$ | 16.71 |

Mean $(||(\u2207u\u0302)RR||F)$ | 16.79 |

Mean $(|||(\u2207u\xaf)RR||F\u2212||(\u2207u\u0302)RR||F|)$ | 2.96 |

Mean $(||(\u2207u\xaf)RR||F)$ | 16.71 |

Mean $(||(\u2207u\u0302)RR||F)$ | 16.79 |

Mean $(|||(\u2207u\xaf)RR||F\u2212||(\u2207u\u0302)RR||F|)$ | 2.96 |

The difference only appears for the rigid body rotation. For both strain and shear, both Schur forms lead to identical magnitudes of the respective measure in the Frobenius norm.

**Remark I:** While the two real Schur forms described above each have a reasonable geometric interpretation, the analogous interpretation of the BRF resulting from the complex Schur form is less clear. As noted earlier, the complex Schur form makes *T* triangular rather than quasi-triangular, with the eigenvalues on the diagonal, and with the strictly upper triangular part [*R* in Eq. (11)] representing shear. Strain in this decomposition, again being the real part of the eigenvalues, is clearly the same as when using the real Schur forms. The non-normal part of the matrix and thereby the shear magnitude in the Frobenius norm is also preserved. Since the normal part of the matrix is concentrated to the diagonal however, we end up with yet another possible measure of rotation. The Schur form preserves the eigenvalues of the velocity gradient, and consequently, the rotation part in the complex Schur form is

where *λ _{c}* is one of the complex eigenvalues.

**Remark II:** The rotation from the first real Schur form structure ($\u2207u\xaf$, with the real eigenvalue top left) is also equal to the rotation of the second real Schur form structure applied to $\u2207uT$, and vice versa. Liu *et al.* implicitly used this to align the rotation axis with the $x\xaf3$ axis instead of the $x\xaf1$ axis by transposing once before computing the Schur form, and transposing again after the computation.^{8}

The original minimization problem proposed by Kolář is solved by all of the above structures of the Schur form. In fact, there are even more potential BRF structures from Kolář's proposed optimization approach to solve the problem. The real Schur form is restricted to the $x\xaf1$ axis as the rotation axis, or the $x\u03021x\u03022$ plane as the rotation plane, but Kolář's algorithm has no such restrictions. Strain is always preserved regardless of which BRF is used, and so is shear, but rotation is not.

To conclude, based on the discussions of the standardized Schur form and the different coordinate system alignments with the flow, our suggestion for computing the triple decomposition is as follows: (i) a standardized Schur form should always be used to make the Givens rotation around the flow rotation axis after the Schur computation obsolete; (ii) based on Euler's rotation theorem and the subsequent definition of a rotation axis, we advocate sorting the eigenvalues on the quasi-diagonal of the Schur form with the real eigenvalue top left, thereby aligning the coordinate system with the rotation axis. Based on this, our proposed algorithm for computing the triple decomposition is as follows:

Compute a standardized real Schur form, sorted with a real eigenvalue top left: $\u2207u=QT\u2207u\xaf\u2009Q$

Subtract shear, strain, and rotation from each other by Eqs. (19)–(21)

The algorithm is implemented below in Python, with the magnitudes of shear, strain, and rotation in the Frobenius norm computed directly:

import scipy.linalg.lapack as la

import numpy as np

def triple_decomposition(grad_u):

def dselect(arg1,arg2): return (arg2==0)

T=la.dgees(dselect,grad_u,sort_t = 1)[0]

sh=np.linalg.norm([T[0,1],T[0,2],T[1,2]+T[2,1]])

el=np.linalg.norm(np.diag(T))

rr=np.sqrt(2*min(abs(T[1,2]),abs(T[2,1]))**2)

return sh, el, rr

## ACKNOWLEDGMENTS

This work was supported by the Swedish Research Council (Grant No. 2018-04854). Access to computing infrastructure was provided by the Swedish National Infrastructure for Computing (Grant No. SNIC 2022/5-23). The authors would like to thank Dr. Václav Kolář for a helpful discussion regarding the optimization problem definition of the triple decomposition.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Joel Kronborg:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Johan Hoffman:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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