3.1. Normalised pressure fluctuations and second moments

A new variable, the normalised pressure fluctuation (denoted as $\hat {\wp }=\textrm{i}\hat {p}/(\gamma \bar {p})$ henceforth), is introduced to simplify the upcoming derivations for the second moments and to remove dependence on imaginary terms. The evolution equation for $\hat {\wp }$ follows (3.7):

(3.10) \begin{align} \frac {{\rm D}\hat {\wp }}{{\rm D}t} =-\kappa _k\hat {u}_k. \end{align}

Next, the second moments must be determined. The velocity spectrum tensor $\Phi _{ij}$ is used to find the one-point Reynolds stress $R_{ij}$ , but its evolution equation is unclosed. We now define the single-point pressure–velocity spectrum vector ( $\beta _i$ ) and pressure spectrum scalar ( $\Gamma$ ) to help close the $\Phi _{ij}$ evolution equation. These additional spectra are defined in a manner similar to that of the velocity spectrum tensor. The symmetric spectra are defined as

(3.11) \begin{align} \Phi _{ij} =\frac {1}{2} \left (\overline { \hat {u}_i\hat {u}^{+}_{\!j}}+\overline { \hat {u}^{+}_i\hat {u}_{\!j}}\right ), \end{align}

(3.12) \begin{align} \beta _{i}=\frac {1}{2}\left (\overline {\hat {\wp }^{+}\hat {u}_i} + \overline {\hat {\wp }\hat {u}_i^{+}}\right ), \end{align}

(3.13) \begin{align} \Gamma =\overline {\hat {\wp }^{+}\hat {\wp }}, \end{align}

where the $+$ superscript represents the complex conjugate. The velocity spectrum tensor $\Phi _{ij}$ is defined as the real components of the ensemble average of the Fourier-transformed velocity fluctuations (Durbin & Pettersson-Reif Reference Durbin and Pettersson-Reif2011). The final first-moment equations are then established by substituting in the normalised pressure fluctuation:

(3.14) \begin{align} \frac {{\rm D}\hat {u}_i}{{\rm D}t}=-\frac {\partial \bar {u}_i}{\partial x_k}\hat {u}_k+\kappa _ia^{2}\hat {\wp }, \end{align}

(3.15) \begin{align} \frac {{\rm D}\kappa _i}{{\rm D}t}=-\kappa _k\frac {\partial \bar {u}_k}{\partial x_i}, \end{align}

(3.16) \begin{align} \frac {{\rm D}\hat {\wp }}{{\rm D}t} =-\kappa _k\hat {u}_k. \end{align}

The appearance of the $\kappa _ia^{2}\hat {\wp }$ term is associated with the dispersion frequency of acoustic waves (Sagaut & Cambon Reference Sagaut and Cambon2018). This introduces new difficulties with the integration of spectral higher-order moments to physical space. With the first-moment evolution equations established, the second-moment equations are derived by substituting in the time derivatives of the first moments and differentiating with respect to time using the product rule:

(3.17) \begin{align} \frac {\rm D}{{\rm D}t}(\Phi _{ij})=-\frac {\partial \bar {u}_i}{\partial x_k}\Phi _{kj}-\frac {\partial \bar {u}_{\!j}}{\partial x_k}\Phi _{ki}+\kappa _ia^{2}\beta _{\!j}+\kappa _{\!j}a^{2}\beta _i, \end{align}

(3.18) \begin{align} \frac {\rm D}{{\rm D}t}(\beta _i)=-\kappa _k\Phi _{ki}-\frac {\partial \bar {u}_i}{\partial x_k}\beta _{k}+\kappa _ia^{2}\Gamma, \end{align}

(3.19) \begin{align} \frac {\rm D}{{\rm D}t}(\Gamma )=-2\kappa _i\beta _i. \end{align}

These final spectral second-moment equations are fully closed and require no new modelling assumption aside from the assumptions made in RDT.

3.2. Transformation to physical space

The final step is to take an inverse Fourier transform to obtain statistics in physical space. This is an integration over all wavenumber components, which is simplified by transforming into spherical coordinates. For example, the evolution equation for the Reynolds stress becomes

(3.20) \begin{align} \frac {\rm D}{{\rm D}t}(R_{ij})= \!\int ^{2\pi}_0 \!\int ^{\pi} _0 \!\int ^{\infty }_{-\infty } \left (-\frac {\partial \bar {u}_i}{\partial x_k}\Phi _{kj}-\frac {\partial \bar {u}_{\!j}}{\partial x_k}\Phi _{ki}+\kappa _ia^{2}\beta _{\!j}+\kappa _{\!j}a^{2}\beta _i\right ) \kappa ^{2}\sin {\theta }\,{\rm d}\kappa \,{\rm d}\theta \,{\rm d}\phi, \end{align}

where $\kappa =\sqrt {\kappa _i\kappa _i}$ and $\kappa _i$ can be decomposed into product of wavevector magnitude $\kappa$ and orientation $\theta, \phi$ . Note that the appearance of $\kappa _i$ means we cannot analytically integrate with respect to the wavevector magnitude as done by Kassinos & Reynolds (Reference Kassinos and Reynolds1994). The solenoidal velocity spectrum tensor, however, can be integrated this way. The latter expression can be solved using a fast Fourier transform to get the most accurate results. However, for this work, we chose to approach the integration through stochastic methods and numerical integration.

3.3. Higher-order tensors

The transformation into Fourier space allows for relatively easier computation of several higher-order tensors. These tensors are used to obtain more information about the structure of turbulence. The structure tensors as defined by Kassinos & Reynolds (Reference Kassinos and Reynolds1996) are of particular interest. The dimensionality tensor is based on the kinetic energy spectrum:

(3.21) \begin{align} D_{ij}=\iiint ^{\infty }_{-\infty }\frac {\kappa _i\kappa _{\!j}}{\kappa ^{2}}\Phi _{mm}\,{\rm d}^{3}{\boldsymbol{\kappa }}, \end{align}

and gives information on whether the turbulence is one- or two-dimensional. The dimensionality gives information on the type of eddies present in the flow. The eddy types are the one-dimensional jettal eddy, the two-dimensional vorticial eddy and the helical eddy (Kassinos & Reynolds Reference Kassinos and Reynolds1994). The structure tensor information is also required to accurately represent rotating flows (Reynolds Reference Reynolds1989). Similarly, the circulicity tensor, which describes the vorticial field, is also found easily in Fourier space:

(3.22) \begin{align} F_{ij}=\epsilon _{inm}\epsilon _{jts}\iiint ^{\infty }_{-\infty }\frac {\kappa _n\kappa _t}{\kappa ^{2}}\Phi _{ms}\,{\rm d}^{3}{\boldsymbol{\kappa }}. \end{align}

In the former equation, $\epsilon _{ijk}$ is the alternating tensor. Finally, higher-rank tensors are also found without requiring additional transport equations. An important fourth-rank tensor,

(3.23) \begin{align} M_{ijpq}=\iiint ^{\infty }_{-\infty }\frac {\kappa _p\kappa _q}{\kappa ^{2}}\Phi _{ij}\,{\rm d}^{3}{\boldsymbol{\kappa }}, \end{align}

gives information regarding the fluctuating pressure field and is often used to represent non-local information as it contains information on two-point statistics (Kassinos & Reynolds Reference Kassinos and Reynolds1994).

3.4. Limitations

The three main assumptions of the current formulation imply limitation to a certain subset of flows. The most substantial assumptions are those of linearisation about a base flow, small density fluctuations compared with the mean and the isentropic assumption. Linearisation is valid only for turbulence subject to rapid deformation. This simplification allows the closure of the equations for the first and second moments without any phenomenological modelling. However, it implicitly prohibits the application of the formulations to slower deformations, as experienced in many engineering flows of interest. We take linearisation as a first step towards more general models, where it can be viewed as a limiting behaviour with the formulation matching predictions from RDT when subject to increasingly large deformations. However, the formulation is expected to gradually become imprecise as deformation decreases.

