2.1 The governing equations

We solve the incompressible LES momentum and continuity equations,

where $\overline{u}_{i}$ are the resolved (large-scale) velocities, and $\overline{p}$ is a modified resolved kinematic pressure that includes the diagonal part of the SG stress tensor. The eddy viscosity is written as $\unicode[STIX]{x1D708}_{t}=l_{S}^{2}|\overline{\unicode[STIX]{x1D70E}}|$ , in terms of a static length scale $l_{S}$ and of the local resolved strain rate (Smagorinsky Reference Smagorinsky1963), where $|\overline{\unicode[STIX]{x1D70E}}|^{2}=2\overline{\unicode[STIX]{x1D70E}}_{ij}\overline{\unicode[STIX]{x1D70E}}_{ij}$ and $\overline{\unicode[STIX]{x1D70E}}_{ij}=(\unicode[STIX]{x2202}_{j}\overline{u}_{i}+\unicode[STIX]{x2202}_{i}\overline{u}_{j})/2$ is the strain-rate tensor of the filtered velocity. We use $(x,y,z)$ to represent the streamwise, vertical (cross-shear) and spanwise coordinates, respectively, and $(u,v,w)$ to denote the respective velocity components. Whenever convenient, as in (2.1), this notation is substituted by subindices, in which case repeated indices imply summation. This is always the case for the vorticity components $\overline{\unicode[STIX]{x1D714}}_{i}$ . There is no explicit filtering in our code, and the smallest flow scales are controlled either by the grid or by the eddy viscosity.

It is usual in most LES to include the molecular viscosity $\unicode[STIX]{x1D708}$ as part of the right-hand side of (2.1), and to write the length $l_{S}$ in terms of some convenient grid spacing $\unicode[STIX]{x1D6E5}_{g}$ and of a ‘universal’ Smagorinsky constant $C_{S}$ . None of those devices are needed here, and will not be used. The main role of the molecular viscosity in LES is to enforce the boundary conditions when $\unicode[STIX]{x1D708}_{t}$ vanishes at the wall. Since there are no walls in our problem, we set $\unicode[STIX]{x1D708}=0$ , and all our simulations run at infinite molecular Reynolds number.

Similarly, the role of $\unicode[STIX]{x1D6E5}_{g}$ in SG models is to accommodate variable grid spacings, assuming that most of the spectral content of the velocity gradients is concentrated at the grid scale. Again, the statistical homogeneity of the flow implies that our grids are uniform, and makes this complication unnecessary.

In fact, the essential role of the eddy viscosity is to introduce a minimum length scale that takes the role of the Kolmogorov length, $\unicode[STIX]{x1D702}$ , independently of the grid. If we estimate the energy dissipation by the SG model as $\unicode[STIX]{x1D700}=\langle \unicode[STIX]{x1D708}_{t}|\overline{\unicode[STIX]{x1D70E}}|^{2}\rangle =l_{S}^{2}\langle |\overline{\unicode[STIX]{x1D70E}}|^{3}\rangle$ , and define an effective Kolmogorov length as $\unicode[STIX]{x1D702}_{t}=(\langle \unicode[STIX]{x1D708}_{t}\rangle ^{3}/\unicode[STIX]{x1D700})^{1/4}$ , it follows from the definition of $\unicode[STIX]{x1D708}_{t}$ that $\unicode[STIX]{x1D702}_{t}\approx l_{S}$ . Here $\langle \cdot \rangle$ represents time and space average, and will occasionally be replaced by $\langle \cdot \rangle _{xz}$ to represent the $y$ -dependent average over horizontal planes, or by $\langle \cdot \rangle _{c}$ to represent $\langle \cdot \rangle _{xz}$ particularised to the central plane $y=0$ . From now on, capital letters are reserved for mean quantities and lower-case ones for fluctuations with respect to that mean, as in $U=\langle \overline{u}\rangle _{xz}$ and $\overline{u}=U+u$ . Primes denote root-mean-squared (r.m.s.) intensities of the fluctuations.

Since most of the computational effort in DNS is dedicated to resolving the smallest dissipative scales, the introduction of a cut-off length in LES drastically reduces the size of the linear systems that have to be solved to obtain invariant solutions.

For the parallel flows that are the subject of this paper, (2.1) can be integrated to

where $u_{\unicode[STIX]{x1D70F}}$ is independent of $y$ and can be used to scale the velocities. Variables in this scaling are denoted by a ‘+’ superscript.

2.2 LES with sinuous symmetry

