Triple Decomposition of Velocity Gradient Tensor in Compressible Turbulence

: The decomposition of the local motion of a ﬂuid into straining, shearing, and rigid-body rotation is examined in this work for a compressible isotropic turbulence by means of direct numerical simulations. The triple decomposition is closely associated with a basic reference frame (BRF), in which the extraction of the biasing effect of shear is maximized. In this study, a new computational and inexpensive procedure is proposed to identify the BRF for a three-dimensional ﬂow ﬁeld. In addition, the inﬂuence of compressibility effects on some statistical properties of the turbulent structures is addressed. The direct numerical simulations are carried out with a Reynolds number that is based on the Taylor micro-scale of Re λ = 100 for various turbulent Mach numbers that range from Ma t = 0.12 to Ma t = 0.89. The DNS database is generated with an improved seventh-order accurate weighted essentially non-oscillatory scheme to discretize the non-linear advective terms, and an eighth-order accurate centered ﬁnite difference scheme is retained for the diffusive terms. One of the major ﬁndings of this analysis is that regions featuring strong rigid-body rotations or straining motions are highly spatially intermittent, while most of the ﬂow regions exhibit moderately strong shearing motions in the absence of rigid-body rotations and straining motions. The majority of compressibility effects can be estimated if the scaling laws in the case of compressible turbulence are rescaled by only considering the solenoidal contributions.


