<figure><img src="image_1.png" /></figure>

HAL
open science

# A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors

Vincent Mons, Claude Cambon, Pierre Sagaut

► To cite this version:

Vincent Mons, Claude Cambon, Pierre Sagaut. A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors. Journal of Fluid Mechanics, 2016, 788, pp.147-182. 10.1017/jfm.2015.705. hal-01276637

HAL Id: hal-01276637

https://hal.science/hal-01276637v1

Submitted on 19 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A spectral model for homogeneous shear-driven anisotropic turbulence in terms of spherically averaged descriptors

Vincent Mons¹,², Claude Cambon³ and Pierre Sagaut⁴

¹Sorbonne Univ, UPMC Univ Paris 06, UMR 7190, Inst Jean Le Rond d'Alembert, F-75005, Paris, France

²CNRS, UMR 7190, Inst Jean Le Rond d'Alembert, F-75005, Paris, France

³Laboratoire de Mécanique des Fluides et d'Acoustique, Université de Lyon, École Centrale de Lyon and CNRS, UMR 5509, Écully, France

⁴Aix Marseille Universite, CNRS, Centrale Marseille, M2P2 UMR 7340, 13451 Marseille, France

A nonlinear spectral model in terms of spherically averaged descriptors is derived for the prediction of homogeneous turbulence dynamics in the presence of arbitrary mean-velocity gradients. The governing equations for the tensor $\hat{R}_{ij}(k, t)$, the Fourier transform of the two-point second-order correlation tensor, are first closed by an anisotropic eddy-damped quasinormal Markovian procedure. This closure is restricted to turbulent flows where linear effects induced by mean-flow gradients have no essential qualitative effects on the dynamics of triple correlations compared with the induced production effects in the equations for second-order correlations. Truncation at the first relevant order of spectral angular dependence allows us to derive from these equations in vector $k$ our final model equations in terms of the wavenumber modulus $k$ only. Analytical spherical integration results in a significant decrease in computational cost. Besides, the model remains consistent with the decomposition in terms of directional anisotropy and polarization anisotropy, with a spherically averaged anisotropic spectral tensor for each contribution. Restriction of anisotropy to spherically averaged descriptors, however, entails a loss of information, and realizability conditions are considered to quantify the upper boundary of anisotropy that can be investigated with the proposed model. Several flow configurations are considered to assess the validity of the present model. Satisfactory agreement with experiments on grid-generated turbulence subjected to successive plane strains is observed, which confirms the capability of the model to account for production of anisotropy by mean-flow gradients. The nonlinear transfer terms of the model are further tested by considering the return to isotropy (RTI) of different turbulent flows. Different RTI rates for directional anisotropy and polarization anisotropy allow us to correctly predict the apparent delayed RTI shown after axisymmetric expansion. The last test case deals with homogeneous turbulence subjected to a constant pure plane shear. The interplay between linear and nonlinear effects is reproduced, yielding the eventual exponential growth of the turbulent kinetic energy.

**Key words:** homogeneous turbulence, turbulence modelling, turbulent flows
# 1. Introduction

Homogeneous anisotropic turbulence is a very important topic in turbulence theory, since it allows for a detailed analysis of linear and nonlinear effects of the mean-flow gradient on turbulence dynamics. A key point is the description of anisotropy, and the derivation of evolution equations that retain most of the structural and dynamical information on the flow evolution. The usual starting point for Fourier-space description is to consider the spectral tensor $\hat{R}_{ij}(k, t)$, which is defined as the Fourier transform of the two-point second-order correlation tensor $R_{ij}(r, t) = \langle u_i(x, t) u_j(x+r, t) \rangle$. Evolution equations can be found for this tensor, which require some closures for nonlinear cubic terms and pressure-velocity correlation terms. A few closed forms have been proposed based on Heisenberg's transfer models, e.g. Canuto & Dubovikov (1996a,b,c). These models do not give a detailed insight into anisotropy, since they rely on a strictly isotropic transfer term as in homogeneous isotropic turbulence. A more general model can be obtained using higher-order closures in place of the quasilocal energy flux, as illustrated in Weinstock (2013), who used the eddy-damped quasinormal Markovian (EDQNM) closure for the isotropic turbulence of Orszag (1970) to derive expressions for the nonlinear transfer term and the nonlinear pressure-strain contribution. This work also relies on an exact treatment of the linear operators induced by mean-velocity gradients in the governing equation for the spectral tensor $\hat{R}_{ij}(k, t)$, which is permitted by the $k$-space level of description of the proposed model. The latter was used to perform an exhaustive analytical study of homogeneous turbulence subjected to a constant pure plane shear without limitations on time or wavenumber.

As in Weinstock (2013), this level of description in $k$-space is the first step of the present model, which entirely includes the nonlinear closure for transfer terms. First, instead of considering the spectral tensor with all of its components, we use its decomposition in terms of directional anisotropy and polarization anisotropy (Cambon & Rubinstein 2006). A more general representation, say $(\mathcal{E}, Z, \mathcal{H})(k, t)$, first introduced by Cambon & Jacquin (1989), also includes a helicity contribution generated by a helicity spectrum $\mathcal{H}(k, t)$. The helicity spectrum, which is associated with the imaginary part of $\hat{R}_{ij}(k, t)$, is not considered here because it remains zero in homogeneous turbulence if it is not initialized or forced, usually in an unphysical way (see also the trivial helicity equation in Cambon & Jacquin (1989) and Cambon et al. (2013)). In addition to their mathematical origin, the recourse to the two different terms $\mathcal{E}(k, t)$ and $Z(k, t)$ has physical meaning. For instance, the structure-based single-point modelling of Kassinos, Reynolds & Rogers (2001) can be related to these spectra, the dimensionality tensor being derived from the angular dependence of $\mathcal{E}(k, t)$, and the stropholysis tensor from $Z(k, t)$. This structure-based modelling is in contrast to other single-point models where anisotropy is characterized by the sole deviatoric tensor associated with the Reynolds stress tensor which is used to express the 'rapid' and 'slow' parts of the pressure-strain rate tensor. The linkage between our spectral approach based on the directional-polarization anisotropy decomposition and the structure-based single-point modelling of Kassinos et al. (2001) is detailed in appendix A of the present paper, along with discussions concerning other single-point models.

The compact decomposition of the spectral tensor discussed above simplifies, without any loss of generality, the derivation of dynamical equations, here in Line equations for $\mathcal{E}(k, t)$ and $Z(k, t)$, and that of closure relations. On the other hand, the numerical treatment of the angular dependence of $\mathcal{E}(k, t)$ and $Z(k, t)$
remains a challenge. As an example, the numerical cost of axisymmetric EDQNM
models (Godeferd & Cambon 1994; Bellet et al. 2006; Favier et al. 2011) is much
higher than that of the basic EDQNM for isotropic turbulence. These models were
developed, from the first case of rotating turbulence (Cambon & Jacquin 1989), for
flows dominated by interacting dispersive waves and without explicit linear effect in
the governing equation for the energy density $\mathcal{E}(k, t)$, requiring the use of the most
complicated EDQNM version to handle the effects induced by the mean flow on triple
correlations, such as the inhibition of the energy cascade in the case of rotating flows.
This configuration is no longer considered here, and we focus in the present paper
on turbulent flows interacting with mean-velocity gradients with a symmetric part,
hence leading to energy production. At least for the flow configurations considered
in this paper (turbulent flows in the presence of mean strain/shear), these production
effects induced by the mean flow, which originate from the linear contributions in
the equations for second-order correlations, prevail in the evolution of anisotropy
compared with the ones that affect triple correlations (Sagaut & Cambon 2008). A
self-consistent nonlinear model will be built accordingly, based on fully tensorial
EDQNM closure for an arbitrary anisotropic second-order spectral tensor but with no
explicit mean-gradient effect in the expression of the nonlinear transfer terms.

Our second step for deriving our final model equations does not include further
closure assumptions, it is purely technical but essential to derive a tractable model,
especially in terms of numerical treatment. The computational cost of the models
in 3D Fourier space discussed above can be significantly reduced by deriving a 1D
problem through integration over spheres with arbitrary radius in Fourier space, $S_k$.
This yields expressions of governing equations for the spherically averaged tensor
$\varphi_{ij}(k, t) = \int \int_{S_k} \hat{R}_{ij}(k, t) \, d^2k$ that depend on the wavevector modulus, $k$, instead of
$k$. This was done by Clark & Zemach (1995) in their diffusion-approximation-based
model, where the prescription of a scale-dependent relaxation rate leads to satisfactory
comparisons with experimental data and to a good description of partial return to
isotropy, which is confirmed in, e.g., Chasnov (1995). However, discussions in Choi
& Lumley (2001), Kassinos et al. (2001) and Zusi & Perot (2013, 2014) about
possible nonlinear return to isotropy can eventually motivate the use of anisotropic
EDQMN closure and the distinction between directional anisotropy and polarization
anisotropy. One can also mention the model of Cambon, Jeandel & Mathieu (1981),
where a closure for nonlinear terms using anisotropic EDQNM was proposed, but
without a separation of directional anisotropy and polarization anisotropy. Spherical
integration requires a parametrization of the second-order spectral tensor to restore at
least a part of its angular dependence, leading to a mathematically consistent model
reduced to spherically averaged descriptors. This will be done in our final model
by using truncated expansions in terms of angular harmonics of the second-order
spectral tensor. A similar approach was considered by Herring (1974) in the context
of the direct-interaction approximation (DIA) of Kraichnan (1959). In this work,
the question of the number of spherical harmonics necessary to accurately describe
the return to isotropy of axisymmetric turbulence is investigated numerically, and
a scale-by-scale relaxation rate is identified. This truncation of spherical harmonic
series lies at the origin of a loss of information that restricts the present model to
moderate anisotropy. Realizability conditions will be considered in order to quantify
the upper boundary of anisotropy intensity that can be investigated with the proposed
model. A complementary approach to the models described above can be found in
Kassinos & Akylas (2012), where explicit angular dependence is preserved, whereas
scale information is lost after integration along rays of fixed orientation in Fourier space.

The key spherical descriptors of our spectral model are now introduced. Starting from a trace–deviator splitting of the real part of the spectral tensor $\hat{R}_{ij}(\mathbf{k}, t)$ restricted to the plane normal to the wavevector $\mathbf{k}$ by virtue of incompressibility,

$$ \hat{R}_{ij}(\mathbf{k}, t) = \frac{1}{2}\hat{R}_{nn}(\mathbf{k}, t)P_{ij}(\mathbf{k}) + \text{Re}\left(\hat{R}_{ij}(\mathbf{k}, t) - \frac{1}{2}\hat{R}_{nn}(\mathbf{k}, t)P_{ij}(\mathbf{k})\right), \quad P_{ij}(\mathbf{k}) = \delta_{ij} - \frac{k_j k_i}{k^2}, \qquad (1.1a,b) $$

the following threefold decomposition in terms of isotropic, directional and polarization parts is obtained:

$$ \hat{R}_{ij}(\mathbf{k}, t) = \underbrace{\frac{E(k, t)}{4\pi k^2} P_{ij}(\mathbf{k})}_{\hat{R}_{ij}^{(iso)}(\mathbf{k}, t)} + \underbrace{\mathcal{E}^{(dir)}(\mathbf{k}, t) P_{ij}(\mathbf{k})}_{\hat{R}_{ij}^{(dir)}(\mathbf{k}, t)} + \hat{R}_{ij}^{(pol)}(\mathbf{k}, t), \qquad (1.2) $$

by separating the trace of $\hat{R}_{ij}(\mathbf{k}, t)$ into purely isotropic and directional anisotropy contributions, according to

$$ \mathcal{E}(\mathbf{k}, t) = \frac{1}{2}\hat{R}_{nn}(\mathbf{k}, t), \quad \mathcal{E}^{(dir)}(\mathbf{k}, t) = \mathcal{E}(\mathbf{k}, t) - \frac{E(k, t)}{4\pi k^2}, \quad E(k, t) = \iint_{S_k} \mathcal{E}(\mathbf{k}, t) \, d^2\mathbf{k}, \qquad (1.3a,b) $$

where $E(k, t)$ corresponds to the kinetic energy spectrum. The polarization part $\hat{R}_{ij}^{(pol)}(\mathbf{k}, t)$ in (1.2) can be expressed in terms of a complex-valued scalar $Z(\mathbf{k}, t)$. This decomposition yields the following splitting of the spherically averaged tensor $\varphi_{ij}(\mathbf{k}, t)$:

$$ \varphi_{ij}(\mathbf{k}, t) = \iint_{S_k} \hat{R}_{ij}(\mathbf{k}, t) \, d^2\mathbf{k} = 2E(k, t) \left( \frac{1}{3}\delta_{ij} + H_{ij}^{(dir)}(\mathbf{k}, t) + H_{ij}^{(pol)}(\mathbf{k}, t) \right). \qquad (1.4) $$

The dimensionless deviatoric part of the spherically averaged spectral tensor is therefore $H_{ij}(\mathbf{k}, t) = H_{ij}^{(dir)}(\mathbf{k}, t) + H_{ij}^{(pol)}(\mathbf{k}, t)$, with

$$ 2E(k, t)H_{ij}^{(dir)}(\mathbf{k}, t) = \iint_{S_k} \hat{R}_{ij}^{(dir)}(\mathbf{k}, t) \, d^2\mathbf{k}, \quad 2E(k, t)H_{ij}^{(pol)}(\mathbf{k}, t) = \iint_{S_k} \hat{R}_{ij}^{(pol)}(\mathbf{k}, t) \, d^2\mathbf{k}. \qquad (1.5a,b) $$

The aim of this paper is to propose a model based on the decomposition (1.4) of the averaged tensor $\varphi_{ij}(\mathbf{k}, t)$, yielding a system of 11 equations for the global state vector $(E, EH_{ij}^{(dir)}, EH_{ij}^{(pol)})(\mathbf{k}, t)$. This system is obtained in two steps. The first step is based on an anisotropic EDQNM closure for the equations that govern the whole spectral tensor $\hat{R}_{ij}(\mathbf{k}, t)$, or equivalently the scalar spectra $\mathcal{E}(\mathbf{k}, t)$ and $Z(\mathbf{k}, t)$ from the directional–polarization anisotropy decomposition. This closure is restricted to flows where linear effects induced by mean-velocity gradients lie at the origin of energy production which prevails in the evolution of anisotropy compared with the induced effects on the dynamics of triple correlations. The mean flow is taken into account
through the linear operators, which reflect rapid distortion theory (RDT) if considered alone, in the governing equations for second-order moments, as done in Weinstock (2013). The second step consists in parametrizing the angular dependence of $\mathcal{E}(k, t)$ and $Z(k, t)$ in terms of the spherically averaged spectral descriptors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$ in order to derive a more tractable model in $k$-space. Even if spherical integration entails a loss of information that restricts the present model to moderate anisotropy and prevents a complete resolution of the flow in 3D Fourier space as in Weinstock (2013), the model developed in the present paper can be used to calculate homogeneous anisotropic flows in a wide range of configurations, without any restriction to a particular symmetry, including when non-stationary mean-velocity gradients are considered. This point will be illustrated by the different applications of the model considered in this paper.

The paper is organized as follows. In § 2, the governing equations for the spectral tensor $\hat{R}_{ij}(k, t)$ are recalled along with those for the scalar spectra $\mathcal{E}(k, t)$ and $Z(k, t)$. The latter are then closed via an anisotropic EDQNM procedure. The spectral model in terms of spherically averaged descriptors is derived in § 3. Section 4 is devoted to the validation of the present model. To this end, several flow configurations are considered. First, multiple straining processes are addressed and comparisons with various experiments are performed. Several cases of return to isotropy (RTI) are then considered, to assess the capability of the model to account for slow pressure terms and possible lack of RTI with separate investigation of directional anisotropy and polarization anisotropy. The last test case deals with homogeneous turbulence subjected to a constant pure plane shear, to check the capability of the model to recover the asymptotic stage of exponential growth of kinetic energy predicted by, e.g., Weinstock (2013). The linkage between single-point models and the model of the present paper is discussed in appendix A.