We study a flow with a nominally uniform mean shear $S=\text{d}U/\text{d}y$ , which is a given constant, in a parallelepiped-shaped computational domain that is periodic in the $(x,z)$ directions, and periodic between points of the lower and upper boundaries that are uniformly shifted in time by the shear (Gerz, Schumann & Elghobashi Reference Gerz, Schumann and Elghobashi1989). The mean flow for turbulence statistics is $Sy$ plus the vertically periodic correction, which is an approximately linear profile as evaluated later in this section. The numerical code for DNS is described in detail in Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016). Equations (2.1) and (2.2) for LES are formulated in terms of the Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016). Equations (2.1) and (2.2) for LES are formulated in terms of the vertical vorticity and of the Laplacian of the vertical velocity (Kim, Moin & Moser Reference Kim, Moin and Moser.1987; Hughes, Oberai & Mazzei Reference Hughes, Oberai and Mazzei2001). The discretisation uses $2/3$ -dealiased Fourier expansions in $(x,z)$ , and sixth-order spectral-like compact finite differences in $y$ , with the shear-periodic boundary conditions embedded in the compact finite-difference matrices for each Fourier mode. As explained in Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016), this avoids recurrent remeshing and the resulting secular loss of enstrophy over long integration times.

There are three dimensionless parameters: a Reynolds number, and two box aspect ratios, $A_{xz}=L_{x}/L_{z}$ and $A_{yz}=L_{y}/L_{z}$ , where the $L_{j}$ are the dimensions of the computational domain along the three coordinate directions. It was shown in Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016) that the correct scales for the large-scale length and velocity of SS-HST are based on the spanwise box dimension, $L_{z}$ and $SL_{z}$ , so that the Reynolds number is $Re_{z}\equiv SL_{z}^{2}/\unicode[STIX]{x1D708}$ in DNS. We will extend this definition to LES using the effective mean eddy viscosity, $\langle \unicode[STIX]{x1D708}_{t}\rangle$ , replacing the molecular viscosity $\unicode[STIX]{x1D708}$ . However, it is more convenient to use the length-scale ratio $R_{S}=L_{z}/l_{S}$ in LES. For example, Piomelli, Rouhi & Geurts (Reference Piomelli, Rouhi and Geurts2015) recently developed a grid-independent LES using as parameter the ratio of the small (effective Kolmogorov) and integral scales, which is essentially the inverse of $R_{S}$ . It can be extended to DNS by substituting $\unicode[STIX]{x1D702}$ for $l_{S}$ ( $L_{z}/\unicode[STIX]{x1D702}\sim Re_{z}^{3/4}$ ) and it turns out that, in DNS of SS-HST, $Re_{z}\approx 4R_{S}^{4/3}$ (Sekimoto et al. Reference Sekimoto, Dong and Jiménez2016). We will also use the integral scale $L_{0}=u_{0}^{3}/\unicode[STIX]{x1D700}$ , where $u_{0}^{2}=\langle u_{i}u_{i}\rangle /3$ .

We study in this paper LES turbulence and invariant solutions in a subspace defined by enforcing two ‘sinuous’ symmetries: (I) a reflection with respect to the plane of $z=0$ followed by a streamwise shift by $L_{x}/2$ ,

and (II) a rotation by $\unicode[STIX]{x03C0}$ around $x=y=0$ followed by a spanwise shift by $L_{z}/2$ ,

Note that no translational symmetries are allowed in this subspace, so that travelling waves are excluded. Moreover, $(\text{I})$ and $(\text{II})$ together with the boundary conditions enforce that the instantaneous plane-averaged streamwise velocities at the top and bottom of the box are $U(\pm L_{y}/2)=\pm SL_{y}/2$ .

Before proceeding to seek equilibrium solutions, the effect of the symmetry restriction was tested in several LES of symmetric SS-HST, which are summarised in table 1. Table 1 also includes two reference unconstrained DNS from Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016), and shows that the length and velocity scales found in DNS also work well in the symmetric LES. Although not included in table 1, the effective Kolmogorov scale in the LES is found to be $\unicode[STIX]{x1D702}_{t}/l_{S}\approx 0.9$ , in approximate agreement with our analysis in § 2.1. Interestingly, both $\unicode[STIX]{x1D702}_{t}/l_{S}\approx 0.9$ –1.0 and $L_{0}/L_{z}\approx 0.4$ –0.7 remain good approximations in the localised equilibrium solutions described below. Note that the vertical translational invariance is broken by the symmetries (2.4) and (2.5), so that a uniform mean shear is no longer guaranteed. Even so, the LES of symmetric turbulence retain an approximately linear mean profile $(|U-Sy|\approx 0.01SL_{z})$ and approximately constant fluctuation profiles. We will see below that the same is not generally true for the profiles of the equilibrium solutions.

