2.1 Solution algorithm in a deforming, anisotropic domain

Rogallo (Reference Rogallo1981) introduced the coordinate transformation $\unicode[STIX]{x1D709}_{i}=\unicode[STIX]{x1D609}_{ij}(t)x_{j}$ , in which the deforming coordinate system $\unicode[STIX]{x1D709}_{i}$ is related to the laboratory coordinates $x_{i}$ through the metric tensor $\unicode[STIX]{x1D609}_{ij}(t)$ . Summation over repeated Latin indices is implied throughout the text. The deforming coordinates are required to move with the mean flow, according to

where $\langle U_{i}\rangle$ is the mean velocity. In this coordinate system the turbulence is homogeneous, and periodic boundary conditions are applicable provided the mean velocity gradients are uniform in space. The continuity and momentum equations become

where $u_{i}$ is the fluctuating velocity, $p$ is the fluctuating pressure, $\unicode[STIX]{x1D70C}$ is the density and $\unicode[STIX]{x1D708}$ is the kinematic viscosity. Nonlinear terms are evaluated using a Fourier pseudo-spectral method, and the numerical solution in Fourier space is advanced in time using a second-order predictor–corrector scheme. Aliasing errors are mitigated by using truncation and grid shifting in wavenumber space (Rogallo Reference Rogallo1981).

For irrotational axisymmetric contraction, the mean deformation tensor is given by (Lee & Reynolds Reference Lee and Reynolds1985)

with $S(t)\geqslant 0$ . The solution domain is orthogonal at all times, such that only the diagonal elements of $\unicode[STIX]{x1D609}_{ij}(t)$ are non-zero. Solving (2.1) leads to (for $\unicode[STIX]{x1D6FC}=1,2,3$ , no summation)

where $\unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}^{0}\equiv \unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}(0)$ and $\exp [f_{\unicode[STIX]{x1D6FC}}(t)]$ is the total strain in each direction with

Figure 1 shows a schematic of how the domain is deformed during the application of strain. At any time $t$ , the length of the domain on each side given by $L_{\unicode[STIX]{x1D6FC}}=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}(t)$ varies as $L_{\unicode[STIX]{x1D6FC}}(t)=L_{\unicode[STIX]{x1D6FC}}^{0}\exp [f_{\unicode[STIX]{x1D6FC}}(t)]$ where $L_{\unicode[STIX]{x1D6FC}}^{0}$ is the initial length. While the number of grid points is fixed, the grid spacings and hence also resolution in each direction vary. Correspondingly, the wavenumbers represented in each direction of the domain are given by $k_{\unicode[STIX]{x1D6FC}}(n_{\unicode[STIX]{x1D6FC}},t)=n_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}(t)$ , where $n_{\unicode[STIX]{x1D6FC}}=-N_{\unicode[STIX]{x1D6FC}}/2+1,\cdots \,,N_{\unicode[STIX]{x1D6FC}}/2$ , and $N_{\unicode[STIX]{x1D6FC}}$ is the number of grid points in the $x_{\unicode[STIX]{x1D6FC}}$ direction. The maximum resolved wavenumber in each direction after truncation is $k_{max,\unicode[STIX]{x1D6FC}}(t)=\sqrt{2}N_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}(t)/3$ . As the simulation proceeds and the domain becomes deformed from its original shape, it is necessary to check that the large scales continue to be well contained within the domain in all directions, while the small scales likewise remain well resolved (Gualtieri & Meneveau Reference Gualtieri and Meneveau2010).

Figure 1. Grid deformation under axisymmetric contraction with a total elongation of 4 in the $x_{1}$ direction. The grid metrics in the $x_{2}$ and $x_{3}$ directions are equal.

2.2 A time-dependent strain as a model for experiment

In most wind tunnel experiments the change in cross-section is gradual, which results in a smooth variation of the mean velocity with distance downstream along the wind tunnel centreline. The mean strain rate along the wind tunnel centreline is thus readily established as a function of the spatial coordinates. However, as a fluid element passes through the wind tunnel, the spatially varying mean strain rate is experienced in a time-dependent manner. To model the mean strain rates from experimental wind tunnels in our time-dependent DNS, we apply a spatially uniform, but time-dependent mean strain rate that corresponds to the local value of the mean strain rate that a fluid element experiences as it passes through the wind tunnel (Pearson Reference Pearson1959). This is accomplished by introducing a convective time $t$ as the time taken for a fluid element travelling with the mean flow to reach a position $x$ from a reference location $x_{a}$ . (Here $x$ alone, or $x_{a}$ etc. without tensor subscripts shall refer to distances measured along the wind tunnel centreline.) We write

According to (2.5)–(2.6), by this time $t$ the grid metric in the extensional direction will be

The integral in (2.8) can be evaluated through a change of variables, $\text{d}\unicode[STIX]{x1D70F}=\text{d}x/\langle U_{1}(x)\rangle$ (suggested by (2.7)). After some simplifications, we obtain

Since the applied mean strain is volume preserving, the product of the three diagonal metric factors $\unicode[STIX]{x1D609}_{11}\unicode[STIX]{x1D609}_{22}\unicode[STIX]{x1D609}_{33}$ is fixed, while axisymmetry requires $\unicode[STIX]{x1D609}_{22}=\unicode[STIX]{x1D609}_{33}$ . As a result, the change in the transverse grid metrics is given by

Since in the simulation the turbulence is advanced in time instead of space, at each time step from $t_{n}$ to $t_{n}+\unicode[STIX]{x0394}t$ , it is necessary to solve (2.7) for a new spatial location $x$ , so that the new metric factors can be evaluated according to (2.9)–(2.10). In addition, at every time step the metric factors are used to form the viscous integrating factors (Rogallo Reference Rogallo1981). The solution to (2.7) and the calculation of the integrating factors are carried out using QUADPACK (Piessens et al. Reference Piessens, de Doncker-Kapenga, Überhuber and Kahaner1983).