Introduction
Turbulence in the context of fluid dynamics often refers to the chaotic nature of fluid motion in terms of velocity and pressure. This kind of definition suggests deterministic chaotic and turbulent characteristics [1], which must be distinguished from random behaviour. However, also for deterministic turbulent flows, some statistical laws apply in the same manner as for linear systems [2], whereas their statistics may widely differ; they usually are of Lévy-type and show long tail distributions that are produced by rare large-scale fluctuations and clustering phenomena of anomalous eddy diffusion. However, different than for random systems in phase space, for example, turbulent ones may exhibit highly ordered structures, called strange attractors [1], which relate to nonlinear, nonloca,l and fractional features of turbulence [3]. Additionally, nonequilibrium [4], nonextensive [5], and universal behaviour [6] characterize turbulent flows. In this article, we focus on universality, and, by this, we mean the independence of statistical properties of small-scale eddies from the overall nature of the fluid (forcing of the fluid, geometry, and size of the basin, etc.). One of the first illustrations of this universal character dates back to 1941, when the Kolmogorov theory emerged. The latter states that, regardless of a fluid's properties, the turbulent flow exhibits a scale-invariant structure at the inertial region. Based on this finding, Kolmogorov was able to characterize this zone as a power law of the energy spectrum with an exponent that is equal to −5/3 [7], which for flows of Reynolds numbers smaller than infinity, is modified by an intermittency correction [8,9].
Subsequently, many areas of turbulence research were driven by the fundamental challenge of identifying generalizable turbulence properties and turbulent flow similarities with the aim of investigating and defining turbulent flow universal patterns [10]. In particular, because fluid mixing is mainly driven by the smallest diffusive scales, many studies have been devoted to the assessment of the fundamental and intrinsic topological aspects of these flow scales. In this regard, the pure kinematic approach provides a general and useful tool to describe the small scales of fluid motions. One of these quantities, which is often of significant interest, is the velocity gradient tensor (VGT), ∇U ≡ A, which can be split into the sum of the strain-rate tensor, S, and the rotation-rate, W tensor, i.e., (see also Ref. [11]) A = S + W.
For many decades, this decomposition has been an essential tool for the turbulence community [12] due to the simplicity of its formulation and the physical insight that it provides. It has also led to the construction of various subgrid-scale models [13,14] and the development of flow visualization tools that are based on simulation realizations [15]. However, recent studies have shown the importance of considering the contribution of shear in both the identification of vortex structure regions and strain rate [15]. Indeed, vorticity cannot distinguish between pure shearing motion and vortical motion, and the strain rate cannot distinguish between straining motion and shearing motion. By considering two-dimensional DNS of a subsonic mixing layer behind a flat plate, [16] showed the importance of shear regions in the identification of vortex structures. Most of the discussions on the structural aspects of turbulence adopt Eulerian vortex detection methods, like the Q-criterion, the ∆-criterion, and its closely related swirling strength criterion (λ ci -criterion [17]), or the λ 2 -criterion [18]. In all of these criteria, a scalar field, which is derived from the decomposition (1) and/or eigenvalues of the VGT, is used as a "marker" to indicate whether a given point in the flow field belongs to a vortex or not [19]. Therefore, vortices are identified as the connected regions that are mapped by such scalar fields. Recently, [20] formulated a Liutex-based VGT (previously called Rortex) decomposition for locally fluid-rotational points (the VGT has complex eigenvalues) in a turbulent flow field. This method [20] of separating the rigid-body-rotation "Liutex" from shear in vorticity is more computationally viable and it has been employed in some studies [17] for the investigation of coherent vortex structures in turbulent flows. In another study, [21] presented a decomposition of the VGT into normal and non-normal tensors, such that the non-normal counterpart represents the local effects of shear. Aside from these methods, [15] originally introduced a robust triple decomposition method, from which a residual vorticity could be found to represent the pure rigid-body rotation of a fluid element. The cornerstone of this method is to decompose A as where N , H and R denote normal-straining, pure-shearing, and rigid-body-rotation components of the VGT, respectively. According to [15], the three components in Equation (2) can be determined while using a suitable reference frame, called the basic reference frame (BRF). The BRF refers to a special local coordinate system in which the effect of shear is maximal, so that the decomposition is unique.
Most of the studies that targeted this triple decomposition method (TDM) were carried out in two-dimensional configurations [22,23]. Under such flows, the identification of the basic reference can be analytically obtained and the extraction of the three motions is straightforward. Conversely, the analytical identification of the BRF is not possible in the general case where turbulence is comprised of three-dimensional fluid motions, and the extraction of vortices and internal shear layers are of considerable interest. In the general case, one can identify the BRF from a finite set of discrete frame representations by performing a series of discrete rotations of the reference frame [15]. In this context, [24] proposed a numerical procedure and then successfully applied it to the framework of three-dimensional incompressible homogeneous isotropic turbulence at Taylor Reynolds numbers, Re λ = 27 and 140. One of the major findings of [24] is that some of the strong shear regions, identified using the TDM and visualized as sheet-like structures, could not have been detected using the classical double decomposition. Indeed, TDM analysis characterizes these regions as high intensity shearing components, whereas, from the double decomposition, they appear as regions where the effects of the strain-rate tensor, S, and the rotation-rate tensor, W, are comparable. With this in mind, and given that only incompressible flows have been considered so far, we must question the accuracy of the TDM for more complex flows that are characterized by high turbulent Mach numbers, Ma t . For instance, unlike the incompressible case, at high Ma t , both the strong nonlinear coupling between the velocity and thermal fields and the occurrence of localized shocklets in the flow field lead to an additional and rapid conversion of the kinetic energy dissipation mechanism. The analysis of the flow topology in compressible flows has been investigated before; see, for example, the works of [25][26][27][28]. A general conclusion from these studies is that the appearances of the intense events tend to increase with Ma t . These studies also indicate that, as Ma t increases, strong compressions are more likely to occur than equally strong expansions.
The present work aims at extending the previous studies to the compressible high Ma t regime. Furthermore, we suggest a new numerical approach to determine the BRF with the aim of achieving high accuracy at a significantly reduced computational cost. We then apply the TDM and the new numerical procedure in order to define the BRF to a new DNS database of isotropic compressible turbulence with a turbulent Mach number ranging from Ma t = 0.12 (i.e., almost incompressible regime) to Ma t = 0.89 (i.e., highly compressible regime).
The remainder of this paper is organized, as follows. Section 2 gives a brief overview of the governing equations with the methods and algorithms that were used in their numerical resolution. Section 3 provides the general statistics of the generated DNS database. The invariants of the VGT are briefly recalled in Section 4. Next, the methodology of A triple decomposition, as well as the procedure of determining the basic reference frame in Section 5, are presented. By exploring its dependency on the turbulent Mach number, the VGT composition of a turbulent flow field is examined in detail in Section 6. Finally, Section 7 draws the key findings and conclusion of this analysis.