Figure 1 displays the longitudinal streamwise-velocity spectrum of the symmetric LES turbulence, compared to that of the reference DNS. Given the different Reynolds numbers and techniques, the agreement is good, showing that symmetry does not influence turbulence greatly. Note that the three LES in figure 1(b) collapse well in terms of $l_{S}$ , and that the resolution of the simulations is fine enough to resolve the smallest LES scales.

Table 1.

Parameters for the turbulence runs. L32 and M32 are reference DNS of SS-HST from Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016). LES $_{s,m,t1,t2,t3}$ are the present LES in the symmetric subspace defined by (I)  $+$  (II) in (2.4) and (2.5). The effective Reynolds number, $Re_{z}$ , and the Kolmogorov viscous scale, $\unicode[STIX]{x1D702}$ , are computed with the molecular viscosity in DNS, and with the averaged eddy viscosity in LES. The scale ratio $R_{S}$ is $L_{z}/l_{S}$ in LES and $L_{z}/\unicode[STIX]{x1D702}$ in DNS. The Reynolds number $Re_{\unicode[STIX]{x1D706}}\equiv u_{0}\unicode[STIX]{x1D706}/\langle \unicode[STIX]{x1D708}_{t}\rangle$ is that based on the Taylor microscale $\unicode[STIX]{x1D706}=\sqrt{15}u_{0}/\langle \unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}\rangle ^{1/2}$ .

Figure 1.

Streamwise velocity spectra $E_{uu}(k_{x})$ . (a) Large-scale normalisation, $E_{uu}/(u^{\prime 2}L_{z})$ , as a function of $k_{x}L_{z}$ . (b) Small-scale normalisation, $E_{uu}^{\ast }\equiv \unicode[STIX]{x1D700}^{-2/3}l^{-5/3}E_{uu}(k_{x}l)$ , with $l=l_{S}$ for LES and $l=\unicode[STIX]{x1D702}$ for DNS. Curves/symbols: – –  $\triangledown$ – – (grey), DNS (L32); —— (grey), DNS (M32); – –  $\bullet$ – – (blue),  $R_{S}=52.6$ ; — $\circ$ — (red),  $R_{S}=101.6$ ; — $\rhd$ — (green),  $R_{S}=203$ . The slope of the chain-dotted diagonal is $-5/3$ .

The size of the collocation grids for the turbulence LES and for typical equilibrium solutions are given in table 2. It follows from the definition of the different quantities that the collocation resolution is $\unicode[STIX]{x0394}x/l_{S}=R_{S}A_{xz}/N_{x}$ , $\unicode[STIX]{x0394}y/l_{S}=R_{S}A_{yz}/N_{y}$ and $\unicode[STIX]{x0394}z/l_{S}=R_{S}/N_{z}$ , where the $N_{j}$ are the collocation points along the three coordinate directions. For example, in the lowest-Reynolds-number LES $_{s}$ in tables 1 and 2, this results in $(\unicode[STIX]{x0394}x,\unicode[STIX]{x0394}y,\unicode[STIX]{x0394}z)=(3.6,1.4,2.7)l_{S}$ . If our SG model is interpreted in terms of the usual Smagorinsky formula $\unicode[STIX]{x1D708}_{t}=(C_{S}\unicode[STIX]{x1D6E5})^{2}|\overline{\unicode[STIX]{x1D70E}}|$ , with $\unicode[STIX]{x1D6E5}$ a representative grid spacing, it follows that $C_{S}=l_{S}/\unicode[STIX]{x1D6E5}=0.42$ , which is a relatively high value for usual LES practice. In this sense, all our LES are overdamped. Their dynamics are controlled not by the grid, but by the SG model, and can be expected to be approximately independent of the resolution. They also only resolved the integral scales of the flow, and very little of its inertial range. The effective Smagorinsky constant for the different simulations is included in table 2.

Table 2.

Grid information for the sinuous-symmetric equilibrium and turbulence LES. EQ(L) and EQ(U) represent lower- and upper-branch equilibrium solutions, respectively, and the turbulence LES are defined in table 1. The effective Smagorinsky constant $C_{S}$ is discussed in § 2.2, with the filter scale defined as $\unicode[STIX]{x1D6E5}_{g}\equiv \sqrt[3]{\unicode[STIX]{x0394}x\unicode[STIX]{x0394}y\unicode[STIX]{x0394}z}$ .