In the experiments of AW, the mean streamwise velocity $\langle U_{1}(x)\rangle$ along the centreline of the wind tunnel is well described by an error function, which yields a Gaussian profile for the mean extensional strain (see figure 2 of AW). We consider the functional form

where the five parameters $a$ , $b$ , $c$ , $x_{a}$ and $x_{b}$ are chosen so the mean strain closely matches the experiments in a non-dimensional sense, and the origin of the $x$ -axis has been placed at the location of maximum mean velocity gradient. In the experiments, some appreciable strain also existed immediately upstream and downstream of the physical beginning and ending locations of the contraction. To incorporate this feature in our DNS, we specify two intermediate locations $x_{i}$ and $x_{f}$ satisfying $x_{a}<x_{i}<x_{f}<x_{b}$ ( $x_{i}$ and $x_{f}$ correspond to the vertical lines in figure 2 of AW). Since these two locations are approximately equally far from the location of maximum strain rate, we take $x_{i}=-x_{f}$ . The straining ‘period’ in the DNS is considered to be the entire phase of the simulation during which there is non-zero mean strain, corresponding to the physical distance between $x_{a}$ and $x_{b}$ . Figure 2(a) presents a sketch of an error function mean velocity profile with these locations marked for later reference. While $x_{i}$ and $x_{f}$ are used to specify the mean velocity variation, we report pre- and post-contraction statistics at $x_{a}$ and $x_{b}$ , respectively.

Figure 2. (a) Schematic of mean velocity profile in experiments, and (b) examples of non-dimensional strain rate $S\unicode[STIX]{x1D70F}_{a}$ as function of convective time in experiments and DNS. In (b): ○ (black) is for the AW $R_{\unicode[STIX]{x1D706}}=260$ experiment, ▫ (blue) for AW $R_{\unicode[STIX]{x1D706}}=40$ experiment, ▴ (magenta) for $S_{0}^{\ast }=25$ numerical, ♦ (green) for $S_{0}^{\ast }=20$ numerical and ▾ (red) for $S_{0}^{\ast }=15$ numerical.

For a systematic procedure for choosing the five curve-fit parameters noted above, we need five constraints to control both the spatial spread of the profile and its maximum velocity gradient (which is proportional to the product $ab$ ). First, since the value of $c$ does not affect the strain rate, it is arbitrary provided it is large enough to ensure that $\langle U_{1}(x)\rangle >0$ for all $x$ . Second, we specify the velocity ratio $\langle U_{1}(x_{f})\rangle /\langle U_{1}(x_{i})\rangle$ at a value close to the experiment. Third, we require that the strain rate $S$ at $x=0$ , denoted by $S_{0}$ , corresponds to a non-dimensional strain $S_{0}^{\ast }$ formed from $S_{0}$ and a large-eddy time scale of the turbulence before strain is applied. Specifically, we define $S_{0}^{\ast }=S_{0}\unicode[STIX]{x1D70F}_{a}$ , where $\unicode[STIX]{x1D70F}_{a}=(2K/\langle \unicode[STIX]{x1D716}\rangle )_{a}$ , with $K$ and $\unicode[STIX]{x1D716}$ being the turbulence kinetic energy and the dissipation rate. Fourth, the starting location $x_{a}$ is chosen to give a desired non-dimensional strain rate at $x_{a}$ based on the experiments. Finally, the ending location $x_{b}$ is chosen to achieve a desired velocity ratio $\langle U_{1}(x_{b})\rangle /\langle U_{1}(x_{a})\rangle$ , i.e. deformation, for the entire straining period.

Figure 2(b) shows a qualitative comparison of the non-dimensional strain rate profiles obtained numerically from three values of $S_{0}^{\ast }$ with experimental data derived from AW at two values of the pre-strain Reynolds number. Experimental data in this figure are obtained first by applying curve fits of a form similar to (2.11) and then differentiating. The non-dimensional strain rates are plotted at convective times $t$ in the laboratory facility obtained using (2.7) and normalized by large-eddy turnover times obtained from table 1 of AW. The numerical strain rates are obtained using the aforementioned procedure with $c=10$ , $bx_{f}=0.9$ , a velocity ratio $\langle U_{1}(x_{f})\rangle /\langle U_{1}(x_{i})\rangle =2.85$ , a starting non-dimensional strain rate of $0.3$ at $x_{a}$ and an overall velocity ratio $\langle U_{1}(x_{b})\rangle /\langle U_{1}(x_{a})\rangle =4$ . The resulting strain rate profiles are similar to those in the experiments, suggesting that the procedures described here model the mean velocity variation in the wind tunnel contraction well.

2.3 Spectral evolution and rapid-distortion theory

To understand the dynamical processes underlying the change in shape of the energy spectrum under applied strain, we need to compute various terms in the spectral evolution equations. In isotropic turbulence, computation of the nonlinear transfer spectrum of the kinetic energy (e.g. Domaradzki & Rogallo Reference Domaradzki and Rogallo1990; Yeung, Brasseur & Wang Reference Yeung, Brasseur and Wang1995) in spherical wavenumber shells is sufficient. However, in this work we also have to consider pressure-strain correlations and systematic anisotropy, which requires individual components of spectra in one- and two-dimensional partitions of wavenumber space. In addition, the mean strain distorts the wavenumbers on a deforming domain, leading to transport of the spectrum in wavenumber space (Pope Reference Pope2000).

In Fourier space, the fluctuating velocity on a deforming periodic domain evolves by