To represent slower deformations, additional modelling approximations must be added to the formulation and will be explored in upcoming work. The assumption of small density fluctuations also causes discrepancies with application to problems with large compression ratios. The rationale is that as the mean density changes, the magnitude of the fluctuations can become significant. While we do assume that the density fluctuation is small compared with the mean, it is weakly enforced. The simplification is used in the derivation of (3.1) and (3.2), where it is assumed that terms scaled by $\rho '/\bar {\rho }$ are small enough to be ignored when compared in magnitude with the other terms. However, $\rho '/\bar {\rho }$ is not assumed to be zero when considering its rate of change. The isentropic relation uses ${\rm D}(\rho '/\bar {\rho })/{\rm D}t$ to derive the pressure equation, and thus the influence of the density fluctuation with respect to the mean is accounted for in the pressure terms. This weak assumption was similarly applied by Jacquin et al. (Reference Jacquin, Cambon and Blin1993) and Cambon et al. (Reference Cambon, Coleman and Mansour1993) in their studies of compressed homogeneous turbulence.

Finally, the isentropic relation limits the current formulation to adiabatic and entropy-conserving flows. While many compressible flows have non-entropy-conserving features such as shocks, we choose to invoke isentropy due to the significant reduction in equations. If isentropy was not invoked, the pressure equation would need to be derived in a similar manner to that of Yu & Girimaji (Reference Yu and Girimaji2007), where the ideal gas law and energy equation are used to relate temperature and density to pressure evolution. This would result in four more transport equations; two first-rank tensors and two scalars. This amounts to eight more transport equations, which adds to the computational cost. This work seeks a simple formulation which is applicable to most regions where compressible turbulence is important. Non-isentropic parts of the flow usually require additional treatment. For example, turbulence models for flows with shocks require special shock-capturing methods or variable-Prandtl-number ad hoc modifications (Roy et al. Reference Roy, Pathak and Sinha2018) to accurately predict the shock, and non-adiabatic flows such as hypersonic boundary layers require special treatment to account for the inhomogeneity from the wall. We believe that the treatment for various issues that appear in non-isentropic flows can be combined with the required modifications to the current formulation, rather than incorporated in the base formulation.

4.1. Stochastic particle properties and identities

To complete the Monte Carlo integration, the velocity spectrum tensor must be sampled. This provides the motivation for representing the velocity spectrum as a stochastic variable from which we treat samples as Lagrangian particles. We call these stochastic samples particles to follow the terminology of Kassinos & Reynolds (Reference Kassinos and Reynolds1994).

Figure 1. Particle properties associated with a cluster of samples.

Under the assumptions previously stated, the stochastic representation is valid for the turbulence considered in this study. We can now formally define the stochastic variables. Suppose there exist multivariate distributions in ${\boldsymbol{\kappa }}$ -space for $\Phi _{ij}({\boldsymbol{\kappa }},t),\ \beta _i({\boldsymbol{\kappa }},t)$ and $\Gamma ({\boldsymbol{\kappa }},t)$ such that ${\mathcal{K}} \subset {{\mathbb{R}}^{3}}$ is the state space of the wavevector. We can obtain samples from ${{\boldsymbol{\kappa }}^{(1)}},{{\boldsymbol{\kappa}}^{(2)}},\ldots, {{\boldsymbol{\kappa}}^{(n)}} \in {\mathcal{K}}$ , where the $(n)$ superscript denotes sample number. The second moments defined previously correspond to realisations of the underlying distributions. Now, we sample (for a range of ${\boldsymbol{\kappa }}$ vectors) from the initial random distributions for the three second-moment probability density functions (PDFs) and advance their solutions forward in time using the evolution equations for the second moments. Note that the stochastic behaviour comes from the initialisation rather than from the evolution equations. Each realisation for $\Phi _{ij}, \beta _i$ , $\kappa _i$ and $\Gamma$ is evolved in time separately.

Particles are assigned properties which are evolved through time and are used to compute second-moment statistics. The assigned properties are the particle velocity vector $\textbf{u}^{*}$ (which can be decomposed into solenoidal and dilatational components), the particle wavevector ${\boldsymbol{\kappa }}^{*}$ and the particle normalised pressure $\wp ^{*}$ . These stochastic particles represent a one-dimensional one-component flow, a basic building block for the turbulence structure. The particle velocity vector, wavevector and normalised pressure can be thought of as stochastic samples of the respective Fourier modes. Note that compared with the incompressible PRM, the particle wavevector and velocity vector are not orthogonal since the velocity field is not divergence-free. However, the solenoidal velocity vector $\textbf{u}^{*s}$ is orthogonal to ${\boldsymbol{\kappa }}^{*}$ and the dilatational velocity vector $\textbf{u}^{*d}$ is parallel to ${\boldsymbol{\kappa }}^{*}$ such that $u_i^{*s}\kappa _i=0$ and $u^{*}_i=u^{*s}_i+u^{*d}_i$ . This is illustrated in figure 1. In addition, the particle wavevector is similar but not identical to the gradient vector defined by Kassinos & Reynolds (Reference Kassinos and Reynolds1994). Since the stochastic particles are analogous to the Fourier modes of the flow variables, the stochastic governing equations are also analogous to the RDT equations. The first-moment equations are

(4.6) \begin{align} \frac {{\rm D}u^{*}_i}{{\rm D}t}=-\frac {\partial \bar {u}_i}{\partial x_k}u^{*}_k+\kappa ^{*}_ia^2\wp ^{*}, \end{align}

(4.7) \begin{align} \frac {{\rm D}\kappa ^{*}_i}{{\rm D}t}=-\kappa ^{*}_k\frac {\partial \bar {u}_k}{\partial x_i}, \end{align}

(4.8) \begin{align} \frac {{\rm D}\wp ^{*}}{{\rm D}t} =-\kappa ^{*}_ku^{*}_k. \\[-12pt] \nonumber \end{align}

The closed set of second-moment stochastic equations are

(4.9) \begin{align} \frac {\rm D}{{\rm D}t}(\Phi ^{*}_{ij})=-\frac {\partial \bar {u}_i}{\partial x_k}\Phi ^{*}_{kj}-\frac {\partial \bar {u}_{\!j}}{\partial x_k}\Phi ^{*}_{ki}+\kappa ^{*}_ia^2\beta ^{*}_{\!j}+\kappa ^{*}_{\!j}a^2\beta ^{*}_i, \end{align}

(4.10) \begin{align} \frac {\rm D}{{\rm D}t}(\beta ^{*}_i)=-\kappa ^{*}_k\Phi ^{*}_{ki}-\frac {\partial \bar {u}_i}{\partial x_k}\beta ^{*}_{k}+\kappa ^{*}_ia^2\Gamma ^{*}, \end{align}

(4.11) \begin{align} \frac {\rm D}{{\rm D}t}(\Gamma ^{*})= -2\kappa ^{*}_i\beta ^{*}_i. \\[-12pt] \nonumber \end{align}

4.2. Equivalency with spectral representation

This section presents a more rigorous justification to show equivalency with spectral and stochastic representations. We turn to the PDF formulation of compressible RDT. While PDF methods are typically used to model the one-point, one-time Navier–Stokes equations (Pope Reference Pope2001), they also provide a more grounded justification for using a stochastic point of view for turbulence. In this section, the marginal PDF is used to relate the random stochastic variables to the Fourier spectra. To show equivalency and consistency between the PDF and traditional formulations, the stochastic velocity spectrum tensor must have the same evolution equation form as the Fourier velocity spectrum:

(4.12) \begin{align} \overline {u'_iu'_{\!j}}=\iiint ^{\infty }_{-\infty }\Phi _{ij}({\boldsymbol{\kappa }}){\rm d}^{3}{\boldsymbol{\kappa }}\equiv \iiint ^{\infty }_{-\infty }\Phi ^{*}_{ij}({\boldsymbol{\mathcal{K}}}){\rm d}^{3}{\boldsymbol{\mathcal{K}}}=\overline {u^{*}_iu^{*}_{\!j}}. \end{align}

The stochastic velocity spectrum is defined similarly to that in Van Slooten & Pope (Reference Van Slooten and Pope1997), where it is equal to the integral of the product of the velocity state space and joint PDF $f^{*}({\textbf{V}},{\boldsymbol{\mathcal{K}}})$ . However, since the velocity also depends on the normalised pressure fluctuation $\wp$ , the full state space of the joint PDF contains a total of nine variables (two vectors and a scalar), $f^{*}({\textbf{V}},{\boldsymbol{\mathcal{K}}},{\mathcal{P}})$ , for which the normalised pressure state space ${\mathcal{P}}$ must first be integrated out to give the joint PDF dependent only on velocity and wavevector states. The stochastic pressure–velocity and pressure spectra are defined in a similar fashion. The stochastic second-order moments are

(4.13) \begin{align} \Phi ^{*}_{ij}({\boldsymbol{\mathcal{K}}})=\iiiint ^{\infty }_{-\infty } V_iV_{\!j}f^{*}({\textbf{V}},{\boldsymbol{\mathcal{K}}},{\mathcal{P}}){\rm d}{\mathcal{P}}\,{\rm d}^{3}{\textbf{V}}=\langle u_i^{*}u_{\!j}^{*}|{\boldsymbol{\kappa }}^{*}={\boldsymbol{\mathcal{K}}}\rangle f^{*}({\boldsymbol{\mathcal{K}}}), \end{align}

(4.14) \begin{align} \beta ^{*}_{i}({\boldsymbol{\mathcal{K}}})=\iiiint ^{\infty }_{-\infty } V_i{\mathcal{P}}f^{*}({\textbf{V}},{\boldsymbol{\mathcal{K}}},{\mathcal{P}}){\rm d}{\mathcal{P}}\, {\rm d}^{3}{\textbf{V}}=\langle u_i^{*}\wp ^{*}|{\boldsymbol{\kappa }}^{*}={\boldsymbol{\mathcal{K}}}\rangle f^{*}({\boldsymbol{\mathcal{K}}}), \end{align}

(4.15) \begin{align} \Gamma ^{*}_{i}({\boldsymbol{\mathcal{K}}})=\iiiint ^{\infty }_{-\infty } {\mathcal{P}}^2f^{*}({\textbf{V}},{\boldsymbol{\mathcal{K}}},{\mathcal{P}})\,{\rm d}\mathcal{P}\,{\rm d}^{3}{\textbf{V}} =\langle \wp ^{*}\wp ^{*}|{\boldsymbol{\kappa }}^{*}={\boldsymbol{\mathcal{K}}}\rangle f^{*}({\boldsymbol{\mathcal{K}}}). \end{align}

The joint PDF $f^{*}({\textbf{V}},{\boldsymbol{\mathcal{K}}},{\mathcal{P}})$ is derived using traditional methods (Heinz Reference Heinz2003) and is based on the stochastic first-moment equations (equations (4.6)–(4.8)):

(4.16) \begin{align} \frac {\partial f^{*}}{\partial t}=\frac {\partial \bar {u}_k}{\partial x_n}\frac {\partial }{\partial \mathcal{K}_n}(\mathcal{K}_kf^{*})+\frac {\partial \bar {u}_n}{\partial x_k}\frac {\partial }{\partial V_n}(V_kf^{*})+\frac {\partial }{\partial V_n}(a^2 \mathcal{K}_n {\mathcal{P}} f^{*})-\frac {\partial }{\partial {\mathcal{P}}}(\mathcal{K}_kV_kf^{*}) .\end{align}

Since an Eulerian PDF is used, the stochastic velocity spectrum evolution equation must match the non-conservative Eulerian velocity spectrum evolution equation (transformed from equation (3.17)):

(4.17) \begin{align} \frac {\partial \Phi _{ij}}{\partial t}=\kappa _m\frac {\partial \bar {u}_m}{\partial x_n}\frac {\partial \Phi _{ij}}{\partial \kappa _n}+\frac {\partial \bar {u}_k}{\partial x_k}\Phi _{ij}-\frac {\partial \bar {u}_{\!j}}{\partial x_k}\Phi _{ki}-\frac {\partial \bar {u}_i}{\partial x_k}\Phi _{kj}-\kappa _ia^2\beta _{\!j}-\kappa _{\!j}a^2\beta _i. \end{align}

Now, the derivative of (4.13) with respect to time is taken to find the stochastic evolution equation. Note that the state space variables ${\textbf{V}}, {\mathcal{P}},{\boldsymbol{\mathcal{K}}}$ are constant in time, but the random variables sampled from them are not. Substituting in (4.16) into the time derivative of (4.13) and integrating each term yields the stochastic velocity spectrum transport equation. Integration by parts is used along with linearity to simplify the solution for each term. The first two terms only depend on integrating with respect to ${\textbf{V}}$ since no $\mathcal{P}$ terms appear:

(4.18) \begin{align} \iiint ^{\infty }_{-\infty } V_iV_{\!j}\frac {\partial \bar {u}_k}{\partial x_n}\frac {\partial }{\partial \mathcal{K}_n}(\mathcal{K}_kf^{*}){\rm d}^{3}{\textbf{V}}=\frac {\partial \bar {u}_k}{\partial x_n}\frac {\partial \Phi ^{*}_{ij}}{\partial \mathcal{K}_n}\mathcal{K}_k+\frac {\partial \bar {u}_k}{\partial x_k}\Phi ^{*}_{ij}, \end{align}

(4.19) \begin{align} \iiint ^{\infty }_{-\infty } V_iV_{\!j}\frac {\partial \bar {u}_n}{\partial x_k}\frac {\partial }{\partial V_n}(V_kf^{*}){\rm d}^{3}{\textbf{V}}=-\frac {\partial \bar {u}_i}{\partial x_k}\Phi ^{*}_{kj}-\frac {\partial \bar {u}_{\!j}}{\partial x_k}\Phi ^{*}_{ki}. \\[-2pt] \nonumber \end{align}

The third term requires four integration steps. Integration is first carried out on $\mathcal{P}$ and then on ${\textbf{V}}$ . Terms are simplified using the stochastic second-moment definitions (equations (4.13)–(4.15)):

(4.20) \begin{eqnarray} &&\def\lumina\hspace*{-1.8pc} \iiiint ^{\infty }_{-\infty } V_iV_{\!j}\frac {\partial }{\partial V_n}(a^2 \mathcal{K}_n {\mathcal{P}} f^{*}){\rm d}\mathcal{P} \,{\rm d}^{3}{\textbf{V}} \nonumber\\[4pt]&& \qquad =V_iV_{\!j}a^2 \mathcal{K}_n {\mathcal{P}} f^{*}-\iiiint ^{\infty }_{-\infty } (\delta _{ni}V_{\!j}+\delta _{nj}V_i)a^2 \mathcal{K}_n {\mathcal{P}} f^{*}{\rm d}\mathcal{P}\, {\rm d}^{3}{\textbf{V}} \nonumber \\[4pt]&& \qquad =-a^2 \mathcal{K}_i \beta ^{*}_{\!j} -a^2 \mathcal{K}_{\!j}\beta ^{*}_i. \end{eqnarray}

Finally, the fourth term also requires four integration steps. Due to linearity and independence of $\mathcal{P}$ and ${\textbf{V}}$ , the last term is zero:

(4.21) \begin{eqnarray} &&\def\lumina\hspace*{-1.8pc} \iiiint ^{\infty }_{-\infty } -V_iV_{\!j}\frac {\partial }{\partial {\mathcal{P}}}(\mathcal{K}_kV_kf^{*}){\rm d}\mathcal{P}\, {\rm d}^{3}{\textbf{V}} =-\iiint ^{\infty }_{-\infty } V_iV_{\!j}\mathcal{K}_kV_kf^{*}\,{\rm d}^{3}{\textbf{V}} \nonumber\\[4pt]&&\qquad + \iiiint ^{\infty }_{-\infty } \frac {\partial }{\partial {\mathcal{P}}}(V_iV_{\!j}) \mathcal{K}_kV_kf^{*}{\rm d}\mathcal{P} \,{\rm d}^{3}{\textbf{V}}=0.\end{eqnarray}

Putting solutions (4.18)–(4.21) together gives the evolution equation for the stochastic velocity spectrum tensor from an Eulerian frame of reference:

(4.22) \begin{align} \frac {\partial \Phi ^{*}_{ij}}{\partial t}=\mathcal{K}_m\frac {\partial \bar {u}_m}{\partial x_n}\frac {\partial \Phi ^{*}_{ij}}{\partial \mathcal{K}_n}+\frac {\partial \bar {u}_k}{\partial x_k}\Phi ^{*}_{ij}-\frac {\partial \bar {u}_{\!j}}{\partial x_k}\Phi ^{*}_{ki}-\frac {\partial \bar {u}_i}{\partial x_k}\Phi ^{*}_{kj}-\mathcal{K}_ia^2\beta ^{*}_{\!j}-\mathcal{K}_{\!j}a^2\beta ^{*}_i. \end{align}

Equations (4.22) and (4.17) show that the two equations are identical. Since the two are identical, the stochastic representation is consistent with the Fourier representation and thus gives the same results when computing turbulent statistics in physical space.

4.3. Particle clusters

The integration with respect to wavevector magnitude presents additional difficulties. Consider a type of particle cluster where particles are grouped together such that they all have the same orientation but different wavevector magnitudes and spectra, as shown in figure 1. Our idea of clusters is nearly identical to that of Kassinos & Reynolds (Reference Kassinos and Reynolds1994), but rather than possessing random directional spectra, the clusters have a specified range of wavevector magnitudes associated with them, which directly changes the spectra. The motivation for clustering is to reduce the number of samples required to compute accurate statistics and separate the particle information on wavevector magnitude and orientation. Clustering also enables the use of non-stochastic methods for the inner integral:

(4.23) \begin{align} R_{ij}\approx \frac {1}{N_c}\sum ^{N_c}_{n=1}4\pi \int _0^{\infty } \big(\kappa ^{(sh))}\big)^2\Phi _{ij}\big(\kappa ^{(sh)},\theta ^{(c)},\phi ^{(c)}\big){\rm d}\kappa, \end{align}