At box aspect ratios $(A_{xz},A_{yz})=(3,1.33)$ , symmetric self-sustaining LES turbulence exists for $R_{S}\gtrsim 65$ . There is a transitional range, $45\lesssim R_{S}\lesssim 65$ , in which the kinetic energy of the flow increases and decreases several times before decaying to laminar after $St=O(1000)$ . Both the transitional range and the decay time depend somewhat on $A_{yz}$ . We will see in § 4 that LES turbulence in taller boxes, $A_{yz}\approx 2$ –3, survives in the range $50\lesssim R_{S}\lesssim 60$ , but collapses intermittently to vertically localised turbulence.

2.3 Searching for invariant solutions in LES

Our main interest is to characterise invariant solutions in LES of SS-HST. Strictly equilibrium solutions are technically impossible in this system. The shear-periodic boundary condition slides the upper computational boundary with a velocity $SL_{y}$ with respect to the lower one, and the numerical configuration only repeats itself after a ‘box’ period $T_{s}\equiv L_{x}/(SL_{y})$ . It was shown in Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016) that this periodic forcing does not interfere strongly with the turbulent solutions as long as the aspect ratios are kept in the range $2<A_{xz}\lesssim 5$ and $A_{xz}\lesssim 2A_{yz}$ . We will approximately respect these constraints here, and find solutions that are numerically indistinguishable from fixed points. Vertically localised solutions can be equilibria at the limit of $A_{yz}\rightarrow \infty$ , since they are independent of the shear-periodic boundary condition, but it should be remembered that all of them are conceptually periodic orbits in a finite computational domain. All of the statistics discussed below are averages over a box period.

Solutions are computed using the Newton–Krylov–hookstep method (Viswanath Reference Viswanath2007) on $\boldsymbol{f}^{T}(\boldsymbol{x})-\boldsymbol{x}=0$ , where $\boldsymbol{x}$ is the vector of independent variables, and $\boldsymbol{f}^{T}:\boldsymbol{x}(0)\rightarrow \boldsymbol{x}(T)$ is the integration of (2.1) over time $T$ , using the evolution code described in the previous section. Because of the periodicity mentioned above, the search time is always an integer multiple of the box period, $T_{m}=mT_{s}$ . The convergence criterion for the relative error of the Newton method is generally taken as $\Vert \boldsymbol{f}^{T_{m}}(\boldsymbol{x})-\boldsymbol{x}\Vert /\Vert \boldsymbol{x}\Vert <10^{-5}$ , where $\Vert \cdot \Vert$ is the $L_{2}$ norm.

The initial condition for the search was taken from snapshots of the LES symmetric turbulence in the above-mentioned marginal range in which turbulence eventually decays to laminar, $R_{S}<50$ , $A_{xz}=3$ and $A_{yz}=1.33$ . In order to find a solution, we fix the aspect ratios and the time period $T_{m}\approx 10$ ( $m=5$ ), then change $R_{S}$ . For $R_{S}\approx 35$ , the Newton search always converges to a trivial laminar state. The first non-trivial solution is obtained spontaneously at $R_{S}=38.6$ , and a trial with $m=7$ converges to the same solution. It is, therefore, a periodic orbit with the period of $T_{m}=T_{s}$ ( $m=1$ ) close to the bifurcation point that marks the lowest $R_{S}$ of our family of solutions. The solutions are always unstable and some examples of the linear stability analysis are shown in the later section and in the Appendix.

From this initial solution, lower and upper branches are continued along the parameters $A_{xz},~A_{yz}$ and $R_{S}$ , overcoming the turning points with a pseudo-arclength method. We have seen that all these solutions are periodic orbits in which the period is $T_{s}$ . The periodicity is especially noticeable for the upper-branch solutions in flat boxes, $A_{yz}=1.33$ –1.64 and $R_{S}\approx 45$ , where the temporal oscillation of the integrated kinetic energy is of order 1 %. They asymptotically approach a fixed point in the limit of $A_{yz}\rightarrow \infty$ , and the oscillations become negligible, $O(10^{-5})$ for $A_{yz}>2$ , for the vertically localised solutions described later in taller boxes.

4.1 The effect of box aspect ratios and characterisation of equilibria

Figure 5(a) shows the dependence on $A_{yz}$ of the Reynolds stress at the central plane for $R_{S}=38.9$ and $A_{xz}=3$ , close to the initial bifurcation. Solutions exist only for $A_{yz}>1.1$ . The Reynolds stress decreases as $A_{yz}$ increases, and may approach a non-zero constant in the limit $A_{yz}\rightarrow \infty$ . The same is true for the velocity fluctuations (not shown). Note that, since the grid becomes relatively coarser as $A_{yz}$ increases, a finer $y$ grid is used for taller boxes (dashed lines), and that the good agreement of the results whenever the two grids overlap confirms numerical convergence.