Governing Equations and Numerical Methods
In this study, the nondimensional form of the three-dimensional compressible Navier-Stokes equations, as described in Samtaney et al. [25], is considered ∂ρ ∂t

∂E ∂t
where Re, Ma, and Pr represent the Reynolds number, the Mach number, and the Prandtl number, respectively. The parameter γ is the ratio of the specific heat at constant pressure and the specific heat at constant volume. In the present simulations, the values of γ and Pr are set to 1.4 and 0.7, respectively. Furthermore, α is defined as α ≡ Pr Re(γ − 1)Ma 2 . The primary variables are the density, ρ, the velocity components, U i , the temperature, T , and the pressure, P. The viscous stress, σ ij , and the total energy per unit volume, E , are defined by where θ = ∂U k /∂x k refers to the velocity divergence or dilatation. The variable θ quantifies the local rate of rarefaction (θ > 0) or compression (θ < 0). Sutherland's law is employed to specify the temperature-dependent and dimensionless dynamical viscosity µ and thermal conductivity κ, as follows: and where c 1 = 1.4042 and T c = 0.40417.

Numerical Methods
The numerical simulation of compressible flows requires the use of highly accurate numerical schemes that are capable of precisely capturing shock waves. Consequently, the viscous fluxes in Equations (3)-(5) are obtained through the application of an eighthorder accurate centered finite difference scheme, while the inviscid fluxes are approximated using a seventh-order accurate weighted essentially non-oscillatory (WENO-Z) scheme of Don and Borges [29].
It is worth recognizing the difficulties in conducting DNS of compressible flows featuring shock waves due to the conflicting requirements of: (i) resolving the broadband turbulence and (ii) capturing shocks. This makes the obtained solution highly dependent of the considered numerical procedures. In practice, the present numerical solver uses the ses the optimal seventh-order accurate flux reconstruction, and the application of the WENO-Z scheme is conditioned to a smoothness criterion that involves the local values of the normalized spatial variations of both pressure and density [30]. This method is equivalent to the superposition of the four candidate stencils of the WENO-Z scheme when each stencil is evaluated with its optimal weight. As a result, the introduction of such a hybridization makes the used method minimally dissipative, thanks to the use of a central scheme in the majority of the domain, and it yields a reasonable predicted spectrum up to rather high wavenumbers.
Regarding the temporal advancement, a semi-discrete system of ordinary differential equations stems from the discretization of the spatial derivatives where c = [ρ, ρU 1 , ρU 2 , ρU 3 , ρE ] is the vector of the conservative variables and Res the vector of the residuals. The time advancement of Equation (10) is performed while using a third-order accurate total variation diminishing (TVD) low-storage Runge-Kutta scheme [31].