where the $(c)$ superscript denotes the cluster index. Generally, the wavevector magnitude evolution depends on the wavevector orientation. However, because of the wide magnitude integration limits, the evolution of the wavevector magnitude sometimes plays a negligible role. Therefore, by splitting up the integration steps, we can explicitly ignore the magnitude evolution and only need to consider the evolution of the wavevector orientations. The evolution equation for the unit wavevector $\kappa _i/k=e_i$ is

(4.24) \begin{align} \frac {{\rm D}e_i}{{\rm D}t}=-e_k\frac {\partial \bar {u}_k}{\partial x_i}+\frac {\partial \bar {u}_k}{\partial x_r}e_ke_re_i. \end{align}

By only considering the orientation evolution, we are in essence resampling the stochastic particles at each time step such that they fall within the predefined wavevector magnitude range. The integral in equation (4.23) can be computed in various ways, such as trapezoidal, quadrature, etc. In this work, we use the trapezoidal and Gauss–Laguerre quadrature integration methods. Note that this resampling technique only works for some cases. For example, it was found that the pure shear case requires the evolution magnitude for stability reasons.

5.1. Initialisation

The initialisation method will change depending on whether the velocity spectrum tensor or directional velocity spectrum tensor is being evolved. The simpler initialisation is the velocity spectrum tensor. Turbulence is initialised from homogeneous isotropic turbulence (HIT), unless otherwise stated. The definition for the velocity spectrum tensor in HIT is given by Durbin & Pettersson-Reif (Reference Durbin and Pettersson-Reif2011) in terms of the energy spectrum:

(5.1) \begin{align} \Phi _{ij}\left({\boldsymbol{\kappa }^{(n)}}\right)=\frac {E\left(\kappa ^{(n)}\right)}{4\pi \left(\kappa ^{(n)}\right)^2}\left (\delta _{ij}-\frac {\kappa ^{(n)}_i\kappa ^{(n)}_{\!j}}{(\kappa ^{(n)})^2}\right ), \end{align}

where $n$ denotes the particle number and $E(K)$ is the energy spectrum. Note that is initialisation is only valid for solenoidal turbulence. The velocity spectrum tensor is assigned to each $n{\rm th}$ particle given its wavevector orientation and magnitude. The initial energy spectrum is specified based on desired initial conditions. In the most general cases, the Von Kármán energy spectrum is used. It is given as (Durbin & Pettersson-Reif Reference Durbin and Pettersson-Reif2011)

(5.2) \begin{align} E(\kappa )=q^2LC_{vK}\frac {(\kappa L)^4}{[1+(\kappa L)^2]^p}, \end{align}

where $L$ is the length scale required for dimensional consistency, $q^2=R_{ii}$ and $p=17/6$ to match the $-5/3$ law. Here $C_{vK}$ is found by enforcing

(5.3) \begin{align} \int _{-\infty }^{\infty }E(k){\rm d}\kappa =\frac {1}{2}q^2, \end{align}

such that $C_{vK}\approx 0.4843$ . When comparing results with RDT/DNS data from Simone et al. (Reference Simone, Coleman and Cambon1997) or Cambon et al. (Reference Cambon, Coleman and Mansour1993), we used the same initial energy spectrum defined in those papers:

(5.4) \begin{align} E(\kappa )=q^2A\kappa ^4\exp {\frac {-2\kappa ^2}{\kappa _p^2}}, \end{align}

where $\kappa _p=8$ is the wavevector magnitude where the energy spectrum peaks and $A$ is a constant such that equation (5.3) is satisfied. To evolve the velocity spectrum tensor, the unit wavevector is initialised randomly about a unit sphere for a specified range of $k$ magnitudes. The number of unique orientations is equal to the cluster sample count, while the number of $\kappa$ discretisations is based on the integration method (logarithmic spacing with $\Delta \kappa$ total range for trapezoidal, for example). To initialise the directional velocity spectrum tensor, only the unit vectors in a unit sphere must be initialised. This is also done randomly. Each unit vector has a unique directional velocity spectrum tensor value such that

(5.5) \begin{align} \frac {1}{N}\sum ^N_{c=1}4\pi \Phi ^{\nearrow \kappa }_{kk}\big(\textbf{e}^{(c)}\big)=\frac {1}{2}q^2. \end{align}

This means that the diagonal values of the Reynolds stress are equal to $({1}/{3})q^2$ . The HIT initialisation for the directional velocity spectrum tensor is then

(5.6) \begin{align} 4\pi \Phi ^{\nearrow \kappa }_{ij}\big({\boldsymbol{\kappa }}^{(c)}/\kappa \big)=\left (\delta _{ij}-\frac {\kappa ^{(c)}_i\kappa ^{(c)}_{\!j}}{\kappa ^2}\right )\int ^\infty _{-\infty }E\big(\kappa ^{(c)}\big){\rm d}\kappa . \end{align}

The energy spectrum integral is numerically evaluated using midpoint integration approximation.