where the metric factors in (2.3) are expressed through time dependence of the wave vector $\boldsymbol{k}$ , $\hat{p}(\boldsymbol{k})$ is the Fourier coefficient of the fluctuating pressure, $G_{i}(\boldsymbol{k})=\text{i}k_{j}\widehat{u_{i}u_{j}}$ and $\text{i}=\sqrt{-1}$ . It is important to distinguish between rapid pressure and slow pressure, which respond to the mean flow and the velocity fluctuations respectively; we write $\hat{p}=\hat{p}^{r}+\hat{p}^{s}$ . Both are obtained by solving well-known Poisson equations. The spectral covariance $\unicode[STIX]{x1D60C}_{ij}(\boldsymbol{k})\equiv \langle \hat{u} _{i}^{\ast }(\boldsymbol{k})\hat{u} _{j}(\boldsymbol{k})\rangle$ (where superscript $\ast$ denotes complex conjugates) is in general complex valued, although because mean shear is absent we are only concerned with the diagonal components, which are real. We can write

where terms on the right-hand side represent, respectively, production due to the mean velocity gradient, redistribution due to rapid pressure and slow pressure, spectral transfer due to nonlinear terms and dissipation due to viscosity. Specifically, these terms are:

The angled brackets in these equations represent averaging over multiple realizations, i.e. ensemble averaging, which we perform for the majority of our simulations detailed later. When axisymmetric contraction is applied the production term is negative for $\unicode[STIX]{x1D60C}_{11}(\boldsymbol{k})$ , but positive for $\unicode[STIX]{x1D60C}_{22}(\boldsymbol{k})$ and $\unicode[STIX]{x1D60C}_{33}(\boldsymbol{k})$ , thus causing anisotropy directly, especially at the large scales. The pressure-strain term exchanges energy among the diagonal components $\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}(\boldsymbol{k})$ and is traceless due to incompressibility, for both rapid and slow pressures. The rapid pressure is present only while strain is applied, while slow pressure is also important during relaxation. Since nonlinear transfer redistributes energy in Fourier space, each component of $T_{ij}(\boldsymbol{k})$ integrates to zero over wavenumber space. Integrating (2.13) over all wavenumbers gives the Reynolds stress transport equation for homogeneous turbulence, which is

where $s_{ij}$ is the fluctuating strain rate.

In (2.13), only the production and rapid pressure-strain terms depend directly on the mean strain rate. As a result, if the strain rate is very strong, i.e. very rapid compared to the time scales of the turbulence, this equation can be simplified by retaining only $P_{ij}(\boldsymbol{k})$ and $\unicode[STIX]{x1D6F1}_{ij}^{r}(\boldsymbol{k})$ , while neglecting viscous dissipation and the nonlinear effects of slow pressure and spectral transfer. The behaviour of the energy spectrum tensor can then be described by inviscid RDT (Townsend Reference Townsend1976; Savill Reference Savill1987). Although in practice the strain rates in experiments and simulations are necessarily finite, comparisons with RDT theory are still useful for a better understanding.

The use of a deforming domain, and hence a time-dependent set of wave vectors in our simulations, requires some care when interpreting spectra. It would be useful to relate the 1-D spectrum $\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}(k_{\unicode[STIX]{x1D6FD}})$ to its equivalent representation $\unicode[STIX]{x1D60C}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}^{0}(k_{\unicode[STIX]{x1D6FD}}^{0})$ as a function of the pre-strain wavenumbers $k_{\unicode[STIX]{x1D6FD}}^{0}$ . This would clearly show how a set of modes (those perpendicular to an initial $k_{1}$ ) are affected by the strain. The integrals of these spectra over the corresponding wavenumbers $k_{\unicode[STIX]{x1D6FD}}$ and $k_{\unicode[STIX]{x1D6FD}}^{0}$ are both equal to $\langle u_{\unicode[STIX]{x1D6FC}}^{2}\rangle$ . At any time $t$ during the straining period, the current and pre-strain wavenumbers are related by

where $\unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}^{0}/\unicode[STIX]{x1D609}_{\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FD}}(t)$ is the total deformation in the $x_{\unicode[STIX]{x1D6FD}}$ direction. A change of variables between the integrals in $k_{\unicode[STIX]{x1D6FD}}$ and $k_{\unicode[STIX]{x1D6FD}}^{0}$ then gives the relation

Following AW, we present most of our 1-D spectral results in terms of the pre-strain wavenumbers during the application of strain, and in terms of the post-strain wavenumbers (which are fixed) during the relaxation period.

Figure 3. Ring decomposition of wavenumber space for axisymmetric turbulence. The axis of axisymmetry $\unicode[STIX]{x1D740}$ is in the $\boldsymbol{k}_{1}$ direction, and the ring is perpendicular to $\unicode[STIX]{x1D740}$ . The longitude with respect to $\boldsymbol{k}_{2}$ is $\unicode[STIX]{x1D703}$ and the colatitude with respect to $\unicode[STIX]{x1D740}$ is $\unicode[STIX]{x1D719}$ .

Statistics in axisymmetric turbulence are rotationally symmetric about a preferred direction $\unicode[STIX]{x1D740}$ (Batchelor Reference Batchelor1946), which in this study is the extensional direction. Figure 3 shows a decomposition of Fourier space motivated by this rotational symmetry, such that spectral quantities are expressed as functions of $k_{1}$ and the wavenumber magnitude in the cross-sectional plane $k_{r}=\sqrt{k_{2}^{2}+k_{3}^{2}}$ . We can also write $k_{r}=k\sin \unicode[STIX]{x1D719}$ , where $k=\sqrt{\boldsymbol{k}\cdot \boldsymbol{k}}$ and $0\leqslant \unicode[STIX]{x1D719}\leqslant \unicode[STIX]{x03C0}$ is the colatitude with respect to $\unicode[STIX]{x1D740}$ . The spectral covariance $\unicode[STIX]{x1D60C}_{ij}(\boldsymbol{k})$ defined earlier is a discrete version of the velocity spectrum tensor $\unicode[STIX]{x1D6F7}_{ij}(\boldsymbol{k})$ (whose integral in wavenumber space gives $\langle u_{i}u_{j}\rangle$ ). We define the axisymmetric spectrum tensor as

