1 Introduction
The instability of shear flows and their subsequent turbulent dynamics is a fundamental problem of fluid mechanics (see e.g. Schmid & Henningson Reference Schmid and Henningson2000; Drazin & Reid Reference Drazin and Reid2004). The roles of stratification and rotation in aiding or suppressing instabilities and modifying the turbulence have been extensively studied owing to their importance in geophysics and astrophysics. These investigations have utilised theoretical, computational and experimental techniques to shed light on these problems and it is now fair to say that the Taylor–Couette system, the rotating Poiseuille system and the rotating Couette system constitute important paradigm problems for understanding the interaction of rotation both with instabilities and with fully developed turbulence.
Since the pioneering work of Taylor (Reference Taylor1923) much experimental and computational resource has been allocated to Taylor–Couette (TC) flow (see e.g. Koschmieder Reference Koschmieder1993). It is therefore interesting that comparatively little attention has been focussed on rotating plane Couette (RPC) flow, which can be shown to be equivalent to TC flow in the narrow gap limit (Faisst & Eckhardt Reference Faisst and Eckhardt2000), despite it being the simplest model problem in which rotation interacts with shear instabilities. This is no doubt due to the difficulty in designing and constructing experiments for the RPC flow. Despite these difficulties, pioneering experiments have been performed by the Royal Institute of Technology in Stockholm (KTH) group (Tillmark & Alfredsson Reference Tillmark and Alfredsson1992; Hiwatashi et al. Reference Hiwatashi, Alfredsson, Tillmark and Nagata2007; Tsukahara, Tillmark & Alfredsson Reference Tsukahara, Tillmark and Alfredsson2010), who have demonstrated a bewildering array of dynamics (similar to that exhibited in TC flow), with 17 different flow regimes being identified.
Much attention has focussed on the anticyclonic case, where rotation is destabilising, as this leads to interesting interactions and dynamics. In particular anticyclonic rotation drives secondary flows that compete with the turbulence in transporting momentum via Reynolds stresses. This competition has been investigated numerically (Bech & Andersson Reference Bech and Andersson1996, Reference Bech and Andersson1997) to a much lesser extent than TC flow. Recently Salewski & Eckhardt (Reference Salewski and Eckhardt2015) have extended this pioneering numerical work to more turbulent flows and have demonstrated that, for a given value of the Reynolds number, the wall shear stress is a non-monotonic function of rotation and that the momentum flux has contributions from both two-dimensional and three-dimensional structures. We also note here that the RPC system has also been utilised to latch onto fully nonlinear states in non-rotating Couette flow in the pioneering work of Nagata (Reference Nagata1998).
Here, rather than investigate the transitions that occur in RPC flow, we shall utilise this system as a paradigm problem for the interaction of wall-bounded shear flows and rotation to test the effectiveness of the generalised quasilinear (GQL) approximation for three-dimensional systems. Previous evaluations of the GQL approximation have been performed for two-dimensional systems. First Marston, Chini & Tobias (Reference Marston, Chini and Tobias2016) examined GQL in the context of two-dimensional driven turbulence on a spherical surface and
$\unicode[STIX]{x1D6FD}$
-plane and showed it to be more effective than the regular QL approximation in reproducing both the dynamics and statistics of these flows away from equilibrium. This remains true even if only a single extra mode is retained in the large scales. Child et al. (Reference Child, Hollerbach, Marston and Tobias2016) also demonstrated the effectiveness of the GQL approximation for the case of axisymmetric (two-dimensional) helical magnetorotational instability (HMRI) at low magnetic Reynolds number. They showed that the GQL approximations performed significantly better in the quasilinear approximation in reproducing the statistics of this turbulent magnetohydrodynamic flow even for a fairly modest set of modes (3–5 modes) retained in the large scales.
The testing of the effectiveness of the GQL approximation represents a vital piece in the formulation of direct statistical simulation (DSS), a programme for simulating the low-order statistics of turbulent flows directly. Many fluid systems in geophysics, astrophysics and experiments are highly turbulent and involve the interaction of non-trivial mean flows with turbulent fluctuations. It has been argued that for such flows a statistical representation of the behaviour may in certain circumstances be preferable to one based on the calculation of detailed dynamics (Lorenz Reference Lorenz1967). In many situations such statistical representations are based either implicitly or explicitly on the quasilinear approximation, which can be shown to be inadequate as the system moves away from statistical equilibrium (see e.g. Tobias & Marston Reference Tobias and Marston2013). By extending the quasilinear approximation, accurate application in a more extensive parameter regime is possible, providing insight into key processes such as wave – mean-flow interactions.
Another reason for developing statistical representations is that such methods may provide a route to the construction of subgrid models for unresolved dynamics that can be systematically improved. As the short wavelength modes in the GQL approximation interact only with the long wavelength modes, and not amongst themselves, they can be replaced exactly by a generalised second-order closure (Marston et al. Reference Marston, Chini and Tobias2016) called GCE2 (or generalised cumulant expansion at second order). A test of the accuracy of GQL is therefore also a test of the accuracy of the GCE2 form of DSS. The extension and testing of the approximation to three-dimensional flows is crucial to the program of DSS; not only are the three-dimensional flows that we discuss here extremely anisotropic in the streamwise and spanwise direction, they are also characterised by the presence of a forward turbulent cascade in contrast with the two-dimensional models that have been investigated previously.
It is also interesting to relate the GQL approximation to a class of restricted nonlinear (RNL) models that have been utilised in the description of (non-rotating) plane Couette flow (Thomas et al. Reference Thomas, Lieu, Jovanović, Farrell, Ioannou and Gayme2014, Reference Thomas, Farrell, Ioannou and Gayme2015; Bretheim, Meneveau & Gayme Reference Bretheim, Meneveau and Gayme2015). These models have been shown to reproduce the mean profiles of plane Couette flow, including a log layer. RNL corresponds to a particular limit of GQL that is quasilinear in one of the translationally invariant directions and fully nonlinear in the other.
In the next section we present the formulation of the problem and introduce the GQL approximation for three-dimensional systems. In § 3 we describe the dynamics and statistics of weakly and moderately rotating Couette flow as revealed by direct numerical simulation and proceed to evaluate how well the QL and GQL approximation can capture the low-order statistics of this evolution. We conclude with a discussion in § 4.
2 Formulation
2.1 Model, equations and numerical method
Rotating Couette set-up.
Separation in wave vector space of modes into ‘low’ and ‘high’ components. (a) The separation employed in this paper. (b) An alternative separation scheme.
We consider rotating plane Couette flow between two parallel plates as shown in figure 1. The parallel plates at
$z=\pm h$
translate with velocity
$\pm U_{0}$
in the
$x$
-direction and the system rotates about the
$y$
-axis with rotation rate
$\unicode[STIX]{x1D6FA}$
. The flow is described by the non-dimensional three-dimensional incompressible Navier–Stokes equations, where we non-dimensionalise lengths with half the gap-width
$h$
and times with
$h/U_{0}$
. This yields
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+\frac{1}{Ro}\hat{\boldsymbol{y}}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}$
is the velocity field,
$p$
is the pressure and the non-dimensional parameters are the Reynolds number
$Re=U_{0}h/\unicode[STIX]{x1D708}$
and the Rossby number
$Ro=U_{0}/(2\unicode[STIX]{x1D6FA}h)$
. (Note that many authors discussing RPC flow utilise the non-dimensional ‘rotation number’ which is the inverse Rossby number (and also termed
$Ro$
just to add to the confusion).) No-slip boundary conditions are applied at
$z=\pm 1$
in non-dimensional co-ordinates – so we set
$$\begin{eqnarray}\displaystyle \boldsymbol{u}=(u,v,w)=(\pm 1,0,0)\quad \text{at }z=\pm 1, & & \displaystyle\end{eqnarray}$$
whilst the domain is assumed to be periodic in the two horizontal directions at
$x=0,L_{x}$
and
$y=0,L_{y}$
. Equation (2.1) together with boundary conditions (2.2) are solved using a pseudo-spectral code, which employs Fourier expansions in the two horizontal periodic directions and uses Chebychev expansions in the
$z$
-direction. Time stepping is performed via a third-order semi-implicit backward differentiation formula (BDF) scheme – the code utilises the flexible Dedalus partial differential equation (PDE) solving framework (Burns et al.
Reference Burns, Vasil, Brown, Lecoanet and Oishi2016). We use a resolution of
$128\times 128\times 64$
for the two lower Reynolds numbers and
$256\times 256\times 128$
for the
$Re=3000$
case. These numerical resolutions were verified as sufficient by comparing with shorter runs at higher resolution. All calculations are run until a statistically steady state is reached.
2.2 The GQL approximation
The generalised quasilinear approximation was introduced by Marston et al. (Reference Marston, Chini and Tobias2016) to extend the often-utilised quasilinear approximation to include self-consistent interactions of large-scale modes; also see Child et al. (Reference Child, Hollerbach, Marston and Tobias2016) and Constantinou, Farrell & Ioannou (Reference Constantinou, Farrell and Ioannou2016). Here we extend the approximation to the three-dimensional case for which there is translational symmetry in two of the three dimensions. Specifically we perform a idempotent decomposition of the velocity and pressure into large-scale and small-scale modes, i.e. we set, say,
$\boldsymbol{u}(x,y,z)=\boldsymbol{u}_{l}+\boldsymbol{u}_{h}$
, where
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{l}(x,y,z)=\mathop{\sum }_{k_{x}=-\unicode[STIX]{x1D6EC}_{x}}^{\unicode[STIX]{x1D6EC}_{x}}\mathop{\sum }_{k_{y}=-\unicode[STIX]{x1D6EC}_{y}}^{\unicode[STIX]{x1D6EC}_{y}}\boldsymbol{u}_{\boldsymbol{ k}}(z)\text{e}^{\text{i}2\unicode[STIX]{x03C0}k_{x}x/L_{x}+\text{i}2\unicode[STIX]{x03C0}k_{y}y/L_{y}},\quad \boldsymbol{u}_{h}=\boldsymbol{u}-\boldsymbol{u}_{l}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}_{l}$
and
$\boldsymbol{u}_{h}$
are the ‘low’ and ‘high’ wavenumber modes respectively. Similarly
$p(x,y,z)=p_{l}+p_{h}$
. The separation between high and low components in wave vector space is shown in figure 2(a). We note here that when
$\unicode[STIX]{x1D6EC}_{x}=\unicode[STIX]{x1D6EC}_{y}=0$
the low modes are the horizontally averaged (mean) modes and the high modes are the fluctuations about that mean (sometimes termed eddies). Other divisions into ‘low’ and ‘high’ are possible; see for instance figure 2(b). Here we adopt the choice of figure 2(a) as, for this choice, no modes with short length scales (and consequently small time scales) are included in the low modes. This is useful as the statistics of the high modes may then be viewed as a subgrid model for the low modes, and there is no mixture of time scales between the evolving low modes and the subgrid high modes.
Set of three triadic interactions retained in the GQL approximation, out of six in total. The other three triads are discarded (Marston et al. Reference Marston, Chini and Tobias2016). Note that outgoing arrows are identical to incoming arrows with negative wavenumber.
The QL approximation is then characterised by the suppression of certain mode interactions in the dynamics. In particular, though mean/eddy
$\rightarrow$
eddy and eddy/eddy
$\rightarrow$
mean interactions are kept, eddy/eddy
$\rightarrow$
eddy interactions are discarded (as for example done in Herring Reference Herring1963; Lindzen & Holton Reference Lindzen and Holton1968; Dickinson Reference Dickinson1969; Tobias, Dagon & Marston Reference Tobias, Dagon and Marston2011). When
$\unicode[STIX]{x1D6EC}_{x},\unicode[STIX]{x1D6EC}_{y}\neq 0$
the generalised quasilinear approximation differs from QL – the triad interactions that are retained and discarded are selected so as to enable closure and preserve the relevant conservation laws (Marston et al.
Reference Marston, Chini and Tobias2016). Here the interactions low/low
$\rightarrow$
low, high/high
$\rightarrow$
low and low/high
$\rightarrow$
high are retained with all other interactions discarded, as shown in figure 3. Again, retaining only this set of interactions, with both
$\unicode[STIX]{x1D6EC}_{x}=0$
and
$\unicode[STIX]{x1D6EC}_{y}=0$
, is QL. If one of
$\unicode[STIX]{x1D6EC}_{x}$
or
$\unicode[STIX]{x1D6EC}_{y}=0$
then the system is quasilinear in that direction, but may include more triad interactions in the other direction. Furthermore as
$\unicode[STIX]{x1D6EC}_{x},\unicode[STIX]{x1D6EC}_{y}\rightarrow \infty$
the GQL system consists solely of fully interacting low modes and returns to fully nonlinear direct numerical simulation. Finally, if the system has one of
$\unicode[STIX]{x1D6EC}_{x}$
or
$\unicode[STIX]{x1D6EC}_{y}=0$
with the other cutoff set so that all the modes are retained then the system is quasilinear in one direction of symmetry and fully nonlinear in the other; the RNL systems described in the introduction with
$\unicode[STIX]{x1D6EC}_{x}=0$
and all modes kept in the
$y$
-direction are in this class of approximation. Thus GQL is a controlled approximation that systematically improves upon QL.
To be clear, when the GQL approximations is implemented numerically the same spectral resolution is employed, but certain triad interactions are suppressed in the streamwise and spanwise directions (all modes and interactions are retained in the
$z$
-direction). There are a number of ways of suppressing the nonlinear interactions in a spectral code, however we choose to derive separate equation sets to be satisfied by the low and high modes; these take the form
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}_{l}+\langle (\boldsymbol{u}_{l}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{l}+(\boldsymbol{u}_{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{h}\rangle _{l}+\frac{1}{Ro}\hat{\boldsymbol{y}}\times \boldsymbol{u}_{l}=-\unicode[STIX]{x1D735}p_{l}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{l}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}_{h}+\langle (\boldsymbol{u}_{l}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{h}+(\boldsymbol{u}_{h}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}_{l}\rangle _{h}+\frac{1}{Ro}\hat{\boldsymbol{y}}\times \boldsymbol{u}_{h}=-\unicode[STIX]{x1D735}p_{h}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{h}, & \displaystyle\end{eqnarray}$$
$\langle \,\rangle _{l}$
and
$\langle \,\rangle _{h}$
represent the projection of the nonlinear terms onto the low and high wavenumbers respectively. We note that this is not the most efficient way of implementing GQL and we stress that this formalism does not therefore currently offer any saving over direct numerical simulation (DNS). Rather we are using this approach to determine the effectiveness of the GQL approximation. If GQL can be demonstrated to yield and accurate description of the statistics then a statistical description of the high modes via direct statistical simulation will offer a significant reduction in computational cost.
Enstrophy
$\int _{V}\unicode[STIX]{x1D714}^{2}\,\text{d}V$
as a function of cutoff wavenumbers
$\unicode[STIX]{x1D6EC}_{x}$
and
$\unicode[STIX]{x1D6EC}_{y}$
.
3 Results
3.1 DNS
We begin the results section by describing the behaviour of the RPC system as revealed by DNS. For this geometry the non-dimensional mean absolute vorticity is given by
$1/Ro+\text{d}\overline{U}/\text{d}z$
(where
$U(z)$
is the mean flow) and so we expect the instability to be enhanced for negative Rossby number and suppressed for positive
$Ro$
(we note in this coordinate system this is opposite to that for Bech & Andersson (Reference Bech and Andersson1996) and subsequent papers). Much of the nonlinear behaviour can – and has – been interpreted in terms of turbulence and secondary flows mixing the mean absolute vorticity in the core of the flow to suppress the instability. For a thorough discussion of the competition between turbulent transport and that provided by the secondary flows see Bech & Andersson (Reference Bech and Andersson1997), Suryadi, Segalini & Alfredsson (Reference Suryadi, Segalini and Alfredsson2014) and the references therein.
In this paper, we set
$L_{x}=10\unicode[STIX]{x03C0}$
,
$L_{y}=4\unicode[STIX]{x03C0}$
, and vary
$Ro$
and
$Re$
as shown in table 1. The fiducial run is weakly rotating
$Ro=-100$
and at moderate
$Re=1300$
. For this choice of parameters the solution is particularly complicated and so presents a significant test for GQL. The solution has three constituent parts; a mean flow, a secondary flow and a turbulent part (Bech & Andersson Reference Bech and Andersson1996). The secondary flow takes the form of (initially three and then two) time-dependent turbulent sinuous rolls largely aligned with the
$x$
-axis (as shown in figure 4
a). Figure 5 shows the mean flow
$U(z)$
and mean shear
$\text{d}U/\text{d}z$
as a function of depth. Here the averaging has been performed over horizontal planes and the last third of the calculation (
$500\leqslant t\leqslant 750$
); the solutions are run up to
$t=750$
from rest. All statistical quantities such as spectra and second cumulants are also averaged in time over this interval. These rolls and the anisotropic turbulence act so as to redistribute the vorticity, forming a core with a weak shear with two matching layers of strong shear (see the blue line of figure 5). The statistically steady state is achieved with the energy input through the boundaries being matched on average by the viscous dissipation, given by
$D_{\unicode[STIX]{x1D708}}=(1/Re)\int \unicode[STIX]{x1D714}^{2}\,\text{d}V$
; we shall use this dissipation rate as a measure to evaluate the effectiveness of the GQL approximation (see table 1). For moderate rotation (
$Ro=-10.0$
) the solution is less turbulent and the (less time-dependent) secondary flow has three relatively straight rolls (figure 4
b). These are more efficient at transporting the mean flow (green line of figure 5) and so the core has a gradient closer to zero with steeper shears in the wings. Finally we note that for more turbulent, weakly rotating flow (
$Re=3000$
,
$Ro=-100$
) both the secondary flow, which again has (on average two) rolls largely aligned with the
$x$
-axis (see figure 4
c) and the turbulence are acting strongly to transport the momentum (red lines of figure 5).
We conclude this section by noting that, as discussed by Diamond et al. (Reference Diamond, Itoh, Itoh and Hahm2005), a measure of the departure from quasilinearity in turbulent fluid systems is often given by the Kubo number
$R=u_{rms}\unicode[STIX]{x1D70F}_{c}/l_{c}$
where
$\unicode[STIX]{x1D70F}_{c}$
is the correlation time of the turbulence and
$l_{c}$
its correlation length. The Kubo number is therefore a measure of the ratio of the advective (nonlinear) term to the time-derivative term in the equation for the fluctuating velocity. For short correlation time turbulence
$R$
is asymptotically small, the nonlinear term may be formally neglected and the system is well represented by a quasilinear approximation. We calculate the correlation time using the first zero of the temporal auto-correlation function and the correlation length via a typical length scale of variation of the two-point (in space) correlation function. For these flows we note that
$R\sim O(0.5)$
, with the Kubo number being smallest for case B. Thus for these systems there is no a priori reason why the quasilinear approximation should hold.
Instantaneous snapshots of the streamwise velocity
$u(x,y)$
at
$z=0.1$
for
$Re=1300$
and (a)
$Ro=-100$
, (b)
$Ro=-10$
and (c)
$Re=3000$
,
$Ro=-100$
.
(a) Mean flows
$U(z)$
and (b) mean shears
$\text{d}U/\text{d}z$
as a function of
$z$
(vertical axis) for
$Re=1300$
,
$Ro=-100$
(blue, dot-dot-dashed line),
$Re=1300$
,
$Ro=-10$
(red dashed line) and
$Re=3000$
,
$Ro=-100$
(green dot-dashed line). We note that the red and green curves lie nearly on top of each other.
3.2 Evaluation of the GQL approximation
Instantaneous snapshots of the evolution of the flow for case A. The snapshots show the streamwise velocity
$u(x,y)$
at
$z=0.1$
and
$t=750$
. The flows are intermittent; therefore it is more meaningful to compare statistics than snapshots.
Having established the dynamics of the RPC system via DNS, we evaluate the GQL approximation for various truncations
$(\unicode[STIX]{x1D6EC}_{x},\unicode[STIX]{x1D6EC}_{y})$
in the streamwise and spanwise directions. Clearly the presence of mean flows and rotation make the turbulence highly anisotropic and so there is no a priori reason why
$\unicode[STIX]{x1D6EC}_{x}$
should be set equal to
$\unicode[STIX]{x1D6EC}_{y}$
. This yields a huge range of choice for these cutoffs and we present only a representative few below. Previous evaluations of GQL have only considered two-dimensional models, with only one direction of splitting into low and high modes. We shall evaluate the degree of success of the approximation by determining how well they reproduce the low-order statistics of the system; mean flows, two-dimensional spectra, two-point correlation functions and total enstrophy.
Mean flow
$U(z)$
as a function of
$z$
(vertical axis) for
$Re=1300$
, (a)
$Ro=-100$
, (b)
$Ro=-10.0$
, (c)
$Re=3000$
,
$Ro=-100$
. Shown are DNS (dashed black line),
$(\unicode[STIX]{x1D6EC}_{x},\unicode[STIX]{x1D6EC}_{y})=(0,0)$
(navy),
$(0,64)$
(red),
$(64,0)$
(cyan),
$(3,3)$
(green),
$(5,5)$
(purple),
$(3,10)$
(pink),
$(6,6)$
(grey), and
$(10,3)$
(brown – top panel only). Note only
$-0.2\leqslant U(z)\leqslant 0.2$
is shown to give a close up of the region of reduced shear.
Although we are evaluating the effectiveness of GQL in describing the statistics of the rotating Couette system it is of interest to determine how the varying degrees of approximation perform in reproducing some of the dynamics of the system. In particular, for case A the system evolves significantly in time and produces coherent structures in the form of rolls largely aligned with the
$x$
-axis. At short time the solution has three largely straight rolls, which can be understood via the linear instability of the modified shear flow (Bech & Andersson Reference Bech and Andersson1997). For longer times these rolls evolve nonlinearly, interacting with the turbulence until the coherent rolls take the form of two wavy rolls as shown in figure 6(a), which is a snapshot of the streamwise velocity (measured on the plane
$z=0.1$
) at the end of the simulation (
$t=750$
). Figure 6(b–h) show how well severe truncations of GQL perform in capturing the dynamics of the system. The quasilinear approximation in figure 6(b) appears to capture the initial instability of the modified shear state to three rolls, but not the subsequent evolution to two wavy rolls. GQL performs better, even at severe truncations. We draw attention to the snapshots in figure 6(e–h) which suggest that these truncations are capable of capturing the long-time nonlinear evolution including the merging of the coherent structures. This will be quantified by measuring the statistics of the flows.
Figure 7 shows how well the various degrees of approximation perform in reproducing the mean flow (or first cumulant –
$U(z)$
) giving a close up of the core of the mean flow. For the weakly rotating case – in a – the suppression of the shear is greatly overestimated by the quasi-linear approximation and fairly well approximated by GQL at all truncations. Note though that GQL performs poorly when either one of
$\unicode[STIX]{x1D6EC}_{x}$
or
$\unicode[STIX]{x1D6EC}_{y}=0$
, with the suppression of the shear in the core being overestimated, for reasons we shall discuss later. For the moderately rotating flows (case B, shown in figure 7
b) the shear in the core of the flow is well reproduced by all the approximations, although that near the boundaries is significantly overestimated by QL. We believe that the improvement over that for case A in all the methods arises because the increase in rotation leads to a suppression of the role of the turbulence and a diminishment of the variability of the secondary flows. The secondary flow rolls are relatively straight, are persistent and their role in transporting momentum is more easily captured by all the approximations. For this flow, although the secondary flow is still present, it is more ordered. Finally figure 7(c) shows the mean flows for the more turbulent, weakly rotating case (
$Re=3000$
,
$Ro=-100$
) with the various degrees of approximation. Again here the QL approximations significantly overestimate the degree of mixing in the core, and consequently overestimate the shear in the boundary layers. This deficiency is again shared by the calculations that are quasilinear in either the streamwise or spanwise direction – with the full GQL approximation calculations performing better.
A quadratic statistic that is often of interest is the energy dissipation rate, which is proportional to the enstrophy. Table 1 gives the total enstrophy for the three runs at various degrees of approximation. The first thing to note is that the enstrophy increases with either increasing rotation rate or increased Reynolds number. This is to be expected as in both of these cases the small scales are more prevalent than for the moderately turbulent and weakly rotating case (case A); see the discussion of the two-dimensional spectra to follow. The dissipation must, of course, in a statistically steady state balance the energy input through the boundaries, which is given by the wall shear stress. Reference back to figure 5(b) shows how the wall shear stress increases from case A to B to C in the same ratio as the enstrophy in table 1. Flows under the QL approximation significantly overestimate the enstrophy; the overestimation of the mixing in the core of the flow leads to strong gradients at the boundaries. These in turn lead to greater wall shear stress required to satisfy the boundary conditions. This leads to a consequent increase in the dissipation via an increase in the mean enstrophy. Moreover flows that are quasilinear in one direction (i.e. either
$\unicode[STIX]{x1D6EC}_{x}=0$
or
$\unicode[STIX]{x1D6EC}_{y}=0$
) also tend to perform poorly, whilst GQL runs with
$\unicode[STIX]{x1D6EC}_{x}$
,
$\unicode[STIX]{x1D6EC}_{y}\geqslant 3$
perform better for all runs. We note at this point that the approach to the DNS values of the enstrophy is non-monotonic with increase in
$\unicode[STIX]{x1D6EC}_{x}$
and
$\unicode[STIX]{x1D6EC}_{y}$
; this is discussed further below.
The reason for this improvement of GQL over QL can be ascertained from an examination of the two-dimensional spectra of
$u$
; an example for the DNS of case A is given in figure 8(a). These are calculated for a fixed horizontal plane
$z=0.1$
just above the midplane and averaged over time. These show the anisotropy of the flow in the horizontal plane for the full solutions; the secondary flow apparent in the snapshot of figure 4(a) dominates this spectrum at
$(k_{x},k_{y})=(0,2)$
and
$(2,2)$
, though there is much energy in the surrounding wavenumbers. Energy cascades from the mean flows and secondary flows to smaller scales in an anisotropic fashion, with the energy falling more rapidly in the streamwise (
$k_{x}$
) direction than in the spanwise (
$k_{y}$
) direction. The shear inhibits the spectral transfer from large scales anisotropically, as expected.
Two-dimensional spectra for
$u$
in the horizontal plane
$z=0.1$
. Case A is shown (table 1). The energy is shown on a (base 10) logarithmic scale given by the colour bar, with values below the minimum of the colour bar not assigned a colour (so they appear white).
The corresponding spectra for the QL and GQL approximations are shown in figure 8(b–h). What is clear from figure 8(b), which shows the spectrum for QL, is that under this approximation the energy is concentrated at large scales in both directions. Because of the quasilinear approximation, these wavenumbers (which of course are considered ‘high’) must be those to which energy feeds in directly from an instability of the mean flow. This instability may have occurred at any stage of the evolution to the statistically steady state. Once these high modes have reached finite amplitude they are sustained by the saturated shear and act to alter the shear flow via self-interaction (though of course some modes that arise from an instability of an intermediate shear may not be sustained by the final shear profile – and these will have decayed away). The ‘high’ modes do not interact with each other though and so are not able to redistribute energy among themselves. It is clear therefore that the mean flow is not directly unstable to modes at small spatial scales and so there is no mechanism for these modes to receive energy.
Figure 8(c,d) demonstrate that if the truncation is quasilinear in one direction then the cascade is severely inhibited in that direction, though transfer can occur in the orthogonal direction, because we have chosen to keep all the modes in this orthogonal direction as ‘low’ modes. Thus the energy transfer is artificially anisotropic on utilising such models. We note again here that ‘minimal models’ (Thomas et al.
Reference Thomas, Farrell, Ioannou and Gayme2015) have been proposed for describing the evolution of shear instabilities that keep all spanwise modes and only the mean streamwise mode (equivalent to
$\unicode[STIX]{x1D6EC}_{x}=0$
,
$\unicode[STIX]{x1D6EC}_{y}=$
all modes for GQL). What we have shown here is that, at least for this problem, this model performs poorly. One would expect the truncation that is quasilinear in the spanwise direction to perform poorly, as indeed it does.
However, provided
$\unicode[STIX]{x1D6EC}_{x}$
,
$\unicode[STIX]{x1D6EC}_{y}$
are both
$\neq 0$
(shown in figure 8
e–h) then energy can be transported via non-local interactions to modes at small spatial scales, with the consequent improvement of the spectra in comparison with the direct numerical simulation. This energy transport involves the successive scattering of high mode energy off the low modes into high modes with a different wavenumber. For example, if
$\unicode[STIX]{x1D6EC}_{x}=\unicode[STIX]{x1D6EC}_{y}=3$
then a high mode with wavenumber
$(6,7)$
can interact with a low mode with wavenumber
$(1,1)$
to put energy into a high mode at wavenumbers
$(7,8)$
. In this non-local (in wavenumber space) interaction energy can be transferred to high wavenumbers to which the evolving shear flow is not directly unstable. Indeed successive interactions can redistribute energy among all the high modes as shown by the spectra. This is, we believe, the key advantage of the GQL approximation over the QL approximation.
Figure 9 shows how well the various approximations perform in reproducing a second cumulant (or two-point correlation functions) of the RPC system for case A. The second cumulant for the streamwise velocity
$u$
is defined in terms of a spatial average over the two directions of translational invariance (
$x$
and
$y$
),
$$\begin{eqnarray}\displaystyle c_{uu}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702},z_{1},z_{2})=\frac{1}{L_{x}L_{y}}\int _{0}^{L_{x}}\int _{0}^{L_{y}}u^{\prime }(x,y,z_{1})u^{\prime }(x+\unicode[STIX]{x1D709},y+\unicode[STIX]{x1D702},z_{2})\,\text{d}x\,\text{d}y, & & \displaystyle\end{eqnarray}$$
where
$u^{\prime }=u-U(z)$
. In figure 9 we plot, as a function of
$\unicode[STIX]{x1D709}$
and
$\unicode[STIX]{x1D702}$
, a two-dimensional slice at
$z_{1}=z_{2}=0.1$
of this four-dimensional object. Information contained in this piece of the second cumulant contributes to the power spectra of figure 8 – however the (logarithmically plotted) spectral power emphasises the transfer of energy to small scales, while the second cumulant makes the ordering wavenumber readily apparent. We note that the full second cumulant, being four-dimensional, also contains information that is not present in the full three-dimensional spectral power. Recall from the discussion of figure 4(a) that the major contribution to the streamwise velocity comes from the secondary flow, which, after a long transient with strong
$(0,3)$
and
$(2,3)$
components, settles into a dynamic
$(0,2)$
and
$(2,2)$
component. This pattern is selected in the second cumulant from DNS shown in figure 9(a). Interestingly the second cumulant from the QL approximation is dominated by the
$(0,3)$
mode, as are some of the cumulants using GQL (see e.g. figure 9(f) for
$\unicode[STIX]{x1D6EC}_{x}=3$
,
$\unicode[STIX]{x1D6EC}_{y}=10$
). As stated earlier, this mode arises through the linear instability of the saturated shear flow, which is captured by QL. QL is however not able to capture the subsequent nonlinear transition to the wavy roll state. The GQL integrations that (at least for this particular parameter choice) perform best in reproducing this second cumulant are shown in figure 9(d,e,h) – interestingly those that include the most modes in the streamwise direction. This behaviour is non-monotonic for small truncation levels (the transition between the straight and wavy rolls is subtle and not always captured by severe truncations). However at moderate cutoffs GQL performs very well – see for example figure 9(i).
(a) Mean shears
$\text{d}U/\text{d}z$
as a function of
$z$
(vertical axis) for
$Re=1300$
,
$Ro=-100$
(red), and for the underresolved case. Note how the underresolved case underestimates the wall shear stress and hence the enstrophy. (b) Second cumulant for underresolved case which should be compared with the fully resolved case in figure 9(a) and GQL with small cutoffs in figure 9(e,h).
We conclude this section by stressing that the GQL approximations performs significantly better than a poorly resolved direct numerical simulation. For most values of the truncation utilised here, a DNS of the same resolution performed pseudo-spectrally will blow up owing to lack of resolution. We note that a
$16\times 16\times 64$
DNS is numerically stable for case A, but massively underresolves the dissipation (yielding an enstrophy of
$9.0$
). Figure 10 shows the shear and the second cumulant from the underresolved run. Clearly at this resolution the scheme is incapable of reproducing both the mean shear near the boundaries and the spatial structure of the two-point correlation function. Thus, as expected, it performs significantly worse than a more severely truncated GQL run.
4 Discussion and conclusions
In this paper we have demonstrated that the generalised quasilinear approximation introduced by Marston et al. (Reference Marston, Chini and Tobias2016) performs significantly better than the quasilinear approximation for the problem of turbulent rotating Couette flow at moderate Reynolds number. We have assessed the efficacy of the approximation in reproducing mean flows, spectra, dissipation rates and two-point correlation functions. The approximation performs better than proposed minimal models that are quasilinear in the streamwise direction, even if only a few modes are kept as ‘low modes’ in the streamwise and spanwise directions. That GQL shows a substantial and systematic improvement over QL for a three-dimensional system with a forward cascade of energy from large to small scales is interesting, because nonlinear interactions that transfer energy only between high modes are dropped in the GQL approximation. The result lays the foundation for direct statistical simulation of three-dimensional wall-bounded flows because the high modes can be closed statistically at second order, yielding a new class of subgrid models that can be systematically improved with increasing cutoff
$\unicode[STIX]{x1D6EC}$
. Future work includes assessing the applicability of GQL for other wall-bounded systems such as Taylor–Couette and pipe flow, and completion of the program of direct statistical simulation by implementation of the realizable GCE2 closure that is exact for GQL.
Acknowledgements
The authors would like to thank G. Chini and J. Oishi. The computations were performed on ARC1 and ARC2, part of the High Performance Computing facilities at the University of Leeds, UK. This work was supported in part by the US NSF under grant no. DMR-1306806 (J.B.M.).