Initialisation of the pressure–velocity and pressure spectra will require knowing their respective energy spectra. Generally, these are not known. For the examples in this work, the velocity–pressure spectrum is set to zero for HIT. However, following the initialisation of the solenoidal/dilatational modes in Simone et al. (Reference Simone, Coleman and Cambon1997), the strong form of acoustic equilibrium (Sarkar et al. Reference Sarkar, Erlebacher, Hussaini and Kreiss1991) is enforced when the pressure spectrum is defined by the dilatational component of the energy spectrum. The ratio of solenoidal to dilatational turbulent kinetic energy $q^2_s/q^2_d$ must first be specified and used to scale the energy spectrum $E(k)$ , such that $E_s(\kappa )+E_d(k)=E(\kappa )$ (where $E_s(\kappa )/E_d(\kappa )=q^2_s/q^2_d$ ), to obtain the following initialisation:

(5.7) \begin{align} \Gamma \big({\boldsymbol{\kappa }^{(n)}}\big)=\frac {E_d\left(\kappa ^{(n)}\right)}{2\pi \left(\kappa ^{(n)}\right)^2\! a_0^2}. \end{align}

The influence of solenoidal and dilatational effects depends on the distortion Mach number $M_d$ (Durbin & Zeman Reference Durbin and Zeman1992). For each case, an initial distortion Mach number $M_d=M_tSq^2/\epsilon$ is specified along with a turbulent Mach number $M_t=q/a$ and turbulent Reynolds number $Re_{t}=q^4/(\nu \epsilon )$ . These three conditions along with the relationship between the total turbulent dissipation rate and energy spectrum,

(5.8) \begin{align} \epsilon =2\nu \int _0^{\infty }\kappa ^2E(\kappa ){\rm d}\kappa, \end{align}

are used to initialise the stochastic spectra. Finally, it is noted that for the pure shear and shear compression case, $\gamma =1.4$ is used whereas $\gamma =5/3$ is used for the axisymmetric compression. All cases use a constant initial mean temperature of $\bar {T}=300$ . The temperature is not assumed to remain constant throughout the simulation, and its change is captured through the change in pressure and speed of sound. The density and temperature fluctuations do not explicitly appear in the evolution equations and thus are not required to be explicitly initialised. However, they can be initialised through the isentropic relations between density, temperature and pressure.

5.2. Algorithm for stochastic representation

All the background necessary for a functional algorithm to evolve the turbulent statistics of compressible turbulence has been discussed in the prior sections. The complete algorithm is now presented in the following pseudo-code to further solidify methodology. Depending on the case, the following variables are known/chosen: the final non-dimensional time $St_f$ , the number of particle clusters $n_{p}$ , the number of wavevector modulus discretisations $n_k$ , the number of time steps $n_t$ , the initial mean speed of sound $a_0$ , the initial turbulent Mach number $M_{t0}$ , the initial distortion Mach number $M_{d0}$ , the initial turbulent Reynolds number $Re_{t0}$ , the initial solenoidal-to-dilatational turbulent kinetic energy ratio $q_s^2/q_d^2$ , the initial shear-to-compression ratio $S_0/C_0$ and the maximum and minimum wavevector magnitudes $\kappa _{{min}}, \kappa _{{max}}$ . The governing equations are solved numerically using finite time discretisation to advance the solution forward in time. We used an explicit fourth-order Runge–Kutta scheme.

Algorithm 1 Compressible PRM

5.3. Axial compression

The axial compression considered in this case is defined by

(5.9) \begin{align} G_{ij}=C\delta _{i1}\delta _{j1}, \quad C(t)=\frac {C_0}{1+C_0t}=\frac {C_0\rho }{\rho _0}, \end{align}

where the compression axis is in the $x_1$ direction (as shown in figure 2 a) and $C$ represents the compression rate. Due to the change in volume/density (varying from $\rho /\rho _0=1$ to $\rho /\rho _0=5$ ), the speed of sound also changes. The relationship between the speed of sound and density ratio is derived from isentropic relations (Mahesh Reference Mahesh1996, Jacquin et al. Reference Jacquin, Cambon and Blin1993):

Figure 2. Schematic of the mean velocity deformations for (a) axial compression, (b) pure shear and (c) sheared compression.

(5.10) \begin{align} a=a_0\sqrt {\frac {1+C_0 t}{(1+C_0t)^\gamma }}. \end{align}

As mentioned previously, the distortion Mach number $M_d$ is equal to the ratio between the initial acoustic time scale and the mean distortion or compression time scale. So if $M_d\ll 1$ , the acoustic time scale is much larger than the mean distortion time scale and a pure acoustic/solenoidal regime is obtained (Durbin & Zeman Reference Durbin and Zeman1992). In this solenoidal regime, the dilatational mode is balanced by the pressure gradient. The pressure is found with the pressure Poisson equation and consequentially describes the dilatational modes (Simone et al. Reference Simone, Coleman and Cambon1997). The energy is bounded in this regime and limits are well established (Durbin & Zeman Reference Durbin and Zeman1992; Jacquin et al. Reference Jacquin, Cambon and Blin1993; Simone et al. Reference Simone, Coleman and Cambon1997). On the other hand, if $M_d\gg 1$ , the dilatational mode is no longer constrained by the pressure (hence the name pressure-released) and is fed directly by compression. For the $M_d\gg 1$ pressure-released limit, the turbulent kinetic energy amplification can be expressed as (Jacquin et al. Reference Jacquin, Cambon and Blin1993)

(5.11) \begin{align} \frac {q^2(t)}{q^2(0)}=\frac {2+(\rho (t)/\rho (0))^2}{3}. \end{align}

For the solenoidal limit $M_d\ll 1$ , turbulent kinetic energy amplification is also found analytically using (Simone et al. Reference Simone, Coleman and Cambon1997)

(5.12) \begin{align} \frac {q^2(t)}{q^2(0)}=\frac {1}{2}\left (1+(\rho (t)/\rho (0))^2\frac {\tan ^{-1}\left (\sqrt {(\rho (t)/\rho (0))^2-1}\right )}{\sqrt {(\rho (t)/\rho (0))^2-1}}\right ). \end{align}

Results for the axial compression cases are compared with DNS data of Cambon et al. (Reference Cambon, Coleman and Mansour1993). The initial conditions used for the three intermediate cases are listed in table 1. Note that the distortion Mach number is defined with the negative compression rate for this problem, $M_d=-M_tCq^2/\epsilon$ . Unless otherwise stated, results presented use 800 particles and 60 discretisations of the wavevector magnitude for a range of $10^{-3}$ to $10^{3}$ . The trapezoidal integration with respect to magnitude was used alongside Monte Carlo integration with respect to orientation. Figure 3 shows the histories of the turbulent kinetic energy amplification for cases A–C with different cluster amounts, varying from 200 to 800. With increasing cluster count, predictions tend to improve and agree closer with the DNS data. The optimal cluster count varies based on the distortion Mach number. We recommend a starting point of 500 clusters for most cases. Figure 3 also shows that increased compressibility ( $M_d$ ) leads to an increase in the kinetic energy amplification. Both solenoidal and dilatational contributions to the kinetic energy are shown in figure 4 and compared with DNS results. The solenoidal mode is nearly unaffected by the increased compression, which corroborates with the theory in Cambon et al. (Reference Cambon, Coleman and Mansour1993). A strong increase of the dilatational energy at the end of the compression is also captured by the CPRM.

Table 1. Initial conditions for axial compression cases.

Figure 3. Histories of turbulent kinetic energy amplification for three cases of varying distortion Mach number in axial compression: $M_{d0}=5$ , black; $M_{d0}=29$ , blue; $M_{d0}=87$ , red; DNS of Cambon et al. (Reference Cambon, Coleman and Mansour1993), circles. Multiple predictions are shown for varying particle cluster counts varying from 200 to 800 clusters.