from which the 1-D spectra can be recovered by

where the factor of 2 is used to collect contributions from both positive and negative values of $k_{1}$ . Axisymmetric representations of the spectra are useful in studies of rotating (Clark di Leoni et al. Reference Clark di Leoni, Cobelli, Mininni and Dmitruk2014) and stably stratified turbulent flows (Godeferd & Staquet Reference Godeferd and Staquet2003). Evolution equations for the axisymmetric spectra and 1-D spectra are readily obtained by integration of (2.13) over rings and planes in wavenumber space, respectively. For isotropic turbulence, the contours of the axisymmetric energy spectrum $E_{A}(k_{1},k_{r})\equiv {\textstyle \frac{1}{2}}A_{ii}(k_{1},k_{r})$ multiplied by $1/\sin \unicode[STIX]{x1D719}$ are circles in the $(k_{1},k_{r})$ plane (Mininni et al. Reference Mininni, Rosenberg and Pouquet2012). In practice, the axisymmetric spectra are formed by summing over Fourier modes residing in rings of finite thickness. Because the computational space is Cartesian, the distribution of modes in rings at small $k_{r}$ is significantly uneven, which can lead to some noise in contours of the axisymmetric spectra. A node-density correction factor is applied in a manner similar to those used for 3-D spectra collected into discrete spherical shells in wavenumber space (Eswaran & Pope Reference Eswaran and Pope1988).

4.1 Single-point moments

We first provide data in table 3 that summarizes the post-contraction state of each run, to be followed by figures that show the evolution of important statistics over the complete straining period. The post-contraction shape of the domains is as shown in the right of figure 1, with a decrease in $\unicode[STIX]{x1D609}_{11}$ and increase in $\unicode[STIX]{x1D609}_{22}$ and $\unicode[STIX]{x1D609}_{33}$ relative to their pre-contraction values in table 1. As the domain deforms and the physical length scales of the turbulence evolve, large-scale sampling and small-scale resolution (second block in table 3) change. The post-contraction domains are still large compared with the integral length scales, indicating that large-scale sampling is still adequate. It is worth noting that the ratio $L_{1}/\ell _{21}$ has dropped compared to its pre-contraction value, despite a factor of 4 increase in $L_{1}$ . This implies a substantial increase in the transverse integral length scale $\ell _{21}$ , which is consistent with the formation of coherent longitudinal vortices in the extensional direction (Rogers & Moin Reference Rogers and Moin1987). As expected, small-scale resolution worsens in the direction where the domain is stretched ( $x_{1}$ ), but improves in the directions where the domain is compressed ( $x_{2}$ and $x_{3}$ ).

Continuing in table 3, the turbulence kinetic energy is amplified in all cases, with the amplification ratio nearly independent of the Reynolds number. The dissipation rate also increases, but less so at higher Reynolds number. It is clear that the turbulence becomes anisotropic, as velocity fluctuations in the extensional and compressive directions are suppressed and amplified, respectively. Anisotropy at the small scales is also seen in the statistics of velocity gradient fluctuations, such as the ratio of transverse to longitudinal velocity gradient variances, which differ from the isotropic value of 2. Apparently, in the higher Reynolds number runs, the ratio $\langle u_{2,1}^{2}\rangle /\langle u_{1,1}^{2}\rangle$ becomes less anisotropic but $\langle u_{3,2}^{2}\rangle /\langle u_{2,2}^{2}\rangle$ is nearly fixed at close to 3, which (as we discuss further later in this subsection) is an indication of quasi two-dimensionality in the cross-sectional plane.

Third and fourth moments of longitudinal velocity gradients are also included in table 3. In 3-D incompressible isotropic turbulence, the skewness of the longitudinal velocity gradient is negative, and related to the phenomenon of vortex stretching (Batchelor Reference Batchelor1953; Tavoularis, Bennett & Corrsin Reference Tavoularis, Bennett and Corrsin1978). However, as reported by AW and others (Mills & Corrsin Reference Mills and Corrsin1959; Sjögren & Johansson Reference Sjögren and Johansson1998), we find that axisymmetric contraction causes the skewness of $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ to undergo a change in sign, with its magnitude increasing with the Reynolds number. In contrast, the skewness of $\unicode[STIX]{x2202}u_{3}/\unicode[STIX]{x2202}x_{3}$ remains negative but its magnitude is reduced compared to that observed in isotropic turbulence. At the same time, we observe an increase in the flatness of $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ and a decrease in the flatness of $\unicode[STIX]{x2202}u_{3}/\unicode[STIX]{x2202}x_{3}$ during the contraction. The flatness of $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ also increases with Reynolds number, which is consistent with AW. A comparison of the post-contraction flatness of $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ between runs with different resolution (Run 2 and 8 versus 1 and 7, respectively) suggests that higher-order derivative statistics in this work are affected by finite resolution in the extensional direction. However, we are mainly interested in lower-order quantities such as spectra, for which the lower-resolution simulations are adequate. We show in the Appendix that the major results from this paper are not affected.