## 2. Closed equations for the second-order spectral tensor in sheared turbulence

We consider incompressible homogeneous turbulence. In the presence of a mean-velocity gradient, the Navier–Stokes equation for the fluctuating velocity $u_i(\mathbf{x}, t)$ includes additional advection and deformation terms linked to the mean – or large-scale – velocity field $U_i(\mathbf{x}, t)$:

$$ \left( \frac{\partial}{\partial t} + u_j \frac{\partial}{\partial x_j} \right) u_i + U_j \frac{\partial u_i}{\partial x_j} + \frac{\partial U_i}{\partial x_j} u_j = - \frac{\partial p}{\partial x_i} + \nu \nabla^2 u_i. \quad (2.1) $$

In (2.1), $p(\mathbf{x}, t)$ is the pressure divided by a reference density, and $\nu$ is the kinematic viscosity of the fluid. In § 2.1, we give the governing equation for the second-order spectral tensor $\hat{R}_{ij}(k, t)$. The $(\mathcal{E}, Z)$ decomposition is applied to the latter tensor in § 2.2 and to its governing equation in § 2.3. Transfer terms appear from the latter decomposition, which are closed via the EDQNM procedure described in § 2.4.

### 2.1. Craya's equations

The mean flow $U(\mathbf{x}, t)$ is characterized by a space-uniform gradient $\lambda_{ij}(t) = (\partial U_i / \partial x_j)(t)$ in accordance with homogeneity for the fluctuations. For the sake of readability, and without loss of generality, we will omit the time dependence of $\lambda_{ij}$ in the following. In anisotropic homogeneous turbulence, all information about two-point second-order correlations is provided by the second-order spectral tensor
$\hat{R}_{ij}(\mathbf{k}, t)$, which is the Fourier transform of the two-point second-order correlation tensor $R_{ij}(\mathbf{r}, t) = \langle u_i(\mathbf{x}, t)u_j(\mathbf{x} + \mathbf{r}, t) \rangle$, with $\mathbf{r}$ the vector separating the two points in physical space. This tensor is defined as

$$ \langle \hat{u}_i^*(\mathbf{p}, t) \hat{u}_j(\mathbf{k}, t) \rangle = \delta(\mathbf{k} - \mathbf{p}) \hat{R}_{ij}(\mathbf{k}, t), \qquad (2.2) $$

where $\hat{u}_i(\mathbf{k}, t)$ is the Fourier transform of the fluctuating velocity $u_i(\mathbf{x}, t)$, and the operators * and $\langle \cdot \rangle$ denote complex conjugate and ensemble average respectively. From the counterpart of (2.1) in Fourier space and using definition (2.2) one derives the governing equation for the tensor $\hat{R}_{ij}(\mathbf{k}, t)$ (see, e.g., Sagaut & Cambon (2008)),

$$ \left( \frac{\partial}{\partial t} - \lambda_{ln} k_l \frac{\partial}{\partial k_n} + 2\nu k^2 \right) \hat{R}_{ij}(\mathbf{k}, t) + M_{in}(\mathbf{k}) \hat{R}_{nj}(\mathbf{k}, t) + M_{jn}(\mathbf{k}) \hat{R}_{ni}(\mathbf{k}, t) = T_{ij}(\mathbf{k}, t), \quad (2.3) $$

with the linear operator induced by mean-velocity gradients

$$ M_{ij}(\mathbf{k}) = \left( \delta_{in} - 2 \frac{k_i k_n}{k^2} \right) \lambda_{nj}. \qquad (2.4) $$

The transfer tensor $T_{ij}(\mathbf{k}, t)$ in (2.3) accounts for triadic interactions between vectors $\mathbf{k}$, $\mathbf{p}$ and $\mathbf{q}$ so that they form a triangle. It is possible to disentangle contributions from the tensorial transfer term, with zero integral over $\mathbf{k}$, and contributions from the fluctuating pressure $W_{ij}(\mathbf{k}, t)$ as follows:

$$ T_{ij}(\mathbf{k}, t) = P_{in}(\mathbf{k}) \tau_{nj}(\mathbf{k}, t) + P_{jn}(\mathbf{k}) \tau_{ni}^*(\mathbf{k}, t) = \tau_{ij}(\mathbf{k}, t) + \tau_{ji}^*(\mathbf{k}, t) - \underbrace{\frac{k_i k_n}{k^2} \tau_{nj}(\mathbf{k}, t) - \frac{k_j k_n}{k^2} \tau_{ni}^*(\mathbf{k}, t)}_{W_{ij}(\mathbf{k}, t)}. \qquad (2.5) $$

The tensor $W_{ij}(\mathbf{k}, t)$ contains a possible RTI mechanism, its integral over $\mathbf{k}$ gives the nonlinear – so-called slow – pressure–strain rate tensor, and both the $T_{ij}(\mathbf{k}, t)$ and $W_{ij}(\mathbf{k}, t)$ terms originate from the same tensor $\tau_{ij}(\mathbf{k}, t)$. In the same way, the tensor $\tau_{ij}(\mathbf{k}, t)$,

$$ \tau_{ij}(\mathbf{k}, t) = k_n \iiint S_{ijn}(\mathbf{k}, \mathbf{p}, t) \mathrm{d}^3\mathbf{p}, \qquad (2.6) $$

is given from the third-order three-point spectral tensor $S_{ijn}(\mathbf{k}, \mathbf{p}, t)$, defined by

$$ i \langle \hat{u}_i(\mathbf{q}, t) \hat{u}_j(\mathbf{k}, t) \hat{u}_n(\mathbf{p}, t) \rangle = \delta(\mathbf{k} + \mathbf{p} + \mathbf{q}) S_{ijn}(\mathbf{k}, \mathbf{p}, t), \qquad (2.7) $$

and the closure is applied to the equation that governs the latter tensor, as shown in § 2.4.

## 2.2. The ($\mathcal{E}$, $Z$) decomposition

A general decomposition of the second-order spectral tensor $\hat{R}_{ij}(\mathbf{k}, t)$, for arbitrary anisotropy, results from a trace–deviator splitting, restricted to the plane normal to the wavevector $\mathbf{k}$ by virtue of incompressibility. It is consistent with the threefold decomposition in terms of the isotropic contribution, directional anisotropy and polarization anisotropy, from (1.2) (Cambon & Jacquin 1989; Cambon & Rubinstein 2006; Sagaut & Cambon 2008). The first term in the decomposition (1.2) of $\hat{R}_{ij}(\mathbf{k}, t)$ corresponds to its isotropic part. The second term characterizes ‘directional anisotropy’
via the scalar $\mathcal{E}^{(\text{dir})}(\mathbf{k}, t) = \mathcal{E}(\mathbf{k}, t) - E(\mathbf{k}, t)/(4\pi k^2)$, which corresponds to the difference between the energy density $\mathcal{E}(\mathbf{k}, t)$ and its spherical average. The third term, which is generated by the scalar $Z(\mathbf{k}, t)$, characterizes the 'polarization anisotropy', or tensorial anisotropy, at a given wavevector.

Following Cambon & Jacquin (1989), the contribution from the polarization anisotropy is generated by the single complex-valued pseudoscalar term $Z(\mathbf{k}, t)$ as follows:

$$ \hat{\mathbf{R}}_{ij}^{(\text{pol})}(\mathbf{k}, t) = \text{Re}(Z(\mathbf{k}, t) N_i(\mathbf{k}) N_j(\mathbf{k})), \quad (2.8) $$

with

$$ Z(\mathbf{k}, t) = \frac{1}{2} \hat{\mathbf{R}}_{mn}(\mathbf{k}, t) N_m(\mathbf{k})^* N_n(\mathbf{k})^*, \quad (2.9) $$

where the vector $N(\mathbf{k})$, introduced by Cambon & Jacquin (1989) and recovered independently by Waleffe (1992) as the helical mode, is defined by

$$ N(\mathbf{k}) = \mathbf{e}^{(2)}(\mathbf{k}) - i\mathbf{e}^{(1)}(\mathbf{k}), \quad \mathbf{e}^{(1)}(\mathbf{k}) = \frac{\mathbf{k} \times \mathbf{n}}{|\mathbf{k} \times \mathbf{n}|}, \quad \mathbf{e}^{(2)}(\mathbf{k}) = \mathbf{e}^{(3)}(\mathbf{k}) \times \mathbf{e}^{(1)}(\mathbf{k}), \quad \mathbf{e}^{(3)}(\mathbf{k}) = \frac{\mathbf{k}}{\mathbf{k} \cdot \mathbf{n}}. \quad (2.10a-d) $$

Here, $(\mathbf{e}^{(1)}(\mathbf{k}), \mathbf{e}^{(2)}(\mathbf{k}), \mathbf{e}^{(3)}(\mathbf{k}))$ is an orthonormal right-handed frame of reference associated with a privileged direction $\mathbf{n}$, which is often referred to as the Craya-Herring frame (Herring 1974; Sagaut & Cambon 2008). The realizability condition, or condition for the Hermitian covariance matrix $\hat{\mathbf{R}}(\mathbf{k}, t)$ to be definite-positive, can be written as (Cambon, Mansour & Godeferd 1997)

$$ |Z(\mathbf{k}, t)| \le \mathcal{E}(\mathbf{k}, t) \quad \forall \mathbf{k}, t. \quad (2.11) $$

## 2.3. Lin equations for the ($\mathcal{E}$, $Z$) decomposition

When taking into account the decomposition (1.2) of the second-order spectral tensor $\hat{\mathbf{R}}_{ij}(\mathbf{k}, t)$ with equation (2.8), (2.3) is equivalent to the following set of two equations in terms of $\mathcal{E}(\mathbf{k}, t)$ and $Z(\mathbf{k}, t)$:

$$ \left( \frac{\partial}{\partial t} - \lambda_{ln} k_l \frac{\partial}{\partial k_n} + 2\nu k^2 \right) \mathcal{E}(\mathbf{k}, t) - \mathcal{E}(\mathbf{k}, t) S_{ij} \alpha_i \alpha_j + \text{Re}(Z(\mathbf{k}, t) S_{ij} N_i(\mathbf{k}) N_j(\mathbf{k})) = T^{(\mathcal{E})}(\mathbf{k}, t), \quad (2.12) $$

$$ \left( \frac{\partial}{\partial t} - \lambda_{ln} k_l \frac{\partial}{\partial k_n} + 2\nu k^2 \right) Z(\mathbf{k}, t) - Z(\mathbf{k}, t) S_{ij} \alpha_i \alpha_j + \mathcal{E}(\mathbf{k}, t) S_{ij} N_i^*(\mathbf{k}) N_j^*(\mathbf{k}) - 2iZ(\mathbf{k}, t) \left( \frac{W_l}{2} \alpha_l - \Omega_E \right) = T^{(Z)}(\mathbf{k}, t), \quad (2.13) $$

where $\alpha_i = k_i/k$, $S_{ij} = (\lambda_{ij} + \lambda_{ji})/2$ is the symmetric part of the mean-velocity gradient and $W_i = \epsilon_{ijn}\lambda_{nj}$ refers to its antisymmetric part (mean vorticity). The rotation vector component $\Omega_E$ expresses the solid-body motion of the local Craya frame with respect to a fixed frame of reference, following characteristic lines. The expression of $\Omega_E$ is given by

$$ \Omega_E = - \frac{k}{|\mathbf{k} \times \mathbf{n}|} \lambda_{ln} n_l e_n^{(1)} - \lambda_{ln} e_l^{(2)} e_n^{(1)}. \quad (2.14) $$

The derivation of the above expressions may be found in Cambon et al. (2013), but with an error of sign in front of the rotation terms in (2.13). The nonlinear transfer
terms on the right-hand sides of (2.12) and (2.13) are obtained by applying the ($\mathcal{E}$, Z) decomposition to the transfer term $T_{ij}(\mathbf{k}, t)$ in (2.3):

$$T^{(\mathcal{E})}(\mathbf{k}, t) = \frac{1}{2} T_{ii}(\mathbf{k}, t) = \frac{1}{2}(\tau_{ii}(\mathbf{k}, t) + \tau_{ii}^{*}(\mathbf{k}, t)), \quad (2.15)$$

$$T^{(Z)}(\mathbf{k}, t) = \frac{1}{2} T_{ij}(\mathbf{k}, t) N_i^*(\mathbf{k}) N_j^*(\mathbf{k}) = \frac{1}{2} (\tau_{ij}(\mathbf{k}, t) + \tau_{ji}^*(\mathbf{k}, t)) N_i^*(\mathbf{k}) N_j^*(\mathbf{k}), \quad (2.16)$$

where the tensor $\tau_{ij}(\mathbf{k}, t)$ is defined by equations (2.5) and (2.6). As mentioned earlier, $T_{ij}(\mathbf{k}, t)$ includes both the ‘true’ transfer tensor, with zero integral, and the contribution $W_{ij}(\mathbf{k}, t)$ involved in the RTI effect. The latter tensor can be generated from a scalar transfer term $T^{(RTI)}(\mathbf{k}, t)$ according to

$$W_{ij}(\mathbf{k}, t) = -\mathrm{Re}(T^{(RTI)}(\mathbf{k}, t)(\alpha_i N_j(\mathbf{k}) + \alpha_j N_i(\mathbf{k}))), \quad (2.17)$$

consistently with $\tau_{ij}(\mathbf{k}, t)k_j = 0$, $\tau_{ij}(\mathbf{k}, t)k_i \neq 0$ and

$$T^{(RTI)}(\mathbf{k}, t) = \alpha_i (\tau_{ij}(\mathbf{k}, t) + \tau_{ji}^*(\mathbf{k}, t)) N_j^*(\mathbf{k}) = \alpha_i \tau_{ij}(\mathbf{k}, t) N_j^*(\mathbf{k}). \quad (2.18)$$

## 2.4. The EDQNM closure for transfer terms

In this section, a ‘triadic’ closure is applied to the equations governing the third-order spectral tensor $S_{ijn}(\mathbf{k}, \mathbf{p}, t)$ defined by (2.7), from which the term $\tau_{ij}(\mathbf{k}, t)$ in (2.6) is derived. From the counterpart of (2.1) in Fourier space and (2.7), its dynamics is obtained via

$$\left( \frac{\partial}{\partial t} + \nu(k^2 + p^2 + q^2) - \lambda_{lm} \left( k_l \frac{\partial}{\partial k_m} + p_l \frac{\partial}{\partial p_m} \right) \right) S_{ijn}(\mathbf{k}, \mathbf{p}, t) + M_{im}(\mathbf{q}) S_{mjn}(\mathbf{k}, \mathbf{p}, t) \\ + M_{jm}(\mathbf{k}) S_{imn}(\mathbf{k}, \mathbf{p}, t) + M_{nm}(\mathbf{p}) S_{ijm}(\mathbf{k}, \mathbf{p}, t) = T_{ijn}(\mathbf{k}, \mathbf{p}, t), \quad (2.19)$$

where $\mathbf{k} + \mathbf{p} + \mathbf{q} = \mathbf{0}$ and $T_{ijn}(\mathbf{k}, \mathbf{p}, t)$ is expressed in terms of a fourth-order spectral tensor

$$T_{ijn}(\mathbf{k}, \mathbf{p}, t) = P_{imp}(\mathbf{q}) \iiint S_{mpjn}(\mathbf{r}, \mathbf{k}, \mathbf{p}, t) \mathrm{d}^3\mathbf{r} + P_{jmp}(\mathbf{k}) \iiint S_{mpin}(\mathbf{r}, \mathbf{q}, \mathbf{p}, t) \mathrm{d}^3\mathbf{r} \\ + P_{nmp}(\mathbf{p}) \iiint S_{mpij}(\mathbf{r}, \mathbf{q}, \mathbf{k}, t) \mathrm{d}^3\mathbf{r}, \quad (2.20)$$