Figure 4. Evolution of solenoidal and dilatational components of turbulent kinetic energy amplification for three cases of varying distortion Mach number in axial compression: $M_{d0}=5$ , black; $M_{d0}=29$ , blue; $M_{d0}=87$ , red; DNS of Cambon et al. (Reference Cambon, Coleman and Mansour1993), circles.

Figure 5. Evolution of solenoidal and dilatational components of the $D_{11}/q^2$ structure tensor component for three cases of varying distortion Mach number in axial compression: $M_{d0}=5$ , black; $M_{d0}=29$ , blue; $M_{d0}=87$ , red; DNS of Cambon et al. (Reference Cambon, Coleman and Mansour1993), circles.

The dimensionality structure tensor $D_{11}$ component is also computed and compared with DNS data. Figure 5 shows the histories of normalised solenoidal and dilatational components of $D_{11}$ and matches the trends of the solenoidal/dilatational components of the DNS data. Note that for axial compression, $D^d_{11}=R_{11}^d$ (Cambon et al. Reference Cambon, Coleman and Mansour1993). The dimensionality provides information on the angular distribution of the spectral energy. In the pressure-released case, the alignment of the wavevector with the compression direction increases $D_{11}/q^2$ towards one, while the angular distribution of the spectral energy remains unchanged. Meanwhile, in the solenoidal limit, the angular distribution of spectral energy in the plane normal to the compression direction opposes the wavevector motion so that a slower positive increase of $D_{11}/q^2$ is displayed. The $M_{d0}=5$ case has poor $D_{11}^d/q^2_d$ prediction, showing that the CPRM does not capture this spectral distribution opposing the wavevector motion. Reasons for this disagreement include wavevector discretisation errors and the resampling technique inaccuracy for the wavevector magnitude, which remains constant due to the infinite integration. In contrast, when evolving the wavevector magnitude, the dimensionality results are over-predicted, but the case near the solenoidal limit was much closer.

While the CPRM is expected to compare well with RDT, the results show some discrepancies between the high-fidelity data. One reason for this is that the comparing data in this section are from the DNS of Cambon et al. (Reference Cambon, Coleman and Mansour1993). We do not expect RDT to match exactly with DNS, and thus small disagreement is reasonable. Cambon et al. (Reference Cambon, Coleman and Mansour1993) also stated that the solenoidal ratio of $D_{11}/q^2$ predicted by DNS was found to be slightly lower than the RDT analytical prediction. This finding is also supported by the results for the solenoidal CPRM $D_{11}/q^2$ predictions. Finally, the CPRM requires different information for initialisation from what is given for DNS. The pressure–velocity spectrum $\beta _i$ was initialised to zero for this case, while the pressure spectrum $\Gamma$ was initialised with the dilatational energy. These initialisations are approximate as the energy spectra for the pressure–velocity and pressure spectra are not known. Therefore, since the evolution of the turbulence statistics depends on initial conditions, the inexact initial conditions contribute to the error seen in the results.

It is worth noting the differences and analogies of the current CPRM formulation to the traditional RDT used by Jacquin et al. (Reference Jacquin, Cambon and Blin1993). The CPRM uses the same set of governing equations and assumptions as Jacquin et al. (Reference Jacquin, Cambon and Blin1993), but ends up with a different and reduced set of evolution equations. This is because the CPRM only used deforming coordinates/Lagrangian perspective as the coordinate transformation whereas the RDT of Jacquin et al. (Reference Jacquin, Cambon and Blin1993) used deforming coordinates and a Craya–Herring transformation. The CPRM also solves for the full-field variables and does not require the Helmholtz decomposition as compared with the RDT which was developed using this decomposition. The normalised pressure variable $\hat {\wp }$ removes the explicit appearance of i in the evolution equation and so the CPRM can be solved without any need to handle imaginary components, whereas the simplified RDT equations presented by Jacquin et al. (Reference Jacquin, Cambon and Blin1993) did have an explicit imaginary number appear in the pressure terms for the general (equations (30) and (32)) but not in the reduced equations for the axial compression case. The last difference is the method of numerical integration. Both the CPRM and RDT of Jacquin et al. (Reference Jacquin, Cambon and Blin1993) require numerical integration, but the CPRM is solved using stochastic methods while the numerical integration method is not specified by Jacquin et al. (Reference Jacquin, Cambon and Blin1993), but likely used some form of fast Fourier transforms.

5.3.1. Axial compression: a surrogate for shock–turbulence interaction

Axial compression can be used as a simplified model for shock–turbulence interaction. Jacquin et al. (Reference Jacquin, Cambon and Blin1993) used this model problem to compare solenoidal and pressure-released limits against the LIA theory of Ribner (Reference Ribner1954). Homogeneous axially compressed turbulence is only a model of shock–turbulence interaction because true shock–turbulence interaction does not meet the quasi-homogeneous requirement that $\ell \ll \delta$ , where $\delta$ is the mean scale (the thickness of the shock). For further evaluation of the CPRM, we follow the canonical problem of Jacquin et al. (Reference Jacquin, Cambon and Blin1993), where an initial flow at Mach number $M_0$ impinges a shock of strength $\Delta \bar {u}/\bar {u}_0=(\rho -\rho _0)/\rho$ , such that the density ratio is given by

(5.13) \begin{align} \frac {\rho }{\rho _0}=\frac {(\gamma +1)M_0^2}{2+(\gamma -1)M_0^2}. \end{align}

Figure 6 shows the turbulent kinetic energy amplification for varying $M_0$ . The disagreement with LIA (computed by Jacquin et al. (Reference Jacquin, Cambon and Blin1993)) is significant, and follows the same results obtained by Jacquin et al. (Reference Jacquin, Cambon and Blin1993). The CPRM performance in figure 6 is identical to figure 3 because the problem dynamics is the same, just with increased compression rate $C$ (where increasing $C$ approaches the pressure-released solution). The CPRM does match the far-field LIA at very short times/low $M_0$ , but does not match the near-field LIA predictions.

Figure 6. Evolution of turbulent kinetic energy amplification for three intermediate cases in axial compression compared with LIA and limiting RDT: $M_{d0}=5$ , black; $M_{d0}=29$ , blue; $M_{d0}=87$ , red.

5.4. Pure shear

The deformation for pure shear is shown in figure 2(b). Pure shear is defined as

(5.14) \begin{align} G_{ij}=S\delta _{i1}\delta _{j2}, \quad\; S(t)=S_0, \end{align}

where $S_0$ is the initial shear rate, which remains constant. No volumetric compression is experienced ( $\rho /\rho _0=1$ ) and thus the speed of sound and mean pressure/density are constant in time: $a=a_0$ , $\bar {p}(t)=\bar {p}_0$ , $\bar {\rho }(t)=\bar {\rho }$ . For anisotropies $b_{ij}=R_{ij}/q^2-\delta _{ij}/3$ , the pressure-released limit is found by neglecting the pressure terms in the evolution equation for $\Phi ^{*}_{ij}$ . The evolution of the turbulent kinetic energy amplification in this limit is, however, found analytically when the initial turbulence conditions are isotropic (Simone et al. Reference Simone, Coleman and Cambon1997):

(5.15) \begin{align} \frac {q^2(t)}{q^2(0)}=1+\frac {(St)^2}{3}. \end{align}

Figure 7. Turbulent kinetic energy amplification evolution in pure shear for varying distortion Mach numbers prescribed in table 2. Solid line given by equation (5.15). Results are comparable with those of Bertsch et al. (Reference Bertsch, Suman and Girimaji2012) and Lavin et al. (Reference Lavin, Girimaji, Suman and Yu2012).