It is also useful to compare DNS, which uses a finite strain rate, to RDT. Following the pre-simulation, we evolve the Fourier coefficients of the velocity field according to well-known relations for inviscid RDT (Townsend Reference Townsend1976; Lee & Reynolds Reference Lee and Reynolds1985). Each velocity field is subjected to a 4 : 1 area ratio axisymmetric contraction. We then ensemble-average statistics over all realizations for each run, e.g. over the 16 simulations comprising Run 1. RDT predicts $K_{b}/K_{a}=2.1$ and $\langle \unicode[STIX]{x1D716}\rangle _{b}/\langle \unicode[STIX]{x1D716}\rangle _{a}=5.5$ for all runs, which are considerably different from the DNS results presented in table 3. The extent of large-scale anisotropy is predicted well by RDT; it gives $\unicode[STIX]{x1D623}_{11}^{b}=-0.3$ and $\unicode[STIX]{x1D623}_{22}^{b}=0.15$ for all runs. There is more discrepancy between the DNS and RDT when examining velocity derivative statistics. For example, RDT estimates $\langle u_{2,1}^{2}\rangle _{b}/\langle u_{1,1}^{2}\rangle _{b}=8$ , which is markedly different from the DNS results. The velocity derivative statistic $\langle u_{3,2}^{2}\rangle _{b}/\langle u_{2,2}^{2}\rangle _{b}$ , however, takes a post-contraction value of 3 with RDT, which is very similar to the DNS results. Table 4 shows the RDT predictions for higher-order velocity derivative statistics. While RDT fails to predict a positive value for the skewness of $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ , a small negative value of the skewness of $\unicode[STIX]{x2202}u_{3}/\unicode[STIX]{x2202}x_{3}$ is consistent with the DNS at lower Reynolds numbers. The large increase in the flatness of $\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}$ at high Reynolds numbers, and the decrease in the flatness of $\unicode[STIX]{x2202}u_{3}/\unicode[STIX]{x2202}x_{3}$ observed in the DNS are not predicted by RDT.

Figure 4. Evolution of (a) $K$ and $\langle \unicode[STIX]{x1D716}\rangle$ normalized by initial values for Runs 4, 6 and 8, (b) budget for $\text{d}K/\text{d}t$ normalized by $\langle \unicode[STIX]{x1D716}\rangle _{a}$ for Runs 4 and 8, and (c) non-dimensional strain rates for Runs 4, 6 and 8. In (a), solid curves (red) for $K(t)/K_{a}$ for the three runs are almost coincident, and dashed curves (blue) are for $\langle \unicode[STIX]{x1D716}(t)\rangle /\langle \unicode[STIX]{x1D716}\rangle _{a}$ , with $R_{\unicode[STIX]{x1D706}}$ increasing in the direction of the arrow. In (b), dashed curves with open symbols are for Run 4, and solid curves with filled symbols for Run 8: ▪ and ▫ (blue) for production, ● and ○ (red) for minus the dissipation, and ▴ and ▵ (black) for overall rate of change. In (c), mean strain rate normalized by large-eddy turnover time $\unicode[STIX]{x1D70F}$ (upper red curves) and Kolmogorov time scale $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ (lower blue curves), with $R_{\unicode[STIX]{x1D706}}$ increasing in the directions of the arrows.

For information on the evolution of the turbulence during the period of axisymmetric contraction, and to facilitate comparison with previous studies, it is useful to show results against the total deformation at a given time, rather than time itself. Figure 4 shows, as a function of $L_{1}(t)/L_{1}^{0}$ (which ranges from 1 to 4 for a 4 : 1 contraction ratio), in (a) the evolution of turbulence kinetic energy and dissipation rate, in (b) various terms in the turbulence kinetic energy budget and in (c) the non-dimensional mean strain rates. To focus on the effects of the Reynolds number we have selected data from Runs 4, 6 and 8 (see table 1). For homogeneous turbulence the turbulence kinetic energy budget is governed by production and dissipation, as

It can be seen that both $K$ and $\langle \unicode[STIX]{x1D716}\rangle$ initially decay since it takes finite time for production (initially zero) to grow to exceed the dissipation. Subsequently, as both variables grow, the evolution of $K$ is almost independent of the Reynolds number, while dissipation grows less rapidly at high Reynolds number. Because of the form of the strain rate profile (figure 2), both $K$ and $\langle \unicode[STIX]{x1D716}\rangle$ decrease towards the end of the straining period when the strain rate is weak. The production term in frame (b) is driven by the mean flow and is essentially the same for all Reynolds numbers simulated. It is also dominant over dissipation during the straining period, which explains why the kinetic energy evolution is nearly identical for all runs. The non-dimensional strain rates in frame (c) indicate that the strain rates are strong compared to the time scales of the large-scale motions, but actually weak from the perspective of the small scales.

The production of turbulence kinetic energy by mean strain requires anisotropy in the Reynolds stresses, which is expressed by the anisotropy tensor $\unicode[STIX]{x1D623}_{ij}=\langle u_{i}u_{j}\rangle /(2K)-\unicode[STIX]{x1D6FF}_{ij}/3$ , and its second and third coordinate-frame invariants. We use the definitions

where here $\unicode[STIX]{x1D702}$ is not to be confused with the Kolmogorov scale. For isotropic turbulence subjected to axisymmetric contraction the invariants satisfy

Figure 5 shows for Runs 4 and 8, (a) the evolution of the Reynolds stresses, (b) the evolution of the anisotropy tensor elements $\unicode[STIX]{x1D623}_{11}$ and $\unicode[STIX]{x1D623}_{22}$ and (c) the invariants $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D702}$ . The results suggest that the large-scale anisotropy is completely determined by the total deformation at any time $t$ , but is independent of the Reynolds number. Since the large-scale anisotropy depends only on the total deformation, the anisotropy development is the same as observed in simulations with constant strain rates (Lee & Reynolds Reference Lee and Reynolds1985), and can be predicted almost exactly by RDT. Although (4.3) is in principle exact, in numerical results it is not guaranteed if the large scales contributing the most to the Reynolds stress tensor are not sampled well in a domain of finite size. Our present results show that the large scales are sufficiently well sampled.