The dependence on the streamwise aspect ratio $A_{xz}$ is shown in figure 5(b). Solutions exist for $1.58\lesssim A_{xz}\lesssim 3.29$ , which covers the range of box aspect ratios ( $A_{xz}\approx 3$ ) identified by Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016) to be good models for wall-bounded shear turbulence. It is interesting that the minimum aspect ratio for steady Nagata equilibrium solutions in plane Couette flow is $A_{xz}\approx 1.62$ at low Reynolds numbers (Jiménez et al. Reference Jiménez, Kawahara, Simens, Nagata and Shiba2005), and that Deguchi (Reference Deguchi2015) showed that they exist only for $A_{xz}\gtrsim 1.5$ , even at high Reynolds numbers. On the other hand, Kawahara & Kida (Reference Kawahara and Kida2001) found unstable periodic orbits for somewhat shorter boxes, $A_{xz}=1.46$ .

Figure 6 shows the continuation along $R_{S}$ for $A_{xz}=3$ and several $A_{yz}$ . The box-averaged statistics of LES (symmetric) turbulence are also shown for comparison. For $R_{S}\lesssim 60$ , those turbulence simulations occasionally become vertically localised around $y=0$ , but they spread again to fill the whole domain (see § 4.3). The box-averaged turbulence statistics in figure 6 include such locally quiescent regions, leading to weaker $v^{\prime }$ and $w^{\prime }$ fluctuations than would otherwise be obtained at the central plane. When $R_{S}\lesssim 50$ , LES turbulence often decays to laminar after these localisation events.

The velocity fluctuations for the upper-branch solutions are quite large in flat boxes and at low Reynolds number, $A_{yz}=1.5$ and $R_{S}\approx 45$ , but that behaviour disappears for taller boxes and saturates beyond $A_{yz}\approx 3$ . These large fluctuations are thus probably an effect of the shear-periodic boundary condition. This is also the range in which the shear periodicity results in the strongest temporal oscillations of the kinetic energy, of the order of 1 %, but we tested that the fluctuations in figure 6(a,c) are not due to the temporal variability. They are also present in the spacial average of instantaneous snapshots. The temporal oscillation of solutions with $A_{yz}>2$ is of the order of $10^{-5}$ .

The relatively poor scaling of the velocity fluctuations of the LES equilibria with $SL_{z}$ can be traced to the poor scaling of $|\unicode[STIX]{x1D70E}|/S$ . Figure 6(d) shows that the dimensionless strain rate of the upper-branch equilibria is roughly unity for $R_{S}=50$ –100, while that of turbulence is in the range of 2–5. It was shown by Sekimoto et al. (Reference Sekimoto, Dong and Jiménez2016) that, although $SL_{z}$ is the natural velocity scale for ‘good’ DNS boxes, the fluctuations in non-optimal boxes scale better with the friction velocity obtained from their total measured stress. In essence, the Reynolds stress and the velocity fluctuations scale with each other.

The same is true for LES and for equilibrium solutions. Figure 7(a,c,d) shows that the fluctuations collapse to a common curve in these ‘wall’ units. The vorticity components also scale well with $u_{\unicode[STIX]{x1D70F}}/l_{S}$ (not shown). In this normalisation, the velocity fluctuations of symmetric LES turbulence at $R_{S}=100$ are $u^{\prime +}\approx 2$ , $v^{\prime +}\approx 1.2$ and $w^{\prime +}\approx 1.4$ , which agree well with those at the top of the logarithmic layer of turbulent channels. When $R_{S}<60$ , the velocity fluctuations of the LES equilibria are not very different from those of turbulence, but we have seen that the lower-branch solutions tend to get concentrated around the critical layer as $R_{S}$ increases. Their statistics then become very different from turbulence.

Figure 7(b) is an enlargement of figure 7(a), showing the dependence on $A_{xz}$ of the minimum bifurcation Reynolds number of equilibria with $A_{yz}=3$ . It turns out that the bifurcation is more dependent on $A_{xz}$ than on $A_{yz}$ . At a fixed $A_{xz}=3$ , the minimum $R_{S}$ is $(38.11,38.00,37.94,37.91,37.89,37.90)$ for $A_{yz}=(1.33,1.5,1.64,1.8,2.0,3.0)$ , with a maximum scatter of 0.6 %. On the other hand, when $A_{yz}$ is fixed at 3, $R_{S}$ changes by 8 %, from 35.05 to 37.90, as $A_{xz}$ changes from 2 to 3.

Figure 7.

