1. Introduction
The components of the velocity gradient tensor (VGT)
$\unicode{x1D63C}$
are defined as
$A_{\textit{ij}}=\partial u_{i}/\partial x_{j}$
, where
$u_i$
is the fluctuating velocity in the
$i$
th direction, and it provides fundamental information about the properties of small-scale structures (Sreenivasan & Antonia Reference Sreenivasan and Antonia1997; Meneveau Reference Meneveau2011; Davidson Reference Davidson2015). The tensor appears in the Taylor expansion of the velocity field around an arbitrary point in a flow, and therefore encapsulates local information around the point. The tensor is sensitive to the non-Gaussian characteristics of turbulent fields due to internal intermittency (Tsinober Reference Tsinober2009). Furthermore, utilising critical point theory, local topological properties can be analysed using the VGT (Perry & Chong Reference Perry and Chong1987; Ooi et al. Reference Ooi, Martin, Soria and Chong1999).
Extensive research on the VGT and its invariants has been conducted in a number of canonical flows, such as homogeneous isotropic turbulence (HIT) (Kerr Reference Kerr1985; Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987; Cantwell Reference Cantwell1993), channel flow (Blackburn et al. Reference Blackburn, Mansour and Cantwell1996; Wallace Reference Wallace2009; Wu et al. Reference Wu, Zhang, Wang, Zou and Chen2020), boundary layers (Ganapathisubramani et al. Reference Ganapathisubramani, Longmire, Marusic and Pothos2005; Bechlars & Sandberg Reference Bechlars and Sandberg2017a ), jet flows (Buxton & Ganapathisubramani Reference Buxton and Ganapathisubramani2010; Hao et al. Reference Hao, Tian and Zhou2021) and mixing layers (Soria et al. Reference Soria, Sondergaard, Cantwell, Chong and Perry1994; Buxton et al. Reference Buxton, de Kat and Ganapathisubramani2013), as well as in more complex flows, for example around single or fractal square grids (Gomes-Fernandes et al. Reference Gomes-Fernandes, Ganapathisubramani and Vassilicos2014; Zhou et al. Reference Zhou, Nagata, Sakai, Suzuki, Ito, Terashima and Hayase2014b ; Paul et al. Reference Paul, Papadakis and Vassilicos2017).
Traditionally, the VGT
$\unicode{x1D63C}$
is decomposed into a strain rate tensor
$\unicode{x1D64E}_{A}$
and a vorticity (or rotation rate) tensor
$\mathsf{\boldsymbol{\Omega}}_{A}$
, which are symmetric and skew-symmetric, respectively. However, this decomposition has important flaws. The vorticity tensor cannot disambiguate between pure shear and rigid-body motion (see Kolář (Reference Kolář2007) and Das & Girimaji (Reference Das and Girimaji2020) for a visualisation of these motions). The strain rate tensor also cannot disambiguate between irrotational normal strain and pure shear.
These limitations have prompted efforts to formulate a triple decomposition of
$\unicode{x1D63C}$
in terms of irrotational normal-strain rate, pure shear and rigid-body rotation. Kolář (Reference Kolář2007) proposed such a triple decomposition that requires an orthogonal transformation matrix, which is obtained from the solution of an optimisation problem. When
$\unicode{x1D63C}$
has one real and two complex conjugate eigenvalues, Gao & Liu (Reference Gao and Liu2019) developed the ‘Rortex’ decomposition method, effectively splitting the rotation rate tensor into solid-body rotation and pure shear. The magnitudes of the two components were visualised in a transitional boundary layer in a flat plate. The work of Liu et al. (Reference Liu, Yu and Gao2022) is in the same direction.
Another decomposition was proposed by Keylock (Reference Keylock2018), which splits
$\unicode{x1D63C}$
into two tensors, a normal and a non-normal tensor, using the complex form of the Schur decomposition. Das & Girimaji (Reference Das and Girimaji2020) also derived a triple decomposition by combining the approaches of Keylock (Reference Keylock2018) and Gao & Liu (Reference Gao and Liu2019). They applied the former when
$\unicode{x1D63C}$
has two complex and one real eigenvalue, and the latter when it has three real eigenvalues. Kronborg & Hoffman (Reference Kronborg and Hoffman2023) proved that the standardised real Schur form is a solution to the optimisation problem posed by Kolář (Reference Kolář2007). In the non-rotational case (all eigenvalues are real), the normal tensor of Keylock (Reference Keylock2018) represents irrotational straining motions, and the non-normal tensor represents pure shear at three mutually perpendicular planes (see details in § 2).
The triple decomposition has provided new physical insights. For example, in forced isotropic turbulence, Das & Girimaji (Reference Das and Girimaji2020) found that at all Reynolds numbers, pure shear contributes the most and rigid-body rotation the least in the Frobenius norm of the VGT. They also noticed that enstrophy intermittency is due to the pure shear component of vorticity (rather than solid-body motion). Recently, Arun & Colonius (Reference Arun and Colonius2024) applied the triple decomposition to head-on collision between vortex rings. They found that during the turbulent decay phase, the contributions of the different components of deformation to the velocity gradient strength (which is not stationary) are roughly invariant in time, suggesting an equilibrium partitioning of the VGT.
The complex Schur decomposition of Keylock (Reference Keylock2018) has also been applied to different flows. Beaumard et al. (Reference Beaumard, Buxton and Keylock2019) studied the characteristics of the VGT in a near wake, two axisymmetric jets and a planar mixing layer. They noticed that the effect of non-normality on the turbulent dynamics is more pronounced than what is typically observed in HIT, attributing this to the role of the pressure Hessian. Yu & Lu (Reference Yu and Lu2020) analysed the filtered VGT in compressible mixing layers and identified the normal effect as the primary factor influencing compressibility, whereas the behaviour of subgrid-scale energy dissipation was found to be primarily linked to the non-normal part. Yu et al. (Reference Yu, Zhao and Lu2021) found that the non-normal component has a substantial impact on the behaviour of the VGT in wall-bounded compressible turbulent flows, where non-normal effects influence both enstrophy and total strain. Keylock (Reference Keylock2022) studied the connection between the dynamics of the VGT and ejection-sweep events in a channel flow. The study revealed that the non-normal terms above the viscous sublayer had a significantly diminished impact on the dynamics of the sweeps.
The aforementioned decompositions of the VGT can also be used to study the relative contributions of the components of the production terms of enstrophy and turbulent dissipation. The two main production terms are due to vortex stretching and strain self-amplification. The relevant expressions were derived in Keylock (Reference Keylock2018) for the normal/non-normal decomposition of the fluctuating VGT. However, mean flow inhomogeneity, i.e. mean strain and vorticity, also contributes to the generation of enstrophy and dissipation. Currently, it is not known how the normal/non-normal decompositions of mean and fluctuating VGTs interact to contribute to the production of these two quantities.
The properties of the normal/non-normal decomposition of the VGT allow us to gain additional insight into the production terms, especially vortex stretching. More specifically, as is shown in § 2, the normal part carries the eigenvalues of the VGT. When all eigenvalues are real, which occurs in regions 4 and 6 of the QR diagram (where
$Q$
and
$R$
are the second and third invariants of the VGT, respectively), the normal part contains only irrotational strain along the three perpendicular eigendirections, and the normal vorticity component vanishes. In this case, production due to vortex stretching is only due to the non-normal vorticity. Region 6, however, is the most populated region, but only the non-normal vorticity produces stretching. This motivates the investigation of the contributions of different QR regions (and the normal and non-normal strain and vorticity components within each region) to the production of enstrophy and turbulent dissipation.
Within the aforementioned context, in the present work, we investigate the properties of small-scale statistics (enstrophy and dissipation), as well as the VGT and its Schur decomposition in a three-dimensional, separated and transitioning shear flow. Our motivation stems from a lack of understanding regarding the effect of strong flow inhomogeneity on the dynamics of inertial and dissipative ranges.
The velocity field is obtained from a direct numerical simulation of the flow around a NACA 0018 wing with a square wingtip profile at
$10^\circ$
angle of attack. This is the same configuration as that examined in Bilbao-Ludena & Papadakis (Reference Bilbao-Ludena and Papadakis2023), but in the present work, the Reynolds number is lower (
$Re_c=5000$
) and the resolution is finer. Under these conditions, the flow separates from the leading edge, forming a transitioning shear layer that does not reattach to the wing surface.
As mentioned above, the mean flow inhomogeneity introduces production terms in the transport equations of enstrophy and total strain (proportional to dissipation) that involve both fluctuating and mean strain and vorticity. We derive new expressions for the decomposition of these terms into normal, non-normal and mixed (interaction) components. Production terms, conditionally averaged on the QR regions with the largest contributions, are also shown.
Our work considers the following questions. (i) How does the presence of mean shear alter the generation of enstrophy and total strain compared with HIT across the recirculation zone and along the transitioning shear layer? (ii) What is the contribution of each region of the QR plane to the production of enstrophy and total strain? (iii) What is the effect of non-normality in the VGT on vortex stretching and strain self-amplification in this flow? What are the dominant components in the Schur decomposition of these terms in areas of the flow that deviate strongly from HIT? (iv) What is the role of non-normality in the production terms of enstrophy and total strain due to mean flow inhomogeneity?
This paper is organised as follows. The traditional and Schur decompositions of the VGT are detailed in § 2. The decomposition of the production terms of the enstrophy and total strain transport equations for inhomogeneous flows into normal, non-normal and mixed components is derived in § 3. Section 4 briefly describes the flow setting and computational aspects of the simulation. Results are presented and discussed in § 5, and we conclude in § 6.
2. Decompositions of the VGT
In the following, the velocity components in the three directions
$(x_1,x_2,x_3)$
are denoted by the vector
$U=(U_1,U_2,U_3)$
. Time averages are designated with angle brackets, and fluctuations are represented with lowercase variables. For example,
$\langle U_{i} \rangle$
and
$u_i$
are the mean and fluctuating velocities, respectively, in the
$i$
th direction
$x_i$
, and
$U_i=\langle U_{i} \rangle +u_i$
. The Reynolds decomposition is also applied to all other variables; for example,
$P=\langle P\rangle +p$
for static pressure.
For incompressible flow, the Navier–Stokes equations take the form
$$\frac {\partial U_i}{\partial t}+ \frac {\partial U_i U_j}{\partial x_j} =- \frac {1}{\rho } \frac {\partial P}{\partial x_i}+ \nu \frac {\partial ^{2} U_{i}}{\partial ^2 x_{j}},$$
$$\frac {\partial U_{i}}{\partial x_{i}} =0,$$
where
$\rho$
is the fluid density and
$\nu$
is its kinematic viscosity. Taking the time average and subtracting it from set (2.1), we get the governing equations for the fluctuating velocity field that read
$$\frac {\partial u_i}{\partial t}+ \frac {\partial \langle U_i \rangle u_j}{\partial x_j}+ \frac {\partial u_i \langle U_j \rangle }{\partial x_j} =- \frac {1}{\rho } \frac {\partial p}{\partial x_i}+\nu \frac {\partial ^{2} u_{i}}{\partial ^2 x_{j}}-\frac {\partial \left (u_iu_j- \langle u_i u_j \rangle \right )}{\partial x_j},$$
$$\frac {\partial u_{i}}{\partial x_{i}} =0.$$
We consider first the VGT of the fluctuating velocity field,
$\unicode{x1D63C}$
, with components
$A_{\textit{ij}}=\partial u_{i}/\partial x_{j}$
. The characteristic equation of
$\unicode{x1D63C}$
is
\begin{equation} \lambda ^{3}_{i}+P_{A}\lambda ^{2}_{i}+Q_{A}\lambda _{i}+R_{A}=0, \end{equation}
where
$\lambda _i$
are the eigenvalues, and the tensor invariants
$P_A$
,
$Q_A$
and
$R_A$
are defined as
\begin{equation} P_{A}=\mathrm {Tr}(\unicode{x1D63C})=\lambda _{1}+\lambda _{2}+\lambda _{3} =0, \end{equation}
due to the incompressibility condition (2.2b ), and
$$Q_{A}=\frac {1}{2} \underbrace{\left (\mathrm{Tr}(\unicode{x1D63C})\right )^2}_{=0} -\frac {1}{2} \mathrm {Tr} (\unicode{x1D63C}^{2}) =-\frac {1}{2} \mathrm {Tr} (\unicode{x1D63C}^{2})= \lambda _{1}\lambda _{2}+\lambda _{2}\lambda _{3}+\lambda _{3}\lambda _{1},$$
$$R_{A}=-\frac {1}{3}\mathrm {Tr} (\unicode{x1D63C}^{3})=-(\lambda _{1}\lambda _{2}\lambda _{3}).$$
The discriminant
\begin{equation} \varDelta _{A}=Q^{3}_{A}+\frac {27}{4} R^{2}_{A} \end{equation}
characterises the eigenvalues. For
$\varDelta _{A}\lt 0$
, all eigenvalues are real, whereas for
$\varDelta _{A}\gt 0$
, there is one real and two complex eigenvalues (that come in a conjugate pair). Depending on the signs of
$Q_A$
,
$R_A$
and
$\varDelta _A$
, six regions can be defined in the
$Q_A R_A$
plane (traditionally known as the
$ QR$
diagram); see table 1 discussed later in § 5.3.
Occupancy fractions F(I) of the six QR regions for HIT at
$Re_{\zeta_{i}}=433$
(data from Keylock (Reference Keylock2018)).
Traditionally, the VGT is decomposed into a symmetric tensor
$\unicode{x1D64E}_{A}=1/2 (\unicode{x1D63C}+\unicode{x1D63C}^\top )$
and a skew-symmetric tensor
$\mathsf{\boldsymbol{\Omega}}_{A}=1/2 (\unicode{x1D63C}-\unicode{x1D63C}^\top )$
. That is,
$\unicode{x1D63C}={\unicode{x1D64E}_{A}+\mathsf{\boldsymbol{\Omega}}_{A}}$
, where
$\unicode{x1D64E}_{A}$
and
$\mathsf{\boldsymbol{\Omega}}_{A}$
are known as the strain rate and rotation rate tensors, respectively. Recently, Keylock (Reference Keylock2018) proposed an alternative decomposition of
$\unicode{x1D63C}$
into a normal tensor,
$\unicode{x1D63D}$
, and a non-normal tensor,
$\unicode{x1D63E}$
:
\begin{equation} \unicode{x1D63C}=\unicode{x1D63D}+\unicode{x1D63E}\,. \end{equation}
This is accomplished by performing the complex Schur transform of the VGT:
\begin{equation} \unicode{x1D63C}=\unicode{x1D650}\unicode{x1D64F}\unicode{x1D650}^{\,\ast} , \end{equation}
where
$\unicode{x1D650}$
is a unitary matrix, i.e.
$\unicode{x1D650}\unicode{x1D650}^{\,\ast} =\unicode{x1D650}^{\,\ast} \unicode{x1D650}=\unicode{x1D644}$
, and the asterisk denotes conjugate transpose (we employ complex variables because the eigenvalues of
$\unicode{x1D63C}$
can be complex). Matrix
$\unicode{x1D64F}$
is upper triangular and has the form
\begin{equation} \unicode{x1D64F} = \begin{bmatrix} \lambda _{1} & N_{1,2} & N_{1,3}\\ 0 & \lambda _{2} & N_{2,3} \\ 0 & 0 & \lambda _{3} \\ \end{bmatrix} = \underbrace { \begin{bmatrix} \lambda _{1} & 0 & 0 \\ 0 & \lambda _{2} & 0 \\ 0 & 0 & \lambda _{3} \\ \end{bmatrix} }_{ \equiv \unicode{x1D647}}+\underbrace { \begin{bmatrix} 0 & N_{1,2} & N_{1,3}\\ 0 & 0 & N_{2,3} \\ 0 & 0 & 0 \\ \end{bmatrix}}_{ \equiv \unicode{x1D649}}=\unicode{x1D647}+\unicode{x1D649} , \end{equation}
where
$\unicode{x1D647}$
is a diagonal matrix that contains the eigenvalues of
$\unicode{x1D63C}$
and
$\unicode{x1D649}$
is strictly upper triangular. Using this decomposition of
$\unicode{x1D64F}$
,
$\unicode{x1D63D}$
is defined as
\begin{equation} \unicode{x1D63D}=\unicode{x1D650}\unicode{x1D647} \unicode{x1D650}^{\,\ast} , \end{equation}
and is clearly a normal tensor, i.e.
$\unicode{x1D63D} \unicode{x1D63D}^\ast = \unicode{x1D63D}^\ast \unicode{x1D63D}$
. Tensor
$\unicode{x1D63E}$
is defined as
\begin{equation} \unicode{x1D63E}=\unicode{x1D650}\unicode{x1D649} \unicode{x1D650}^{\,\ast} , \end{equation}
and it is non-normal.
The normal tensor
$\unicode{x1D63D}$
carries the eigenvalues of
$\unicode{x1D63C}$
and can (exclusively) determine the isotropic part of the pressure Hessian tensor
$H_{\textit{ij}}=\partial ^{2} p/\partial x_{i}\partial x_{j}$
,
$\unicode{x1D643}^{\,iso}= (1/3)\mathrm {Tr}(\unicode{x1D643}\,)\unicode{x1D644}$
, where
$\unicode{x1D644}$
is the identity matrix. Indeed, it can be easily shown (see Meneveau Reference Meneveau2011) that
\begin{equation} \unicode{x1D643}^{\,iso}= \frac {1}{3}\mathrm {Tr}(\unicode{x1D643}\,)\unicode{x1D644}=\frac {1}{3}\nabla ^{2} p\unicode{x1D644}=-\frac {1}{3}\mathrm {Tr}\big(\unicode{x1D63C}^{2}\big)\unicode{x1D644}=\frac {2}{3} Q_{A}\unicode{x1D644}. \end{equation}
On the other hand, the deviatoric part of the pressure Hessian
$\unicode{x1D643}^{\,dev}$
requires knowledge of the velocity field in the whole domain (in order to extract pressure from the Poisson equation from which
$\unicode{x1D643}^{\,dev}$
is computed; see also Ohkitani & Kishiba Reference Ohkitani and Kishiba1995). Keeping the isotropic part of the pressure Hessian (and ignoring
$\unicode{x1D643}^{\,dev}$
) in the exact transport equation of the VGT leads to the restricted Euler equations, which can be written in terms of the two invariants
$Q_A$
and
$R_A$
. Since these can be computed directly from the eigenvalues, matrix
$\unicode{x1D63D}$
is sufficient to describe the restricted Euler equations.
When all the eigenvalues
$\lambda _i$
are real, there is a straightforward geometrical interpretation of the tensors
$\unicode{x1D63D}$
and
$\unicode{x1D63E}$
. Tensor
$\unicode{x1D63D}$
is diagonal in the directions defined by the columns of matrix
$\unicode{x1D650}$
; see (2.11). Because
$\unicode{x1D650}$
is unitary, these directions are perpendicular to each other. So,
$\unicode{x1D63D}$
describes normal straining in these directions, and
$\unicode{x1D63E}$
represents pure shear (see also Das & Girimaji Reference Das and Girimaji2020). Note that
$\unicode{x1D649}$
is strictly upper triangular and its eigenvalues are 0 (as are the eigenvalues of
$\unicode{x1D63E}$
). On the other hand, in the eigenvalue decomposition of the VGT,
$\unicode{x1D63C}=\boldsymbol {\mathit {\Phi }} \unicode{x1D647} \boldsymbol {\mathit {\Phi }}^{-1}$
, the columns of
$\boldsymbol {\mathit {\Phi }}$
(i.e. the eigenvectors) are not perpendicular to each other, so the non-normal effects are included multiplicatively in the (skewed) eigenvectors. The benefit of the Schur transform is that it decomposes the normal and non-normal effects additively. When there are two complex eigenvalues, the decomposition cannot be interpreted in terms of pure shear, rigid-body rotation and irrotational strain because the transformation involves complex numbers.
A ratio can be defined to quantify the degree of non-normality of the VGT:
\begin{equation} k_{BC}= \frac { \lVert \unicode{x1D63D} \lVert -\lVert \unicode{x1D63E} \lVert }{ \lVert \unicode{x1D63D}\lVert +\lVert \unicode{x1D63E} \lVert }, \end{equation}
where
$\lVert \cdot \lVert$
denotes the Frobenius norm of a tensor, for example,
$\lVert \unicode{x1D63D} \lVert = \sqrt { \mathrm {Tr}(\unicode{x1D63D}\unicode{x1D63D}^\ast )}$
. Positive values of
$k_{BC}$
indicate that the dynamics is dominated by normal effects, whereas negative values indicate enhanced non-normal dynamics.
We can decompose
$\unicode{x1D63D}$
and
$\unicode{x1D63E}$
into strain and rotation rate tensors, i.e.
\begin{equation} \unicode{x1D63D}=\unicode{x1D64E}_B+\mathsf{\boldsymbol{\Omega}}_{B}, \quad \quad \unicode{x1D63E}=\unicode{x1D64E}_C+\mathsf{\boldsymbol{\Omega}}_{C}, \end{equation}
and since
$\unicode{x1D63C}=\unicode{x1D63D}+\unicode{x1D63E}$
, we have
\begin{equation} \unicode{x1D64E}_A=\unicode{x1D64E}_B+\unicode{x1D64E}_C, \quad \quad \mathsf{\boldsymbol{\Omega}}_{A}=\mathsf{\boldsymbol{\Omega}}_{B}+\mathsf{\boldsymbol{\Omega}}_{C}. \end{equation}
Note that since
$\unicode{x1D63D}$
carries the eigenvalues of
$\unicode{x1D63C}$
, when all
$\lambda _i$
are real (in regions 4 and 6 of the
$Q_AR_A$
plane), the normal tensor
$\unicode{x1D63D}$
contains only irrotational strain and the vorticity component
$\mathsf{\boldsymbol{\Omega}}_{B}$
vanishes.
It can be easily shown that
$$\quad\lVert \unicode{x1D64E}_{A} \lVert ^{2} =\lVert \unicode{x1D64E}_{B} + \unicode{x1D64E}_{C} \lVert ^{2}=\lVert \unicode{x1D64E}_{B}\lVert ^{2} + \lVert \unicode{x1D64E}_{C} \lVert ^{2},$$
$$\lVert \mathsf{\boldsymbol{\Omega}}_{A} \lVert ^{2} =\lVert \mathsf{\boldsymbol{\Omega}}_{B} + \mathsf{\boldsymbol{\Omega}}_{C} \lVert ^{2}=\lVert \mathsf{\boldsymbol{\Omega}}_{B}\lVert ^{2} + \lVert \mathsf{\boldsymbol{\Omega}}_{C} \lVert ^{2},$$
and
$\lVert \unicode{x1D64E}_{C} \lVert ^{2}=\lVert \mathsf{\boldsymbol{\Omega}}_{C} \lVert ^{2}$
(see Appendix A for the proof).
The second invariant
$Q_A$
can be written as
\begin{equation} Q_{A}=-\frac {1}{2} \mathrm {Tr} \big(\unicode{x1D63C}^{2}\big) = -\frac {1}{2} \mathrm {Tr} \left ( \left ( \unicode{x1D64E}_{A}+\mathsf{\boldsymbol{\Omega}}_{A}\right ) ^{2}\right ) =-\frac {1}{2} \left [ \mathrm {Tr} \left ( \unicode{x1D64E}_{A} ^{2}\right ) + \mathrm {Tr} \left ( \mathsf{\boldsymbol{\Omega}}_{A}^{2}\right ) \right ], \end{equation}
because
$\mathrm {Tr} (\unicode{x1D64E}_{A}\mathsf{\boldsymbol{\Omega}}_{A} )=\mathrm {Tr} (\mathsf{\boldsymbol{\Omega}}_{A}\unicode{x1D64E}_{A})=0$
. For the third invariant
$R_A$
, we have
\begin{equation} R_{A}=-\frac {1}{3}\mathrm {Tr}\big(\unicode{x1D63C}^{3}\big)=-\frac {1}{3}\mathrm {Tr}\big( \unicode{x1D63D}^{3}\big)=R_{B}, \end{equation}
because
$\unicode{x1D63C}$
and
$\unicode{x1D63D}$
have the same eigenvalues. Thus,
\begin{equation} R_{A} = -\frac {1}{3}\mathrm {Tr}\big(\unicode{x1D64E}_{A}^{3}\big)-\mathrm {Tr}\big (\mathsf{\boldsymbol{\Omega}}_{A}^{2}\unicode{x1D64E}_{A}\big )=R_{B}=-\frac {1}{3}\mathrm {Tr}\big(\unicode{x1D64E}_{B}^{2}\big)-\mathrm {Tr}\big (\mathsf{\boldsymbol{\Omega}}_{B}^{2}\unicode{x1D64E}_{B}\big ). \end{equation}
Also
$R_C=-{1}/{3} \mathrm {Tr} (\unicode{x1D63E}^3 )=-{1}/{3} \mathrm {Tr} (\unicode{x1D650} \unicode{x1D649}^3 \unicode{x1D650}^{\,\ast} )=-{1}/{3}\mathrm {Tr} (\unicode{x1D649}^3 )=0$
, because
$\unicode{x1D649}^3$
has 0 entries in the main diagonal and the trace is invariant under a similarity transformation. Therefore we have
\begin{align} R_C &=-\frac {1}{3}\mathrm {Tr}\big(\unicode{x1D64E}_{C}^{3}\big)-\mathrm {Tr}\big (\mathsf{\boldsymbol{\Omega}}_{C}^{2}\unicode{x1D64E}_{C}\big )=0 \Rightarrow \mathrm {Tr}\big(\unicode{x1D64E}_{C}^{3}\big) = - 3 \mathrm {Tr}\big (\mathsf{\boldsymbol{\Omega}}_{C}^{2}\unicode{x1D64E}_{C}\big ) \Rightarrow \det (\unicode{x1D64E}_{C}) \notag\\&=-\mathrm {Tr}\big (\mathsf{\boldsymbol{\Omega}}_{C}^{2}\unicode{x1D64E}_{C}\big ), \end{align}
because
$\mathrm {Tr} (\unicode{x1D64E}_{C}^{3} )=3 \det (\unicode{x1D64E}_{C})$
.
We can also define the VGT of the time-average velocity field,
$\langle A_{\textit{ij}} \rangle =\partial \langle U_{i} \rangle /\partial x_{j}$
, which can be also decomposed into normal and non-normal tensors,
$\langle \unicode{x1D63C} \rangle =\langle \unicode{x1D63D} \rangle + \langle \unicode{x1D63E} \rangle$
. As before, these tensors can be further decomposed into the strain and vorticity rate tensors, i.e.
$\langle \unicode{x1D63D} \rangle =\unicode{x1D64E}^{M}_{B} + \mathsf{\boldsymbol{\Omega}}^{M}_{B}$
and
$\langle \unicode{x1D63E}\rangle =\unicode{x1D64E}^{M}_{C} + \mathsf{\boldsymbol{\Omega}}^{M}_{C}$
, where the superscript
$M$
indicates mean flow.
3. Production terms of enstrophy and dissipation and their Schur decomposition in inhomogeneous flows
The transport equation for enstrophy
$\langle \omega _{i}\omega _{i} \rangle$
is (Tennekes & Lumley Reference Tennekes and Lumley1972)
\begin{align} \frac {1}{2} \frac {\partial \langle \omega _{i}\omega _{i} \rangle }{\partial t}={}&\underbrace { - \langle U_{j} \rangle \frac {\partial }{\partial x_{j}} \frac {1 }{2} \langle \omega _{i}\omega _{i} \rangle }_{\mathcal {C}_{\omega }} \underbrace {- \langle u_{j} \omega _{i} \rangle \frac {\partial \langle \Omega _{i} \rangle }{\partial x_{j}}}_{\mathcal {P}_{\omega }} \underbrace {- \frac {1 }{2} \frac {\partial }{\partial x_{j} } \langle u_{j} \omega _{i} \omega _{i} \rangle }_{\mathcal {T}_{\omega }} \underbrace {+ \langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle }_{\mathcal {EF}_{\omega }} \notag \\& \underbrace {+ \langle \omega _{i} \omega _{j} \rangle \langle S_{\textit{ij}}\rangle }_{\mathcal {EM}_{\omega }} \underbrace {+ \langle \omega _{i} s_{\textit{ij}}\rangle \langle \Omega _{j} \rangle }_{\mathcal {M}_{\omega }} \underbrace {+ \nu \frac {\partial ^{2}}{\partial x_{j}\partial x_{j}} \left\langle \frac {1}{2} \omega _{i}\omega _{i}\right\rangle }_{\mathcal {D}_{\omega }} \underbrace {- \nu \left\langle \frac {\partial \omega _{i}}{\partial x_{j}} \frac {\partial \omega _{i} }{\partial x_{j}} \right\rangle }_{\epsilon _{\omega }}\,, \end{align}
where
$\mathcal {C}_{\omega }$
is convection by mean flow,
$\mathcal {P}_{\omega }$
represents production by mean vorticity field,
$\mathcal {T}_{\omega }$
is transport by turbulent fluctuations,
$\mathcal {EF}_{\omega }$
is production via vortex stretching due to fluctuating strain field,
$\mathcal {EM}_{\omega }$
is production due to mean strain rate field,
$\mathcal {M}_{\omega }$
is the mixed production term and
$\mathcal {D}_{\omega }$
and
$\mathcal {\epsilon }_{\omega }$
represent the viscous diffusion and dissipation, respectively. We assume steady state, so the transient term on the left-hand side of the equation is taken to be 0.
All the production terms can be expressed using the normal and non-normal strain and rotation tensors
$\unicode{x1D64E}_B, \unicode{x1D64E}_C$
and
$\mathsf{\boldsymbol{\Omega}}_{B}, \mathsf{\boldsymbol{\Omega}}_{C}$
, respectively. Starting with
$\mathcal {EF}_{\omega }$
, it can easily shown that
$\mathcal {EF}_{\omega }=4 \mathrm {Tr}\langle \mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A}\rangle$
, so we have
\begin{align} \mathrm {Tr}\big(\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A}\big) &=\mathrm {Tr} \big( (\mathsf{\boldsymbol{\Omega}}_{B}+\mathsf{\boldsymbol{\Omega}}_{C})^2(\unicode{x1D64E}_{B}+\unicode{x1D64E}_{C})\big)\notag\\ &= \mathrm {Tr}\left (\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{B}\right )+2\mathrm {Tr}\left (\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{B}\right )+\mathrm {Tr}\left (\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}\right )\notag\\ &\quad +\mathrm {Tr}\left (\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{C}\right )+2\mathrm {Tr}\left (\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{C}\right )+\mathrm {Tr}\left (\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}\unicode{x1D64E}_{C}\right ), \end{align}
but
$\mathrm {Tr} (\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{B} )=\mathrm {Tr} (\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{C} )=\mathrm {Tr} (\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{B} )=0$
(see Appendix B) and using (2.22) we get
\begin{equation} \mathrm {Tr}\big(\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A}\big)=\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B}\unicode{x1D64E}_{B})-\det (\unicode{x1D64E}_{C})+\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}\unicode{x1D64E}_{B} ). \end{equation}
Keylock (Reference Keylock2018) proved this expression ((2.21) in his paper) from the transport equations of the invariants
$Q_B^S$
,
$Q_B^{\Omega}$
,
$Q_C^S$
and
$Q_C^{\Omega}$
, but here we have adopted a derivation from first principles, i.e. starting from the decompositions (2.16) of
$\unicode{x1D64E}_A$
and
$\mathsf{\boldsymbol{\Omega}}_{A}$
and substituting directly into
$\mathrm {Tr} (\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A} )$
. This approach allows us to decompose also the mixed terms, such as
$\mathcal {EM}_{\omega }=\langle \omega _{i} \omega _{j} \rangle \langle S_{\textit{ij}}\rangle$
, that represent enstrophy production due to vortex stretching from the mean strain rate field. To this end, we have
$\mathcal {EM}_{\omega }=4 \mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{A}^{2}\rangle \unicode{x1D64E}_{A}^M )$
, where the superscript
$M$
denotes the mean field, as already mentioned. Similarly to (3.2), we get
\begin{align} \mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}^{M}_{A}\right ) &= \mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}^{M}_{B}\right )+2\mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B}\unicode{x1D64E}^{M}_{B}\right )+\mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}^{M}_{B}\right )\notag\\ &\quad+\mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}^{M}_{C}\right )+2\mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}^{M}_{C}\right )+\mathrm {Tr}\left ( \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}^{M}_{C}\right ), \end{align}
but now none of the terms is identically zero. The reason for this is explained in Appendix B and is due to the fact that the unitary matrices arising from the Schur decomposition of the VGT for the mean and fluctuating fields,
$\unicode{x1D650}^M$
and
$\unicode{x1D650}$
, respectively, are not equal. Therefore,
$\unicode{x1D650}^{\,\ast} \unicode{x1D650}^M \ne \unicode{x1D644}$
. To the best of our knowledge, (3.4) has not appeared in the literature before. A similar expression can be easily derived for the other mixed production term,
$\mathcal {M}_{\omega }=\langle \omega _{i} s_{\textit{ij}}\rangle \langle \Omega _{j} \rangle$
.
Next we examine the transport equation of total strain
$\langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
that reads (Tsinober Reference Tsinober2009)
\begin{align} \frac {1}{2} \frac {\partial \langle s_{\textit{ij}} s_{\textit{ij}} \rangle } {\partial t}={}&\underbrace {- \langle U_{k} \rangle \frac {\partial }{\partial x_{k}} \frac {1 }{2} \langle s_{\textit{ij}} s_{\textit{ij}} \rangle }_{\mathcal {C}_{s}} \underbrace {- \langle u_{k} s_{\textit{ij}} \rangle \frac {\partial \langle S_{\textit{ij}} \rangle }{\partial x_{k}}}_{\mathcal {P}_{s}} \underbrace {- \frac {1}{2} \frac {\partial }{\partial x_{k} } \langle u_{k} s_{\textit{ij}} s_{\textit{ij}} \rangle }_{\mathcal {T}_{s}} \notag\\& \underbrace {- { \frac {1 }{4} \langle \omega _{i} \omega _{j} s_{\textit{ij}}\rangle }}_{\mathcal {EF}_{s}} \underbrace {-\langle s_{\textit{ij}} s_{jk} s_{ki}\rangle }_{\mathcal {SF}_{s}} \underbrace {- 2 \langle s_{\textit{ij}} s_{ik} \rangle \langle S_{kj} \rangle }_{\mathcal {SM}_{s}} \underbrace {+\frac {1}{2} \langle \omega _{i} s_{\textit{ij}} \rangle \langle \Omega _{j} \rangle }_{\mathcal {M}_{s}} \notag\\& \underbrace {- \left\langle s_{\textit{ij}} \frac {\partial ^{2} p}{\partial x_{i} \partial x_{j} } \right\rangle }_{\mathcal {H}_{s}} \underbrace {+ \nu \left\langle s_{\textit{ij}} \nabla ^{2}s_{\textit{ij}} \right\rangle }_{\mathcal {D}_{s}}\,, \end{align}
where
$\mathcal {C}_{s}$
is convection by mean flow,
$\mathcal {P}_{s}$
represents production by mean strain rate gradient,
$\mathcal {T}_{s}$
is transport by turbulent fluctuations,
$\mathcal {EF}_{s}$
is production via vortex stretching due to fluctuating strain rate,
$\mathcal {SF}_{s}$
is strain self-amplification,
$\mathcal {SM}_{s}$
is strain amplification by mean strain rate,
$\mathcal {M}_{s}$
denotes mixed production and
$\mathcal {H}_{S}$
and
$\mathcal {D}_{s}$
are the pressure Hessian and viscous dissipation terms, respectively. Again, we assume steady state, so the transient term is assumed to be 0. It can be easily shown that
$ ( \mathcal {EF}_{s}+\mathcal {SF}_{s} ) - 1/2 \mathcal {EF}_{\omega }=3 \langle R_A \rangle$
, and thus the sign of
$ \langle R_A \rangle$
physically expresses the excess (or deficit) of the production of
$1/2\langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
with respect to the production of
$1/4 \langle \omega _i \omega _i \rangle$
. The previous equation is also valid in the instantaneous sense; thus
$3R_A$
represents the production of
$1/2 s_{\textit{ij}} s_{\textit{ij}}- 1/4 \omega _i \omega _i$
.
Production term
$\mathcal {EF}_{s}=-1/4\mathcal {EF}_{\omega }=-\mathrm {Tr}\langle \mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A}\rangle$
and its decomposition are given by (3.3). Term
$\mathcal {SF}_{s}=-\langle s_{\textit{ij}} s_{jk} s_{ki}\rangle =-\mathrm {Tr}(\unicode{x1D64E}_{A}\unicode{x1D64E}_A \unicode{x1D64E}_A)=-3 \det (\unicode{x1D64E}_{A})$
, and we have
\begin{align} \mathrm {Tr}(\unicode{x1D64E}_{A}\unicode{x1D64E}_A\unicode{x1D64E}_A)=\mathrm {Tr}\big((\unicode{x1D64E}_{B}+\unicode{x1D64E}_C)^3\big)= \mathrm {Tr}\big (\unicode{x1D64E}_{B}^3\big )+3\mathrm {Tr}\big (\unicode{x1D64E}_{B}^2\unicode{x1D64E}_{C}\big )+3\mathrm {Tr}\big (\unicode{x1D64E}_{B}\unicode{x1D64E}_{C}^2\big )+\mathrm {Tr}\big(\unicode{x1D64E}_{C}^3\big).\notag\\[2pt] \end{align}
However,
$\mathrm {Tr} (\unicode{x1D64E}_{B}^3 )=3\det (\unicode{x1D64E}_{B})$
,
$\mathrm {Tr} (\unicode{x1D64E}_{B}^2\unicode{x1D64E}_C )=0$
(see Appendix B),
$\mathrm {Tr} (\unicode{x1D64E}_{B}\unicode{x1D64E}_{C}^2 )=-\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C} \mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_B)$
(see Appendix A) and
$\mathrm {Tr} (\unicode{x1D64E}_{C}^3 )=3\det (\unicode{x1D64E}_{C})$
. Substituting in (3.6), and dividing by
$(-3)$
, we get
\begin{equation} -\det (\unicode{x1D64E}_{A})=-\det (\unicode{x1D64E}_{B})-\det (\unicode{x1D64E}_{C})+\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}\unicode{x1D64E}_{B}). \end{equation}
This is equation (2.20) in Keylock (Reference Keylock2018) proved again through direct substitution of the first of (2.16) into
$\mathrm {Tr}(\unicode{x1D64E}_{A}\unicode{x1D64E}_A\unicode{x1D64E}_A)$
.
We now turn our attention to
$\mathcal {SM}_{s}=-2 \langle s_{\textit{ij}} s_{ik} \rangle \langle S_{kj} \rangle$
, which involves the mean strain rate field. We have
$\mathcal {SM}_{s}=-2 \mathrm {Tr} ( \langle \unicode{x1D64E}_A \unicode{x1D64E}_A \rangle \unicode{x1D64E}_A^M )$
, and applying derivation from first principles, we get the following expression:
\begin{align} \mathrm {Tr}\left (\unicode{x1D64E}_{A} \unicode{x1D64E}_{A} \unicode{x1D64E}^{M}_{A}\right ) &= \mathrm {Tr}\left ( \unicode{x1D64E}_{B}\unicode{x1D64E}_{B} \unicode{x1D64E}^{M}_{B}\right )+2\mathrm {Tr}\left ( \unicode{x1D64E}_{C}\unicode{x1D64E}_{B}\unicode{x1D64E}^{M}_{B}\right )+\mathrm {Tr}\left ( \unicode{x1D64E}_{C}\unicode{x1D64E}_{C} \unicode{x1D64E}^{M}_{B}\right )\notag\\ &\quad+\mathrm {Tr}\left ( \unicode{x1D64E}_{B}\unicode{x1D64E}_{B} \unicode{x1D64E}^{M}_{C}\right )+2\mathrm {Tr}\left ( \unicode{x1D64E}_{C}\unicode{x1D64E}_{B} \unicode{x1D64E}^{M}_{C}\right )+\mathrm {Tr}\left ( \unicode{x1D64E}_{C}\unicode{x1D64E}_{C} \unicode{x1D64E}^{M}_{C}\right ), \end{align}
which is similar to (3.4) and again none of the terms is identically 0. This equation has not appeared before in the literature either. The other mixed production term
$\mathcal {M}_{s}=1/2 \langle \omega _{i} s_{\textit{ij}} \rangle \langle \Omega _{j} \rangle =1/2 \mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{A}\unicode{x1D64E}_A \rangle \mathsf{\boldsymbol{\Omega}}_{A}^M )$
can also be decomposed in the same way. In the following, we apply the previously derived decompositions into normal and non-normal components to the three-dimensional field around a finite wing.
4. Flow configuration and computational details
We consider the flow around a finite NACA 0018 wing with a square wingtip profile at
$10^\circ$
angle of attack. The aspect ratio of the wing is equal to 2, and in order to reduce the computational cost, only half of the wing is simulated. The geometry is the same as in Bilbao-Ludena & Papadakis (Reference Bilbao-Ludena and Papadakis2023), but in the present work, the Reynolds number, based on the chord
$C$
and the free-stream velocity
$U_\infty$
, is
$Re_c=5000$
(in the aforementioned reference, it was
$10^4$
).
The number of finite volume cells is approximately
$80$
million, and the flow was simulated using 1024 cores. It took about
$60$
time units (defined as
$C/U_\infty$
) for the flow to reach a statistically steady state. The simulation was then restarted and continued for another
$200$
time units to collect statistics. Details about the code, validation against reference data from the literature and grid convergence tests for
$Re_c=10^{4}$
are provided in Bilbao-Ludena & Papadakis (Reference Bilbao-Ludena and Papadakis2023). In the following,
$x,y$
and
$z$
represent the streamwise, cross-stream and spanwise directions, respectively. The wingtip is located at
$z/C=-1$
and the symmetry plane at
$z/C=0$
. The velocity components in the three directions (normalised with
$U_\infty$
) are denoted as
$U_1,U_2,U_3$
, respectively. In the present and following sections, the analysis focuses on the
$x{-}y$
plane,
$z/C=-0.2$
, which is close to the symmetry plane.
Contours of the ratio of the grid size to Kolmogorov length scale,
$V^{1/3}/ \eta$
, in the
$x{-}y$
plane at
$z/C=-0.2$
. The dashed line A (located at
$x/C=1.2$
) goes through the recirculating zone and line B (located at
$y/C=0.28$
) marks the separated, transitioning shear layer.
Figure 1 shows contours of the ratio of the grid size (defined as the cubic root of the cell volume,
$V^{ 1/3}$
) to the Kolmogorov length scale; the maximum value is approximately 1.2 at the tip of the recirculation zone. Such fine resolution and long-running time are required to balance the enstrophy and total strain transport equations, (3.1) and (3.5), which contain third-order correlations that can take positive and negative values. Two lines are shown in figure 1: line A crosses the recirculation zone while line B marks the separated and transitioning shear layer. Figure 2(a) indicates sub-Kolmogorov resolution along line A, while figure 2(b) shows that along line B; for
$x/C\approx 1.3{-}1.8$
, the resolution is very close to Kolmogorov.
We are interested in the small-scale dynamics in the recirculating and transitional flow regions. A total of 30 and 90 probes were placed along lines A and B, respectively, to obtain velocity and pressure data for further processing.
Variation of the ratio of the grid size to the Kolmogorov length scale,
$V^{1/3}/ \eta$
along (a) line A and (b) line B. The vertical dashed lines in (a) denote the boundaries of the recirculation zone.
5. Results and discussion
5.1. Mean flow characteristics
We first present and discuss the characteristics of the time-averaged flow. Figure 3(a) shows contours of mean streamwise velocity,
$\langle U_{1} \rangle$
. The flow separates from the leading edge but does not reattach on the suction side due to the high angle of attack. The boundary of the recirculation bubble is marked with a solid white line (contour value
$\langle U_{1} \rangle =0$
), and it extends up to approximately 1.4
$C$
from the leading edge. Downstream of the bubble, figure 4(b) shows that the free-stream velocity is slowly recovered.
Contours of (a) streamwise velocity
$\langle U_{1} \rangle$
, (b) TKE
$k$
and (c) TKE production rate
$-\langle u_iu_j\rangle \langle S_{\textit{ij}} \rangle$
, in the
$x{-}y$
plane at
$z/C=-0.2$
.
Figure 3(b) shows contours of the turbulent kinetic energy (TKE),
$k$
. Inside the recirculating bubble, the values of
$k$
are finite despite the fact that production is negligible (as evidenced by figure 3
c). Instead, turbulence is sustained by turbulent diffusion from the top and bottom shear layers. Higher values of
$k$
are found outside the bubble and extend over a wider region compared with production. In the transitional region,
$k$
values peak at
$ x/C\approx 1.4$
(see also figure 4
b), which is an indication that the transition process at this point has completed. Figure 3(c) shows contours of the TKE production rate,
$-\langle u_iu_j\rangle \langle S_{\textit{ij}} \rangle$
. Regions of strong production are found above and below the recirculation zone due to the presence of the separating and attached shear layers emanating from the suction and pressure sides, respectively.
Distribution of
$\langle U_{1} \rangle$
(black line) and
$k$
(blue line) along (a) line A and (b) line B. The dashed vertical lines in (a) indicate the boundaries of the recirculation zone.
The Taylor microscale in the
$i$
th direction is defined as
\begin{equation} \zeta _{i}=\sqrt {\frac {\langle u_{i} u_{i} \rangle }{ \left\langle \frac {\partial u_{i}} {\partial x_{i}} \right\rangle ^2 }}, \end{equation}
and the corresponding Reynolds number is
\begin{equation} Re_{\zeta _{i}}=\frac { { \sqrt { u_{i}u_{i} }} \zeta _{i}} {\nu }. \end{equation}
On average, the Taylor microscales are about
$5\,\%$
of the chord length. Figure 5 shows the variation of
$Re_{\zeta _{i}}$
along lines A and B, with maximum values of approximately
$80$
reached. This range corresponds to the lower end of HIT studies (Iyer et al. Reference Iyer, Sreenivasan and Yeung2022), but is higher compared with other direct numerical simulation studies of complex, three-dimensional flows that focus on small scales; for example, in the wake of square grid elements where
$Re_{\zeta _{i}}=30{-}40$
(Zhou et al. Reference Zhou, Nagata, Sakai, Suzuki, Ito, Terashima and Hayase2014a
, Reference Zhou, Nagata, Sakai, Ito and Hayase2016; Paul et al. Reference Paul, Papadakis and Vassilicos2017). In the present case, there are no homogeneous directions, and very fine resolution as well as long integration times are needed in order to balance the transport equations of enstrophy and total strain; this necessitates smaller values of
$Re_{\zeta _{i}}$
compared with HIT. Figure 5(a) shows the largest deviations between the three Reynolds numbers are outside the recirculation zone. Note that the high values of
$Re_{\zeta _{1}}$
upstream of transition are unrealistic, as
$k$
values are close to zero (see figure 4
b); we mark this region with dashed lines.
Reynolds number based on the Taylor microscale (a) along line A and (b) along line B. Dashed lines indicate the region upstream of transition.
Figure 6 shows contours of
$0.5 \langle \omega _{i}\omega _{i} \rangle$
and
$ \langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
. The two quantities are maximised at the tip of the recirculation zone and have similar values. Figure 7 shows their distribution along lines A and B, and it can be seen that
$0.5 \langle \omega _{i}\omega _{i} \rangle$
has two peaks above and below the boundary of the recirculation zone. Using (2.19), it can be easily shown that the time average of the second invariant of the VGT is
$\langle Q_A \rangle =0.5 (0.5 \langle \omega _{i}\omega _{i} \rangle -\langle s_{\textit{ij}} s_{\textit{ij}} \rangle )$
. For homogeneous (but not necessarily isotropic) turbulence (HT),
$Q_A$
can be written in divergence form (see Davidson Reference Davidson2015); thus,
$\langle Q^{HT}_{A}\rangle =0$
, so
$0.5 \langle \omega _{i}\omega _{i} \rangle ^{HT}=\langle s_{\textit{ij}} s_{\textit{ij}} \rangle ^{HT}$
. This is one of the Betchov (Reference Betchov1956) relations that requires only homogeneity and incompressibility. Therefore, figure 7(a) is the first indication that deviations from HT appear close to the boundaries of the recirculation zone. In figure 7(b), three sections are demarcated along line B, namely growth, peak and decay. The largest deviation from HIT (albeit small) is between the growth and peak sections. In the decay section,
$0.5 \langle \omega _{i}\omega _{i} \rangle \approx \langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
, which agrees well with HT.
Contours of (a)
$0.5 \langle \omega _{i}\omega _{i} \rangle$
and (b)
$ \langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
in the
$x{-}y$
plane at
$z/C=-0.2$
.
Plot of
$0.5 \langle \omega _{i}\omega _{i} \rangle$
(blue line) and
$ \langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
(black line) distribution (a) along line A and (b) along line B.
Figure 8 presents instantaneous snapshots of
$0.5 \omega _{i}\omega _{i}$
and
$ s_{\textit{ij}} s_{\textit{ij}}$
, and visualises why the time averages peak at the tip of the recirculation zone. The vortical structures identified by
$Q_A=2$
are marked with green contour lines. In figure 8(a,c,e),
$0.5 \omega _{i}\omega _{i}$
is pronounced inside the vortical structures; see, for example,
$E1$
,
$E2$
and
$E3$
. The dense areas of
$0.5\omega _{i}\omega _{i}$
generally overlap with high values of
$ s_{\textit{ij}} s_{\textit{ij}}$
, as seen in figure 8(b,d,f). However, this is not always the case; see vortex
$E4$
, for example. At
$t=194$
,
$E1$
and
$E2$
are Kelvin–Helmholtz vortices and are clustered around the downstream end of the bubble. As these eddies are advected downstream, they are strained and interact with vortical structures that originate from the trailing edge, like
$E3$
, generating complex vorticity and strain fields (
$t=194{-}194.4$
) that peak around the tip of the recirculation zone. These large structures break up into smaller structures with a wide range of sizes, as is the case for
$E3$
during
$t=194.2{-}194.4$
. Vortex
$E4$
is an isolated spiky blob of vorticity with a high value of
$0.5 \omega _{i}\omega _{i}$
. Further downstream, the observed structures are more randomly aligned and less densely packed, as exemplified by
$E5$
in figure 8(a).
Sequence of instantaneous contours of (a,c,e)
$0.5 \omega _{i}\omega _{i}$
at (a)
$t=194$
, (c)
$t=194.2$
, (e)
$t=194.4$
and (b,d,f)
$ s_{\textit{ij}} s_{\textit{ij}}$
at (b)
$t=194$
, (d)
$t=194.2$
, (f)
$t=194.4$
. The light green lines mark vortical structures with
$Q_A=2$
, while the white line marks the recirculation zone boundary.
5.2. Budgets of enstrophy and total strain transport equations
We now turn our attention to the enstrophy transport (3.1). Figure 9 shows isosurfaces of the two production terms arising from vortex stretching due to mean strain,
$\mathcal {EM}_{\omega }=\langle \omega _{i} \omega _{j} \rangle \langle S_{\textit{ij}}\rangle$
, and fluctuating strain,
$\mathcal {EF}_{\omega }= \langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle$
(light green and pink colour, respectively). In the background, vortical structures (with
$Q_A=2$
) are shown. Enstrophy production starts right after transition occurs. It can be observed that
$\mathcal {EM}_{\omega }$
appears upstream of
$\mathcal {EF}_{\omega }$
. The two production terms overlap at
$x/C \approx 1.4$
, i.e. around the tip of the recirculation zone.
Isosurfaces of enstrophy production due to vortex stretching by the mean strain rate
$\mathcal {EM}_{\omega }=\langle \omega _{i} \omega _{j} \rangle \langle S_{\textit{ij}}\rangle$
and strain rate fluctuations
$\mathcal {EF}_{\omega }= \langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle$
(light green and pink colour, respectively, with a value of 300). Vortical structures are represented by the background grey isosurface with
$Q_A=2$
.
Budgets of the transport equation of
$0.5\langle \omega _{i}\omega _{i} \rangle$
along (a) line A and (b) line B. The dashed vertical lines in (a) denote the boundaries of the recirculation zone and in (b) demarcate the growth, peak and decay sections (see also figure 7
b).
Budgets of the transport equation of
$\langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
along (a) line A and (b) line B. The dashed vertical lines in (a) denote the boundaries of the recirculation zone and in (b) demarcate the growth, peak and decay regions.
Balancing the enstrophy equation is challenging for the flow examined due to the lack of homogeneous directions. In order to increase the sample size, we include data over a small spanwise distance,
$0.08C$
, around
$z/C=-0.2$
. The budgets of the transport equation for enstrophy along lines A and B are shown in figure 10. Outside the recirculation zone (see figure 10
a), there are two peaks of production due to vortex stretching by the mean flow,
$\mathcal {EM}_{\omega }$
, but within the bubble, vortex stretching due to small-scale strain,
$\mathcal {EF}_{\omega }$
, becomes dominant, and it is balanced by dissipation. The turbulent transport term,
$\mathcal {T}_{\omega }$
, changes sign outside the recirculation zone, while inside this zone it has negligible values. The convection term,
$\mathcal {C}_{\omega }$
, also plays a crucial role in balancing the equation outside the recirculation zone.
In the transition region (see figure 10
b), the primary contributor to the rapid increase in enstrophy in the growth section is the mixed production term,
$\mathcal {EM}{\omega }$
. This term reaches its maximum before the peak of enstrophy, while
$\mathcal {EF}{\omega }$
continues to grow at a much slower pace, reaching its peak further downstream. Term
$\mathcal {EM}{\omega }$
reaches values slightly higher than the maximum of
$\mathcal {EF}{\omega }$
. This picture is consistent with the instantaneous snapshot of figure 9. In the decay section,
$\mathcal {EM}_{\omega }$
is reduced smoothly due to the lack of mean shear, and
$\mathcal {EF}{\omega }$
also decays, but at a slower rate. The balancing error (indicated by the solid black line) along line A is approximately
$2\,\%{-}4\,\%$
and reaches approximately
$13\,\%{-}14\,\%$
with respect to the production peak along B. This is not perfect, but it is within acceptable limits. Our choices of low Reynolds number, fine resolution and long integration time have all been instrumental in achieving this.
We now focus on the transport equation for
$\langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
, (3.5). Again we average over a small spanwise distance
$0.08C$
around
$z/C=-0.2$
to increase the sample size. The budgets along lines A and B are shown in figure 11. Inside the recirculation bubble (see figure 11
a), production is sustained by the fluctuating strain self-amplification term
$\mathcal {SF}_{s}$
and balanced mainly by dissipation (and with vortex stretching
$\mathcal {EF}_{s}$
having smaller contribution). There are two peaks of
$\mathcal {SM}_{s}$
(production by mean strain) outside the recirculation zone that bracket the peaks of
$\mathcal {SF}_{s}$
close to the recirculation boundary. The convection term
$\mathcal {C}_{s}$
plays an important role outside the recirculation zone. Along line B (see figure 11
b)
$\mathcal {SF}_{s}$
and
$\mathcal {SM}_{s}$
are the main production terms. The peak of
$\mathcal {SF}_{s}$
is delayed with respect to that of
$\mathcal {SM}_{s}$
but the value is significantly higher. Vortex stretching
$\mathcal {EF}_{s}$
and dissipation
$\mathcal {D}_{s}$
act as the main sink terms. The other terms, such as pressure Hessian
$\mathcal {H}_{s}$
and convection
$\mathcal {C}_{s}$
play a secondary role. The balancing error along line A is
$\approx 8\,\%$
, and along line B the maximum value is
$\approx 11\,\%$
of the peak value of
$\mathcal {SF}_{s}$
. Again the balancing is not perfect, yet the percentage error is small, which makes the results acceptable.
5.3. The QR regions and their contributions to vorticity stretching and strain self-amplification
In this section, we turn our attention to the production terms,
$\mathcal {EF}_{\omega }=\langle \omega _i \omega _j s_{\textit{ij}}\rangle$
and
$\mathcal {SF}_{s}= - \langle s_{\textit{ij}} s_{jk} s_{ki}\rangle$
, and consider the contributions of the six regions of the QR diagram to these terms. To simplify the notation, we refer to the joint probability density function of the
$Q_A$
and
$R_A$
invariants as the QR diagram, instead of the
$Q_A R_A$
diagram. As a reference, the flow topology of the six regions and their occupancy fractions (defined as the percentage of realisations) for HIT at
$Re_{\zeta_{i}}=433$
are shown in table 1 (data from Keylock (Reference Keylock2018), who processed the results stored in the John Hopkins database (Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008)). The
$QR$
diagram is traditionally divided into four regions based on the signs of the discriminant
$\varDelta _A$
and of
$R_A$
. However, in the present study, we divide it into six regions, based on the signs of three parameters (
$Q_A$
,
$R_A$
and
$\varDelta _A$
), in order to compare directly with the HIT results of table 1.
Figure 12 shows the joint probability density function in the decay section (segment
$x/C=2.46{-}3.43$
) of line B. At this spatial location, the flow is acquiring a quasi-HIT behaviour, and contours resemble the universal teardrop shape. Similar shapes have been reported in mixing layers (Soria et al. Reference Soria, Sondergaard, Cantwell, Chong and Perry1994), boundary layers (Chacin & Cantwell Reference Chacin and Cantwell2000; Bechlars & Sandberg Reference Bechlars and Sandberg2017b
), turbulent jets (Beaumard et al. Reference Beaumard, Buxton and Keylock2019) and single-square grid-generated turbulence (Paul et al. Reference Paul, Papadakis and Vassilicos2017).
The QR diagram for the segment
$x/C=2.46{-}3.43$
of the decay section of the transitional flow. The six QR regions are marked. The dashed line indicates the locus of points where
$\varDelta _A=0$
. Contours are in logarithmic scale.
Figure 13 shows the occupancy fractions of each QR region along lines A and B. The results were obtained by processing the signals of the instantaneous
$Q_A$
and
$R_A$
invariants at each point and computing the fraction of time the VGT spends in each region. The horizontal dashed lines indicate the HIT values reported in table 1. In figure 13(a), the occupancy of region 2 has two clear peaks at
$y/C\approx 0.05$
and
$y/C\approx 0.25$
, which correspond to the peaks of
$k$
,
$0.5 \langle \omega _{i}\omega _{i} \rangle$
and
$ \langle s_{\textit{ij}} s_{\textit{ij}} \rangle$
, while region 6 has a clear peak in the middle of the bubble, matching the HIT value. Away from the recirculation zone, enstrophy is zero; as a result, only regions 4 and 6 are occupied. In figure 13(b), the dominant region 2 in the growth section gives way to region 6 in the peak section. The latter dominates over half of the decay section. At the end of the domain, the occupancy of region 2 increases, closely matching the value for HIT. However, region 6 has lower occupancy than HIT, and there are also small differences in the other regions. This may be due to the flow inhomogeneity, and also the fact that the local
$Re_{\zeta_{i}}$
is much smaller than the HIT value.
Occupancy fractions of QR regions
$I=1{-}6$
along (a) line A and (b) line B. The horizontal dashed lines correspond to HIT for
$Re_{\zeta_{i}}=433$
(see table 1).
Contributions of different QR regions to (a,b)
$\langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle$
and (c,d)
$-\langle s_{\textit{ij}} s_{jk} s_{ki} \rangle$
along lines A and B.
The contributions of each QR region to the time-average values of vortex stretching,
$\langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle$
, and strain self-amplification,
$-\langle s_{\textit{ij}} s_{jk} s_{ki}\rangle$
, are shown in figure 14. For
$\langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle$
, regions 1 and 2 are the leading contributors across the recirculation zone, while the other regions play a secondary role. The dominance of region 2 can be partly explained by the higher occupancy, as indicated in figure 13. Region 1, where
$Q_A\gt 0$
,
$R_A\gt 0$
, acts as a sink. The contribution of the rest of the regions is very small, even though they have similar occupancy to region 1. In figure 14(b), region 2 is also dominant. Chu et al. (Reference Chu, Wang and Lu2014) discovered that the small-scale structures belonging to the SF/S topology (regions 2 and 3) in the outer layer of compressible turbulent boundary layers play the most important role in enstrophy production. At the peak section, there is a mild contribution also from the other regions, especially regions 1 and 6. For
$-\langle s_{\textit{ij}} s_{jk} s_{ki}\rangle$
(see figure 14
c,d), region 6 is the main contributor (approximately one half) along both lines A and B. The rest of the regions collectively make up approximately the other half of the time-average value.
5.4. The normal and non-normal components of the VGT and their contribution to the production of enstrophy and dissipation
The spatial variation of the non-normality index,
$k_{BC}$
, defined in (2.14), is shown in figure 15(a,b). Inside the recirculation zone, and for
$x/C\approx 3{-}3.5$
in the transitional flow section, the values of
$ k_{BC}$
match approximately those of HIT. In the growth area (
$x/C=1.2{-}1.3$
), velocity gradients are highly non-normal, whereas downstream of the peak (
$x/C=1.8{-}2.5$
), they smoothly return to normal (
$ k_{BC}\gt 0$
); see figure 15(b). In the growth section, Kelvin–Helmholtz instabilities grow by extracting energy from the separated shear-layer velocity profile; this process is absent in HIT. In figure 15(c,d), the profiles of the occupancy fraction of non-normality (across all QR regions) are shown. Inside the recirculation zone, at
$y/C\approx 0.1{-}0.2$
, and at the end of the decay regions,
$x/C\approx 3{-}3.5$
, the occupancy fraction of non-normality is
$\approx 50\,\%$
, close to the HIT value. The areas with the most frequent non-normality are found outside the recirculating bubble and in the growth section. Notice the inverse trend between
$k_{BC}$
and
$F(I)|_{k\lt 0}$
in figure 15(b,d); the highest
$k_{BC}$
is correlated with the smallest occupancy fraction.
Variation of the non-normality index
$k_{BC}$
along (a) line A and (b) line B, and occupancy fraction
$F(I)$
of non-normality (
$k_{BC}\lt 0$
) across all QR regions (
$I=1{-}6$
) along (c) line A and (d) line B. The dashed pink lines correspond to HIT values for
$Re_{\zeta {i}}=433$
(see Keylock Reference Keylock2018).
Figure 16 shows the variation of
$\langle \| \mathsf{\boldsymbol{\Omega}}_{A} \|^{2} \rangle$
(figure 16
a,b) and
$\langle \| \unicode{x1D64E}_{A} \|^{2} \rangle$
(figure 16
c,d), along with their normal and non-normal components in the two sections of the flow; refer to (2.17) and (2.18). Inside the recirculating region,
$\langle \| \mathsf{\boldsymbol{\Omega}}_{B} \|^{2} \rangle \approx \langle \| \mathsf{\boldsymbol{\Omega}}_{C} \|^{2} \rangle$
(see figure 16
a), but the non-normal component,
$\langle \| \mathsf{\boldsymbol{\Omega}}_{C} \|^{2} \rangle$
, is slightly higher outside of the bubble. In the transitional zone (see figure 16
b), the two components are of similar magnitude. Figure 16(c,d) shows
$\langle \| \unicode{x1D64E}_{B} \|^{2} \rangle \approx \langle \| \unicode{x1D64E}_{C} \|^{2} \rangle$
in the middle of the recirculating bubble, whereas in the transitional region, the two components are equal from
$x/C\approx 1.6$
until the end of the domain. The results are consistent with the
$k_{BC}$
profiles shown in figure 15(a,b), where
$k_{BC}\approx 0$
. Outside the recirculating zone, however, there is a clear dominance of the non-normal component,
$\langle \| \unicode{x1D64E}_{C} \|^{2} \rangle$
. In the transitional region, as turbulence starts to develop,
$\langle \| \unicode{x1D64E}_{C} \|^{2} \rangle$
grows faster than
$\langle \| \unicode{x1D64E}_{B} \|^{2} \rangle$
, again in agreement with figure 15.
Distribution of the different components of
$\langle \| \mathsf{\boldsymbol{\Omega}}_{A} \|^{2} \rangle$
(a,b) and
$\langle \| \unicode{x1D64E}_{A} \|^{2} \rangle$
(c,d) along line A (a,c) and line B (b,d).
We proceed now with the analysis of the Schur decomposition of the production terms due to strain self-amplification,
$\mathcal {SF}_{s}=-\langle s_{\textit{ij}} s_{jk} s_{ki}\rangle$
, and vortex stretching,
$\mathcal {EF}_{\omega }= \langle \omega _{i} \omega _{j} s_{\textit{ij}} \rangle$
, (3.7) and (3.3), respectively. Figure 17 shows the decomposition of
$-\langle \det (\unicode{x1D64E}_{A}) \rangle$
(equal to
$-1/3 \langle s_{\textit{ij}} s_{jk} s_{ki}\rangle$
). When all QR regions are considered (figure 17
a,b), the interaction (mixed) term,
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}) \rangle$
, is the main contributor to
$-\langle \det (\unicode{x1D64E}_{A}) \rangle$
, followed by the normal,
$-\langle \det (\unicode{x1D64E}_{B}) \rangle$
, component. We also compute the terms conditioned on QR region 6, the dominant region for
$\mathcal {SF}_{s}$
, as shown in figure 14(c,d). The results are shown in figure 17(c,d). Now,
$-\langle \det (\unicode{x1D64E}_{B}) \rangle$
is the dominant term, followed by
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}) \rangle$
. The component
$-\langle \det (\unicode{x1D64E}_{C}) \rangle$
has small values in both averaging approaches (it is close to zero when conditionally averaging on region 6). It has been reported that the median values of this term are close to zero for other flow cases (see Beaumard et al. Reference Beaumard, Buxton and Keylock2019). This indicates that when averaging in region 6, almost all the strain is carried by
$\unicode{x1D64E}_{B}$
, while vorticity is carried by
$\mathsf{\boldsymbol{\Omega}}_{C}$
only, as already mentioned. Since in region 6
$Q_A\lt 0$
, strain is in excess of vorticity, and this explains why the term
$-\langle \det (\unicode{x1D64E}_{B}) \rangle$
dominates over
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}) \rangle$
in this region (for both the recirculating as well as the transitional sections of the flow).
Distribution of
$-\langle \det (\unicode{x1D64E}_{A}) \rangle$
(a,b) and conditionally averaged
$-\langle \det (\unicode{x1D64E}_{A})|_{I=6} \rangle$
(c,d) and their components along line A (a,c) and line B (b,d).
The variation of
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A}) \rangle$
(equal to
$1/4 \langle \omega _{i} \omega _{j} s_{\textit{ij}}\rangle$
) and its Schur components, given by (3.3), is shown in figure 18. Figure 18(a,b) shows the results from averaging over all QR regions, while figures 18(c,d) and 18(e,f) show conditionally averaged results over QR regions 1 and 2, respectively; these are the two dominant regions for this quantity (see figure 14
a,b). When considering all regions (figure 18
a,b), the mixed term,
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}) \rangle$
, is dominant, followed by the normal term,
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{B}) \rangle$
, for both flow sections. When results are conditioned on QR regions 1 and 2 (figure 18
c–f), the normal term,
$\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B} \unicode{x1D64E}_{B}) \rangle$
, is dominant, but the values in these two regions have opposite signs and tend to cancel out when performing standard time averaging over all QR regions. Term
$-\langle \det (\unicode{x1D64E}_{C}) \rangle$
remains mostly positive in all panels, but with significantly lower values compared to the other terms, as in figure 17.
Distribution of
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A}) \rangle$
(a,b) averaged over all QR regions, conditionally averaged over region 1
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A})|_{I=1}\rangle$
(c,d), conditionally averaged over region 2
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{A}^{2} \unicode{x1D64E}_{A})|_{I=2}\rangle$
(e,f) and their Schur components along line A (a,c,e) and line B (b,d,f).
We now explore the Schur decomposition of the mean VGT. First, we analyse the spatial footprint of the norms
$\lVert \unicode{x1D64E}^{M}_{B} \rVert ^{2}$
and
$ \lVert \unicode{x1D64E}^{M}_{C} \rVert ^{2}$
. Recall that
$\lVert \unicode{x1D64E}_{A}^M \lVert ^{2} =\lVert \unicode{x1D64E}_{B}^M\lVert ^{2} + \lVert \unicode{x1D64E}_{C}^M \lVert ^{2}$
; this is the equivalent of (2.17) for the time-average strain rate (denoted by the superscript M). Figure 19 shows contours of
$\lVert \unicode{x1D64E}^{M}_{B} \rVert ^{2}$
and
$ \lVert \unicode{x1D64E}^{M}_{C} \rVert ^{2}$
. It is observed that the peak values of the normal component,
$\lVert \unicode{x1D64E}^{M}_{B} \rVert ^{2}$
, are clustered at the leading edge of the wing and are an order of magnitude smaller than the non-normal component,
$ \lVert \unicode{x1D64E}^{M}_{C} \rVert ^{2}$
. The latter is maximised in the separated shear layers at the suction side and the attached boundary layer at the pressure side. Term
$\lVert \unicode{x1D64E}^{M}_{B} \rVert ^{2}$
captures the normal strain at the leading edge (compression along the stagnation streamline), while
$ \lVert \unicode{x1D64E}^{M}_{C} \rVert ^{2}$
has a much higher footprint in the shear layers (recall that it is equal to
$ \lVert \boldsymbol {\mathit {\Omega }}^{M}_{C} \rVert ^{2}$
; see Appendix A).
Contours of (a)
$ \lVert \unicode{x1D64E}^{M}_{B} \rVert ^{2}$
and (b)
$ \lVert \unicode{x1D64E}^{M}_{C} \rVert ^{2}$
in the
$x{-}y$
plane at
$z/C=-0.2$
.
We can now turn our attention to the mixed production term
$\mathcal {EM}_{\omega }=\langle \omega _{i} \omega _{j} \rangle \langle S_{\textit{ij}}\rangle$
, which represents the generation of enstrophy due to mean strain. Figure 20 shows the distribution of
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{A}^{2}\rangle \unicode{x1D64E}^{M}_{A} )$
(equal to
$1/4\mathcal {EM}_{\omega }$
) and its Schur components; see (3.4). A clear trend is visible: the terms that involve the non-normal tensor
$\unicode{x1D64E}^{M}_{C}$
– i.e.
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B}\rangle \unicode{x1D64E}^{M}_{C} )$
,
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}\rangle \unicode{x1D64E}^{M}_{C} )$
and
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B}\rangle \unicode{x1D64E}^{M}_{C} )$
– are the main contributors at the peaks for both flow sections, whereas the terms involving the normal strain rate tensor,
$\unicode{x1D64E}^{M}_{B}$
, have significantly lower values. The fact that these components are weak is consistent with figure 19, where negligible values of
$\unicode{x1D64E}^{M}_{B}$
are found across the recirculating and transitional regions.
Distribution of
$\mathrm {Tr}(\langle \mathsf{\boldsymbol{\Omega}}_{A}^{2}\rangle \unicode{x1D64E}^{M}_{A})$
and its Schur components along (a) line A and (b) line B.
Finally, we examine the mixed production term
$\mathcal {SM}_{s}=-2 \langle s_{\textit{ij}} s_{ik} \rangle \langle S_{kj} \rangle$
, which represents the generation of strain square due to the mean strain rate field. Figure 21 presents the distribution of
$- \mathrm {Tr} ( \langle \unicode{x1D64E}_A \unicode{x1D64E}_A \rangle \unicode{x1D64E}_A^M )$
(equal to
$1/2 \mathcal {SM}_{s}$
) and its Schur components (see (3.8)) along lines A and B. A similar pattern to that of
$\mathrm {Tr}(\langle \mathsf{\boldsymbol{\Omega}}_{A}^{2}\rangle \unicode{x1D64E}^{M}_{A})$
is observed. The three main contributors are
$-2\mathrm {Tr} (\langle \unicode{x1D64E}_{C}\unicode{x1D64E}_{B}\rangle \unicode{x1D64E}^{M}_{C} )$
,
$-\mathrm {Tr} (\langle \unicode{x1D64E}_{C}\unicode{x1D64E}_{C} \rangle \unicode{x1D64E}^{M}_{C} )$
and
$-\mathrm {Tr} ( \langle \unicode{x1D64E}_{B}\unicode{x1D64E}_{B}\rangle \unicode{x1D64E}^{M}_{C} )$
; they all involve the non-normal strain rate,
$\unicode{x1D64E}^{M}_{C}$
.
We conclude that both mixed production terms are maximised outside the recirculation zone and between the growth and peak regions of the transition section. Both terms involve the mean strain rate,
$\unicode{x1D64E}^{M}$
, and the Schur decomposition showed that the dominant contributions arise from the non-normal component,
$\unicode{x1D64E}^{M}_{C}$
.
Distribution of
$-\mathrm {Tr}(\langle \unicode{x1D64E}_{A}\unicode{x1D64E}_{A}\rangle \unicode{x1D64E}^{M}_{A})$
and its Schur components along (a) line A and (b) line B.
6. Conclusions
We performed direct numerical simulation of the recirculating shear flow around a NACA 0018 wing with a square wingtip profile and investigated the dynamics of the small-scale statistics. The simulation was carried out at a Reynolds number of
$5000$
and
$10^\circ$
angle of attack. The low Reynolds number allowed us to reach Kolmogorov grid resolution and, in certain areas of the flow, sub-Kolmogorov grid resolution. Such fine resolution is required to balance the transport equations of enstrophy and dissipation. Statistical analysis was conducted along two sections of the flow, one that crosses the recirculation zone (line A) and the other along the transitioning shear layer (line B).
We characterised the role of mean shear in the generation of enstrophy and dissipation, a mechanism absent in HIT. Vortex stretching by mean shear (
$\mathcal {EM}{\omega }$
) grows fast in the separated shear layer and peaks earlier than vortex stretching by fluctuating shear (
$\mathcal {EF}{\omega }$
). We computed the occupancy fractions of the six regions of the
$QR$
plane along lines A and B. Areas that approximate HIT were found in the middle of the recirculating flow. Similarly, in the turbulent shear flow along line B, the results (partially) match HIT in the decay section.
In order to investigate the impact of non-normality, we also performed Schur decomposition of the VGT. Within the recirculation zone, the value of the non-normality index
$k_{BC}$
is close to 0, similar to HIT. Outside the recirculation zone, however, the values are (weakly) negative, indicating the presence of non-normal effects. The strongest deviation from HIT, with strong negative
$k_{BC}$
, was detected in the growth section of the separating shear layer.
We also identified the role of normal and non-normal effects on the production of enstrophy and total strain. The dominant component of enstrophy production due to vortex stretching from fluctuating strain is the interaction (mixed) term
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}) \rangle$
. For strain self-amplification, the leading term is also
$\langle \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_{B}) \rangle$
, followed closely by the normal term
$-\langle \det (\unicode{x1D64E}_{B}) \rangle$
. We also applied the Schur decomposition to the mean VGT and derived new decompositions for the mixed production terms that arise due to mean flow inhomogeneity (and thus are absent in HIT). The decomposition has led to the emergence of new terms compared with the fluctuating counterparts. It was found that generation of enstrophy due to vortex stretching from mean strain is mainly due to
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B}\rangle \unicode{x1D64E}^{M}_{C} )$
,
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}\rangle \unicode{x1D64E}^{M}_{C} )$
and
$\mathrm {Tr} (\langle \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{B}\rangle \unicode{x1D64E}^{M}_{C} )$
.
We close the paper with some thoughts on possible future research directions. The incorporation of mean flow inhomogeneity into the Schur decomposition opens up the application to many other inhomogeneous turbulent flows, such as channel flow, boundary layers, impinging jets, forward-facing steps, etc. In the present paper we focused on the analysis of the two production terms, vortex stretching and strain self-amplification. It is well known that vortex stretching always involves a change in length scale (see Tennekes & Lumley Reference Tennekes and Lumley1972) and thus is one of the mechanisms that mediates the cascade of energy from large to small scales. This motives the analysis of the energy cascade through the lens of Schur decomposition. What is the role of the normal/non-normal parts of the VGT in the interscale energy transfer? A mathematically rigorous approach to analyse this process in inhomogeneous flows is the transport equation of the second-order structure function, known as the Karman–Howarth equation (see Hill (Reference Hill2002) for the derivation and Yao et al. (Reference Yao, Mollicone and Papadakis2022) and Yao & Papadakis (Reference Yao and Papadakis2023) for application to a transitioning boundary layer). One can write the nonlinear interscale transfer term (that involves derivatives of triple products of velocity differences) as a function of the non/non-normal tensors at the two points that define the structure function and analyse the different contributions. That would provide a fresh, and hopefully useful, perspective to study this energy transfer process.
Acknowledgements
The authors would like to thank the anonymous referees for bringing to our attention relevant literature and for the many pertinent and insightful comments that have significantly improved the paper.
Funding
We thank also the UK Turbulence Consortium for providing computational time at the UK supercomputing facility ARCHER2 via EPSRC grant no. EP/R029326/1. Additional computational resources were provided by the UK Materials and Molecular Modelling Hub, which is partially funded by EPSRC (EP/P020194/1, EP/W032260/1 and EP/T022213/1) and Imperial College London. J.C.B.-L. would like to acknowledge the financial support by the National Program of Scholarships and Educational Credit (PRONABEC) of Peru as part of the Beca Generación Bicentenario PhD scholarship. G.P. is supported by EPSRC (EP/X017273/1 and EP/W001748/1) and NERC (NE/Z503861/1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Proof that
$\lVert \unicode{x1D64E}_{C} \lVert ^{2}=\lVert \mathsf{\boldsymbol{\Omega}}_{C} \lVert ^{2}$
and
$\mathrm {Tr}(\unicode{x1D64E}_{C} \unicode{x1D64E}_C \unicode{x1D64E}_B)=-\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C} \mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_B)$
We have
\begin{align} \unicode{x1D64E}_{C}\unicode{x1D64E}_{C}^{\,\ast} &=\frac {1}2(\unicode{x1D63E}+\unicode{x1D63E}^{\,\ast} )\frac {1}2(\unicode{x1D63E}^{\,\ast} +\unicode{x1D63E}) =\frac {1}{4}(\unicode{x1D63E}\unicode{x1D63E}^{\,\ast} +\unicode{x1D63E}\unicode{x1D63E}+\unicode{x1D63E}^{\,\ast} \unicode{x1D63E}^{\,\ast} +\unicode{x1D63E}^{\,\ast} \unicode{x1D63E})\notag \\ &=\frac {1}{4}(\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}\unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}\unicode{x1D650}^{\,\ast} ). \end{align}
The products of two strictly upper (lower) triangular matrices, denoted as
$\unicode{x1D649}\unicode{x1D649}$
and
$\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast}$
, have zero entries on their main diagonal; therefore,
$\mathrm {Tr}(\unicode{x1D649}\unicode{x1D649})=\mathrm {Tr}(\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} )=0$
. Since the trace remains invariant under a similarity transformation,
$\mathrm {Tr}(\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}\unicode{x1D650}^{\,\ast} )=\mathrm {Tr}(\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} )=0$
. However, the matrices
$\unicode{x1D649}\unicode{x1D649}^{\,\ast}$
and
$\unicode{x1D649}^{\,\ast} \unicode{x1D649}$
have non-zero trace; therefore,
\begin{equation} \lVert \unicode{x1D64E}_{C} \lVert ^2= \mathrm {Tr}(\unicode{x1D64E}_{C}\unicode{x1D64E}_{C}^{\,\ast} )= \frac {1}{4} \mathrm {Tr}{(\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}\unicode{x1D650}^{\,\ast} )}. \end{equation}
Tensor
$\mathsf{\boldsymbol{\Omega}}_{C}$
is a skew-Hermitian tensor, i.e.
$-\mathsf{\boldsymbol{\Omega}}_{C}=\mathsf{\boldsymbol{\Omega}}_{C}^\ast$
. Thus,
\begin{align} \mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}^\ast &=\frac {1}2(\unicode{x1D63E}-\unicode{x1D63E}^{\,\ast} ) (-\frac {1}{2}(\unicode{x1D63E}-\unicode{x1D63E}^{\,\ast} )) =\frac {1}{4}(-\unicode{x1D63E}\unicode{x1D63E}+\unicode{x1D63E}\unicode{x1D63E}^{\,\ast} +\unicode{x1D63E}^{\,\ast} \unicode{x1D63E}-\unicode{x1D63E}^{\,\ast} \unicode{x1D63E}^{\,\ast} )\notag \\ &=\frac {1}{4}(-\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}\unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}\unicode{x1D650}^{\,\ast} -\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} ). \end{align}
Since
$\mathrm {Tr}(\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}\unicode{x1D650}^{\,\ast} )=\mathrm {Tr}(\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} )=0$
, we get
\begin{equation} \lVert \mathsf{\boldsymbol{\Omega}}_{C} \lVert ^ 2= \mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}^\ast )=\frac {1}{4} \mathrm {Tr}{(\unicode{x1D650}\unicode{x1D649}\unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D649}\unicode{x1D650}^{\,\ast} )}=\lVert \unicode{x1D64E}_{C} \lVert ^{2}. \end{equation}
Since
$\unicode{x1D64E}_C=\unicode{x1D64E}_C^\ast$
and
$-\mathsf{\boldsymbol{\Omega}}_{C}=\mathsf{\boldsymbol{\Omega}}_{C}^\ast$
, from (A1) and (A3), we have
\begin{equation} \unicode{x1D64E}_{C}\unicode{x1D64E}_{C}+\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}= \unicode{x1D64E}_{C}\unicode{x1D64E}_{C}^{\,\ast} -\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C}^\ast = \frac {1}{2}\unicode{x1D650}(\unicode{x1D649}\unicode{x1D649}+\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} ) \unicode{x1D650}^{\,\ast} =\frac {1}{2}\unicode{x1D650} \unicode{x1D63F} \unicode{x1D650}^{\,\ast} , \end{equation}
where
$\unicode{x1D63F}=\unicode{x1D649}\unicode{x1D649}+\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast}$
, and therefore
\begin{equation} (\unicode{x1D64E}_{C}\unicode{x1D64E}_{C}+\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C})\unicode{x1D64E}_B=\frac {1}{2}\unicode{x1D650}\unicode{x1D63F}\unicode{x1D650}^{\,\ast} \frac {1}{2} (\unicode{x1D63D}+\unicode{x1D63D}^\ast )=\frac {1}{4}\unicode{x1D650} (\unicode{x1D63F}\unicode{x1D647}+\unicode{x1D63F}\unicode{x1D647}^\ast ) \unicode{x1D650}^{\,\ast} . \end{equation}
Tensor
$\unicode{x1D649}\unicode{x1D649}$
$(\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} )$
is strictly upper (lower) diagonal, and
$\unicode{x1D647}$
is a diagonal matrix. Therefore the products
$\unicode{x1D63F}\unicode{x1D647} = \unicode{x1D649}\unicode{x1D649}\unicode{x1D647} + \unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} \unicode{x1D647}$
and
$\unicode{x1D63F}\unicode{x1D647}^\ast =\unicode{x1D649}\unicode{x1D649}\unicode{x1D647}^\ast +\unicode{x1D649}^{\,\ast} \unicode{x1D649}^{\,\ast} \unicode{x1D647}^\ast$
have 0 diagonal entries, thus
$\mathrm {Tr}(\unicode{x1D63F}\unicode{x1D647})=\mathrm {Tr}(\unicode{x1D63F}\unicode{x1D647}^\ast )=0$
. This means that
$\mathrm {Tr} ((\unicode{x1D64E}_{C}\unicode{x1D64E}_{C}+\mathsf{\boldsymbol{\Omega}}_{C}\mathsf{\boldsymbol{\Omega}}_{C})\unicode{x1D64E}_B )=0$
, because trace remains invariant under a similarity transformation; so,
$\mathrm {Tr}(\unicode{x1D64E}_{C} \unicode{x1D64E}_C \unicode{x1D64E}_B)=-\mathrm {Tr}(\mathsf{\boldsymbol{\Omega}}_{C} \mathsf{\boldsymbol{\Omega}}_{C} \unicode{x1D64E}_B)$
.
Appendix B. Triple product traces (with or without mean flow)
We now consider the triple products of two normal tensors and one non-normal tensor. If the tensors involve only the fluctuating field, the traces are identically zero. As an example, let us consider
$\mathrm {Tr} (\unicode{x1D64E}_{B}\unicode{x1D64E}_{B}\unicode{x1D64E}_{C})$
. We can write
\begin{align} \unicode{x1D64E}_{B} \unicode{x1D64E}_{B} &=\frac {1}{2}(\unicode{x1D63D}+\unicode{x1D63D}^\ast )\frac {1}{2}(\unicode{x1D63D}+\unicode{x1D63D}^\ast ) =\frac {1}{4}(\unicode{x1D63D}\unicode{x1D63D}+\unicode{x1D63D}\unicode{x1D63D}^\ast +\unicode{x1D63D}^\ast \unicode{x1D63D}+\unicode{x1D63D}^\ast \unicode{x1D63D}^\ast )\notag \\ &=\frac {1}{4}(\unicode{x1D650}\unicode{x1D647}\unicode{x1D647}\unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D647}\unicode{x1D647}^\ast \unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D647}^\ast \unicode{x1D647}\unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D647}^\ast \unicode{x1D647}^\ast \unicode{x1D650}^{\,\ast} )=\frac {1}{4}\unicode{x1D650}\unicode{x1D63F}\unicode{x1D650}^{\,\ast} , \end{align}
where
$\unicode{x1D63F}=\unicode{x1D647}\unicode{x1D647}+\unicode{x1D647}\unicode{x1D647}^\ast + \unicode{x1D647}^\ast \unicode{x1D647} +\unicode{x1D647}^\ast \unicode{x1D647}^\ast$
is a diagonal matrix. Then
$\unicode{x1D64E}_{C}$
can be written as
\begin{equation} \unicode{x1D64E}_{C}=\frac {1}2(\unicode{x1D63E}+\unicode{x1D63E}^{\,\ast} )=\frac {1}{2}(\unicode{x1D650}\unicode{x1D649}\unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} ). \end{equation}
Multiplying (B1) with (B2) we obtain
\begin{equation} \unicode{x1D64E}_{B} \unicode{x1D64E}_{B} \unicode{x1D64E}_{C}=\frac {1}{4}(\unicode{x1D650}\unicode{x1D63F}\unicode{x1D650}^{\,\ast} )\frac {1}{2}(\unicode{x1D650}\unicode{x1D649}\unicode{x1D650}^{\,\ast} +\unicode{x1D650}\unicode{x1D649}^{\,\ast} \unicode{x1D650}^{\,\ast} ) =\frac {1}{8}\unicode{x1D650}(\unicode{x1D63F}\unicode{x1D649}+\unicode{x1D63F}\unicode{x1D649}^{\,\ast} )\unicode{x1D650}^{\,\ast} . \end{equation}
Products
$\unicode{x1D63F}\unicode{x1D649}$
and
$\unicode{x1D63F}\unicode{x1D649}^{\,\ast}$
have zero diagonal entries, thus
$\mathrm {Tr} (\unicode{x1D63F}\unicode{x1D649})=\mathrm {Tr} (\unicode{x1D63F}\unicode{x1D649}^{\,\ast} )=0$
, and since the trace remains invariant under a similarity transformation,
$\mathrm {Tr} (\unicode{x1D64E}_{B}\unicode{x1D64E}_{B}\unicode{x1D64E}_{C})=0$
. Similarly, we can prove that
$\mathrm {Tr} (\unicode{x1D64E}_{B}\unicode{x1D64E}_{B}\unicode{x1D64E}_{C})=0$
,
$\mathrm {Tr} (\mathsf{\boldsymbol{\Omega}}_{B}\mathsf{\boldsymbol{\Omega}}_{B}\unicode{x1D64E}_{C})=0$
and
$\mathrm {Tr} (\unicode{x1D64E}_{B}\mathsf{\boldsymbol{\Omega}}_{B}\unicode{x1D64E}_{C})=0$
.
On the other hand,
$\mathrm {Tr} (\unicode{x1D64E}_{B} \unicode{x1D64E}_{B} \unicode{x1D64E}_{C}^M)\ne 0$
. This is the same product as before, but the non-normal tensor is now obtained from the mean velocity field (thus the superscript M). The product
$\unicode{x1D64E}_{B} \unicode{x1D64E}_{B}$
is given by (B1), but
$\unicode{x1D64E}_{C}^M$
is written as
\begin{align} \unicode{x1D64E}^{M}_{C} &=\frac {1}{2}\big (\unicode{x1D63E}^{M}+\unicode{x1D63E}^{M^\ast }\big )=\frac {1}{2}\big (\unicode{x1D650}^{M}\unicode{x1D649}^{M}\unicode{x1D650}^{M^\ast }+\unicode{x1D650}^{M}\unicode{x1D649}^{M^\ast }\unicode{x1D650}^{M^\ast }\big ), \end{align}
where
$\unicode{x1D649}^{M}$
and
$\unicode{x1D649}^{M^\ast }$
are the strictly upper and lower triangular tensors corresponding to the time-average velocity field, and likewise for the unitary tensor
$\unicode{x1D650}^{M}$
. Multiplying (B1) with (B4) we obtain
\begin{align} \unicode{x1D64E}_{B} \unicode{x1D64E}_{B} \unicode{x1D64E}^{M}_{C} &=\frac {1}{8} (\unicode{x1D650}\unicode{x1D63F}\unicode{x1D650}^{\,\ast} ) \Big(\unicode{x1D650}^{M}\unicode{x1D649}\unicode{x1D650}^{M^\ast }+\unicode{x1D650}^{M}\unicode{x1D649}^{M^\ast }\unicode{x1D650}^{M^\ast }\Big)\notag \\ &=\frac {1}{8}\Big(\unicode{x1D650}\unicode{x1D63F}\unicode{x1D650}^{\,\ast} \unicode{x1D650}^{M}\unicode{x1D649}\unicode{x1D650}^{M^\ast } +\unicode{x1D650}\unicode{x1D63F}\unicode{x1D650}^{\,\ast} \unicode{x1D650}^{M}\unicode{x1D649}^{M^\ast }\unicode{x1D650}^{M^\ast }\Big). \end{align}
The unitary matrices of the fluctuating and mean velocity fields are different,
$\unicode{x1D650} \ne \unicode{x1D650}^{M}$
; thus
$\unicode{x1D650}^{\,\ast} \unicode{x1D650}^M \ne \unicode{x1D644}$
, and therefore
$\mathrm {Tr}(\unicode{x1D64E}_{B} \unicode{x1D64E}_{B} \unicode{x1D64E}^{M}_{C}) \ne 0$
. This also applies for products that involve the vorticity tensor.
Open access