Figure 5. Evolution of (a) Reynolds stresses normalized by $q_{a}^{2}=2K_{a}$ , (b) components of the Reynolds stress anisotropy tensor and (c) anisotropy tensor invariants for Run 4 (dashed curves with open symbols) and Run 8 (solid curves with filled symbols). Symbols ▪ and ▫ (red) are for $\langle u_{1}^{2}\rangle /q_{a}^{2}$ , $\unicode[STIX]{x1D623}_{11}$ and $\unicode[STIX]{x1D709}$ in each respective figure. Symbols ● and ○ (blue) are for $\langle u_{2}^{2}\rangle /q_{a}^{2}$ , $\unicode[STIX]{x1D623}_{22}$ and $\unicode[STIX]{x1D702}$ in each respective figure. Dashed lines in (c) are for the two-dimensional isotropic limit.

Figure 6. Reynolds stress budgets normalized by initial dissipation rate $\langle \unicode[STIX]{x1D716}\rangle _{a}$ during the application of strain for Run 4 (dashed curves with open symbols) and Run 8 (solid curves with filled symbols). Budget terms for $\text{d}\langle u_{1}^{2}\rangle /\text{d}t$ are shown in (a) and budget terms for $\text{d}\langle u_{2}^{2}\rangle /\text{d}t$ in (b): ▾ and ▿ (red) are for production, ▴ and ▵ (magenta) for rapid pressure strain, ♦ and ♢ (green) for slow pressure strain, ● and ○ (blue) for dissipation, and the sum of all budget terms is marked by ▪ and ▫ (black).

The development of anisotropy in the Reynolds stress tensor can be analysed further using the Reynolds stress transport equation (2.19). Figure 6 presents the terms in the balance equations for (a) $\langle u_{1}^{2}\rangle$ and (b) $\langle u_{2}^{2}\rangle$ , using data from Runs 4 and 8. For the production terms, in this flow $P_{11}<0$ whereas $P_{22}>0$ . The rapid pressure-strain term quickly counteracts the anisotropy generated by the production term. The slow pressure-strain term, on the other hand, becomes significant only at later times, and is more important at higher Reynolds number. In frame (b) the production effect is clearly the dominant term in the Reynolds stress budget for $\langle u_{2}^{2}\rangle$ , being resisted only weakly by the rapid pressure strain at early times and viscous dissipation at later times. In both frames of this figure the production terms (of either sign) reach peak amplitude at intermediate times, which is a consequence of the time-dependent strain rate and is different from results from simulations of constant strain (Lee & Reynolds Reference Lee and Reynolds1985).

Figure 7. Evolution of (a) mean-square vorticities normalized by dissipation rate and (b) velocity derivative statistics for Run 4 (dashed curves with open symbols) and Run 8 (solid curves with filled symbols). In (a), ▪ and ▫ (black) are for $\unicode[STIX]{x1D708}\langle \unicode[STIX]{x1D714}_{1}^{2}\rangle /\langle \unicode[STIX]{x1D716}\rangle _{a}$ , ● and ○ (red) for $\unicode[STIX]{x1D708}\langle \unicode[STIX]{x1D714}_{2}^{2}\rangle /\langle \unicode[STIX]{x1D716}\rangle _{a}$ and ▴ and ▵ (blue) for $\unicode[STIX]{x1D708}\langle \unicode[STIX]{x1D714}_{3}^{2}\rangle /\langle \unicode[STIX]{x1D716}\rangle _{a}$ . Let $u_{i,j}=\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}$ . In (b), ▪ and ▫ (red) are for $\langle u_{2,1}^{2}\rangle /\langle u_{1,1}^{2}\rangle$ , ● and ○ (blue) for $\langle u_{1,2}^{2}\rangle /\langle u_{2,2}^{2}\rangle$ , ♦ and ♢ (magenta) for $-\langle u_{1,1}^{2}\rangle /\langle u_{2,3}u_{3,2}\rangle$ , ▴ and ▵ (green) for $\langle u_{3,2}^{2}\rangle /\langle u_{2,2}^{2}\rangle$ , and ▾ and ▿ (black) for $-\langle u_{2,2}^{2}\rangle /\langle u_{2,3}u_{3,2}\rangle$ . Values for $\langle u_{3,2}^{2}\rangle /\langle u_{2,2}^{2}\rangle$ and $-\langle u_{2,2}^{2}\rangle /\langle u_{2,3}u_{3,2}\rangle$ in two-dimensional isotropic turbulence marked by horizontal dashed lines at 3 and 1, respectively.

Although large-scale statistics show little dependence on the Reynolds number, small-scale statistics, such as the dissipation rate in figure 4(a), show a strong Reynolds number dependence. It was already observed in table 3 that the small scales become anisotropic during the straining period. The geometry of axisymmetric contraction also suggests that vorticity components in different directions ( $\unicode[STIX]{x1D714}_{1}$ versus $\unicode[STIX]{x1D714}_{2}$ and $\unicode[STIX]{x1D714}_{3}$ ) will have different statistics. We therefore investigate the behaviour of both vorticity and velocity gradient fluctuations below, with data given in the two frames of figure 7, respectively.

In homogeneous turbulence we can write

Figure 7(a) shows, for Runs 4 and 8, the three mean-squared vorticities, normalized by the viscosity and the pre-contraction dissipation rate. As expected from the geometry, axisymmetric contraction amplifies and aligns vorticity in the extensional direction while reducing vorticity in the compressive directions (Rogers & Moin Reference Rogers and Moin1987), but the effects are weaker at higher Reynolds number. This dependence on the Reynolds number can be understood by noting that (lower curves in figure 4 c) as the Reynolds number increases (in the order Runs 4, 6, 8) the value of the non-dimensional strain rate $S\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ becomes smaller. In other words, as the range of time scales in the flow widens with increasing Reynolds number, the strain rate becomes weaker with respect to the small scales. Hence, at higher Reynolds numbers, the fluctuating vorticity is less responsive to the applied strain, as reflected in the weaker amplification of $\langle \unicode[STIX]{x1D714}_{1}^{2}\rangle$ in figure 7(a).