(a,c,d) As in figure 6, but scaled by $u_{\unicode[STIX]{x1D70F}}$ at $y=0$ : (a,b)  $u_{c}^{\prime +}$ , (c)  $v_{c}^{\prime +}$ , (d)  $w_{c}^{\prime +}$ . (b) Detail of (a): solid with $\triangle$ , red, $A_{xz}=2$ ; solid with $\Box$ , blue, $A_{xz}=2.5$ ; black with $\circ$ ,  $A_{xz}=3$ .

Figure 8.

Velocity fluctuations for $A_{xz}=A_{yz}=3$ : (a)  $R_{S}=50$ , (b)  $R_{S}=62.5$ ; only half of each plot is shown, using the symmetry in $y$ . Left side is lower branch. Right side is upper branch. Curves: - - - - (black), $u^{\prime }$ ; —— (blue), $v^{\prime }$ ; – - – - – (red), $w^{\prime }$ .

Figure 9.

(a,b) As in figure 8, for the momentum balance: —— (blue), $\langle -uv\rangle _{xz}/u_{\unicode[STIX]{x1D70F}}^{2}$ ; - - - - (red), $\langle 2\unicode[STIX]{x1D708}_{t}\overline{\unicode[STIX]{x1D70E}}_{xy}\rangle _{xz}/u_{\unicode[STIX]{x1D70F}}^{2}$ ; —— (black), total stress.

Figure 10.

Height of the velocity streak, defined by the distance $d_{s}$ between the peaks of $u^{\prime }$ in the core part of the solutions in figure 8. Only peaks in the inner core are considered, and the secondary peaks in the outer one-component layer are ignored. (a) Plot for $A_{xz}=3$ and $A_{yz}=1.33,1.5,1.64,1.8,2,3$ . Lines and symbols as in figure 6. (b) Plot for $A_{xz}=2,2.5,3$ and $A_{yz}=3$ . Lines and symbols as in figure 7(b).

Figure 11.

Vertically localised upper-branch solutions. Plots as in figure 3, for $A_{xz}=3,~A_{yz}=3$ and (a)  $R_{S}=42.7$ and (b)  $R_{S}=62.5$ .

Figure 12.

Isosurfaces of the second invariant of the velocity gradient ( $Q=0.2Q_{max}$ ) coloured by $\unicode[STIX]{x1D714}_{x}$ , with the same colour code as in the streamwise-averaged $\unicode[STIX]{x1D714}_{x}$ maps in the left panels of figure 11(a,b): $\unicode[STIX]{x1D714}_{x}<0$ (dark grey), $\unicode[STIX]{x1D714}_{x}>0$ (light grey). From top to bottom, each panel shows: a three-dimensional view, a side view of the region $0.5\leqslant z/L_{z}\leqslant 1$ , and a cross-section of cross-flow velocity vectors at $x/L_{z}=2.25$ , marked as a chain-dotted vertical line in the side views. The primary streamwise rollers visible in these cross-sections are not always compact enough to appear in the $Q$ isosurfaces and are marked as A in the lateral views. Flow as in figure 11(a,b). Only the active part, $|y|/L_{z}<0.5$ , is shown in all cases.

4.2 The structure of the upper-branch equilibria

We focus next on the flow structure of the vertically localised upper-branch equilibria in a box with $A_{xz}=3$ and $A_{yz}=3$ . The right-hand side of the two panels in figure 8 shows that the velocity fluctuations of these solutions decay exponentially away from $y=0$ , which is a common feature of localised solutions in other shear flows (Schneider, Gibson & Burke Reference Schneider, Gibson and Burke2010; Gibson & Brand Reference Gibson and Brand2014). The tail of $u^{\prime }$ is much stronger than those of $v^{\prime }$ and $w^{\prime }$ . This can be shown to be due to a streamwise-constant ‘streak’ of the streamwise velocity that can be approximated as a roughly sinusoidal spanwise perturbation of the mean profile, $u_{s}\sim \sin (2\unicode[STIX]{x03C0}z/L_{z})$ . This weak perturbation is the far tail of the stronger sinuous streak of $u$ concentrated near $y=0$ and, because it is almost independent of $x$ , is essentially independent of $v$ , $w$ and the shear. The left-hand part of both panels in figure 8 represents lower-branch solutions, which share many of the characteristics of the upper branch. As the Reynolds number increases, all structures become more complex, as seen in figure 8(b) and later in figures 11 and 12, and the core of the structures develop substructures that could perhaps be interpreted as a first indication of a turbulent cascade. The momentum balance of these flows is shown in figure 9. The Reynolds stress is dominant in the core, $|y|/L_{z}<0.5$ , but only the mean shear $\text{d}U/\text{d}y$ contributes to the eddy viscosity in the outer part.