Dns Database
The decaying isotropic turbulence is an important model that is used to investigate the free compressible turbulence. In this context, ensuring that the turbulence is fully developed is important. Prior to the generation of the DNS database, the numerical procedure is first validated by comparing its result against the DNS data of [34]. The validation test case consists of a decaying isotropic turbulence evolving in a three-dimensional periodic box with side length 2π that is uniformly discretized on 128 points per each coordinate direction. The velocity field is initialized by an isotropic state that is prescribed by the energy spectrum where k is the integer non-dimensional wave number and k 0 is the most energetic wave number. A finite density, temperature, and dilatation fluctuations are also generated, as described by [35]. The calculation is performed at a Reynolds number based on the longitudinal Taylor micro-scale Re λ = ρλU /µ = 50, where λ = U 2 1 1/2 / (∂U 1 /∂x) 2 1/2 and U = U i U i /3 1/2 , and at a turbulent Mach number Ma t = U /a = 0.52, where a is the speed of sound. The initial flow-field is allowed to decay temporally for four dimensionless time units, 4τ t , where the time unit, herein, is the eddy turnover time that is defined as τ t = λ/U . This choice allows for the energy spectrum E(k) to develop an inertial range that decays as k −5/3 . In Figure 1, we depict the energy spectra in wave number space at t/τ t = 6.5. These results show a very good agreement with the data of [34]. Moreover, the energy spectrum is marginally consistent with a −5/3 scaling law in a smaller frequency range due to the limit of the spatial resolution in the present DNS (Note that a recently developed approach, based on the extraction of the rigid rotation part from the fluid motion (called Liutex-based vortex), allows following the −5/3 almost exactly over the whole range of frequencies [10]). Following this validation, the DNS database featuring six direct numerical simulations of temporal turbulence decay is conducted on a 512 3 point grid. Table 1 summarizes the one-point statistics of the six simulated flows, obtained by freezing the temporal turbulence decay of periodic boxes, after approximately three eddy turnover times, which ensures that the turbulence is deemed to be developed and realistic. The turbulent Mach number, Ma t , which is varied between 0.12 and 0.89, is the main parameter of interest in this parametric study. The Taylor micro-scale Reynolds number, Re λ , is considered to be common to all the cases and set almost to 100. The wave number k 0 = 4 is selected to ensure that the fully developed homogeneous isotropic turbulence field has a broad range of length scales [36]. The Kolmogorov length scale, which is defined from the kinetic energy dissipation rate per unit volume (η = µ/Re 3 / ε/ρ 1/4 ), is resolved on the retained grid discretization.
Indeed, the resolution parameter ∆x/η is in the range 1.15 ≤ η/∆x ≤ 1.81, where ∆x is the uniform grid spacing in each coordinate direction. This can be also verified in terms of the maximum resolved wave number in the grid (the half of the number of grid points in each direction). In the present DNS database, k max η finds its values in the range 2.22 < k max η < 2.73. This range ensures that the small scales are well resolved in all of the database cases [37]. In order to measure the onset of a realistic turbulence, the mean values of the velocity derivative skewness Sk 3 and flatness F l 3 are used herein as criteria. These are defined, respectively, as As it can be noticed from Table 1, the magnitudes of both Sk 3 and F l 3 increase with increasing Ma t . This is a well-known feature that is due to the onset of shocklets in compressible turbulence [38]. It is noteworthy to point out that, despite the comparability between the shock thickness and the grid length, as well as the fact that the shock thickness is not directly resolved by the WENO-Z scheme, the total amount of dissipation across a shock is independent of numerical viscosity [39]. Indeed, as shown by [40], as long as the shocklet thickness is small, the total amount of dissipation across the shocks is independent of viscosity, and it depends only on the jump conditions across the shocks, which are preserved by the WENO-Z scheme.