To further assess departures from local isotropy using statistics of the velocity gradients, we can compare the simulation data with a number of relations that apply for 3-D incompressible isotropic turbulence, such as (with $\unicode[STIX]{x1D6FC}\neq \unicode[STIX]{x1D6FD}$ )

Figure 7(b) shows the evolution of ratios formed from (4.5) and (4.6) during the straining period. The velocity derivative statistics are initially isotropic, but depart from isotropy as the strain is applied. The statistics are more anisotropic at the lower Reynolds number because the strain rate is more rapid with respect to the small scales. In their experiments, AW measured $\langle u_{2,1}^{2}\rangle /\langle u_{1,1}^{2}\rangle$ for a wide range of Reynolds numbers (see figure 8 in AW). The post-contraction values for this statistic in the DNS are similar to those reported by AW, and also remain closer to the isotropic value as the Reynolds number is increased.

It is worth noting that figure 7(b) shows that during straining, $\langle u_{3,2}^{2}\rangle /\langle u_{2,2}^{2}\rangle$ and $-\langle u_{2,2}^{2}\rangle /\langle u_{2,3}u_{3,2}\rangle$ approach asymptotic limits 3 and 1, respectively. This can be explained by generalizing (4.5) and (4.6) to $n$ -dimensional isotropic turbulence, where $n\geqslant 2$ . We have, for $\unicode[STIX]{x1D6FC}\neq \unicode[STIX]{x1D6FD}$ (Pope Reference Pope2000; Gotoh et al. Reference Gotoh, Watanabe, Shiga, Nakano and Suzuki2007)

Substituting $n=2$ into (4.7) and (4.8), we obtain the limiting values observed in figure 7(b) for $\langle u_{3,2}^{2}\rangle /\langle u_{2,2}^{2}\rangle$ and $-\langle u_{2,2}^{2}\rangle /\langle u_{2,3}u_{3,2}\rangle$ , respectively. This suggests that the small scales under an axisymmetric contraction of sufficient strength tend to asymptotically approach the limiting state of isotropic turbulence in two dimensions. At the same time, a number of relations for velocity gradient statistics conforming to a state of local axisymmetry (George & Hussein Reference George and Hussein1991) are well verified in our DNS data. In other words, the post-contraction state of the small scales is one of local axisymmetry with velocity gradients in the extensional direction becoming asymptotically small compared to gradients in the compressive directions. (This is supported by the line for the ratio $-\langle u_{1,1}^{2}\rangle /\langle u_{2,3}u_{3,2}\rangle$ in the figure approaching almost zero on the scales chosen.)

Figure 8. Axisymmetric energy spectrum for (a–c) Run 4 and (d–f) Run 8 (a,d) before the application of strain, (b,e) half way through straining ( $L_{1}/L_{1}^{0}=2$ ) and (c,f) at the end of the straining period. Contour levels decrease by a factor of 10. Spectra plotted against the instantaneously distorting wavenumbers, and multiplied by 2 to recover $K$ when integrating over $k_{r}$ and non-negative $k_{1}$ . Simulation cutoff wavenumbers marked by outermost black boundary. Spectra normalized by $\sin \unicode[STIX]{x1D719}=k_{r}/k$ to obtain circular contours in isotropic turbulence; data for $k_{r}=0$ omitted from plot.

4.2 Spectral evolution

Results discussed above indicate that anisotropy in the mean flow and the large-scale motions ultimately lead to strong anisotropy in the small scales. This observation implies that isotropy in the flow is scale dependent, which is best explored through spectral quantities in wavenumber space. We consider an axisymmetric representation of the energy spectrum (§ 2.3), 1-D spectra like those measured in the experiments of AW, and various terms in the spectral energy budget (based on (2.13)).

Figure 8 shows contour plots of the axisymmetric energy spectrum $E_{A}(k_{1},k_{r})$ for Runs 4 and 8 before the straining period (a,d), half way through the straining period (b,e), and at the end of the straining period (c,f). If the turbulence is isotropic the energy spectrum would depend only on $k=\sqrt{k_{1}^{2}+k_{r}^{2}}$ , which means isocontours of $E_{A}(k_{1},k_{r})$ (after the normalization discussed in § 2.3) would be circles in the ( $k_{1},k_{r}$ ) plane. This is indeed the case for the pre-contraction spectra in frames (a) and (d), except that the shape of the contours is distorted near the simulation cutoff wavenumbers, which are marked by the outermost black elliptical boundaries in the plots. As the strain is applied (from the left column to the right column), the spectra change significantly. The post-contraction energy spectra in frames (c) and (f) are highly anisotropic at all scales of motion (with non-circular contours), which is consistent with the large-scale and small-scale anisotropy observed for single-point moments in § 4.1. The contours in frame (c) for Run 4 show evidence of stronger anisotropy than those in frame (f) for Run 8. This is in agreement with the results for single-point moments, which showed that the small scales become more anisotropic during the application of strain at lower Reynolds number. During straining, the domain is lengthened in the $x_{1}$ direction but shortened in $x_{2}$ and $x_{3}$ . This results in wavenumber distortion, where the wavenumbers in $k_{1}$ decrease while those in $k_{2}$ and $k_{3}$ increase. Close observation of figure 8 reveals that the cutoff wavenumber boundary for each simulation is distorted as the strain is applied. Wavenumber distortion causes energy to accumulate in long-wavelength (low-wavenumber) modes in the extensional direction, which is consistent with the formation of long coherent vortical structures in the extensional direction (Rogers & Moin Reference Rogers and Moin1987).