Figure 10 presents the height of the velocity streak for the different equilibria, defined as twice the distance to $y=0$ of the first maximum of the $u^{\prime }$ profile of the solutions in figure 8. The height of the lower-branch streak stays roughly constant below $R_{S}\approx 50$ , and approaches zero quickly above that limit, as the solutions collapse to the critical layer. In contrast, the height of upper-branch solutions increases drastically near the bifurcation point and reaches a maximum near $R_{S}\approx 45$ . This is also the range in which the velocity fluctuations become very strong in figure 6, and only appears in flat boxes with $A_{yz}=1.33{-}1.64$ . We argued in the discussion of figure 6 that upper-branch solutions tend to be limited by the height of these flat boxes, and figure 10(a) confirms this interpretation by showing that the maximum height of the streak reaches the box height. On the other hand, this interaction does not take place when $A_{yz}\gtrsim 2$ , confirming the independence of those solutions from the box dimensions.

The flow structures of the upper-branch solutions at $R_{S}=42.7$ and 62.5 are visualised in figures 11 and 12. The streak is represented by the $\overline{u}=0$ isosurface that spans $|y|/L_{z}<0.5$ in figure 11 (see also figure 3 a,b). The corresponding vortical structures are located in the flanks of the streak, as seen in the streamwise-averaged cross-stream maps of $\unicode[STIX]{x1D714}_{x}$ in figure 11. The vortical structures are shown by themselves in figure 12, and are surprisingly complex for an equilibrium solution. This is especially true for the higher-Reynolds-number case in figures 11(b) and 12(b), which appear to include a double structure that is nevertheless in equilibrium. In fact, there is a first indication of two separate scales in these flow fields. The smaller tube-like vortical structures are isosurfaces of the second invariant of the velocity gradient tensor ( $Q$ criterion), and are coloured by the streamwise vorticity. They do not always coincide with the larger-scale structures of the cross-stream velocity, which are visible in the cross-sections in the bottom part of figure 12. These ‘rollers’ are too diffuse to appear in the $Q$ -map, which may actually appear empty in the region in which the roller is dominant (see the region labelled as A in the side view in figure 12 b, and note that, because of the symmetry of the flow, a mirror-symmetric arrangement is present in the upstream half of the box). On the other hand, it is clear from the cross-sections that the roller dominates the velocity field.

The low-Reynolds-number solution in figure 12(a) shows that the $Q$ -vortices are parts of the larger streamwise rollers that have been sheared and stretched by the mean flow. The cross-sections in the bottom part of figure 12 show that these vortices are long enough to spill into neighbouring periodic boxes, so that they appear as double in the cross-section. The inclination angle of these streamwise rollers and vortices is roughly 15 $^{\circ }$ at all Reynolds numbers. In the higher-Reynolds-number case in figure 12(b), the streamwise vortices are strong enough to create new vortices, roughly perpendicular to them, which are labelled as B in the side view in figure 12(b). They rotate in the opposite sense to the primary streamwise rollers and are aligned in the direction of the strain produced by them. A similar generation mechanism of secondary smaller vortices has been investigated in the homogeneous isotropic turbulence (Goto Reference Goto2012).

The velocity spectra in figure 13 show that the upper-branch solutions acquires more small-scale structures as $R_{S}$ increases, approaching the spectrum of turbulence LES at similar Reynolds numbers. Even though the turbulence state has more small scale, the large-scale end of its spectrum is quantitatively similar to the upper-branch solutions.

Figure 13.

Streamwise velocity spectrum $E_{uu}^{\ast }=E_{uu}(k_{x})\unicode[STIX]{x1D700}^{-2/3}l_{S}^{-5/3}$ as a function of $k_{x}l_{S}$ : $\triangledown$  (blue), upper-branch solution at $R_{S}=42.7$ ; $\circ$  (black), upper-branch solution at $R_{S}=62.5$ ; $\triangle$  (red), turbulence LES at $R_{S}=62.5$ . The dash-dotted line represents the inertial theory $E_{uu}^{\ast }(k_{x})=0.6(l_{S}k_{x})^{-5/3}$ . The spectra and dissipation rates of the upper-branch solutions are averaged over $|y|/L_{z}<0.5$ .

Figure 14.