The Velocity Gradient Tensor and the Invariants of Its Characteristic Equation
The velocity-gradient tensor, A, is a 3 × 3 matrix that is given by A ij := ∂U i /∂x j . The characteristic equation for A is where I is the identity matrix and the first three invariants P, Q, and R in Equation (15) can be calculated as where S and W are, respectively, the strain-and rotation-rate tensor components of A, defined as As shown by [41], the expression P, Q and R is of a topological significance, because the P − Q − R space is divided into several regions, and each region represents a particular flow topology. It is more convenient to analyze the flow topology in the Q − R plane due to the spatial complexity of different zones present in the P − Q − R space. To proceed, we reformulate the problem in terms of the anisotropic part of the deformation rate tensor, i.e., A = A − θI /3, where θ = ∂U i /∂x i denotes the divergence of the velocity. Therefore, we observe that, as in incompressible turbulence, by construction P = 0. The use of the anisotropic tensor aims at (i) comparing the statistical properties of compressible turbulent structures and their dynamics with the incompressible case for which P ≡ 0 and (ii) facilitating the discussions of local topological flow structures, since it reduces the space of comparison to the Q − R plane [26,27]. All of the quantities computed from the anisotropic tensor are indicated with the superscript . In the new formulation, the second, Q , and third invariants, R , are given by These two invariants are of great importance from a topological point of view, because the sign of the discriminant function, ∆ A , which characterizes a different streamline flow pattern, depends on their magnitude We can distinguish four non-degenerate topology types in the Q -R space (see Figure 2) • ∆ A > 0, R > 0: compressing towards an unstable focus region (UFC), • ∆ A > 0, R < 0: stretching away from a stable focus region (SFS), • ∆ A < 0, R > 0: two saddles with an unstable node (UNSS), and • ∆ A < 0, R < 0: two saddles with a stable node (SNSS).
For a detailed description of each type, we refer the reader to Ooi et al. [42]. The presence of pure shearing motions and the actual swirling motion of a vortex in both strain-rate and vorticity often limits our understanding of the local streamline topology and the fundamental phenomena in turbulence in general. The velocity gradient dynamics and small-scale behavior of turbulenc are conventionally based on the double decomposition of the VGT in the S and W tensors. To further ease the interpretation of the classical VGT dynamics, the next section suggests introducing (i) an originally additive decomposition that decomposes an arbitrary instantaneous field of the motion into three elementary types of motions A and (ii) a new numerical procedure to efficiently compute the components of this decomposition.