Simple analytical solutions only exist for the pressure-released limit. Solenoidal limits are evaluated by decreasing the distortion Mach number $M_d\ll 1$ and using the solenoidal form of the evolution equations for $\Phi ^s_{ij}$ . Cluster count varied from 500 near the solenoidal limit to 700 near the pressure-released limit. A 20-point trapezoidal integration was used for the wavevector magnitude. Due to instabilities at high $St$ values, the evolution of the wavevector magnitude was taken into account in the evolution equations for the spectra (no resampling used). A very fine time step was required at high $St$ values, reaching as low as $\Delta t/t_f=0.001$ .

Evolution of the turbulent kinetic energy amplification is presented in figure 7. As expected, the turbulent kinetic energy amplification is bounded by the pressure-released limit and increases with initial $M_d$ . Amplification increase is much greater than that of the DNS data from Simone et al. (Reference Simone, Coleman and Cambon1997). The early-time monotonic increase is predicted by the CPRM but the later-time linear behaviour does not match the DNS. This disagreement is likely due to the re-meshing technique used in the DNS that resulted in loss of turbulent kinetic energy and thus lower growth rates. However, the CPRM does predict the three-stage evolution of the turbulent kinetic energy amplification as described by Bertsch et al. (Reference Bertsch, Suman and Girimaji2012), Lavin et al. (Reference Lavin, Girimaji, Suman and Yu2012) and Kumar et al. (Reference Kumar, Bertsch and Girimaji2014). The intermediate flattening of the turbulent kinetic energy evolution at $St \approx 4{-}6$ is seen in figure 7. In the first stage, irrespective of the initial modal Mach number $M_m=S/a|{\boldsymbol{\kappa }}|$ , rapid growth of the turbulent kinetic energy is experienced and matches the rate of the pressure-released solution. This stage is defined by $St/\sqrt {M_m(0)}\lt 1$ , where the fluctuating pressure change is too slow to influence the velocity field, which grows unbounded. Transition from the first to the second stage begins when the pressure starts to influence the velocity field. Through this influence, the kinetic energy evolution begins to depart from the pressure-released solution. We see a sudden reduction/stabilisation in the kinetic energy amplification growth in figure 7 around $St\approx 1.5$ . The flattening of kinetic energy is bounded by $1\lt M_m(t)\lt \sqrt {M_m(0)}$ . The third stage of kinetic energy evolution begins where the kinetic energy grows rapidly again. This occurs when the acoustic time scale $a_0t|\boldsymbol{\kappa _0}|$ becomes smaller than the mean shear time scale $St$ . This occurs at approximately $a_0t|\boldsymbol{\kappa _0}| \approx 3$ . At this stage, the dilatational component of the flow is small and thus the turbulence is predominately governed by the solenoidal dynamics.

Another important feature of compressible pure shear is the ‘crossover point’ displayed in the evolution of the shear anisotropy component $b_{12}$ . A crossover point at $St\approx 4$ , where trends with $M_d$ reverse, is predicted by the CPRM in figure 8. Before the crossover, $b_{12}$ increases with $M_d$ but decreases soon afterwards. The rate of decrease depends on the distortion Mach number $M_d$ . The crossover occurs due to coupling between the solenoidal and dilatational modes (Simone et al. Reference Simone, Coleman and Cambon1997).

Figure 8. Comparison of shear anisotropy evolution in pure shear between the CPRM (a) and RDT (b) for varying distortion Mach numbers prescribed in table 2. The top solid line is the pressure-released limit while the lower solid line is the solenoidal limit.

Predictions in figure 8 do not match as well with the RDT data of Simone et al. (Reference Simone, Coleman and Cambon1997) at large $St$ . It is noted in Simone et al. (Reference Simone, Coleman and Cambon1997) that the ‘waviness’ at large times in the RDT histories is the result of the wavenumber discretisation and resolution at large wavevectors. Since the integration with respect to wavevector magnitude is done numerically, this also impacts the long-time solution of the CPRM and is likely the cause for the mismatch between the CPRM and RDT of Simone et al. (Reference Simone, Coleman and Cambon1997). The histories of the normal $b_{ij}$ components are shown in figure 9. The RDT/DNS data for these quantities are not available, so we simply compare with the pressure-released and solenoidal limits. For most of the range, the predictions stay bounded. Recall that analytical solutions for pressure-released and solenoidal anisotropy limits do not exist, so they are approximated by using large/small $M_d$ values with pressure terms ignored or using the solenoidal model. This approximation is a reason why it does not always bound the lower/higher $M_d$ cases. The oscillation due to wavenumber discretisation is also pronounced and causes the long-time solution to go past the limits.

The equipartition between the pressure fluctuations and dilatational component of the turbulent kinetic energy is also recorded at large mean shear time scales in figure 10. The equipartition variable is defined as

(5.16) \begin{align} \phi _p=\frac {\overline {u'_2u'_2}\bar {p}\gamma \bar {\rho }}{\overline {p'p'}}=\iiint ^{\infty }_{-\infty }\frac {\Phi _{22}}{a^2\Gamma }{\rm d}^{3}{\boldsymbol{\kappa }}. \end{align}

Bertsch et al. (Reference Bertsch, Suman and Girimaji2012) used this definition of equipartition between the pressure and dilatational velocity components because the flow-normal kinetic energy component is nearly equal to the dilatational kinetic energy of the flow field at late times due to the alignment of the wavenumber vector with the flow-normal direction in homogeneous shear. The evolution of the equipartition as predicted by the CPRM shows a similar sharp decay and oscillatory behaviour near $\phi _p=1$ as the RDT analysis done by Bertsch et al. (Reference Bertsch, Suman and Girimaji2012). Since the pressure dilatation is a driving mechanism for distribution of energy between the mechanical and internal modes, it is expected that the equipartition follows a similar trend to that of the three stages of the pressure correlation/turbulent kinetic energy evolution as described previously. Indeed, it is seen that the start of the plateau of turbulent kinetic energy coincides with the equipartition of energy between dilatation and pressure modes ( $St\approx 6$ ). The evolution of the pressure statistics $\overline {p'p'}/\bar {p}^2\gamma ^2$ is also verified by comparing with the RDT data of Lavin et al. (Reference Lavin, Girimaji, Suman and Yu2012) in figure 11.

Figure 9. Evolution of normal anisotropy components in pure shear for varying distortion Mach numbers prescribed in table 2: pressure-released/solenoidal limits, top/bottom solid lines; $b_{11}$ , red; $b_{22}$ , blue; $b_{33}$ , black.

Figure 10. Evolution of the equipartition variable $\phi _p$ for varying higher distortion Mach numbers prescribed in table 2. Relaxation to equipartition is obtained when the turbulent kinetic energy curve reaches the second stage.

Figure 11. Evolution of the normalised pressure statistics $\overline {p'p'}/\bar {p}^2\gamma ^2$ for cases A–J in table 2. Data match the trends seen in Lavin et al. (Reference Lavin, Girimaji, Suman and Yu2012).

Table 2. Homogeneous isotropic turbulence initial conditions for pure-shear cases.

We compare our model results with those of Yu & Girimaji (Reference Yu and Girimaji2007). For this comparison, the histories of $b_{12}$ are compared for the same range of $M_d$ given in table 2. From figure 12, it is clear that the model of Yu & Girimaji (Reference Yu and Girimaji2007) does predict the correct crossover point, but the small details like the decay afterwards or the magnitude of the crossover point peaks do not match the RDT data as well as the CPRM does. This is likely due to the fact that Yu & Girimaji (Reference Yu and Girimaji2007) did not account for the appearance of the wavevector magnitude in the evolution equations which requires one to evolve the velocity spectrum tensor rather than the directional velocity spectrum tensor. The appearance of the $\kappa ^{3}$ factor in the pressure terms after conversion to spherical coordinates renders it impossible to analytically factor out the directional velocity spectrum tensor. The new model also addresses a discrepancy in the past compressible RDT model where the definitions for second moments in Fourier space considered the entire complex expression, which does not align with the symmetric velocity spectrum tensor definition by Durbin & Pettersson-Reif (Reference Durbin and Pettersson-Reif2011) used to compute the Reynolds stress. This resulted in complex evolution equations where real and imaginary components needed to be treated.