Figure 9. Comparison of pre- and post-contraction longitudinal (a,c,e) and transverse (b,d,f) 1-D spectra, for DNS Run 4 at low Reynolds number (a,b), DNS Run 8 at higher Reynolds number (c,d) and the high Reynolds number experiment of AW (e,f), all shown as functions of pre-contraction wavenumbers. Experimental data are reproduced from figure 11 of AW by permission of the authors. Solid lines (black) are for pre-contraction DNS or experiment, dashed (red) for post-contraction DNS or experiment and dashed-dotted (blue) for post-contraction RDT. Insets in each frame show the same spectra but multiplied by the wavenumber (e.g. $k_{1}^{0}\unicode[STIX]{x1D60C}_{22}^{0}(k_{1}^{0})$ ), such that the areas under the curve on log–linear scales give the mean-squared velocities. For the DNS, the transverse spectra $\unicode[STIX]{x1D60C}_{22}^{0}(k_{1}^{0})$ and $\unicode[STIX]{x1D60C}_{33}^{0}(k_{1}^{0})$ are averaged with each other when producing frames (b) and (d). Such averaging is motivated by the axisymmetry of the turbulence.

To characterize the scale-dependent anisotropy between different velocity components, we show in figure 9 the 1-D spectra of the $u_{1}$ and $u_{2}$ velocity fluctuations for low and high Reynolds numbers in the DNS (top and middle rows, respectively), compared with the highest Reynolds number data in the AW experiments (bottom row, from AW figure 11, with permission). The spectra are, based on rationale discussed earlier in § 2.3, all shown as functions of the initial (pre-contraction) wavenumbers. The effect of strain on $\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ is a decrease at low wavenumbers but an increase at high wavenumbers. As the Reynolds number increases, a pronounced rightward shift is seen to develop in the peak of $k_{1}^{0}\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ in both the DNS and experiment (red dashed lines in the insets of frames (c) and (e)). The effect of the mean strain on $\unicode[STIX]{x1D60C}_{22}^{0}(k_{1}^{0})$ is, in contrast, an increase in the spectrum at all values of $k_{1}^{0}$ , with a milder change in the spectrum shape and a weaker Reynolds number dependence. The suppression of $\langle u_{1}^{2}\rangle$ and amplification of $\langle u_{2}^{2}\rangle$ are also indicated by changes in the area under the curves in the inset to each figure.

In figure 9 we have included inviscid RDT results for comparison. While the theory appears to reasonably predict the shape of $\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ at low wavenumbers, it can be seen that the theory fails to capture the rightward shift in $k_{1}^{0}\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ at high Reynolds number as discussed above, while it over-predicts the amplification of $\unicode[STIX]{x1D60C}_{22}^{0}(k_{1}^{0})$ at low wavenumbers. This indicates the physical mechanisms neglected in RDT, namely nonlinear energy transfer, slow pressure strain and viscous dissipation, play a significant role in the evolution of the spectral structure of the Reynolds stresses.

The roles and relative importance of physical processes governing the evolution of the spectra described above can be studied using the spectral budget equation (2.13), which was written for the case of a single Fourier mode. The balance equations for the 1-D spectra parameterized by the initial wavenumbers are obtained by integrating (2.13) over planes perpendicular to $k_{1}$ , and then dividing by the total deformation (like in (2.21)). We then multiply the terms by $k_{1}^{0}$ to study the evolution of the compensated spectra (the insets in figure 9). We are interested in determining what causes the rightward shift in the peak of $k_{1}^{0}\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ at high Reynolds numbers, and what causes the general disagreement between the RDT and DNS spectra at high wavenumbers in all simulations. In figure 10 the balance terms are shown at different times and different Reynolds numbers. As expected for axisymmetric contraction, the production term is negative for the spectrum of $u_{1}$ , but positive for the spectrum of $u_{2}$ . The rapid pressure strain counteracts the generation of anisotropy by taking the opposite sign of the production term, but a net tendency for anisotropy still persists. The nonlinear term is negative for low $k_{1}^{0}$ and positive for high $k_{1}^{0}$ , indicating that there is a forward cascade of energy to high $k_{1}^{0}$ . This forward cascade opposes the tendency for energy to pile up near the $k_{1}^{0}=0$ plane during straining (see figure 8). The overall magnitude of the nonlinear transfer term is greater for $k_{1}^{0}\unicode[STIX]{x1D60C}_{22}^{0}(k_{1}^{0})$ compared to $k_{1}^{0}\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ , which is likely due to the fact that $\langle u_{2}^{2}\rangle$ is amplified during the straining, while $\langle u_{1}^{2}\rangle$ is suppressed.

In our simulations, the magnitude of the strain reaches a maximum and then decreases. To see how this is reflected in the spectral budgets, we present data for Run 8 at two deformations in figure 10(c–f). The production and rapid pressure-strain terms (which depend on the mean strain rate) in the bottom row are smaller than those in the middle row. As the production terms weaken, the relative importance of the nonlinear terms (especially at intermediate and high wavenumbers) increases. This is also the case for the slow pressure-strain term, which peaks at intermediate wavenumbers close to the peak observed for $k_{1}^{0}\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ . We may conclude, therefore, that the slow pressure strain, which is neglected in RDT theory, is the main contributor to rightward shift in $k_{1}^{0}\unicode[STIX]{x1D60C}_{11}^{0}(k_{1}^{0})$ both in our DNS and the experiments of AW. The increasing importance of slow pressure strain at later times noted here is also consistent with a similar feature in figure 6 addressed in § 4.1. On the other hand, the slow pressure-strain term is not as important to the evolution of $k_{1}^{0}\unicode[STIX]{x1D60C}_{22}^{0}(k_{1}^{0})$ , which (besides production) is heavily influenced by nonlinear transfer at intermediate and high wavenumbers.