Triple Decomposition of Velocity Gradient Tensor
In this section, we present a new numerical procedure to compute the triple decomposition of motion (TDM) that was introduced by [15]. The main idea behind this method is to subtract the effective pure shearing tensor from the VGT, and then decompose the residual vorticity (which directly and accurately measures the pure rigid-body rotation of a given fluid element) into symmetric and anti-symmetric parts. Therefore, the 3 × 3 velocity gradient tensor, A, whose components are given by A ij := ∂U i /∂x j , is decomposed into three elementary motions: irrotational straining/elongation tensor, N , pure-shear tensor, H, and rigid-body-rotation-rate tensor, R, see Equation (2). The pure-shear stress tensor, H, is a purely asymmetric tensor that satisfies The residual tensor, defined as A res = A − H, is further decomposed into two components that are associated with the irrotational straining motion (elongation) N and rigid-body rotation R. These transformations are illustrated with some elementary examples in Figure 3. We note that the triple decomposition (2) entails a considerable level of effort and the technique that is outlined here is presented in greater detail in [24]. Figure 3. Two-dimensional sketch of the fluid motion deformation due to (a) normal-strain-rate tensor, (b) shear tensor, and (c) rigid-body-rotation tensor. Adapted from [43].
The main difficulty of the triple decomposition (2) lies in the calculation of the pureshear tensor, H, since the other two tensors are computed straightforwardly by subtracting H from A. The decomposition (2) must be applied in an appropriate reference frame that originates from a specific point in the flow field. This reference frame is called the "basic reference frame" (BRF). In the BRF, an effective pure-shear motion is shown "in a clearly visible manner" set by the condition (20), and the effect of extraction of a shear tensor is maximized. In order to fulfill the latter condition, Kolář [15] introduced an interaction scalar Is that is defined by Then, the BRF is determined from the condition where O(3) is the set of 3 × 3 orthogonal matrices, and Is(Q) is the interaction scalar that is calculated from the velocity gradient tensor A, which is computed under a generic unitary orthogonal transformation Q A tilde quantity represents a quantity that was evaluated in the rotated reference frame. The derivation of the orthogonal matrix (rotation matrix whose yaw, pitch, and roll angles are φ, θ, and θ, respectively) Q applied to a Cartesian coordinate x is provided in [15], and it reads (where c and s denote the cos and sin functions, respectively) where φ ∈ [0, π], θ ∈ [0, π], and ψ ∈ [0, π/2] are the three rotation angles. In the procedure that was proposed by Nagata et al. [24] to compute the BRF, these three rotation angles are varied with a uniform step size ∆, yielding a finite number of reference frame rotations: (l∆, m∆, n∆), where (l, m, n) are three integers satisfying 0 < l ≤ π/∆ , 0 < m ≤ π/∆ , and 0 < n ≤ π/2∆ . Therefore, the determination of the accurate location of BRF requires significant computational effort as a small step size ∆ must be used. Our strategy instead consists in the formulation of the BRF determination as an optimization problem of the three bounded variables φ ∈ [0, π], θ ∈ [0, π], and ψ ∈ [0, π/2], and a fixed velocity gradient tensor, A, associated with the point Here, to solve Equation (25), we use the gradient-based optimizer sequential leastsquare quadratic programming (SLSQP), which is an SQP-like algorithm with a Broyden-Fletcher-Goldfarb-Shanno (BFGS) update the formula for approximating the inverse of the Hessian matrix of the objective function. In practice, a free modern FORTRAN SLSQP optimizer version [44] is coupled with our solver to achieve this purpose. To highlight the behavior of the proposed approach, we carry out a relatively small simulation while using 128 3 grid points.
In Figure 4, we show the probability density function (PDF) of the number of iterations and the number of evaluations of the Jacobian of the objective function performed by the optimizer to locate the BRF. We observe a left-skewed tendency and peak at approximately 10 iterations, which allows for a fast and more accurate location of the BRF when compared to the previous discrete approaches [24,45]. When comparing the difference between the present method and the one of Nagata et al. [24], the difference gets smaller as the step size ∆ becomes smaller. Besides, the method (SLSQP) exhibits no sensitivity to the initial guess of (φ, θ, ψ).  This step results in the velocity gradient tensor A in the BRF, and then can, therefore, be decomposed into the following constituent tensors A = A res + H. Here, the residual tensor A res is expressed as where the sign function of a generic real number x gives -1 if x < 0, 0 if x = 0, and 1 if x > 0. Finally, the tensor A res can be further decomposed into symmetric and anti-symmetric The resulting decomposition is finally expressed in the laboratory reference frame as Φ = Q ΦQ for a generic tensor Φ. The pure-shear tensor is divided into its symmetric H S and anti-symmetric H W counterparts, i.e., The symmetric-shear tensor H S , along with normal-strain-rate tensor N , recovers the strain-rate tensor S, while the anti-symmetric-shear tensor H W m along with rigid-bodyrotation tensor N constitutes the rotation-rate or vorticity tensor, i.e., From Equation (29), we notice that the pure-shear motions contribute to both vorticity and strain-rate. Figure 2 illustrates the shape of the local streamline that is associated with the tensors associated with the triple decomposition in the Q − R plane [43]: • ∆ A < 0: Non-rotational geometries (A = N + H ), the tensor A has only real eigenvalues; and, • ∆ A > 0: Rotational geometries (A = N + R + H ), the tensor A has complex eigenvalues, the flow is locally rotational, and the three tensors N , R , and H are in general non-zero.