While the CPRM does not match the RDT data exactly, its predictions are much improved compared with past models. Again, due to the formulation of the CPRM and the requirement of the $\beta _i$ and $\Gamma$ spectra, the initialisation is not exactly the same as that of the RDT of Simone et al. (Reference Simone, Coleman and Cambon1997) since we do not know the pressure–velocity and pressure energy spectra. This is one reason why the predictions do not match the RDT results exactly, along with possible discretisation errors. Note that the CPRM and the RDT of Simone et al. (Reference Simone, Coleman and Cambon1997) are still very similar methods, sharing the same differences and analogies as the comparisons made between the CPRM and the RDT of Jacquin et al. (Reference Jacquin, Cambon and Blin1993) discussed in the previous section.

The distribution of the directional velocity spectrum over the unit- ${\boldsymbol{\kappa }}$ sphere is visualised in figure 13. These distributions evolve from a homogeneous isotropic (uniform) distribution over the unit sphere. The magnitudes of $\Phi ^{\nearrow \kappa *}_{kk}$ and $\Phi ^{\nearrow \kappa *}_{12}$ tend to decrease with increasing $M_d$ . This matches the behaviour seen in figure 8. Another interesting feature is that the distribution patterns of $\Phi ^{\nearrow \kappa *}_{kk}$ and $\Phi ^{\nearrow \kappa *}_{12}$ are identical for the same case and only vary in magnitude range.

Figure 12. Evolution of shear anisotropy components in pure shear for varying distortion Mach numbers with the new CPRM (black) compared with the model of Yu & Girimaji (Reference Yu and Girimaji2007) (blue).

Figure 13. Visualisation of the directional velocity spectrum trace and shear component in pure shear for varying distortion Mach numbers prescribed in cases (a) A, (b) D, (c) G and (d) J. Results are shown at time $St=10$ .

5.5. Sheared compression

Sheared compression is a combination of the last two cases. It is defined by

(5.17) \begin{align} G_{ij}=C\delta _{i1}\delta _{j1}+S\delta _{i1}\delta _{j2}, \quad C(t)=\frac {C_0}{1+C_0t},\quad S(t)=\frac {S_0}{1+C_0t}. \end{align}

The shear and compression rates change with time to satisfy the compressible Euler equations. To match the results of Mahesh (Reference Mahesh1996), anisotropic initial spectra for the three second moments are required. This is done by applying pure shear for $t_{ps}$ time, defined by the non-dimensional time $\mathcal{S}_0=S_0t_{ps}$ . Note that for this problem, only solenoidal turbulence is considered. Thus the solenoidal model is used to evolve the directional velocity spectrum tensor $\Phi ^{\nearrow \kappa *}_{ij}$ . The initial shear is intended to generate the anisotropic field but is negligible compared with the applied compression, i.e. $S_0/C_0\ll 1$ . Thus a low shear-compression ratio of $S_0/C_0=0.1$ is used (Mahesh Reference Mahesh1996). Only $\mathcal{S}_0$ varies across the different cases, so all other initial variables remain the same. For this case, $M_{t0}=0.025$ , $Re_{t0}=358$ and $M_{d0}=0.5$ . Again, 800 clusters are used, and wavevector magnitude integration is not required since the directional spectrum is evolved. The histories of $b_{12}=R_{12}/q^2$ are plotted in figure 14.

Figure 14. Evolution of solenoidal component of $R_{12}/q^2$ for varying anisotropic initial condition in sheared compression: blue, $\mathcal{S}_0=1$ ; black, $\mathcal{S}_0=2$ ; red, $\mathcal{S}_0=3$ ; circle, RDT of Mahesh (Reference Mahesh1996).

Each curve has different initial total shear, which increases from red to blue. Agreement with the RDT of Mahesh (Reference Mahesh1996) shows that the CPRM solenoidal method works well for anisotropic initial conditions. Slight discrepancies are likely due to imperfect initialisation, noting that the worst results also have the worst initialisation (start off with the most error). The initial anisotropic velocity spectrum was not given by Mahesh (Reference Mahesh1996) and thus the CPRM was used to obtain the initial spectrum by running the pure shear case for $\mathcal{S}_0$ total shear. Therefore, we are unable to match the initial spectrum exactly to that of Mahesh (Reference Mahesh1996). Figure 14 displays important trends of decreasing shear anisotropy and change of sign for the case with sufficiently large total volumetric strain (which correlates with larger initial total shear).

Figure 15. Evolution of solenoidal component of turbulent kinetic energy for varying anisotropic initial condition in sheared compression: blue, $\mathcal{S}_0=1$ ; black, $\mathcal{S}_0=2$ ; red, $\mathcal{S}_0=3$ ; circles, RDT of Mahesh (Reference Mahesh1996).

Figure 16. Evolution of normal components of $R_{ij}/q^2$ in sheared compression for anisotropic initial condition $\mathcal{S}_0=3$ : $R_{11}/q^2$ , blue; $R_{22}/q^2$ , black; $R_{33}/q^2$ , red; RDT of Mahesh (Reference Mahesh1996), circles.

Figure 15 shows the turbulent kinetic energy amplification for the same cases. Again, due to inexact initialisation, the turbulent kinetic energy amplification is also not exactly predicted but still has the same increasing trends for the three different total shear values. Figure 16 shows the normal components of $R_{ij}/q^2$ compared against the RDT of Mahesh (Reference Mahesh1996) and shows nearly exact agreement, where the compression increases the contribution of $u_1'$ to the turbulent kinetic energy and decreases that of $u'_2$ and $u_{3}^{\prime}$ . Only results for the maximum initial total shear $\mathcal{S}_0=3$ are given in figure 16.

5.6. Plane strain and axisymmetric contraction

The CPRM is applicable to any generic mean deformations that are rapid enough such that RDT still applies. To demonstrate this, the CPRM is applied to plane strain and axisymmetric deformations, defined by

(5.18) \begin{align} \frac {\partial \bar {u}_i}{\partial x_{\!j}}_{ps} = \left [ \begin{array}{ccc} S & \ \ 0 \ \ & 0 \\ 0 & \ \ -S \ \ & 0 \\ 0 & \ \ 0 \ \ & 0 \\ \end{array} \right ],\quad \frac {\partial \bar {u}_i}{\partial x_{\!j}}_{axc} = \left [ \begin{array}{ccc} S & \ \ 0 \ \ & 0 \\ 0 & \ \ -S/2 \ \ & 0 \\ 0 & \ \ 0 & \ \ -S/2 \\ \end{array} \right ]. \end{align}

For this problem, DNS/RDT datasets are not available and thus we cannot determine the model’s accuracy. Instead, we simply verify that the model is performing adequately by comparing predictions with the approximate solenoidal and pressure-released limits. For all cases, 800 clusters were used along with a 20-weight Gauss–Laguerre quadrature. Figure 17(a) shows the normal $b_{ij}$ components for plane strain, while figure 17(b) shows $b_{ij}$ for axisymmetric contraction. The results fall within the limits, showing that the the qualitative behaviour of the CPRM is correct. It is important to note that while cases A and J are near the limits, they are not so small/large that we expect them to follow the limiting lines exactly. The results show that the model can be applied to generic rapid deformations, although more high-fidelity data are needed to evaluate the accuracy of the CPRM.

Figure 17. Evolution of normal components of $b_{ij}$ subject to (a) plane strain and (b) axisymmetric contraction for initial conditions A, D, G and J listed in table 2. Same legend as in figure 9.