with

$$\langle \hat{u}_m(\mathbf{r}) \hat{u}_p(\mathbf{r}') \hat{u}_j(\mathbf{k}) \hat{u}_n(\mathbf{p}) \rangle = S_{mpjn}(\mathbf{r}, \mathbf{k}, \mathbf{p}) \delta(\mathbf{r} + \mathbf{r}' + \mathbf{k} + \mathbf{p}). \quad (2.21)$$

So far, the expressions that have been given for the governing equations for the second- and third-order spectral tensors are exact. They strictly reproduce the infinite hierarchy of moments up to $N = 3$, with equations for $N$-order moments having both linear closed terms and contributions from $N + 1$ moments being induced by basic nonlinearity. We now want to break this infinite hierarchy at the order $N = 3$. Equation (2.19) can be rewritten in the following form:

$$\left( \frac{\partial}{\partial t} + \nu(k^2 + p^2 + q^2) \right) S_{ijn}(\mathbf{k}, \mathbf{p}, t) = T_{ijn}(\mathbf{k}, \mathbf{p}, t) + L_{ijn}(\mathbf{k}, \mathbf{p}, t) = R_{ijn}(\mathbf{k}, \mathbf{p}, t), \quad (2.22)$$

where the tensor $R_{ijn}(\mathbf{k}, \mathbf{p}, t)$ gathers the linear operators induced by mean-velocity gradients through $L_{ijn}(\mathbf{k}, \mathbf{p}, t)$ and the fourth-order spectral tensor $T_{ijn}(\mathbf{k}, \mathbf{p}, t)$.
The problem then is to determine a closed expression for $R_{ijn}(\mathbf{k}, \mathbf{p}, t)$. The historical procedure, developed for homogeneous isotropic turbulence, starts from the quasinormal (QN) approximation (Millionschikov 1941; Proudman & Reid 1954), which states that the fluctuating velocity probability distributions are not too far from normal laws, in order to close the nonlinear contributions in $R_{ijn}(\mathbf{k}, \mathbf{p}, t)$. This assumption translates into vanishing fourth-order cumulants and can be written as

$$ \langle \hat{u}_m(\mathbf{r}) \hat{u}_p(\mathbf{r}') \hat{u}_j(\mathbf{k}) \hat{u}_n(\mathbf{p}) \rangle = \langle \hat{u}_m(\mathbf{r}) \hat{u}_p(\mathbf{r}') \rangle \langle \hat{u}_j(\mathbf{k}) \hat{u}_n(\mathbf{p}) \rangle + \langle \hat{u}_m(\mathbf{r}) \hat{u}_j(\mathbf{k}) \rangle \langle \hat{u}_p(\mathbf{r}') \hat{u}_n(\mathbf{p}) \rangle \\ + \langle \hat{u}_m(\mathbf{r}) \hat{u}_n(\mathbf{p}) \rangle \langle \hat{u}_p(\mathbf{r}') \hat{u}_j(\mathbf{k}) \rangle. \qquad (2.23) $$

Injecting (2.23) into (2.20) and using definitions (2.2) and (2.21) leads to the quasinormal contribution of the transfer term $T_{ijn}(\mathbf{k}, \mathbf{p}, t)$:

$$ T_{ijl}^{(QN)}(\mathbf{k}, \mathbf{p}, t) = 2(P_{imn}(\mathbf{q})\hat{R}_{mj}(\mathbf{k}, t)\hat{R}_{nl}(\mathbf{p}, t) + P_{jmn}(\mathbf{k})\hat{R}_{ml}(\mathbf{p}, t)\hat{R}_{ni}(\mathbf{q}, t) \\ + P_{lmn}(\mathbf{p})\hat{R}_{mi}(\mathbf{q}, t)\hat{R}_{nj}(\mathbf{k}, t)), \qquad (2.24) $$

with $P_{imn}(\mathbf{k}) = 1/2(k_m P_{in}(\mathbf{k}) + k_n P_{im}(\mathbf{k}))$; the projector $P_{ij}(\mathbf{k})$ is defined in (1.1). However, it was shown by O'Brien & Francis (1963) and Ogura (1963) that the purely quasinormal approximation fails in decaying isotropic turbulence for long elapsed times, yielding negative energy spectra at small $k$. Orszag (1970) showed that the improper treatment of relaxation effects in the purely quasinormal approximation lies at the origin of this lack of realizability. Consequently, he introduced an eddy viscosity, or eddy-damping (ED), term in the governing equation for third-order correlations. Without any additional assumption, $R_{ijn}(\mathbf{k}, \mathbf{p}, t)$ from (2.22) can be written as

$$ R_{ijn}(\mathbf{k}, \mathbf{p}, t) = T_{ijn}^{(QN)}(\mathbf{k}, \mathbf{p}, t) + T_{ijn}^{(IV)}(\mathbf{k}, \mathbf{p}, t) + L_{ijn}(\mathbf{k}, \mathbf{p}, t). \qquad (2.25) $$

In this equation, the only unknown, and unclosed, term is $T_{ijn}^{(IV)}(\mathbf{k}, \mathbf{p}, t)$, which represents the contribution from fourth-order cumulants. A natural extension of Orszag's introduction of eddy damping is to write

$$ T_{ijn}^{(IV)}(\mathbf{k}, \mathbf{p}, t) = -(\eta(\mathbf{k}, t) + \eta(\mathbf{p}, t) + \eta(\mathbf{q}, t))S_{ijn}(\mathbf{k}, \mathbf{p}, t), \qquad (2.26) $$

while keeping in mind that fourth-order cumulants may act as a linear relaxation of triple correlations, which will reinforce the dissipative operator in (2.22) when added to the purely viscous terms on its left-hand side. The simple relationship (2.26), which is isotropic with a single eddy-damping scalar term, can be discussed with respect to two-time triadic theories as follows. In isotropic turbulence, the EDQNM is possibly derived from the basic DIA (Kraichnan 1959). However, such a derivation of $\eta(\mathbf{k}, t)$ from the two-time response tensor is unsatisfactory and yields an inappropriate 'sweeping' time scale, so that more complicated Lagrangian theories must be called into play, such as Lagrangian-history DIA (LHDIA) (Kraichnan & Herring 1978) and Lagrangian renormalized approximation (LRA) (Kaneda 1981). The identification of a possible anisotropic shear-dependent alternative to (2.26) from these theories is a difficult task because they are essentially non-Markovian, and reduction to a single-time expression was only performed at very weak anisotropy (see, e.g., Yoshida, Ishihara & Kaneda (2003) for LRA). Finally, the test-field model (Kraichnan 1972), as an almost Markovian theory, is fully consistent with a linear relationship such as (2.26), but with a possibly more complex tensorial structure,
too intricate here to be applied to a practical model. Concerning the choice of the relevant time scale to evaluate $\eta(k, t)$, many proposals exist, with 'sweeping', 'straining', linear and nonlinear time scales, see, e.g., Schiestel (1987), Kim & Chen (1989), Rubinstein (1996) and Zhou (2010). In the presence of mean shear, one can introduce a scalar mean straining time scale based on $\sqrt{\lambda_{ij}\lambda_{ij}}$, with $\lambda_{ij}$ the mean-velocity gradient. However, the effects of mean shear/strain are explicit and closed (in terms of third-order correlations) in (2.22) via the operator $L_{ijn}(k, p, t)$. The contribution of the latter tensor was considered and analytically solved in previous EDQNM versions, which are mentioned in the introduction (such as, e.g., the model of Cambon & Jacquin (1989) for rotating flows). However, this approach will not be considered here for three reasons: (i) it was motivated and really needed for the study of anisotropic flows without energy production and in the presence of interacting dispersive waves leading the dynamics of triple correlations, which is not the general case here; (ii) it renders the tensorial structure of the EDQNM model much more complicated, via a threefold product of Green's functions, and explicitly dependent on the type of mean shear, preventing easy further projection on spherical harmonics; (iii) it is not correct when the direct (linear) effect of the mean shear/strain yields exponential growth, with a lack of convergence of the time integral (Cambon & Scott 1999). The eddy-damping coefficient $\eta(k, t)$ is finally chosen as

$$ \eta(k, t) = A \sqrt{\int_{0}^{k} p^2 E(p, t) dp}, \qquad (2.27) $$

following Pouquet et al. (1975), which is an improved variant of Orszag (1970)'s proposal. The constant is fixed at $A = 0.36$ to recover a well-admitted value of the Kolmogorov constant (André & Lesieur 1977). The same equation is used in Weinstock (2013). Equation (2.22) can be now integrated using (2.25), (2.24) and (2.26), (2.27), and neglecting the contribution of $L_{ijn}(k, p, t)$. The solution thus obtained involves time integrals that are further simplified by a Markovianization procedure, which amounts to truncating the proper time memory of triple correlations. Even in the isotropic case, the eddy damping allows a posteriori Markovianization to be justified, arguing that the time variation is much more rapid in the exponential eddy-damping term than in the second-order moments present in the quasinormal term. Further details can be found in Cambon et al. (1981) and Sagaut & Cambon (2008). In the end, this anisotropic EDQNM closure for the tensor $\tau_{ij}(k, t)$ defined by (2.6) amounts to

$$ \tau_{ij}(k, t) = k_l \iiint \theta_{kpq} T^{(QN)}(k, p, t) d^3p, \qquad (2.28) $$

where

$$ \theta_{kpq} = \frac{1 - e^{-\mu_{kpq}t}}{\mu_{kpq}}, \quad \mu_{kpq} = \nu(k^2 + p^2 + q^2) + \eta(k, t) + \eta(p, t) + \eta(q, t). \qquad (2.29a,b) $$

The expression (2.28) of $\tau_{ij}(k, t)$ obtained from the EDQNM approximation is injected into (2.15) and (2.16) in order to obtain closed-form expressions of the transfer terms $T^{(\mathcal{E})}(k, t)$ and $T^{(Z)}(k, t)$:

$$ T^{(\mathcal{E})}(k, t) = \iint\int \theta_{kpq} 2kp [(\mathcal{E}'' + \text{Re} X'')[(xy + z^3)(\mathcal{E}' - \mathcal{E}) - z(1 - z^2)(\text{Re} X' - \text{Re} X)] + \text{Im} X''(1 - z^2)(x \text{Im} X - y \text{Im} X')] d^3p, \qquad (2.30) $$
$$ T^{(Z)}(\mathbf{k}, t) = \iiint \theta_{kpq} 2kpe^{-2i\lambda} [(\mathcal{E}'' + \operatorname{Re} X'') [(xy + z^3)(\operatorname{Re} X' - X) - z(1 - z^2)](\mathcal{E}' - \mathcal{E}) \\ + i(y^2 - z^2)\operatorname{Im} X'] + i \operatorname{Im} X''(1 - z^2)[x(\mathcal{E} + X) - iy \operatorname{Im} X')] d^3\mathbf{p}, \quad (2.31) $$

with $\mathcal{E} = \mathcal{E}(\mathbf{k}, t)$, $\mathcal{E}' = \mathcal{E}(\mathbf{p}, t)$, $\mathcal{E}'' = \mathcal{E}(\mathbf{q}, t)$, $X = Z(\mathbf{k}, t)e^{2i\lambda}$, $X' = Z(\mathbf{p}, t)e^{2i\lambda'}$, $X'' = Z(\mathbf{q}, t)e^{2i\lambda''}$, $x = -\mathbf{p} \cdot \mathbf{q}/(pq)$, $y = -\mathbf{k} \cdot \mathbf{q}/(kq)$ and $z = -\mathbf{k} \cdot \mathbf{p}/(kp)$. Here, $\lambda, \lambda'$ and $\lambda''$ are angles that characterize the rotation of the plane of the triad around $\mathbf{k}, \mathbf{p}$ and $\mathbf{q}$ respectively; $x, y$ and $z$ refer to the cosines of the internal angles of the triangle formed by the triad. The above expressions for the transfer terms $T^{(\mathcal{E})}(\mathbf{k}, t)$ and $T^{(Z)}(\mathbf{k}, t)$ are also derived in Cambon et al. (1997) in a different way. Similarly, one can deduce a closed-form expression of the transfer term $T^{(RTI)}(\mathbf{k}, t)$ defined by (2.18), which generates the 'slow' component of the pressure-strain rate tensor:

$$ T^{(RTI)}(\mathbf{k}, t) = \iiint \theta_{kpq} 2e^{-i\lambda} p(xy + z) \sqrt{1 - z^2}(\mathcal{E}'' + \operatorname{Re} X'') \\ \times [(\mathcal{E} + X)(zk - qx) - k(z(\mathcal{E}' + \operatorname{Re} X') - i \operatorname{Im} X')] d^3\mathbf{p}. \quad (2.32) $$

### 3. Dynamical equations for spherically averaged descriptors

Strictly speaking, the Lin equations (2.12) and (2.13) with the closed-form expressions (2.30) and (2.31) of the transfer terms can be solved. However, important practical difficulties arise from the $\mathbf{k}$ dependence of the second-order spectral tensor $\hat{R}_{ij}(\mathbf{k}, t)$, or equivalently from those of $\mathcal{E}(\mathbf{k}, t)$ and $Z(\mathbf{k}, t)$. In order to circumvent these difficulties, one solution is to integrate analytically the closed Lin equations over a sphere of radius $\mathbf{k}$. This analytical integration requires a representation of the tensor $\hat{R}_{ij}(\mathbf{k}, t)$, which is described in § 3.1. This representation involves spherically averaged descriptors whose governing equations, which are the main result of this paper, are derived in the remainder of this section. Limitations of the model in terms of anisotropy intensity are quantified in § 3.5.

#### 3.1. Representation of the second-order spectral tensor in terms of spherically averaged descriptors

Here, we use for $\hat{R}_{ij}(\mathbf{k}, t)$ the representation proposed by Cambon & Rubinstein (2006). This representation involves spherically averaged descriptors and is obtained by treating directionality and polarization anisotropy separately. It is written as

$$ \hat{R}_{ij}(\mathbf{k}, t) = \underbrace{\frac{E(k, t)}{4\pi k^2} P_{ij}(\mathbf{k})}_{\hat{R}_{ij}^{(iso)}(\mathbf{k}, t)} - \underbrace{\frac{15 E(k, t)}{4\pi k^2} P_{ij}(\mathbf{k}) H_{pq}^{(dir)}(\mathbf{k}, t) \alpha_p \alpha_q}_{\hat{R}_{ij}^{(dir)}(\mathbf{k}, t)} + \underbrace{5 \frac{E(k, t)}{4\pi k^2} \left( P_{ip}(\mathbf{k}) P_{jq}(\mathbf{k}) + \frac{1}{2} P_{ij}(\mathbf{k}) \alpha_p \alpha_q \right) H_{pq}^{(pol)}(\mathbf{k}, t)}_{\hat{R}_{ij}^{(pol)}(\mathbf{k}, t)} \quad (3.1) $$

or equivalently

$$ \mathcal{E}(\mathbf{k}, t) = \frac{E(k, t)}{4\pi k^2} (1 - 15 H_{ij}^{(dir)}(k, t) \alpha_i \alpha_j), \quad Z(\mathbf{k}, t) = \frac{5 E(k, t)}{2 4\pi k^2} H_{ij}^{(pol)}(k, t) N_i^*(\mathbf{k}) N_j^*(\mathbf{k}), \quad (3.2a,b) $$
where the tensor $P_{ij}(k)$ is defined in (1.1) and $\alpha_i = k_i/k$. The tensors $\hat{R}_{ij}^{(iso)}(k, t)$, $\hat{R}_{ij}^{(dir)}(k, t)$ and $\hat{R}_{ij}^{(pol)}(k, t)$ identify the isotropic, directional and polarization parts of $\hat{R}_{ij}(k, t)$ respectively. The representation (3.1) is constructed so that the trace-free tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$, which depend only on $k$, measure the directional anisotropy and polarization anisotropy according to (1.5). Injecting the representation of $\mathcal{E}(k, t)$ and $Z(k, t)$ (3.2) into (2.12) and (2.13) and (2.30) and (2.31) allows us to integrate the latter analytically over a sphere of radius $k$ and to derive a system of governing equations in terms of the spherically averaged descriptors $E(k, t)$, $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$. In view of (3.1), (3.2) and (1.4), the latter completely determine the second-order spectral tensor $\hat{R}_{ij}(k, t)$ and its spherically integrated counterpart $\varphi_{ij}(k, t)$, however restricted to moderate anisotropy. This point is discussed in § 3.5.

## 3.2. Dynamics, final closure

The system of governing equations for the spherically averaged descriptors $E(k, t)$, $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$, which is the main result of this article, is written as follows:

$$ \left( \frac{\partial}{\partial t} + 2\nu k^2 \right) E(k, t) = S^L(k, t) + T(k, t), \qquad (3.3) $$

$$ \left( \frac{\partial}{\partial t} + 2\nu k^2 \right) E(k, t) H_{ij}^{(dir)}(k, t) = S_{ij}^{L(d)}(k, t) + S_{ij}^{NL(d)}(k, t), \qquad (3.4) $$

$$ \left( \frac{\partial}{\partial t} + 2\nu k^2 \right) E(k, t) H_{ij}^{(pol)}(k, t) = S_{ij}^{L(p)}(k, t) + S_{ij}^{NL(p)}(k, t), \qquad (3.5) $$

with

$$ 2 \left( \frac{\delta_{ij}}{3} T(k, t) + S_{ij}^{NL(d)}(k, t) + S_{ij}^{NL(p)}(k, t) \right) = S_{ij}(k, t) + P_{ij}(k, t). \qquad (3.6) $$

The tensors $S^L(k, t)$, $S_{ij}^{L(d)}(k, t)$ and $S_{ij}^{L(p)}(k, t)$ account for the linear terms corresponding to the interactions with the mean flow, whereas $T(k, t)$, $S_{ij}^{NL(d)}(k, t)$ and $S_{ij}^{NL(p)}(k, t)$ correspond to nonlinear transfer terms. The tensor $P_{ij}(k, t)$ is the spherically integrated spectral counterpart of the slow pressure–strain rate tensor, to which an RTI is conventionally attributed. The tensor $S_{ij}(k, t)$ corresponds to a ‘true’ transfer tensor with $\int_0^\infty S_{ij}(k, t) dk = 0 \forall t$. Since the tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$ are symmetric and trace-free, the system (3.3)–(3.5) forms a set of 11 independent equations.

## 3.3. Closure for the linear terms induced by mean-gradient effects

These terms are exact and linear in the equation governing the full spectral tensor, or equivalently $\mathcal{E}(k, t)$ and $Z(k, t)$. Taken independently with zero contribution from third-order correlations, they reflect the RDT limit for the evolution of two-point second-order velocity correlations. In order to obtain the corresponding spherically averaged terms $S^L(k, t)$, $S_{ij}^{L(d)}(k, t)$ and $S_{ij}^{L(p)}(k, t)$, one has to analytically solve the spherical averaging of tensorial products of vectors $\alpha = k/k$. This is done following
the procedure described in Cambon et al. (1981). Performing the spherical integration of the linear terms in (2.3) or that of the linear terms in (2.12), (2.13) with the representation (3.1) or (3.2) leads to

$$S^L(k, t) = -2S_{lm} \frac{\partial}{\partial k} (kE H_{lm}^{(dir)} - 2E S_{lm} (H_{lm}^{(dir)} + H_{lm}^{(pol)}), \quad (3.7)$$

$$S_{ij}^{L(dir)}(k, t) = \frac{2}{15} S_{ij} E - \frac{2}{7} E \left( S_{jl} H_{il}^{(pol)} + S_{il} H_{jl}^{(pol)} - \frac{2}{3} S_{lm} H_{lm}^{(pol)} \delta_{ij} \right) \\ + \frac{2}{7} \left( S_{il} \frac{\partial}{\partial k} (k E H_{lj}^{(dir)}) + S_{lj} \frac{\partial}{\partial k} (k E H_{li}^{(dir)}) - \frac{2}{3} S_{lm} \frac{\partial}{\partial k} (k E H_{lm}^{(dir)}) \delta_{ij} \right) \\ - \frac{1}{7} E \left( S_{jl} H_{li}^{(dir)} + S_{il} H_{lj}^{(dir)} - \frac{2}{3} S_{lm} H_{lm}^{(dir)} \delta_{ij} \right) + E \left( A_{jn} H_{ni}^{(dir)} + A_{in} H_{jn}^{(dir)} \right) \\ - \frac{1}{15} S_{ij} \frac{\partial}{\partial k} (kE), \quad (3.8)$$

$$S_{ij}^{L pol}(k, t) = -\frac{2}{5} E S_{ij} - \frac{12}{7} E \left( S_{lj} H_{li}^{(dir)} + S_{il} H_{lj}^{(dir)} - \frac{2}{3} S_{lm} H_{lm}^{(dir)} \delta_{ij} \right) \\ - \frac{2}{7} \left( S_{jl} \frac{\partial}{\partial k} (k E H_{il}^{(pol)}) + S_{il} \frac{\partial}{\partial k} (k E H_{lj}^{(pol)}) - \frac{2}{3} S_{lm} \frac{\partial}{\partial k} (k E H_{lm}^{(pol)}) \delta_{ij} \right) \\ + \frac{1}{7} E \left( S_{il} H_{lj}^{(pol)} + S_{jl} H_{li}^{(pol)} - \frac{2}{3} S_{lm} H_{lm}^{(pol)} \delta_{ij} \right) \\ - \frac{1}{3} E \left( A_{il} H_{lj}^{(pol)} + A_{jl} H_{li}^{(pol)} \right), \quad (3.9)$$

with $E = E(k, t)$, $H_{ij}^{(dir)} = H_{ij}^{(dir)}(k, t)$, $H_{ij}^{(pol)} = H_{ij}^{(pol)}(k, t)$, $S_{ij} = (\lambda_{ij} + \lambda_{ji})/2$, $A_{ij} = (\lambda_{ij} - \lambda_{ji})/2$; $\lambda_{ij}$ refers to the mean-flow gradient as previously stated.

### 3.4. Closure for the terms mediated by third-order correlations

The transfer terms $T(k, t)$, $S_{ij}^{NL(dir)}(k, t)$ and $S_{ij}^{NL(pol)}(k, t)$ are obtained from the spherical integration of the expressions of the transfer terms $T(\mathcal{E})(\mathbf{k}, t)$ and $T^{(Z)}(\mathbf{k}, t)$ closed by the EDQNM procedure (2.30) and (2.31) and using the representation (3.2) for $\mathcal{E}(\mathbf{k}, t)$ and $Z(\mathbf{k}, t)$. It is consistent to retain only linear contributions from the tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$ in the terms present on the right-hand sides of the system (3.3)–(3.5) in view of the discussion in § 3.5. Their low contribution is checked in § 4.7. In anisotropic triadic closure, the new difficulty is to solve the integral over the orientation of the plane of the triad (Cambon et al. 1997; Sagaut & Cambon 2008), which is performed analytically. The final results are

$$T(k, t) = \iint_{\Delta_k} \theta_{kpq} 16\pi^2 p^2 k^2 q (xy + z^3) \mathcal{E}_0'' (\mathcal{E}_0' - \mathcal{E}_0) dp dq, \quad (3.10)$$

$$S_{ij}^{NL(dir)}(k, t) = \iint_{\Delta_k} \theta_{kpq} 4\pi^2 p^2 k^2 q \mathcal{E}_0'' [(y^2 - 1)(xy + z^3)(\mathcal{E}_0' - \mathcal{E}_0)H_{ij}^{(pol)''} \\ + z(1 - z^2)^2 \mathcal{E}_0' H_{ij}^{(pol)'}] dp dq$$
$$ + \iint_{\Delta_k} \theta_{kpq} 8\pi^2 p^2 k^2 q (xy + z^3) \mathcal{E}_0'' [(3y^2 - 1)(\mathcal{E}_0' - \mathcal{E}_0)H_{ij}^{(dir)''} \\ + (3z^2 - 1)\mathcal{E}_0' H_{ij}^{(dir)'} - 2\mathcal{E}_0 H_{ij}^{(dir)}] \mathrm{d}p \mathrm{d}q, \quad (3.11) $$

$$ \begin{aligned} S_{ij}^{\mathrm{NL}(\mathrm{pol})}(k, t) = & \iint_{\Delta_k} \theta_{kpq} 4\pi^2 p^2 k^2 q \mathcal{E}_0'' [(xy + z^3)((1+z^2)\mathcal{E}_0' H_{ij}^{(\mathrm{pol})'} - 4\mathcal{E}_0 H_{ij}^{(\mathrm{pol})}) \\ & + z(z^2 - 1)(1+y^2)(\mathcal{E}_0' - \mathcal{E}_0)H_{ij}^{(\mathrm{pol})''} + 2z(z^2 - y^2)\mathcal{E}_0' H_{ij}^{(\mathrm{pol})'} \\ & + 2yx(z^2 - 1)\mathcal{E}_0 H_{ij}^{(\mathrm{pol})''}] \mathrm{d}p \mathrm{d}q \\ & + \iint_{\Delta_k} \theta_{kpq} 24\pi^2 p^2 k^2 qz(z^2 - 1)\mathcal{E}_0'' [(y^2 - 1)(\mathcal{E}_0' - \mathcal{E}_0)H_{ij}^{(\mathrm{dir})''} \\ & + (z^2 - 1)\mathcal{E}_0' H_{ij}^{(\mathrm{dir})'}] \mathrm{d}p \mathrm{d}q, \end{aligned} \quad (3.12) $$

$$ P_{ij}(k, t) = \iint_{\Delta_k} \theta_{kpq} 16\pi^2 p^2 k^2 q (yz + x) \mathcal{E}_0'' [\mathcal{E}_0'(y(z^2 - y^2)(6H_{ij}^{(dir)''} + H_{ij}^{(pol)''}) \\ -(xz + y)H_{ij}^{(pol)''}) - y(z^2 - x^2)\mathcal{E}_0(6H_{ij}^{(dir)''} + H_{ij}^{(pol)''})] \mathrm{d}p \mathrm{d}q, \quad (3.13) $$

with $\mathcal{E}_0 = (E(k, t))/(4\pi k^2)$, $\mathcal{E}_0' = (E(p, t))/(4\pi p^2)$, $\mathcal{E}_0'' = (E(q, t))/(4\pi q^2)$, $H_{ij}^{(0)} = H_{ij}^{(0)}(k, t)$, $H_{ij}^{(0)} = H_{ij}^{(0)}(p, t)$ and $H_{ij}^{(0)'} = H_{ij}^{(0)}(q, t)$, where $H_{ij}^{(0)}$ may refer to either $H_{ij}^{(dir)}$ or $H_{ij}^{(pol)}$. The integrals over $p$ and $q$ are performed over the domain $\Delta_k$, so that $k, p$ and $q$ are the lengths of the sides of the triangle formed by $k, p$ and $q$. The expression of the 'true' transfer $S_{ij}(k, t)$ can be deduced from (3.6) and (3.10)–(3.13).

### 3.5. Properties of the model

Spherical averaging of the Lin equations (2.12)–(2.13) allows us to obtain a model for anisotropic turbulence that can be used to calculate anisotropic turbulent flows at both very high and low Reynolds numbers, with good resolution of both large and small scales and over very long evolution times. Besides, (3.3)–(3.5) remain consistent with the decomposition of $\hat{\mathbf{R}}(k, t)$ in terms of directional anisotropy and polarization anisotropy. However, spherical integration also has a drawback since it implies a loss of information, which restricts the present model to moderately anisotropic flows. Looking at expansions of the scalars $\mathcal{E}(k, t)$ and $Z(k, t)$ in terms of powers of $\alpha = k/k$, the first approach by Cambon & Rubinstein (2006) can be translated into the following:

$$ \mathcal{E}(k, t) = \frac{E(k, t)}{4\pi k^2} \left( 1 + U_{ij}^{(dir)2}(k, t)\alpha_i\alpha_j + U_{ijmn}^{(dir)4}(k, t)\alpha_i\alpha_j\alpha_m\alpha_n + \dots \right), \quad (3.14) $$

$$ Z(k, t) = \frac{1}{2} \frac{E(k, t)}{4\pi k^2} \left( U_{ij}^{(pol)2}(k, t) + U_{ijm}^{(pol)3}(k, t)\alpha_m + U_{ijmn}^{(pol)4}(k, t)\alpha_m\alpha_n + \dots \right) N_i^*(k)N_j^*(k). \quad (3.15) $$

The latter form is consistent with both scalar (for (3.14)) and tensor (for (3.15)) spherical harmonic expansion generated by the rotation group SO³ decomposition, in agreement with Rubinstein, Kurien & Cambon (2015). With the identification

$$ U_{ij}^{(dir)2}(k, t) = -15H_{ij}^{(dir)}(k, t), \quad U_{ij}^{(pol)2}(k, t) = 5H_{ij}^{(pol)}(k, t), \quad (3.16a,b) $$
the representation (3.2) is interpreted as the first-order truncation of expansions (3.14) and (3.15). The truncation of (3.14) and (3.15) lies at the origin of the loss of information, which limits the degree of anisotropy that can be investigated with the present model. The degree of anisotropy permitted by the representation (3.1) can be derived from realizability conditions. One may use condition (2.11). In order to derive a simple condition in terms of the tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$, one may also consider the weaker condition $\mathcal{E}(k, t) \ge 0 \forall k, t$, which already proves to be very restrictive. In view of (3.2), this condition is equivalent to

$$ \max_i \Lambda_i(H^{(dir)}(k, t)) \le \frac{1}{15} \forall k, t, \qquad (3.17) $$

where $\Lambda_i(H^{(dir)}(k, t))$ refers to the eigenvalues of $H_{ij}^{(dir)}(k, t)$. Condition (3.17) can help to quantify the upper boundary of anisotropy intensity to ensure that the present model represents correctly the corresponding turbulent flow. In §4.7, criterion (3.17) is computed for the different flow configurations studied in this paper in order to check whether the representation (3.2) and the governing equations (3.3)–(3.5) can describe the corresponding turbulent flow. Since the representation (3.2) is restricted to the description of moderate anisotropy, we discard quadratic contributions from the tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$ which appear when the representation (3.2) is injected into (2.30)–(2.32). It is also important to keep in mind that the nonlinear transfer terms of the system (3.3)–(3.5) are not relevant when the mean flow acts only on triple correlations such as in the case for solid-body rotation, but are restricted to mean flows leading to linear production effects in the equations for second-order correlations. This is due to the specific EDQNM version used in §2.4. In the end, the present model is suited for the study of turbulent flows where the anisotropy is moderate and where linear effects induced by mean-velocity gradients play a negligible role in the dynamics of triple correlations compared with those directly induced in the equations for second-order correlations. Since only flows dominated by production effects induced by a mean strain or shear are considered in §4, we expect the EDQNM closure of §2.4 to be valid in these configurations. Besides, the linear contributions in (3.4) and (3.5) originate from the exact linear terms of (2.12) and (2.13). Accordingly, it is inferred that the main source of possible discrepancies between the predictions obtained with (3.3)–(3.5) and experimental/direct numerical simulation (DNS) results is due to the loss of information following the first-order truncation of (3.14) and (3.15).

# 4. First applications and results: RTI and effect of mean shear

## 4.1. Physical and numerical set-up

The predictions of the model presented in §3 are first compared with data from the experiments of Gence & Mathieu (1979, 1980). In these studies, two successive plane strains with different orientations are applied to grid-generated turbulence, and the RTI of the turbulence thus obtained is investigated in the latter experiment. We then consider the experiment of Chen, Meneveau & Katz (2006), where turbulence is subjected to a non-stationary straining–relaxation–destraining cycle. Both linear and nonlinear phenomena come into play in the response of turbulence to this particular straining. The case of turbulence subjected to an axisymmetric expansion or contraction is then investigated. The set-up of the numerical simulations performed with the present model is similar to that of the DNS of Zusi & Perot (2014). In particular, we want to check the ability of the present model to capture the initial increase of anisotropy when the strain is removed, only in the case of axisymmetric
expansion, which is observed in the above-mentioned study. The next test case deals with an initially isotropic flow which is suddenly subjected to a mean shear and is then released. A comparison with the use of RDT to initialize the anisotropy of the flow will allow us to validate the linear part of the present model in the case of homogeneous shear turbulence. This case will also allow us to confirm the consistency between the present model and the permanence of large eddies (PLE). Finally, we consider the case where the mean shear is maintained during the evolution of the flow. Homogeneous shear turbulence has been studied experimentally (Tavoularis & Corrsin 1981; Tavoularis 1985; Rohr et al. 1988; Tavoularis & Karnik 1989; Shen & Warhaft 2000), numerically (Rogers, Moin & Reynolds 1986; Ishihara, Yoshida & Kaneda 2002; Brethouwer 2005; Isaza & Collins 2009) and theoretically (Lumley 1967; Leslie 1973; Yoshida et al. 2003; Weinstock 2013), and the results obtained with the system (3.3)–(3.5) will be compared with the predictions of these different studies. In the last subsection, it is checked that condition (3.17) is verified for all of the flow configurations discussed above.

As briefly mentioned above, concerning the study of the RTI of shear-released turbulence, the anisotropy of the flow can be initialized in two ways. A first possibility consists in deriving analytically an initial anisotropic condition thanks to RDT. The other option is to start with an isotropic field and to generate the anisotropy with the linear terms of (3.3)–(3.5) that account for the interactions with a mean flow. For $St \ll 1$, where $S$ refers to the characteristic shear rate, these two ways of introducing anisotropy must be equivalent. This will help to validate the model developed in this paper. The complete RDT solution for homogeneous shear turbulence can be found in Townsend (1976), Piquet (2001) and Sagaut & Cambon (2008). The energy spectrum used for the RDT calculations, or directly to initialize the simulations, is written as

$$E(k) = C\epsilon^{2/3}k^{-(5/3)}f(kL)g(k\eta), \quad (4.1)$$

with

$$f(x) = \left( \frac{x}{\left( x^{1.5} + 1.5 - \frac{\sigma}{4} \right)^{1/1.5}} \right)^{(5/3)+\sigma}, \quad g(x) = \exp(-5.2((x^4 + 0.4^4)^{1/4} - 0.4)). \quad (4.2a,b)$$

The functions defined by (4.2) have been proposed by Pope (2000) and Meyers & Meneveau (2008) respectively. Here, $\epsilon$ refers to the dissipation rate, $L$ to the integral length scale, $\eta$ to the Kolmogorov scale, and $\sigma$ is the slope at large scales of the energy spectrum ($E(k \to 0) \propto k^\sigma$). In all simulations we set $\sigma = 2$, which corresponds to Saffman turbulence. The integral scale $L$ and the Reynolds number are prescribed at the beginning of the simulation. As a consequence, $\eta$ is initialized thanks to the relation $\eta = ((3/20)Re_\lambda^2)^{-(3/4)}L$, where $Re_\lambda$ is the Reynolds number based on the Taylor microscale $\lambda = \sqrt{10\kappa v/\epsilon}$. The constant $C$ in (4.1) is adjusted in order to impose the initial value of the turbulent kinetic energy $\kappa$ so that $\kappa(t=0)/\kappa_0 = 1$, where $\kappa_0$ is a reference energy. The spectral mesh used to perform the simulations covers a very broad wavenumber range in the spectral domain. The smallest wavenumber of the mesh $k_0$ (which roughly corresponds to the largest resolved scale) is chosen so that $k_0 = 10^{-8}k_L(t=0)$, where $k_L = 1/L$ is the wavenumber associated with the integral length scale $L$. The largest resolved wavenumber $k_N$ is chosen such that a
resolution of at least one decade at small scales, with respect to the initial value of the Kolmogorov scale $\eta(t=0)$, is granted. The relation $k_N = 10k_\eta(0)$ is thus imposed, where $k_\eta = 1/\eta$. The total number of elements $N+1$ is recovered so that $k_N = r^N k_0$, where $r$ represents the constant aspect ratio between contiguous elements of the mesh. In the present work, $r=1.122$, which means that each decade in the spectral space is discretized by 20 mesh elements. The temporal integration is made by a forward Euler scheme, and the Courant–Friedrichs–Lewy condition is based on the Kolmogorov scale (Lesieur 2008).

The anisotropy at the spectral level is evaluated with the deviatoric tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$ introduced in § 3.1, and $H_{ij}(k, t) = H_{ij}^{(dir)}(k, t) + H_{ij}^{(pol)}(k, t)$. The global anisotropy is quantified via the tensors $b_{ij}^{(dir)}(t)$, $b_{ij}^{(pol)}(t)$ and $b_{ij}(t) = b_{ij}^{(dir)}(t) + b_{ij}^{(pol)}(t)$, which originate from the $(\mathcal{E}, Z)$ decomposition of the Reynolds stress tensor:

$$ \langle u_i u_j \rangle(t) = 2K_{ij}(t) = 2\mathcal{H}(t) \left( \frac{\delta_{ij}}{3} + b_{ij}^{(dir)}(t) + b_{ij}^{(pol)}(t) \right), \quad \mathcal{H}(t) = K_{ii}(t). \quad (4.3) $$

The tensors $b_{ij}^{(dir)}(t)$ and $b_{ij}^{(pol)}(t)$ are calculated from the tensors $H_{ij}^{(dir)}(k, t)$ and $H_{ij}^{(pol)}(k, t)$ according to

$$ b_{ij}^{(dir)}(t) = \int_{0}^{\infty} E(k, t) H_{ij}^{(dir)}(k, t) dk / \mathcal{H}(t), \quad b_{ij}^{(pol)}(t) = \int_{0}^{\infty} E(k, t) H_{ij}^{(pol)}(k, t) dk / \mathcal{H}(t). \quad (4.4a,b) $$

We define the invariants $II(t)$ and $III(t)$ as

$$ II(t) = b_{ij}(t)b_{ji}(t), \quad III(t) = b_{ik}(t)b_{kj}(t)b_{ij}(t). \quad (4.5a,b) $$

In order to quantitatively characterize the RTI process, we introduce the ratio $\rho(t)$ between a characteristic time of turbulence decay and a characteristic time of RTI, namely

$$ \rho(t) = \frac{\mathcal{H}(t) \left( \frac{\mathrm{d}\mathcal{H}}{\mathrm{d}t}(t) \right)}{II(t) \left( \frac{\mathrm{d}II}{\mathrm{d}t}(t) \right)}. \quad (4.6) $$

For the comparison with the experiment of Chen et al. (2006), we use a two-component surrogate of the anisotropic tensor $b_{ij}(t)$, defined as

$$ \tilde{b}_{ij}(t) = \frac{K_{ij}(t)/K_{ij}(t_0)}{K_{11}(t)/K_{11}(t_0) + K_{22}(t)/K_{22}(t_0)} - \frac{1}{2}\delta_{ij} \quad (i, j = 1, 2), \quad (4.7) $$

where $t_0$ refers to the time at which strain starts. The budget terms of the governing equation for the spherically integrated second-order spectral tensor $\varphi_{ij}(k, t)$ defined by (1.4),

$$ \frac{\partial}{\partial t} \varphi_{ij}(k, t) = -2\nu k^2 \varphi_{ij}(k, t) + P_{ij}(k, t) + S_{ij}(k, t) + L_{ij}(k, t), \quad (4.8) $$

are also considered. The tensor $P_{ij}(k, t)$ is the spherically integrated spectral counterpart of the slow pressure–strain rate tensor and $S_{ij}(k, t)$ is a true transfer
tensor with zero integral over $k$, the expressions of these two tensors are given in § 3.
The term $L_{ij}(k, t)$, defined as

$$L_{ij}(k, t) = 2 \left( \frac{\delta_{ij}}{3} S^L(k, t) + S_{ij}^{L(dir)}(k, t) + S_{ij}^{L pol}(k, t) \right), \quad (4.9)$$

includes all contributions of the linear terms of (3.3)–(3.5). The results reported in the figures are often scaled with an initial turbulent characteristic time $\tau_0$, defined as $\tau_0 = \mathcal{K}(0)/\epsilon(0)$. The kinetic energy spectrum $E(k, t)$ is scaled by $k_{max}$ and $E_{max}$, defined by

$$\max_{k} E(k, t) = E_{max}(t) = E(k_{max}(t), t). \quad (4.10)$$

A similar scaling is used for the cross-correlation spectrum $\varphi_{13}(k, t)$.

## 4.2. Comparison with the experiments of Gence and Mathieu

In the experiments of Gence & Mathieu (1979, 1980), a plane strain is first applied to quasi-isotropic grid turbulence. A second strain is then applied, whose principal axes have been rotated by an angle $\alpha$ in the plane of the first strain. The mean-velocity gradients corresponding to the first and second strains are respectively

$$\lambda = \begin{pmatrix} 0 & 0 & 0 \\ 0 & S & 0 \\ 0 & 0 & -S \end{pmatrix}, \quad \lambda = \begin{pmatrix} 0 & 0 & 0 \\ 0 & S \cos(2\alpha) & -S \sin(2\alpha) \\ 0 & -S \sin(2\alpha) & -S \cos(2\alpha) \end{pmatrix}. \quad (4.11a,b)$$

In Gence & Mathieu (1980), the original experimental device is extended in such a way that the turbulence can develop downstream without a mean-velocity gradient. For these experiments, $S \approx 2.9\tau_0^{-1}$ with $\tau_0 = \mathcal{K}(0)/\epsilon(0)$, where the origin corresponds to the entrance of the distorting duct. The simulations are initialized with $Re_\lambda = 60$. No detailed spectral information is available in these works, and we assume that the turbulence is weakly axisymmetrically dilated at the entrance of the distorting duct. The initial condition is obtained from RDT and the degree of anisotropy is adjusted so that it coincides with the first measured values of the invariant $II$ defined by (4.5). Experimental data for the downstream evolution of the invariant $II$ are reported in figure 1 along with numerical results obtained with the system of governing equations (3.3)–(3.5). This figure shows a satisfactory agreement between experimental and numerical results, especially taking into account the uncertainty in the initial condition and a possible homogeneity fault in the experimental device. The first-order truncations in the description of anisotropy (§ 3.5) may also lie at the origin of discrepancies between experimental and numerical results. The system of governing equations (3.3)–(3.5) allows us to correctly capture the evolution of anisotropy, both in the straining regions and during the relaxation phases. Only the period of RTI for the angle $\alpha = \pi/4$ is not fully satisfactory, mainly because the boundary between the straining and relaxation regions in the experiments does not appear to be as clear as in the simulations. The case of straining without rotation in the second part of the distorting duct ($\alpha = 0$), followed by a relaxation phase, is further illustrated in figure 2. The present model properly captures the evolution of the anisotropy indicators $b_{ij}$ and that of the turbulent kinetic energy, both in the region dominated by linear effects and in the purely nonlinear one.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 1. Evolution of the invariant $\mathcal{I}$ versus the position in the distorting duct of length $L_d$ for the experiments of (a) Gence & Mathieu (1979) and (b) Gence & Mathieu (1980). Symbols correspond to experimental data and lines are obtained with the system of governing equations (3.3)–(3.5). Various values of the angle $\alpha$ between the principal axes of the two successive plane strains are investigated: $\alpha = 0$ (□, ----–), $\alpha = \pi/8$ (+, ····), $\alpha = \pi/4$ (○, ——), $\alpha = 3\pi/8$ (△, ----) and $\alpha = \pi/2$ (×, —).</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 2. Evolution of (a) the anisotropy indicators $b_{11}$ (△, ----–), $b_{22}$ (□, ----–) and $b_{33}$ (○, —) and (b) that of the turbulent kinetic energy $\mathcal{K}$ (○, —) versus the position in the distorting duct of length $L_d$ for the experiment in Gence & Mathieu (1980) without rotation in the second part of the distorting duct ($\alpha = 0$). Symbols correspond to experimental data and lines are obtained with the system of governing equations (3.3)–(3.5).</figcaption></figure>

### 4.3. Turbulence subjected to a straining-relaxation-destraining cycle

The present model is further tested by comparing its predictions with the experiment of Chen et al. (2006), where a piston is used to apply plane straining and distraining on turbulence generated by active grids. The mean-velocity gradient in the experiment is of the form

$$ \lambda(t) = \begin{pmatrix} S(t) & 0 & 0 \\ 0 & -S(t) & 0 \\ 0 & 0 & 0 \end{pmatrix}, \qquad (4.12) $$

where the temporal evolution of $S(t)$ is given by figure 3(a). Initially, the mean flow corresponds to plane straining ($S(t) > 0$), until $t/\tau_0 \simeq 0.5$. After a relaxation
<figure><img src="image_7.png" /><figcaption>FIGURE 3. (a) Temporal evolution of the strain $S(t)$ applied to the turbulence; (b) experimental values (○), numerical values obtained with the present model (—) and RDT prediction (---) for the temporal evolution of the anisotropy indicator $\tilde{b}_{11}(t)$ in the experiment of Chen et al. (2006).</figcaption></figure>

phase ($0.5 \le t/\tau_0 \le 0.7$), destraining ($S(t) < 0$) is applied to the turbulence. In this experiment, the Taylor-microscale-based Reynolds number at the beginning of the straining cycle is $Re_\lambda \simeq 400$. The maximum value of the strain $S(t)$ reached in the experiment is $\simeq 9.5 \tau_0^{-1}$. Figure 3(b) illustrates the temporal evolution of the anisotropy indicator $\tilde{b}_{11}(t)$ defined by (4.7). Experimental and numerical values obtained with the present model are reported, along with the RDT prediction, provided by Chen et al. (2006), corresponding to the mean flow defined by (4.12) and figure 3(a). The temporal evolution of $\tilde{b}_{11}(t)$ shows good agreement between the experiment and the present model. From the comparison with RDT results, it appears that nonlinear phenomena are significant on a quantitative level. This is partly due to the presence of a relaxation phase in the straining cycle. Thus, the validity of both linear and nonlinear contributions in the system of governing equations (3.3)–(3.5) can be confirmed by the comparison with this experiment.

### 4.4. Turbulence subjected to axisymmetric expansion or contraction

We now consider the case of turbulence subjected to an axisymmetric expansion or contraction. The corresponding mean-velocity gradients are respectively

$$ \lambda = \begin{pmatrix} S/2 & 0 & 0 \\ 0 & S/2 & 0 \\ 0 & 0 & -S \end{pmatrix}, \quad \lambda = \begin{pmatrix} -S/2 & 0 & 0 \\ 0 & -S/2 & 0 \\ 0 & 0 & S \end{pmatrix}, \qquad (4.13a,b) $$

with $S > 0$. These two configurations have been recently investigated in the DNS of Zusi & Perot (2014), where both effects of axisymmetric expansion and contraction on initially approximately isotropic turbulence and its subsequent relaxation have been studied. We want to check whether the present model is able to reproduce one of the main observation of this paper: the two-stage RTI of turbulence after an axisymmetric expansion. In this configuration, it has been observed that the anisotropy continues to increase after the release of the strain during a short period of time, before decreasing in a more conventional way. The reader has to keep in mind that the microscale-based Reynolds number in this work is moderate ($Re_\lambda \simeq 50$) and that two-stage RTI
<figure><img src="image_8.png" /></figure>

FIGURE 4. Time evolution of the different components of the global anisotropy tensor $b_{ij}(t)$ (----- $b_{11}$, ----- $b_{22}$, --- $b_{33}$) for (a) axisymmetric expansion and (b) contraction. The thick black lines correspond to numerical values obtained with the present model, while the thin black and grey lines refer to the DNS results of Zusi & Perot (2014) for two different initial conditions. The vertical dotted line marks the time at which strain is released.

is observed only after an expansion with sufficiently high strain rate ($S \ge 3\tau_0^{-1}$). A complete parametric study of this two-stage RTI should be performed, but this is beyond the scope of the present paper. Detailed numerical analyses of the RTI of axisymmetric turbulence can be found in Herring (1974), Chasnov (1995), Davidson, Okamoto & Kaneda (2012) and Mons, Meldi & Sagaut (2014).

First, a quantitative comparison between the predictions of the present model and the DNS results of Zusi & Perot (2014) is performed. Zusi & Perot (2014) use two different initial conditions in their work which are difficult to characterize. Besides, these initial fields are anisotropic and non-axisymmetric. Predictions concerning the time evolution of the different components of the tensor $b_{ij}(t)$ obtained with the present model and the DNS of Zusi & Perot (2014) are reported in figure 4. Both expansion and contraction cases are investigated with $S \simeq 3.4\tau_0^{-1}$, and simulations are initialized with $Re_\lambda = 50$. The production periods are well recovered by the model. In the RTI regions, significant discrepancies between the DNS results obtained with the two initial fields can be noticed, and the results obtained with the present model lie between these two predictions. Thus, given the uncertainties in the initial condition and a potential additional source of discrepancies originating from the limitation to moderate anisotropy of the present model (§ 3.5), these comparisons are relatively satisfactory. In particular, the major difference between the expansion and contraction cases is recovered: after the release of the strain, the anisotropy immediately decreases in the contraction case, whereas it slightly increases in the expansion case before decreasing. To further illustrate this observation, complementary runs are performed starting from initially isotropic turbulence, instead of an anisotropic and non-axisymmetric field, in order to get clearer results.

Figures 5 and 6 illustrate numerical results obtained with the present model for initially isotropic turbulence submitted to an axisymmetric expansion or contraction and then released. Simulations are initialized with $Re_\lambda = 50$ and we set $S = 5\tau_0^{-1}$. The axisymmetric expansion or contraction is maintained from $t = 0$ to $St = 1$. We first consider the results for axisymmetric expansion, which are reported in figures 5(a,b) and 6(a). The temporal evolution of the anisotropy indicator $b_{33}(t)$ is given in
<figure><img src="image_9.png" /><figcaption>FIGURE 5. (a,c) Time evolution of the anisotropy indicator $b_{33}(t)$ (—); (b,d) time evolution of $b_{33}^{(dir)}(k, t)$ (—) and $b_{33}^{(pol)}(k, t)$ (---) with that of their time derivatives; (a,b) and (c,d) correspond to the axisymmetric expansion and contraction cases respectively. The vertical dotted line marks the time at which strain is released.</figcaption></figure>

figure 5(a). Since we consider axisymmetric turbulence, the tensor $b_{ij}(t)$ is diagonal, with $b_{33}(t) = -b_{11}(t)/2 = -b_{22}(t)/2$. As observed above, $b_{33}(t)$ slightly increases after the release of the strain before decaying. Further insight into this two-stage RTI is given by figure 5(b), where the temporal evolution of the directivity and polarization components of $b_{33}(t) = b_{33}^{(dir)}(t) + b_{33}^{(pol)}(t)$ is reported, along with their time derivatives. Contrary to $b_{33}(t)$, both $b_{33}^{(dir)}(t)$ and $b_{33}^{(pol)}(t)$ immediately decay in magnitude after the release of the strain. However, since $b_{33}^{(dir)}(t)$ is negative, whereas both $b_{33}^{(pol)}(t)$ and $b_{33}(t)$ are positive, and it first decreases in magnitude at a faster rate than $b_{33}^{(pol)}(t)$, the total anisotropy indicator $b_{33}(t)$ increases. As soon as the polarization component decays at a faster rate than the directivity one, the global anisotropy decreases. The difference in the relaxation rates between directional anisotropy and polarization anisotropy in the present model, which lies at the origin of the increase of the total anisotropy, may be related to the difference in the relaxation rates between the dimensionality and circulicity tensors in Kassinos & Reynolds (1997) and Kassinos et al. (2001) (see appendix A for further details).

The case of RTI after an axisymmetric contraction is illustrated in figures 5(c,d) and 6(b). The anisotropy indicator $b_{33}(t)$ decays in magnitude immediately after the release of the strain, as $b_{33}^{(dir)}(t)$ and $b_{33}^{(pol)}(t)$ do. The polarization component, of the same sign as $b_{33}(t)$, decays in magnitude at a faster rate than the directivity one as soon as the strain is released. The qualitative difference in the RTI between expansion and contraction is even more striking by looking at figure 6(a,b), where the evolution of the invariant $II(t)$ is indicated for the expansion and contraction cases respectively. The temporal evolution of the invariant $III(t)$ and that of the ratio $\rho(t)$ between the
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 6. (a,b) Temporal evolution of the invariant $II(t)$ (—) for the axisymmetric (a) expansion and (b) contraction cases; (c,d) temporal evolution of (c) the invariant $III(t)$ and (d) the ratio $\rho(t)$ between characteristic times of decay and RTI for the expansion (— — — —) and contraction (— —) cases. The vertical dotted line marks the time at which strain is released.</figcaption></figure>

characteristic times of decay and RTI defined by (4.6) are reported in figure 6(c,d). The expansion case is associated with a positive value for $III(t)$ (axial component $\langle u_3 u_3 \rangle(t)$ greater than the transverse ones), whereas $III(t)$ is negative in the contraction case. Figure 6(d) indicates that the RTI process is slower in the expansion case than in the contraction case ($\rho(t)$ takes a greater value if the RTI is faster and is negative during the slight increase of anisotropy in the expansion case). This result is consistent with the discussions in Gence (1983) and Choi & Lumley (2001).

4.5. The RTI of shear-released turbulence

The RTI of initially isotropic (Saffman) turbulence that is suddenly subjected to a mean shear and is then released is investigated. The mean-velocity gradient that corresponds to the application of a shear on the turbulent flow is

$$ \lambda = \begin{pmatrix} 0 & 0 & S \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (4.14) $$

The simulations are initialized with $Re_\lambda = 10^4$. The anisotropy of the flow is introduced by using either the RDT solution for homogeneous shear flow with $St=0.42$ (this value allows us to observe a significant departure from isotropy of the diagonal components of the Reynolds stress tensor) or the linear terms of (3.3)–(3.5)
<figure><img src="image_14.png" /><figcaption>FIGURE 7. Global anisotropy indicators — $b_{ij}$, — $b_{ij}^{(dir)}$ and — — — $b_{ij}^{(pol)}$ in the case of shear-released Saffman turbulence. The grey lines correspond to the use of RDT solutions for the initialization, whereas the black lines correspond to the use of the linear terms in (3.3)–(3.5) for the introduction of anisotropy. The vertical black dotted line marks the limit at $t = 0.042\tau_0$ after which the shear is released when the linear terms are used.</figcaption></figure>

with $S = 10\tau_0^{-1}$ if $t < 0.042\tau_0$ and $S = 0$ if not, starting from an isotropic field. Results are illustrated in figures 7–9.

It appears from figure 7 that solutions initialized with RDT or with the linear terms (3.7)–(3.9) coincide, which confirms the validity of the linear terms of the present model in the case of homogeneous shear turbulence. The evolution of anisotropy is similar to that in the case of axisymmetric Saffman turbulence (Chasnov 1995; Davidson et al. 2012; Mons et al. 2014). After a transient regime which corresponds to the RTI of small scales, an asymptotic anisotropic state is reached due to the fulfilment of the PLE hypothesis ($E(k, t) = E(k, 0)$ for $k \ll k_{max}(t)$ if $1 \le \sigma \le 3$ with $E(k \to 0, t) \propto k^\sigma$). A similar behaviour was also observed with the model of Clark & Zemach (1995) for the relaxation of turbulence after the release of plane strain. Even the off-diagonal component $b_{13}(t)$ converges towards a non-zero value, which means that the cross-correlation $\langle u_1 u_3 \rangle(t)$ is maintained even though the mean shear is released.

The RTI of small scales during the transient regime is illustrated in figure 8 for the components of the spherically integrated second-order spectral tensor defined by (1.4), $\varphi_{13}(k, t)$ and $\varphi_{33}(k, t)$. Both the directional anisotropy and the polarization anisotropy become negligible in the inertial range for $t \ge 10\tau_0$. This result is also valid for the diagonal components $\varphi_{11}(k, t)$ and $\varphi_{22}(k, t)$ that are not illustrated here. Figure 8(b) indicates that the initial cross-correlation spectrum $\varphi_{13}(k, 0)$ given by RDT evolves like $k^{-5/3}$ in the inertial range. Then, due to purely nonlinear processes, this shape
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 8. (a) Energy spectrum $E(k, t)$, (b) cross-correlation spectrum $\varphi_{13}(k, t)$ and spectra (c) $H_{13}^{(\text{dir})}(k, t)$, (d) $H_{13}^{(\text{pol})}(k, t)$, (e) $H_{33}^{(\text{dir})}(k, t)$ and (f) $H_{33}^{(\text{pol})}(k, t)$ for shear-released Saffman turbulence. The RDT is used for the initialization. The curves have been sampled at the normalized time: $\underline{---} t/\tau_0 = 0$; $--- t/\tau_0 = 10^{-3}$; $---- t/\tau_0 = 10^{-2}$; $----- t/\tau_0 = 10^{-1}$; $\underline{---} t/\tau_0 = 1$; $--- t/\tau_0 = 10^{-4}$.</figcaption></figure>

is modified and evolves like $k^{-7/3}$ at $t = 10\tau_0$. This result is consistent with the predictions of Lumley (1967), Leslie (1973), Yoshida et al. (2003) and Weinstock (2013).

The budget terms of the governing equation (4.8) are reported in figure 9 for the spectra $\varphi_{13}(k, t)$ and $\varphi_{33}(k, t)$ at $t = 0.1\tau_0$ and $t = 10\tau_0$. During the transient regime corresponding to the RTI of small scales ($t = 0.1\tau_0$), the nonlinear pressure–strain rate tensor $P_{ij}(k, t)$ has a significant influence in the inertial range. Its contribution is positive for the negative cross-correlation spectrum $\varphi_{13}(k, t)$, which would be identically zero in the isotropic case, and also positive for the spectrum $\varphi_{33}(k, t)$, since this component has been relatively damped by the mean shear via linear effects. Once small scales have returned to isotropy ($t = 10\tau_0$), the contribution of
<figure><img src="image_21.png" /></figure>

FIGURE 9. Budget terms of the equation $(\partial/\partial t)\varphi_{ij}(k, t) = -2\nu k^2 \varphi_{ij}(k, t) + P_{ij}(k, t) + S_{ij}(k, t)$, premultiplied by $k$ (----- $kS_{ij}(k, t)$; ....... $-2\nu k^3 \varphi_{ij}(k, t)$; ---- $kP_{ij}(k, t)$; $k(\partial/\partial t)\varphi_{ij}(k, t)$), and normalized by $\max_k(|k(\partial/\partial t)\varphi_{ij}(k, t)|)$. Shear-released Saffman turbulence is considered and the RDT is used for the initialization; (a,c) refer to the component $\varphi_{13}(k, t)$, whereas (b,d) refer to the component $\varphi_{33}(k, t)$. The budget terms have been sampled at the normalized time $t/\tau_0 = 10^{-1}$ (a,b) and $t/\tau_0 = 10$ (c,d).

$P_{ij}(k, t)$ is negligible $\forall k$ for the component $\varphi_{33}(k, t)$ (as well as for the other diagonal components), even though large scales have not returned to isotropy. The shape of the transfer term $S_{33}(k, t)$ is then virtually identical to that in the isotropic case. Concerning the cross-correlation component, nonlinear terms act predominantly at scales close to the integral length scale. These results are consistent with previous studies dealing with the RTI of Saffman turbulence, where it is observed that small scales quickly return to isotropy after the release of the shear/strain, whereas large scales fully retain anisotropy due to the fulfilment of the PLE hypothesis.

## 4.6. Homogeneous shear turbulence

Finally, we address the case of homogeneous turbulence subjected to a constant maintained mean shear. The corresponding mean-velocity gradient is given by (4.14). Weinstock (2013) performed an exhaustive analytical study of this configuration in 3D Fourier space, without limitations on time or wavenumber. In this subsection, it is checked that the present model is able to recover results established in this prior work, among others. The corresponding simulation realized with the present model is initialized with $Re_\lambda = 50$, the shear rate is fixed at $S = 2\tau_0^{-1}$ and the turbulence is initially isotropic. The temporal evolutions of the components of the deviatoric tensor $b_{ij}(t)$ and that of the kinetic energy $K(t)$ are reported in figure 10. The different components $b_{ij}(t)$ and the ratio $\varepsilon(t)/SK(t)$ reach constant values after $St \approx 20$, as
<figure><img src="image_22.png" /><figcaption>FIGURE 10. Temporal evolutions of (a) the different components of the global anisotropy tensor $b_{ij}(t)$ and the ratio $\varepsilon/S\mathcal{H}$ (----- $b_{11}$, ------ $b_{13}$, ...... $b_{22}$, ------ $b_{33}$ and ------ $\varepsilon/S\mathcal{H}$) and (b) those of the turbulent kinetic energy $\mathcal{K}(t)$ (----), $k_{max}$ (----) and $E_{max}$ (-----) defined in (4.10) and normalized by their initial values for homogeneous shear turbulence. The grey curves are plotted with $\gamma = 0.337$.</figcaption></figure>

observed in, e.g., Weinstock (2013). Considering the governing equation for the turbulent kinetic energy $\mathcal{K}(t)$,

$$ \frac{d}{dt}\mathcal{K}(t) = -2S\mathcal{K}(t)b_{13}(t) - \varepsilon(t), \quad (4.15) $$

one can deduce that for $St \ge 20$

$$ \mathcal{K}(t) \sim e^{\gamma St}, \quad \gamma = \text{const.} = -2b_{13} - \frac{\varepsilon}{\mathcal{K}S}. \quad (4.16) $$

As illustrated in figure 10(b), the results obtained with the present model are consistent with such an exponential growth of the kinetic energy. The model also predicts the temporal evolutions of the peak of the energy spectrum $E_{max}$ and the corresponding wavenumber $k_{max}$ defined in (4.10) as $E_{max}(t) \sim e^{(3/2)\gamma St}$ and $k_{max}(t) \sim e^{-\gamma St/2}$. All of these results are consistent with dimensional analysis and with the study performed in Weinstock (2013), in which $\gamma = 0.115$. The estimated value of the growth rate $\gamma$ in a set of experiments (Tavoularis & Corrsin 1981; Tavoularis & Karnik 1989) ranges between 0.08 and 0.12, whereas values between 0.1 and 0.2 can be found in DNS studies (Rogers et al. 1986; Brethouwer 2005; Isaza & Collins 2009). Although a clear consensus about the value of the growth rate $\gamma$ cannot be found in the literature, and its sensitivity with respect to initial conditions, Reynolds number and shear rate has to be further investigated, the present model overestimates $\gamma$ ($\simeq 0.337$) in comparison with values found in the literature. This lack of quantitative agreement originates from the first-order truncations in the description of anisotropy (§ 3.5). Compared with previous test cases considered in the present paper, the cumulated shear reached in the simulation illustrated in figures 10 and 11 is much higher (see § 4.7), which provides asymptotic information at high Reynolds numbers. Concerning results in spectral space, the kinetic energy spectrum $E(k, t)$ at $St = 50$ (figure 11a) displays a $-5/3$ slope in the inertial range, as reported in experiments or DNS. Figure 11(b) reports the cross-correlation spectrum $\varphi_{13}(k, t)$
<figure><img src="image_23.png" /><figcaption>FIGURE 11. (a) Energy spectrum $E(k, t)$ and (b) cross-correlation spectrum $\varphi_{13}(k, t)$ at $St = 0$ (thick black line) and $St = 15$ (thin black lines). (c–f) Budget terms of the governing equation $(\partial/\partial t)\varphi_{ij}(k, t) = -2\nu k^2 \varphi_{ij}(k, t) + P_{ij}(k, t) + S_{ij}(k, t) + L_{ij}(k, t)$, premultiplied by $k$ (----- $kS_{ij}(k, t)$; ...... $-2\nu k^3 \varphi_{ij}(k, t)$; ------ $kP_{ij}(k, t)$; --- $kL_{ij}(k, t)$; --- $k(\partial/\partial t)\varphi_{ij}(k, t)$), and normalized by $\max_k(|k(\partial/\partial t)\varphi_{ij}(k, t)|)$. The curves are sampled at $St = 50$ for homogeneous shear turbulence and refer to the components (c) $\varphi_{11}(k, t)$, (d) $\varphi_{13}(k, t)$, (e) $\varphi_{22}(k, t)$ and (f) $\varphi_{33}(k, t)$. The position of the shear scale $L_s$ is also shown by vertical grey dashed lines.</figcaption></figure>

at $St = 50$. The latter evolves like $k^{-7/3}$ in the inertial range, as predicted theoretically (Lumley 1967; Leslie 1973; Yoshida et al. 2003; Weinstock 2013) and observed in the experiments of Shen & Warhaft (2000) and the DNS of Ishihara et al. (2002). This result supports a *posteriori* choice of eddy damping (2.26) and (2.27). The budget terms of the governing equation (4.8) are illustrated in figure 11(c-f). The position of the shear scale $L_s = \sqrt{\varepsilon/S^3}$ (Corrsin 1958) is also reported in the figures. For $1/L \ll k \ll 1/L_s$, the flow is expected to be dominated by production terms, whereas for $k \gg 1/L_s$, the contribution of nonlinear transfers should prevail. This interpretation is reasonably well supported by the present results. Figure 11(c-f) can help in visualizing the nonlinear process that lies at the origin of the exponential
<figure><img src="image_24.png" /><figcaption>FIGURE 12. Energy spectrum (a–c) and maximal eigenvalue of the tensor $H_{ij}^{(\text{dir})}(k, t)$ (d–f) for the different test cases investigated in this paper. The value 1/15 from the criterion (3.17) is displayed as dashed lines. (a,d) The case of plane straining (case $\alpha = \pi/8$ of figure 1a); the curves are sampled at the end of the simulation ($St \approx 1.4$). (b,e) The case of axisymmetric expansion of figure 4(a) ($St = 0.5$). (c,f) Results concerning homogeneous shear turbulence (figures 10 and 11); the curves are sampled at $St = 50$.</figcaption></figure>

growth of the turbulent kinetic energy. Since the component $\varphi_{33}(k, t)$ is less than the other diagonal components by the mean shear, the nonlinear pressure–strain rate tensor $P_{33}(k, t)$ is strongly positive and redistributes energy from the components $\varphi_{11}(k, t)$ and $\varphi_{22}(k, t)$ to $\varphi_{33}(k, t)$ ($P_{ii}(k, t) = 0 \forall k, t$). Since the cross-correlation component of the Reynolds stress tensor $2K_{13}(t)$ is fed by the mean shear via the term $-2SK_{33}(t)$, its production is enhanced, allowing an increase of the growth rate of the turbulent kinetic energy $\mathcal{K}(t)$ via the term $-2S\mathcal{K}(t)b_{13}(t) = -2SK_{13}(t)$.

### 4.7. Realizability condition

Due to spherical integration and truncation of expansions (3.14) and (3.15), the present model is limited to moderate anisotropy, and the criterion (3.17) is derived in § 3.5 in order to quantify the upper boundary of anisotropy intensity that can be investigated with the representation (3.2) and the corresponding system of governing equations (3.3)–(3.5). In this last subsection, it is checked that the different test cases considered in this paper can be described using the representation (3.2) and (3.3)–(3.5). Figure 12 reports the criterion (3.17) for typical flow configurations investigated in this paper: plane straining, axisymmetric expansion and pure plane shear. Rather than the absolute value of the typical strain or shear rate $S$, it is the accumulated strain/shear $St$ that gives a good estimation of the anisotropy introduced in the flow. The interplay between linear and nonlinear contributions in (3.3)–(3.5) also affects the evolution of anisotropy in the flow, but this effect must be studied on a case-by-case basis in view of the significant differences that exist between, e.g., rotational and irrotational mean flows (Sagaut & Cambon 2008). It is worth keeping in mind that for the purpose of spherical integration, the same representation (3.2) is used for both linear contributions and transfer terms in (2.12) and (2.13), the
latter being closed by the EDQNM. The different cases of plane straining/destraining investigated in this paper are associated with small values for the accumulated strain. Figure 12(d) illustrates the case $\alpha = \pi/8$ of figure 1(a) from the experiments of Gence and Mathieu. This configuration corresponds to the accumulated strain $St \simeq 1.4$, and the anisotropy of the flow at all scales is significantly lower than the upper bound of (3.17). This result is similar for the other cases of plane straining/destraining in this paper. The case of axisymmetric expansion of figure 4(a) is investigated in figure 12(e). The accumulated strain is $St = 0.5$ and the criterion (3.17) is respected at all scales. The same applies to the axisymmetric contraction case. Finally, the case of homogeneous shear turbulence is illustrated in figure 12(f) at $St = 50$, and criterion (3.17) is verified at all scales. Thus, all of the flow configurations considered in the present paper can be described by the representation (3.2) and the corresponding system of governing equations (3.3)–(3.5).

## 5. Conclusion and perspectives

Modelling of anisotropy in homogeneous flows has been considered at two different levels. The first one is the three-dimensional spectral level, in which a decomposition of the spectral tensor for arbitrary anisotropy lends support to a splitting of anisotropy at any subsequent level in terms of directional anisotropy and polarization anisotropy. In a second step, a model for spherically integrated quantities has been proposed, which is based on 11 coupled equations. It is dedicated to turbulent flows where the anisotropy is moderate and where linear effects induced by mean-velocity gradients play a negligible role in the dynamics of triple correlations compared with the induced production effects in the equations for the second-order correlations. This model is not restricted to a particular symmetry and can be used for a wide range of flow configurations, as illustrated by the different applications of the model considered in this paper. A satisfactory agreement with the experiments of Gence & Mathieu (1979, 1980) has been observed, which confirms the capability of the model to account for production of anisotropy by mean-flow gradients.

Concerning the RTI of initially deformed or sheared turbulence, the model is consistent with the PLE and ensures correct rapid RTI of the smallest scales. Different relaxation rates for directional anisotropy and polarization anisotropy allow us to correctly interpret the apparent delay in the RTI after axisymmetric expansion (Zusi & Perot 2014). In addition, our model fits well the recent experiment of Chen et al. (2006) with a non-stationary straining cycle. For turbulence continuously subjected to a pure plane shear, the model ensures a correct asymptotic regime with constant values for the components of the dimensionless deviatoric tensor $b_{ij}(t)$ associated with the Reynolds stress tensor. In addition, it can reproduce the exponential growth of the turbulent kinetic energy mediated by nonlinear pressure redistribution terms. However, first-order truncations in the description of anisotropy, which are the main sources of possible discrepancies between the present model and DNS/experimental results, prevent a good quantitative agreement with typical values of the growth rate found in the literature. For this configuration, an alternative approach, which is entirely formulated in 3D Fourier space and without limitation in time, can be found in Weinstock (2013).

Beyond moderate anisotropy, a more complex version of the model could possibly combine exact linear operators in $k$-vectors, as in the left-hand side of (2.12)–(2.13), with transfer terms only generated by low-order angular harmonic expansions. Such a model could reproduce the dominant RDT dynamics for the largest scales, and the quasi-isotropic behaviour for scales smaller than a Corrsin scale.
# Appendix A. Relationships between two-point spectral tensors and single-point tensors

The equations addressed in this paper involve several tensors: production by mean-velocity gradients, dissipation, transfer, pressure–strain rate. In the following, we give their expressions systematically in terms of $k$, then in terms of $k$, obtained by spherically averaging the first one, and their final single-point contribution obtained by integrating on $k$, as in full Reynolds stress models (RSMs) with an additional structure-based tensor. Time dependence is implied (not explicitly written below).

## A.1. Production term by space-uniform mean-velocity gradient

The production term by the space-uniform mean-velocity gradient in (2.3) is $\lambda_{il}\hat{R}_{ij}(k) + \lambda_{jl}\hat{R}_{li}(k)$. By using the representation (3.1), its spherically integrated counterpart is

$$ \iint_{S_k} \lambda_{il} \hat{R}_{ij}(k) + \lambda_{jl} \hat{R}_{il}(k) \, d^2 k = \frac{4}{3} S_{ij} E + 2E (\lambda_{il} H_{lj}^{(dir)} + \lambda_{jl} H_{li}^{(dir)}) + 2E (\lambda_{il} H_{lj}^{(pol)} + \lambda_{jl} H_{li}^{(pol)}), \quad (A1) $$

with $E = E(k)$, $H_{ij}^{(dir)} = H_{ij}^{(dir)}(k)$ and $H_{ij}^{(pol)} = H_{ij}^{(pol)}(k)$. The single-point counterpart of the production term is given by

$$ \iiint \lambda_{il} \hat{R}_{ij}(k) + \lambda_{jl} \hat{R}_{il}(k) \, d^3 k = \frac{4}{3} S_{ij} \mathcal{H} + 2 \mathcal{H} (\lambda_{il} b_{lj}^{(dir)} + \lambda_{jl} b_{li}^{(dir)}) \\ + 2 \mathcal{H} (\lambda_{il} b_{lj}^{(pol)} + \lambda_{jl} b_{li}^{(pol)}), \quad (A2) $$

where the tensors $b_{ij}^{(dir)}$ and $b_{ij}^{(pol)}$ are defined by (4.4).

## A.2. Dissipation term

The dissipation term in (2.3) is $2\nu k^2 \hat{R}_{ij}(k)$. Its spherically integrated counterpart is

$$ \iint_{S_k} 2\nu k^2 \hat{R}_{ij}(k) \, d^2 k = 4\nu k^2 E \left( \frac{\delta_{ij}}{3} + H_{ij}^{(dir)} + H_{ij}^{(pol)} \right). \quad (A3) $$

The corresponding single-point contribution is

$$ \iiint 2\nu k^2 \hat{R}_{ij}(k) \, d^3 k = \frac{2}{3} \delta_{ij} \varepsilon + \varepsilon_{ij}^{(dir)} + \varepsilon_{ij}^{(pol)}, \quad (A4) $$

where $\varepsilon_{ij}^{(dir)} = \int_0^\infty 4\nu k^2 E(k) H_{ij}^{(dir)}(k) dk$, with a similar definition for $\varepsilon_{ij}^{(pol)}$.

## A.3. ‘Rapid’ contribution of the pressure–strain rate tensor

The rapid contribution of the pressure–strain rate tensor is given by $2\lambda_{in}\alpha_l(\alpha_i \hat{R}_{nj}(k) + \alpha_j \hat{R}_{ni}(k))$. Its spherically integrated counterpart is

$$ \iint_{S_k} 2\lambda_{in} \alpha_l (\alpha_i \hat{R}_{nj}(k) + \alpha_j \hat{R}_{ni}(k)) \, d^2 k \\ = \frac{4}{5} E S_{ij} - \frac{12}{7} E \left( S_{lj} H_{li}^{(dir)} + S_{li} H_{lj}^{(dir)} - \frac{2}{3} S_{lm} H_{lm}^{(dir)} \delta_{ij} \right) $$
$$ 
\begin{aligned}
&+4E(A_{jl}H_{il}^{(dir)} + A_{il}H_{jl}^{(dir)}) + \frac{12}{7}E(S_{lj}H_{li}^{(pol)} + S_{li}H_{lj}^{(pol)} - \frac{2}{3}S_{lm}H_{lm}^{(pol)}\delta_{ij}) \\
&+\frac{4}{3}E(A_{il}H_{lj}^{(pol)} + A_{jl}H_{li}^{(pol)}).
\end{aligned} 
\quad (A 5) 
$$

Similarly to the transition from (A1) to (A2), the single-point contribution of the rapid component of the pressure–strain rate tensor is obtained from (A5) by replacing $E(k)$, $H_{ij}^{(dir)}(k)$ and $H_{ij}^{(pol)}(k)$ with $\mathcal{K}$, $b_{ij}^{(dir)}$ and $b_{ij}^{(pol)}$ respectively. The closure of the rapid pressure–strain rate tensor in RSMs is generally applied to the tensor $M_{ijpq}$, defined by

$$ M_{ijpq} = \iiint \alpha_p \alpha_q \hat{R}_{ij}(\mathbf{k}) \, d^3\mathbf{k}, \quad (A 6) $$

and made non-dimensional by the turbulent kinetic energy. The closure is usually in terms of the deviatoric tensor $b_{ij}$. Tensorial expansions range from linear (Launder, Reece & Rodi 1975; Lumley 1975) with a single tuned constant, to quadratic (Speziale, Sarkar & Gatski 1991) and even cubic (Craft, Ince & Launder 1996; Craft & Launder 2001) with increase in the number of tuned constants as the degree of nonlinearity increases.

### A.4. Transfer tensor from linear origin

The spherically averaged version of the term $\lambda_{ln}k_l(\partial\hat{R}_{ij}/\partial k_n)(\mathbf{k})$ is

$$ 
\begin{aligned}
\iint_{S_k} \lambda_{ln} k_l \frac{\partial \hat{R}_{ij}(\mathbf{k})}{\partial k_n} d^2\mathbf{k} = & \frac{4}{7} \left( S_{il} \frac{\partial}{\partial k} (k E H_{lj}^{(dir)}) + S_{lj} \frac{\partial}{\partial k} (k E H_{li}^{(dir)}) - 3 S_{lm} \frac{\partial}{\partial k} (k E H_{lm}^{(dir)}) \delta_{ij} \right) \\
& - \frac{4}{7} \left( S_{jl} \frac{\partial}{\partial k} (k E H_{il}^{(pol)}) + S_{il} \frac{\partial}{\partial k} (k E H_{lj}^{(pol)}) - \frac{2}{3} S_{ln} \frac{\partial}{\partial k} (k E H_{ln}^{(pol)}) \delta_{ij} \right) - \frac{2}{15} S_{ij} \frac{\partial}{\partial k} (k E).
\end{aligned} 
\quad (A 7) 
$$

The $k$ integral of the above expression is 0.

### A.5. Transfer tensor from nonlinear origin and 'slow' part of the pressure–strain rate tensor

The expression of the 'true' nonlinear transfer tensor, with zero integral over $k$, is $\tau_{ij}(\mathbf{k}) + \tau_{ji}^*(\mathbf{k})$, and that of the slow part of the pressure–strain rate tensor is $-\alpha_i \alpha_n \tau_{nj}(\mathbf{k}) - \alpha_i \alpha_n \tau_{ni}^*(\mathbf{k})$, with $\tau_{ij}(\mathbf{k})$ defined by (2.6). The spherically integrated counterparts of these tensors are given by (3.6) and (3.10), (3.13). Concerning RSMs, the slow pressure–strain rate tensor is generally closed as $-C\epsilon b_{ij}$ in order to relax the dimensionless deviatoric part of the Reynolds stress tensor $b_{ij}$ weighted by the dissipation rate $\epsilon$. More or less complicated expressions were proposed instead of the constant $C$.

### A.6. Kassinos et al. (2001) structure-based tensors

The dimensionality tensor can be derived from a special index contraction of the tensor $M_{ijpq}$ defined by (A 6), as $D_{ij} = M_{llij}$, whereas the two other contractions give 0 or the Reynolds stress tensor itself. Accordingly, it is found that

$$ D_{ij} = \iiint \alpha_i \alpha_j \hat{R}_{mm}(\mathbf{k}) \, d^3\mathbf{k} = \int_0^\infty 2E(k) \left( \frac{1}{3}\delta_{ij} - 2H_{ij}^{(dir)}(k) \right) dk = 2\mathcal{K} \left( \frac{1}{3}\delta_{ij} - 2b_{ij}^{(dir)} \right). \quad (A 8) $$
The circulicity tensor $F_{ij}$ is not independent and is given by

$$F_{ij} = 2 \mathcal{H} \left( \frac{1}{3} \delta_{ij} + b_{ij}^{(\text{dir})} - b_{ij}^{(\text{pol})} \right). \qquad (A9)$$

The last relevant tensor used by Kassinos et al. (2001) is the ‘stropholysis’ one, derived from $Q_{ijk} = \epsilon_{ipq} M_{jqpk}$, which can be expressed in the present formalism as

$$Q_{ijk} = \mathcal{H} \left( \frac{1}{3} \epsilon_{ikj} - 2 \epsilon_{ipj} b_{kp}^{(\text{dir})} + \frac{4}{3} \epsilon_{ikp} b_{jp}^{(\text{pol})} + \frac{2}{3} \epsilon_{ipj} b_{kp}^{(\text{pol})} \right). \qquad (A10)$$

In a fully symmetrized form, the spectral counterpart (k-vector) of this third tensor is given by contributions from $\alpha_k$ Im $(Z(k)N_i(k)N_j(k))$ (Sagaut & Cambon 2008). The latter expression is related to the term $U_{ijm}^{(\text{pol})3}(k)$ in (3.15), and the contribution from our final model is zero, because of the truncation of the development of $Z(k)$ in terms of second-order angular harmonics. The true ‘stropholysis’ effect, which breaks mirror symmetry but is distinct from helicity, is given by the imaginary part of $Z(k)$: it is dynamically created by the last term on the left-hand side of (2.13) in the presence of rotational mean flows.

### A.7. Towards the $\epsilon$-equation

An equation for $\epsilon_{ij}$ can easily be found, in which $d\epsilon_{ij}/dt$ results from the balance of different terms, obtained by integrating, over $k$ then over $k$, nonlinear transfer terms, purely viscous contributions and linear mean-gradient terms. Only the scalar $\epsilon$-equation is considered now, for comparison with RSMs. It reduces to

$$\frac{d\epsilon}{dt} = 2\nu \int_{0}^{\infty} k^2 T(k) dk - \int_{0}^{\infty} (2\nu k^2)^2 E(k) dk + 2\nu \int_{0}^{\infty} k^2 S^L(k) dk. \qquad (A11)$$

In single-point models, the two first terms on the right-hand side are globally closed as $-C_{\epsilon 2} \epsilon^2 / \mathcal{H}$, whereas the last term is closed as $-C_{\epsilon 1} \lambda_{mn} \langle u_m u_n \rangle \epsilon / \mathcal{H}$. It should be noted that $\int_0^\infty T(k) dk = 0$, but that its integral weighted by $k^2$ is positive, at least in isotropic turbulence, and corresponds to $\langle (\partial u_i / \partial x_j) \omega_i \omega_j \rangle$ or a nonlinear vortex stretching term in physical space, with $\omega_i$ the fluctuating vorticity.

### REFERENCES

* ANDRÉ, J. C. & LESIEUR, M. 1977 Influence of helicity on the evolution of isotropic turbulence at high Reynolds number. *J. Fluid Mech.* **81**, 187–207.

* BELLET, F., GODEFERD, F. S., SCOTT, J. F. & CAMBON, C. 2006 Wave turbulence in rapidly rotating flows. *J. Fluid Mech.* **562**, 83–121.

* BRETHOUWER, G. 2005 The effect of rotation on rapidly sheared homogeneous turbulence and passive scalar transport. Linear theory and direct numerical simulation. *J. Fluid Mech.* **542**, 305–342.

* CAMBON, C., DANAILA, L., GODEFERD, F. S. & SCOTT, J. F. 2013 Third-order statistics and the dynamics of strongly anisotropic turbulent flows. *J. Turbul.* **14** (3), 121–160.

* CAMBON, C. & JACQUIN, L. 1989 Spectral approach to non-isotropic turbulence subjected to rotation. *J. Fluid Mech.* **202**, 295–317.

* CAMBON, C., JEANDEL, D. & MATHIEU, J. 1981 Spectral modelling of homogeneous non-isotropic turbulence. *J. Fluid Mech.* **104**, 247–262.
CAMBON, C., MANSOUR, N. N. & GODEFERD, F. S. 1997 Energy transfer in rotating turbulence. *J. Fluid Mech.* **337**, 303–332.

CAMBON, C. & RUBINSTEIN, R. 2006 Anisotropic developments for homogeneous shear flows. *Phys. Fluids* **18**, 085106.

CAMBON, C. & SCOTT, J. F. 1999 Linear and nonlinear models of anisotropic turbulence. *Annu. Rev. Fluid Mech.* **31**, 1–53.

CANUTO, V. M. & DUBOVIKOV, M. S. 1996a A dynamical model for turbulence. I. General formalism. *Phys. Fluids* **8** (2), 571–586.

CANUTO, V. M. & DUBOVIKOV, M. S. 1996b A dynamical model for turbulence. II. Shear-driven flows. *Phys. Fluids* **8** (2), 587–598.

CANUTO, V. M. & DUBOVIKOV, M. S. 1996c A dynamical model for turbulence. III. Numerical results. *Phys. Fluids* **8** (2), 599–613.

CHASNOV, J. R. 1995 The decay of axisymmetric homogeneous turbulence. *Phys. Fluids* **7** (3), 600–605.

CHEN, J., MENEVEAU, C. & KATZ, J. 2006 Scale interactions of turbulence subjected to a straining–relaxation–destraining cycle. *J. Fluid Mech.* **562**, 123–150.

CHOI, K.-S. & LUMLEY, J. L. 2001 The return to isotropy of homogeneous turbulence. *J. Fluid Mech.* **436**, 59–84.

CLARK, T. T. & ZEMACH, C. 1995 A spectral model applied to homogeneous turbulence. *Phys. Fluids* **7** (7), 1674–1694.

CORRSIN, S. 1958 On local isotropy in turbulent shear flow. *NACA RM 58B11*.

CRAFT, T. J., INCE, N. Z. & LAUNDER, B. E. 1996 Recent developments in second-moment closure for buoyancy-affected flows. *Dyn. Atmos. Oceans* **23**, 99–114.

CRAFT, T. J. & LAUNDER, B. E. 2001 Principles and performance of TCL-based second-moment closures. *Flow Turbul. Combust.* **66**, 355–372.

DAVIDSON, P. A., OKAMOTO, N. & KANEDA, Y. 2012 On freely decaying, anisotropic, axisymmetric Saffman turbulence. *J. Fluid Mech.* **706**, 150–172.

FAVIER, B. F. N., GODEFERD, F. S., CAMBON, C., DELACHE, A. & BOS, W. J. T. 2011 Quasi-static magnetohydrodynamic turbulence at high Reynolds number. *J. Fluid Mech.* **681**, 434–461.

GENCE, J. N. 1983 Homogeneous turbulence. *Annu. Rev. Fluid Mech.* **15**, 201–222.

GENCE, J. N. & MATHIEU, J. 1979 On the application of successive plane strains to grid-generated turbulence. *J. Fluid Mech.* **93**, 501–513.

GENCE, J. N. & MATHIEU, J. 1980 The return to isotropy of an homogeneous turbulence having been submitted to two successive plane strains. *J. Fluid Mech.* **101**, 555–566.

GODEFERD, F. S. & CAMBON, C. 1994 Detailed investigation of energy transfers in homogeneous stratified turbulence. *Phys. Fluids* **6**, 2084–2100.

HERRING, J. R. 1974 Approach of axisymmetric turbulence to isotropy. *Phys. Fluids* **17** (5), 859–872.

ISAZA, J. C. & COLLINS, L. R. 2009 On the asymptotic behaviour of large-scale turbulence in homogeneous shear flow. *J. Fluid Mech.* **637**, 213–239.

ISHIHARA, T., YOSHIDA, K. & KANEDA, Y. 2002 Anisotropic velocity correlation spectrum at small scales in a homogeneous turbulent shear flow. *Phys. Rev. Lett.* **88**, 154501,1–4.

KANEDA, Y. 1981 Renormalized expansions in the theory of turbulence with the use of the Lagrangian position function. *J. Fluid Mech.* **107**, 131–145.

KASSINOS, S. C. & AKYLAS, E. 2012 Advances in particle representation modeling of homogeneous turbulence: from the linear PRM version to the interacting viscoelastic IPRM. In *ERCOFTAC Series*, vol. 18, pp. 81–101.

KASSINOS, S. C. & REYNOLDS, W. C. 1997 Advances in structure-based turbulence modeling. In *Annual Research Briefs – Center for Turbulence Research*, pp. 179–193.

KASSINOS, S. C., REYNOLDS, W. C. & ROGERS, M. M. 2001 One-point turbulence structure tensors. *J. Fluid Mech.* **428**, 231–248.

KIM, S.-W. & CHEN, C.-P. 1989 A multiple-time-scale turbulence model based on variable partitioning of the turbulent kinetic energy spectrum. *Numer. Heat Transfer B* **16**, 193–211.

KRAICHNAN, R. H. 1959 The structure of isotropic turbulence at very high Reynolds numbers. *J. Fluid Mech.* **5**, 497–543.
KRAICHNAN, R. H. 1972 Test-field model for inhomogeneous turbulence. *J. Fluid Mech.* **56**, 287–304.

KRAICHNAN, R. H. & HERRING, J. R. 1978 A strain-based Lagrangian-history turbulence theory. *J. Fluid Mech.* **88**, 355–367.

LAUNDER, B. E., REECE, G. J. & RODI, W. 1975 Progress in the development of a Reynolds-stress turbulent closure. *J. Fluid Mech.* **68**, 537–566.

LESIEUR, M. 2008 *Turbulence in Fluids*, 4th edn. Springer.

LESLIE, D. C. 1973 *Developments in the Theory of Turbulence*. Clarendon.

LUMLEY, J. L. 1967 Similarity and the turbulent energy spectrum. *Phys. Fluids* **10** (4), 855–858.

LUMLEY, J. L. 1975 *Lectures Series*, 76. Von Karman Institute.

MEYERS, J. & MENEVAU, C. 2008 A functional form of the energy spectrum parametrizing bottleneck and intermittency effects. *Phys. Fluids* **20**, 065109.

MILLIONSCHIKOV, M. D. 1941 Theory of homogeneous isotropic turbulence. *Dokl. Akad. Nauk SSSR* **33**, 22–24.

MONS, V., MELDI, M. & SAGAUT, P. 2014 Numerical investigation on the partial return to isotropy of freely decaying homogeneous axisymmetric turbulence. *Phys. Fluids* **26**, 025110.

O'BRIEN, E. F. & FRANCIS, G. C. 1963 A consequence of the zero fourth cumulant approximation. *J. Fluid Mech.* **13**, 369–382.

OGURA, Y. 1963 A consequence of the zero fourth cumulant approximation in the decay of isotropic turbulence. *J. Fluid Mech.* **16**, 33–40.

ORSZAG, S. A. 1970 Analytical theories of turbulence. *J. Fluid Mech.* **41**, 363–386.

PIQUET, J. 2001 *Turbulent Flows – Models and Physics*, 2nd edn. Springer.

POPE, S. B. 2000 *Turbulent Flows*. Cambridge University Press.

POUQUET, A., LESIEUR, M., ANDRÉ, J. C. & BASDEVANT, C. 1975 Evolution of high Reynolds number two-dimensional turbulence. *J. Fluid Mech.* **72**, 305–319.

PROUDMAN, I. & REID, W. H. 1954 On the decay of a normally distributed and homogeneous turbulent velocity field. *Phil. Trans. R. Soc. Lond. A* **297**, 163–189.

ROGERS, M. M., MOIN, P. & REYNOLDS, W. C. 1986 The structure and modelling of the hydrodynamic and passive scalar fields in homogeneous shear turbulence. *Report No. TF-25*. Department of Mechanical Engineering, Stanford University, Stanford, CA.

ROHR, J. J., ITSWEIRE, E. C., HELLAND, K. N. & VAN ATTA, C. W. 1988 An investigation of the growth of turbulence in a uniform-mean-shear flow. *J. Fluid Mech.* **187**, 1–33.

RUBINSTEIN, R. 1996 A relaxation approximation for time-dependent second-order effects in shear turbulence. *Theor. Comput. Fluid Dyn.* **8**, 377–386.

RUBINSTEIN, R., KURIEN, S. & CAMBON, C. 2015 Scalar and tensor spherical harmonics expansion of the velocity correlation in homogeneous anisotropic turbulence. *J. Turbul.* **16** (11), 1058–1075.

SAGAUT, P. & CAMBON, C. 2008 *Homogeneous Turbulence Dynamics*. Cambridge University Press.

SCHIESTEL, R. 1987 Multiple-time-scale modeling of turbulent flows in one-point closures. *Phys. Fluids* **30** (3), 722–731.

SHEN, X. & WARHAFT, Z. 2000 The anisotropy of the small scale structure in high Reynolds number ($R_d \sim 1000$) turbulent shear flow. *Phys. Fluids* **12**, 2976–2989.

SPEZIALE, C. G., SARKAR, S. & GATSKI, T. B. 1991 Modelling the pressure-strain correlation of turbulence: an invariant dynamical systems approach. *J. Fluid Mech.* **227**, 245–272.

TAVOULARIS, S. 1985 Asymptotic laws for transversely homogeneous turbulent shear flows. *Phys. Fluids* **28**, 999–1001.

TAVOULARIS, S. & CORRSIN, S. 1981 Experiments in nearly homogeneous turbulent shear flow with a uniform mean temperature gradient. Part 1. *J. Fluid Mech.* **104**, 311–347.

TAVOULARIS, S. & KARNIK, U. 1989 Further experiments on the evolution of turbulent stresses and scales in uniformly sheared turbulence. *J. Fluid Mech.* **204**, 457–478.

TOWNSEND, A. A. 1976 *The Structure of Turbulent Shear Flow*, 2nd edn. Cambridge University Press.

WALEFFE, F. 1992 The nature of triad interactions in homogeneous turbulence. *Phys. Fluids* **4**, 350–363.

WEINSTOCK, J. 2013 Analytical theory of homogeneous mean shear turbulence. *J. Fluid Mech.* **727**, 256–281.
YOSHIDA, K., ISHIHARA, T. & KANEDA, Y. 2003 Anisotropic spectrum of homogeneous turbulent shear flow in a Lagrangian renormalized approximation. *Phys. Fluids* **15** (8), 2385–2397.

ZHOU, Y. 2010 Renormalization group theory for fluid and plasma turbulence. *Phys. Rep.* **488**, 1–49.

ZUSI, C. J. & PEROT, J. B. 2013 Simulation and modeling of turbulence subjected to a period of uniform plane strain. *Phys. Fluids* **25**, 110819.

ZUSI, C. J. & PEROT, J. B. 2014 Simulation and modeling of turbulence subjected to a period of axisymmetric contraction or expansion. *Phys. Fluids* **26**, 115103.