Results and Discussion
In Figure 5, we show the isocontour lines (of the logarithm) of the second, Q and third, R , invariants for the values of the turbulent Mach number, Ma t , considered in the present database. As expected, all of the joint PDFs exhibit the classical tear-drop shape as in the incompressible case, with a statistical preference in the SFS and UNSS quadrants. Within these quadrants, the statistical points are mostly aligned with the right branch that corresponds to ∆ A = 0. The compression motion tends to reveal a stronger alignment with the right branch, as indicated by a rather sharp joint PDF with the curve ∆ A = 0 [27]. Table 2 provides a quantitative comparison between the conditional data of each case in percentage, distributed over the possible local topologies (the row sum is 100%). As expected, the most frequently occupied regions are the SFS and UNSS regions, which are dominated by positive enstrophy and strain productions. Furthermore, the probability of occurrence of each topology in the almost incompressible case with Ma t = 0.12 is very close to that found in the compressible turbulence. This fact proves that (i) the dilatation has a negligible effect on the occurrence of each topology and (ii) compressibility marginally affects the invariant statistics. The occupancy of the non-rotational region that is characterized by R ≈ 0 is close to the sum of the conditional data of the UNSS and SNSS topologies.   Figure 2 for the velocity gradient tensor, A , expressed as a percentage of realizations. The last column indicates the non-rotational non-degenerate topology, where R is the Frobenius norm of R . For the sake of clarity, from this point onward, we adopt the notations s Φ = 2Φ S ij 2Φ S ij and ω Φ = 2Φ W ij 2Φ W ij to indicate the strain-rate magnitude and vorticity magnitude (or enstrophy), respectively, of a generic tensor Φ with symmetric, Φ S , and anti-symmetric, Φ W , parts. Note that ω H = s H from the construction of H . Therefore, the strain-rate magnitude, ω A , and rotation-rate magnitude, s A , (in the sense of Frobenius norm) are expressed in terms of ω H , ω R , s H , and s N .

Ma
To show the flow structures that are deduced from the triple decomposition of A , in Figure 6 we plot the isosurfaces for the four constituents of vorticity and strain-rate for Ma t = 0.5. A vortex filament is viewed as a vortex tube of variable cross-section for ω R , as shown in Figure 6a. We note that the shearing motion contribution for the vorticity and strain production corresponds to sheet-like structures, which is a marker of high energy dissipation/strain rate regions [46], as shown in Figure 6b. In contrast, the instantaneous snapshot of the s R isosurface shows a highly intermittent turbulence in space.
A possible two-dimensional illustration of the decomposition of vorticity and strain rate magnitudes is presented in Figure 7 for Ma t = 0.12 (an almost incompressible case) and Figure 8 for Ma t = 0.50 (an example of a compressible case). All of the other cases display similar patterns. In the following, the focus is, mainly, put on the qualitative description of the decomposition for the compressible case.
Because the vorticity is a local characteristic of the flow and includes both the solid body rotation and the shearing motion of a fluid element, we can isolate the two effects by comparing Figure 8a,b. Figure 8b indicates that the vortices are frequently located in high shear zones. This implies that the criteria of vortex detection (e.g., the Q-criterion), which simultaneously detects shear layers and vortex regions, sometimes yields results that are difficult to interpret. Indeed, they exclusively locate a region that is to be either dominated by strain or by rotation. Regions of high intermittency are associated with large values of ω R (vortex with a strong-body rotation) and s R (strain with a strong-body rotation), as depicted in Figure 8c-e.
To deepen our investigation of the statistical properties of the triple decomposition, in Figure 9 we report the normalized ω A and s A and their components ω H , ω R , s H , and s N . As in the incompressible case, all of the PDFs exhibit a two-stage evolution of their shapes, where they first increase up to a peak, and then decrease slowly to lower values. The effect of compressibility is slightly significant for Ma t > 0.5. The PDFs of ω A and s A become wider and extend slowly to much higher values. Moreover, the PDF of the strain-rate magnitude has a less wide tail at higher values and it is more intermittent in comparison to the rotation-rate magnitude at smaller and moderate values. This observation, which holds true independently from the compressibility intensity considered within our investigation, is a well-known characteristic of spatial intermittency [47,48]. The same figure shows that the pure-shear magnitudes ω H and s H exhibit a wider tail and a flatness, and they are more intermittent than the other components. As a result, one can conclude that the pure-shear motions have a greater contribution in the heavy-tailed PDF of the VGT magnitude. In particular, the large probability values of ω A / ω A ≈ 0.5 are mostly associated with the shear, rather than the rigid-body rotation. The rigid-body rotation strain magnitude, s N , first reaches a peak at s N / s A ≈ 0.25, and it decreases while the pure-shear contribution increases to reach its maximum value at s H / s A ≈ s A / s A ≈ 0.75, as the strain-rate magnitude. The compressibility seems to slightly affect the fat tail of the PDFs for small values of ω and s.  Figure 10 shows the PDFs of the ratio H /ω H , s H /ω R , s R /ω H , and s R /ω R for various Ma t . These distributions suggest that vortices are most probably located in regions where the shear-strain rate is equal to, or larger than, the maximum vorticity of the vortices. However, the most probable case is the one of low shear, which confirms that the vorticity is not always a good indicator for measuring the strength of the local swirling. Moreover, as the mode of the PDFs of the four quantities portrayed in Figure 10 is likely associated with the intensity of the compressibility effects. The overall mode of the distribution increases by increasing Ma t .
We further address the coupling between the vorticity and strain by looking at the two-dimensional joint PDF of the vorticity magnitude (enstrophy) and its two components against the total strain rate. Figure 11a,d indicate that strong vorticity coexists with strong strain, even though the correlation between ω A and s A seems to be rather weak, and even weaker for lower Ma t [37]. The shearing-motion contribution in the production of enstrophy is likely more correlated with the strain rate than the rigid-body rotation. In order to quantitatively assess this feature, we compute the Pearson's bi-variate coefficients of correlation between each component. This coefficient detects any existence of linear relations between pairs of variables and it reads r(x, y) = cov(x, y) var(x)var(y) , where cov(x, y) ≡ s xy stands for the covariance of the variables x, y, whereas var(x) ≡ s x and var(y) ≡ s y are their respective standard deviations.  We notice that the high correlation between ω A and s A is mostly due to the shearingmotion of the vorticity, as depicted in Table 3, and the correlation between strain rate and component of vorticity is almost constant, regardless of the intensity of the compressibility.