(a) Vertical velocity fluctuation intensity profiles, $\langle v^{2}\rangle _{xz}^{1/2}/(SL_{z})$ , as a function of $y/L_{z}$ and $St$ . The contours are $[0.02:0.02:0.08]$ , and the mean value of $v^{\prime }$ is approximately $0.07SL_{z}$ (see figure 6). Here $A_{xz}=3$ , $A_{yz}=3$ and $R_{S}=50$ . (b) Homogeneous turbulence at $St=1318$ . (c) Localised state at $St=1525$ . The left panels in (b) and (c) show the streamwise-averaged $\unicode[STIX]{x1D714}_{x}$ , and the right ones show isosurfaces of $\unicode[STIX]{x1D714}_{x}$ and $\overline{u}$ , as in figure 3.

4.3 Intermittent visiting of turbulence to vertically localised states

We mentioned in § 4.1 that, when $R_{S}$ is relatively low, LES turbulence collapses intermittently to a localised state around $y=0$ , and that these states persist for a long time when $R_{S}\gtrsim 50$ . These are the statistics plotted as diamonds in figure 6. Figure 14(a) shows the temporal evolution of the profile of local $v^{\prime }$ . Light colours represent active turbulent regions and dark ones are overdamped ‘laminar’ areas. Localised turbulence occurs when laminarisation does not extend over the whole height of the box, such as in $St=1500{-}2000$ .

Note that localisation only takes place in flows with sinuous symmetry, so that the profiles are symmetric with respect to $y=0$ , but that the symmetry allows localisation either at the centre plane, $y=0$ , or at the top or bottom of the box, $y=\pm L_{y}/2$ . The strongest event in figure 14(a) is localised at $y=0$ $(St=1500{-}2000)$ , but there are several minor ones (e.g. at $St\approx 400$ ) localised at the top of the box.

The structure of the flow at a fully turbulent moment is shown in figure 14(b), and shows that vorticity is spread over the whole box. Figure 14(c) represents the flow during a localisation event. Its similarity to the localised lower-branch equilibrium solutions is striking. Since we saw in § 3 that these solutions are edge states whose unstable manifold leads to bursting, it is tempting to interpret localisation as the occasional approach of fully developed turbulence to the localised solution along its stable manifold, and its later spreading to a fully turbulent state along the unstable manifold. A similar behaviour was reported in a plane channel by Itano & Toh (Reference Itano and Toh2001), and Kawahara (Reference Kawahara2005) used the occasional approaches to the edge state to laminarise turbulent Couette flow.

Figure 15.

(a) The local cross-flow velocity fluctuation intensity $u_{\bot }$ defined in (4.1) at: (black) $y=L_{y}/2$ ( $u_{t}$ ); (grey) $y/L_{z}=0$ ( $u_{c}$ ). The black dashed line is the mean value, $\langle u_{\bot }\rangle$ , and the red dashed lines are $\langle u_{\bot }\rangle \pm \unicode[STIX]{x1D70E}$ , where $\unicode[STIX]{x1D70E}$ is the standard deviation. They are defined as averages over the centre and top of the box, as explained in the text. (b) The intermittency factor $(1-\unicode[STIX]{x1D6FE})$ , where $\unicode[STIX]{x1D6FE}$ is defined in (4.2): — $\circ$ —,  $R_{S}=50.0$ ; — $\bullet$ —,  $R_{S}=52.6$ ; — $\triangledown$ — (grey),  $R_{S}=55.5$ ; – –  ${\vartriangle}$ – –,  $R_{S}=62.5$ ; — $\Box$ — (grey),  $R_{S}=101.6$ .

We next quantify the frequency of these localisation events. Figure 15(a) shows the evolution of the transverse velocity fluctuation intensity, defined as

The darker line is the intensity at $y=0$ , and the lighter one is the intensity at $y=\pm L_{y}/2$ . The thin straight solid and the dashed lines are the mean, $\langle u_{\bot }\rangle$ , and deviation, $\unicode[STIX]{x1D70E}$ , which are the averages of these two positions. By defining $u_{c}=u_{\bot }(0)$ and $u_{t}=u_{\bot }(L_{y}/2)$ , fully turbulent profiles should have large $u_{c}\approx u_{t}\approx \langle u_{\bot }\rangle$ , while localised states are characterised by having one of these intensities much weaker than the fully turbulent one. The quantity

where $P(u_{c},u_{t})$ is the joint probability density function of both intensities, and $t_{u}$ is a threshold, measures the probability that the points $y=0$ or $L_{y}/2$ are turbulence according to $u_{c}\geqslant t_{u}$ or $u_{t}\geqslant t_{u}$ . Since the localised state has the overdamped laminar region, which appears intermittently as shown in figure 15(a), $(1-\unicode[STIX]{x1D6FE})$ indicates the intermittency factor of the localised turbulence. Figure 15(b) shows that the frequency of localisation decreases as $R_{S}$ increases.