Conclusions
The triple decomposition of the VGT appears to be a promising candidate for the existing vortex-identification methods. It provides more insights into the normal-strain and pure rotation on important small-scale features than those that were obtained by the more simple decomposition of A into symmetric, S, and anti-symmetric, W, parts. However, this decomposition suffers from the difficult task of identifying the basic reference frame (BRF) because it is based on the natural decomposition of fluid motion into straining, shearing, and rigid-body rotation. The BRF is determined by rotating the local reference frame that is associated with a fluid element and finding the rotation angles maximizing the interaction scalar coefficient Is. The numerical results show that the usual classification of the local streamline topologies can enhance our understanding by analyzing the VGT in terms of the normal strain, the rigid-body rotation, and the pure-shear tensors. The rigid-body rotation negligibly contributes, on average, to both the vorticity and strainrate magnitudes for all of the tested turbulent Mach numbers of the database, while shearing motion is the main cause of the strong intermittency that is exhibited by the vorticity magnitude. Indeed, shear-magnitude increases steeply, while rigid-body-rotation magnitude decreases at extreme vorticity magnitude values. Therefore, it is reasonable to infer that shear dominates in regions of high intermittency. The left-skewed distribution of the vorticity and strain rate magnitudes are most probably attributed to the rigid-body rotation distribution, while the high values are most probably associated with the pureshear motions. The strain-rate magnitude is found to be more correlated with the shear tensor than with the rigid-body rotation tensor. The study of the anisotropic tensor leads to the conclusion that compressibility has a marginal effect on the statistical quantities that were investigated in this work. More turbulence characteristics need to viewed with this new physical intuition introduced by the triple decomposition of the velocity gradient tensor.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.