1 Introduction
The discovery of very-large-scale motions (VLSMs) has attracted significant interest in research into wall-bounded turbulence over the past decade (e.g. Hutchins & Marusic Reference Hutchins and Marusic2007). The VLSM features as a long streaky motion of streamwise turbulent kinetic energy in the outer region, and it is typically very energetic at sufficiently high Reynolds numbers (
$Re_{\unicode[STIX]{x1D70F}}\gtrsim O(10^{3})$
where
$Re_{\unicode[STIX]{x1D70F}}$
is the friction Reynolds number). It was initially proposed that this long streaky structure may be formed by the concatenation of the large-scale vortical structures, known as the large-scale motions (LSMs) (Kovasznay, Kibens & Blackwelder Reference Kovasznay, Kibens and Blackwelder1970), which themselves were speculated to be formed by merger and/or growth of near-wall hairpin vortices via a ‘bottom-up’ process (for further details, the reader may refer to a recent summary on this proposition by Adrian (Reference Adrian2007)). However, there has been a growing body of recent evidence that the outer structures are largely independent of the near-wall process: for instance, disruption of the near-wall process with wall roughness affects the outer statistics very little (Flores, Jiménez & del Álamo Reference Flores, Jiménez and del Álamo2007).
We have recently shown that the coherent structures in the outer region sustain themselves, even in the absence of the motions in the near-wall and logarithmic regions (Hwang & Cossu Reference Hwang and Cossu2010b ; Rawat et al. Reference Rawat, Cossu, Hwang and Rincon2015). The self-sustaining outer structure is composed of two structural elements – a long streak and short quasistreamwise vortices, which respectively correspond to the VLSM and the LSMs (Hwang Reference Hwang2015). The self-sustaining process was also found to be almost identical to that in the near-wall region: (1) the streaky structure (VLSM) is amplified by the vortical structures (LSMs) via the lift-up effect (e.g. Cossu, Pujals & Depardon Reference Cossu, Pujals and Depardon2009; Pujals et al. Reference Pujals, García-Villalba, Cossu and Depardon2009; Hwang & Cossu Reference Hwang and Cossu2010a ; Willis, Hwang & Cossu Reference Willis, Hwang and Cossu2010); (2) the amplified streak undergoes a rapid meandering motion via secondary instability or transient growth (Park, Hwang & Cossu Reference Park, Hwang and Cossu2011); (3) the subsequent nonlinear regeneration of the streamwise vortical structures (Hwang & Bengana Reference Hwang and Bengana2016).
The existence of a self-sustaining process at large scale in the outer region is of particular theoretical importance, because it indicates that the outer structures are probably organised around invariant solutions of the system, often referred to as ‘exact coherent structures’ (e.g. Nagata Reference Nagata1990; Waleffe Reference Waleffe2001; Faisst & Eckhardt Reference Faisst and Eckhardt2003; Wedin & Kerswell Reference Wedin and Kerswell2004; Hall & Sherwin Reference Hall and Sherwin2010; Park & Graham Reference Park and Graham2016, and many others). The simplest non-trivial exact solutions of the Navier–Stokes equations are typically in the form of a stationary or travelling wave, being equilibria or relative equilibria in phase space. Together with unstable periodic and/or relative periodic orbits (e.g. Kawahara & Kida Reference Kawahara and Kida2001), these invariant solutions have been shown to form a skeleton for solution trajectories in phase space (Gibson, Halcrow & Cvitanovic Reference Gibson, Halcrow and Cvitanovic2008; Willis, Cvitanovic & Avila Reference Willis, Cvitanovic and Avila2013, Reference Willis, Cvitanovic and Avila2016), and their understanding from a dynamical systems viewpoint has been at the heart of recent advancement in the understanding of bypass transition and low-Reynolds-number turbulence.
The goal of the present study is to demonstrate that such invariant solutions are also the driving mathematical mechanism of the large-scale outer structures given with the VLSMs and the LSMs in high-Reynolds-number turbulent channel flow. This task, however, has often been understood to be challenging. The principal difficulty lies in the emergence of a huge number of invariant solutions, which significantly hampers identification of which solutions are most relevant to given coherent structures of interest. Furthermore, the computation of the invariant solutions at such a high Reynolds number is often numerically very sensitive and expensive, yielding a substantial technical barrier. To bypass these difficulties, Rawat et al. (Reference Rawat, Cossu, Hwang and Rincon2015) and Rawat, Cossu & Rincon (Reference Rawat, Cossu and Rincon2016) recently computed a set of coherent invariant solutions at Reynolds numbers up to
$Re_{\unicode[STIX]{x1D70F}}=128$
by modelling all the smaller-scale structures around the structure of interest via large-eddy simulations (LESs) with the Smagorinsky model. In these studies, the Smagorinsky constant
$C_{s}$
was taken as a continuation parameter to replace the surrounding unsteady smaller-scale motions with an elevated eddy viscosity, as in Hwang & Cossu (Reference Hwang and Cossu2010b
, Reference Hwang and Cossu2011). Here, we extend this approach to a different regime of much higher Reynolds numbers up to
$Re_{\unicode[STIX]{x1D70F}}\simeq 1000$
, at which the VLSMs and the LSMs emerge energetically in the flow domain, and show that the invariant solutions are directly linked with the formation of the large-scale structures.
2 Numerical method
We consider a turbulent channel, with half-height
$h$
, in which
$x$
,
$y$
and
$z$
denote the streamwise, wall-normal and spanwise direction respectively. The two walls are located at
$y=0$
and
$y=2h$
. We use a Navier–Stokes solver that is well documented in Bewley (Reference Bewley2014). In this solver, the streamwise and spanwise directions are discretised using Fourier series with
$2/3$
dealiasing rule, whereas the wall-normal direction is discretised using second-order central difference. A set of LESs are considered with the static Smagorinsky model, as in previous studies (e.g. Hwang & Cossu Reference Hwang and Cossu2010b
, Reference Hwang and Cossu2011; Rawat et al.
Reference Rawat, Cossu, Hwang and Rincon2015, Reference Rawat, Cossu and Rincon2016), i.e.
$\tilde{\unicode[STIX]{x1D70F}}_{ij}-\unicode[STIX]{x1D6FF}_{ij}/3\tilde{\unicode[STIX]{x1D70F}}_{kk}=-2\unicode[STIX]{x1D708}_{t}\tilde{\unicode[STIX]{x1D61A}}_{ij}$
with
$\unicode[STIX]{x1D708}_{t}=(C_{s}\tilde{\unicode[STIX]{x1D6E5}})^{2}\tilde{\unicode[STIX]{x1D64E}}\mathscr{D}$
, where
$\tilde{\cdot }$
denotes the filtered quantity,
$\unicode[STIX]{x1D61A}_{ij}$
the strain rate tensor,
$C_{s}$
the Smagorinsky constant,
$\tilde{\unicode[STIX]{x1D6E5}}=(\tilde{\unicode[STIX]{x1D6E5}}_{1}\tilde{\unicode[STIX]{x1D6E5}}_{2}\tilde{\unicode[STIX]{x1D6E5}}_{3})^{1/3}$
the nominal filter width,
$\tilde{\unicode[STIX]{x1D64E}}=(2\tilde{\unicode[STIX]{x1D61A}}_{ij}\tilde{\unicode[STIX]{x1D61A}}_{ij})^{1/2}$
the norm of the strain rate tensor, and
$\mathscr{D}=1-\exp [-(y^{+}/A^{+})^{3}]$
with
$A^{+}=25$
is the van Driest damping function. The Smagorinsky constant for turbulent channel flow is typically set in the range between
$C_{s}=0.05$
and
$C_{s}=0.10$
(e.g. Moin & Kim Reference Moin and Kim1982; Härtel & Kleiser Reference Härtel and Kleiser1998). In particular,
$C_{s}=0.05$
has been shown to provide the best performance in terms of accurate generation of first- and second-order turbulent velocity statistics (a posteriori test), while
$C_{s}=0.10$
is known for the model to provide the best turbulent dissipation in comparison to the true value obtained in a direct numerical simulation (DNS) (a priori test) (Härtel & Kleiser Reference Härtel and Kleiser1998; Meneveau & Katz Reference Meneveau and Katz2000). Mason & Cullen (Reference Mason and Cullen1986) showed that the Smagorinsky constant
$C_{s}$
actually acts as the filter width of the LES. Therefore, artificially increasing
$C_{s}$
allows one to damp the small-scale motions without losing actual resolution of the large-scale structures, as also demonstrated in Hwang & Cossu (Reference Hwang and Cossu2010b
). It is also important to note that the static Smagorinsky model prevents any energy transfer from the modelled residual stress to the resolved motions, ensuring that the resolved motions sustain themselves at the increased
$C_{s}$
. All the simulations in this study are performed by imposing constant volume flux across the channel.
In the present study, computation of the invariant solutions is restricted to the minimal unit for the self-sustaining process of the outer structures in Hwang & Cossu (Reference Hwang and Cossu2010b
), i.e.
$L_{x}=3.0h$
and
$L_{z}=1.5h$
where
$L_{x}$
and
$L_{z}$
are the streamwise and spanwise domain sizes, respectively. An important benefit of the minimal unit is that it realises a low-dimensional dynamical set of coherent structures in the form of a VLSM (long outer streak) and LSMs (short outer vortices) without significant distortion of turbulent statistics (Hwang & Cossu Reference Hwang and Cossu2010b
; Rawat Reference Rawat2014; Hwang & Bengana Reference Hwang and Bengana2016). Given the symmetry of turbulent statistics around the channel midplane, we will also focus on seeking invariant solutions with mirror symmetry about
$y=h$
. For this purpose, only the bottom half of the channel is solved by imposing the symmetric boundary condition at the channel midplane (i.e.
$\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y=0$
,
$v=0$
and
$\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}y=0$
at
$y=h$
where
$u$
,
$v$
and
$w$
are the streamwise, wall-normal and spanwise velocities, respectively). Except for this setting, all the simulation parameters, including the number of grid points, at two Reynolds numbers considered,
$Re_{m}=20\,133$
and
$Re_{m}=38\,133$
, are the same as those in Hwang & Cossu (Reference Hwang and Cossu2010b
, Reference Hwang and Cossu2011), as summarised in table 1. The artificially increased
$C_{s}$
values, by which the small-scale structures are replaced with the eddy viscosity, but not the self-sustaining outer motions, are therefore also the same as those in these works (see the parameters of simulations
$S550$
and
$S950$
in table 1). Finally, it should be noted that the two sets of grid points, respectively for
$Re_{m}=20\,133$
and
$Re_{m}=38\,133$
, in the present study are chosen such that the standard LESs of the present study ensure good resolution for the small-scale near-wall motions (Zang Reference Zang1991, see also figure 1), as in our previous studies (Hwang & Cossu Reference Hwang and Cossu2010b
, Reference Hwang and Cossu2011). An examination of the energy spectra reveals that these numbers of grid points are also found to provide good spatial resolution for the computed invariant solutions, at least for the elevated
$C_{s}$
. We also note that the resolution of the invariant solution for
$Re_{m}=38\,133$
is finer than that in Waleffe (Reference Waleffe2001).
Simulation parameters in the present study (before dealiasing). Here,
$Re_{m}=2U_{m}h/\unicode[STIX]{x1D708}$
where
$U_{m}$
is the bulk velocity. Note that
$U_{m}=(2/3)U_{l}$
where
$U_{l}$
is the centreline velocity of the corresponding laminar flow with the same volume flux. The simulations tagged with ‘
$F$
’ indicate full simulations resolving near-wall motions, while those tagged with ‘
$S$
’ are simulations with only self-sustaining outer motions made by increasing the Smagorinsky constant
$C_{s}$
.
A reference simulation is first performed to check turbulence statistics of the half-channel LES simulation at
$Re_{\unicode[STIX]{x1D70F}}\simeq 950$
, with the value
$C_{s}=0.05$
(
$F950$
) known to provide the best statistical fit to the full DNS result (Hwang & Cossu Reference Hwang and Cossu2010b
). In figure 1, its first- and second-order turbulence statistics are compared with those of full-channel LES with the same streamwise and spanwise computational domain (Hwang & Cossu Reference Hwang and Cossu2011) as well as those of DNS by del Álamo et al. (Reference del Álamo, Jiménez, Zandonade and Moser2004) at
$Re_{\unicode[STIX]{x1D70F}}=934$
. Overall, the half-channel simulation generates fairly good turbulence statistics compared with those from the full-channel LES, which itself shows reasonable agreement with the data from full DNS. The only appreciable difference between the half-channel and full-channel simulations appears in the wall-normal velocity fluctuation very near the channel centre, due to the symmetry condition. This indicates that the half-channel simulation does not lose important physical features, except around the very centre of the channel, where mainly dissipation of structures is expected due to the small shear in this region.
(a) Mean velocity and (b) turbulent velocity fluctuations: ——, half-channel simulation with
$C_{s}=0.05$
(
$F950$
); - - - -, full-channel simulation with
$C_{s}=0.05$
(Hwang & Cossu Reference Hwang and Cossu2011); – ⋅ – ⋅ –, DNS at
$Re_{\unicode[STIX]{x1D70F}}=934$
by del Álamo et al. (Reference del Álamo, Jiménez, Zandonade and Moser2004).
3 The invariant structures
3.1 Computation of invariant solutions
Now, we increase
$C_{s}$
such that the simulation contains only the large-scale self-sustaining structures in the given computational domain (simulations
$S550$
and
$S950$
in table 1). By doing so, all the structures smaller than the large-scale structures are removed, while their roles are modelled with the artificially elevated eddy viscosity. It is very important to note that the removal with an appropriate increase in
$C_{s}$
does not significantly distort the statistics and the self-sustaining dynamics of the large-scale structures themselves, as extensively discussed in the previous studies (Hwang & Cossu Reference Hwang and Cossu2010b
; Hwang Reference Hwang2015; Rawat et al.
Reference Rawat, Cossu, Hwang and Rincon2015; Hwang & Bengana Reference Hwang and Bengana2016). Therefore, the computed invariant solutions at the artificially elevated
$C_{s}$
would conceptually represent the structures in the ‘presence of the surrounding small scales’, which is modelled by the eddy viscosity, while enabling us to compute the invariant solutions with the relatively low resolutions at the high Reynolds numbers considered. The invariant solutions of the system with the elevated
$C_{s}$
are sought in the subspace satisfying the so-called shift-reflect symmetry
$$\begin{eqnarray}[u,v,w,p](x,y,z)=[u,v,-w,p](x-L_{x}/2,y,-z),\end{eqnarray}$$
where
$p$
is the pressure, together with the mirror symmetry about
$y=0$
. It should be mentioned that this specific symmetry is intentionally posed to find the invariant solutions containing the ‘sinuous’ mode of streak instability (i.e. streak meandering motions along the streamwise direction). The sinuous mode of streak instability has been consistently found to be the dominant mechanism of the streak breakdown in our previous theoretical analysis (Park et al.
Reference Park, Hwang and Cossu2011) as well as the minimal channel simulation for the outer structures (Hwang & Bengana Reference Hwang and Bengana2016). Indeed, in this study, we have also found that imposing the symmetry (3.1) in the present half-channel simulation with the increased
$C_{s}$
(
$S550$
and
$S950$
) does not engender any significant difference from the simulation without this symmetry (see figure 4). It is worth mentioning that Waleffe (Reference Waleffe2001) also imposed this symmetry for computation of his invariant solutions for the same reason.
To compute the invariant solutions, we have implemented a Newton–Krylov–Hookstep method and applied it to the present LES solver. Details of the method are given in Willis et al. (Reference Willis, Cvitanovic and Avila2013), and it is similar to that proposed by Viswanath (Reference Viswanath2007). This method computes an invariant solution by minimising the relative error between an initial guess for an initial flow field and the same field time-stepped by an interval
$T$
and shifted in the streamwise direction by a distance
$-l_{x}$
. For an equilibrium state the choice of
$T$
is arbitrary, and the phase speed of the equilibrium is
$c=l_{x}/T$
. Throughout this study, the computation of the invariant solutions is carried out with
$T=16.7h/U_{m}$
. All the solutions are computed to a relative-error tolerance of
$10^{-7}{-}10^{-8}$
between the initial and shifted end states. Since some of the invariant solutions are expected to sit on the so-called ‘edge’ state, which refers to the phase-space boundary manifold between the basic and the chaotic states (e.g. Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006), we start by computing the edge state for the
$S550$
simulation in the given subspace, using the standard bisection technique to obtain a good initial guess for the Newton iteration (e.g. Duguet, Willis & Kerswell Reference Duguet, Willis and Kerswell2008; Avila et al.
Reference Avila, Mellibovsky, Roland and Hof2013). Several instantaneous flow fields on the edge state are given for initial guess of the Newton solver, and we found an invariant solution propagating downstream with a constant speed (i.e. a travelling wave) from
$S550$
. Numerical continuation in
$C_{s}$
is subsequently performed, as in Rawat (Reference Rawat2014) and Rawat et al. (Reference Rawat, Cossu, Hwang and Rincon2015).
Figure 2(a) shows the bifurcation of the invariant solutions with
$C_{s}$
by plotting
$Re_{\unicode[STIX]{x1D70F}}$
of each of the solutions. Here, we note that
$Re_{\unicode[STIX]{x1D70F}}$
in this case is given as a measure of the friction of each solution, revealing its relevance to high-Reynolds-number flows, i.e.
$Re_{\unicode[STIX]{x1D70F}}=2u_{\unicode[STIX]{x1D70F}}/U_{m}Re_{m}$
. Therefore,
$Re_{\unicode[STIX]{x1D70F}}$
given here is different for each invariant solution and represents its friction. Continuation reveals that the invariant solutions experience a saddle-node bifurcation as
$C_{s}$
is gradually lowered from a large value. Two invariant solutions are found to emerge at the critical
$C_{s}$
(
$\simeq 0.335$
), as often observed in transitional Reynolds numbers: one has a low drag (lower-branch solution) and the other has a high drag (upper-branch solution). The solutions obtained at
$Re_{m}=20\,133$
are further continued to a higher Reynolds number,
$Re_{m}=38\,133$
. Essentially, the same behaviour with
$C_{s}$
is obtained at this Reynolds number, as shown in figure 2(b).
Bifurcation of invariant solutions with
$C_{s}$
: (a)
$Re_{m}=20\,133$
; (b)
$Re_{m}=38\,133$
. Here, the simulations
$S550$
and
$S950$
respectively correspond to
$C_{s}=0.2$
in (a) and
$C_{s}=0.3$
in (b).
Visualisation of the invariant solutions: (a)
$U950$
; (b)
$L950$
. In both (a) and (b), the blue isosurfaces indicate
$u^{+}=-2.96$
and the green ones
$\unicode[STIX]{x1D708}_{t}=7.2\unicode[STIX]{x1D708}$
. The red isosurfaces in (a) and (b) are
$Q^{+}=3.6\times 10^{-5}$
and
$Q^{+}=2.0\times 10^{-5}$
, respectively.
The wall-normal profile of (a) mean velocity, and (b) streamwise, (c) wall-normal and (d) spanwise velocity fluctuations: ——,
$S950$
; - - - -, the symmetry-constrained
$S950$
; –
$\cdot \cdot$
–
$\cdot \cdot$
–,
$U950$
; – ⋅ – ⋅ –,
$L950$
.
3.2 Spatial structure of the invariant solutions
The computed invariant solutions are visualised in figure 3. Both of the upper- and lower-branch solutions are characterised by a ‘wavy’ streak (blue isosurfaces) and streamwise vortices on the flank (red isosurfaces), clearly reflecting their tight physical link to the self-sustaining process of the outer coherent structures, i.e. streak generation via the lift-up effect with the vortices, and sinuous-mode streak instability with nonlinear feeding of the vortices. The upper-branch solution exhibits a strongly wavy streak and intense streamwise vortices, while the lower-branch solution is composed of a relatively straight streak and weak streamwise vortices (see the levels of the red isosurfaces in figure 3). This feature is consistent with that of the invariant solutions in, e.g., Waleffe (Reference Waleffe2001). The invariant solutions here, however, are obtained by modelling the surrounding small-scale structures with an eddy viscosity. The eddy viscosity is typically found to be quite strong around the streak, where high local shear is expected (the green isosurfaces in figure 3), as also found by Rawat et al. (Reference Rawat, Cossu, Hwang and Rincon2015, Reference Rawat, Cossu and Rincon2016). However, here the related turbulent dissipation is found to dominate over the kinematic viscosity because of the high Reynolds numbers considered – the maximum eddy viscosities of the upper- and the lower-branch solutions are respectively found as
$\unicode[STIX]{x1D708}_{t}=8.48\unicode[STIX]{x1D708}$
and
$\unicode[STIX]{x1D708}_{t}=7.6\unicode[STIX]{x1D708}$
, the values being an order of magnitude larger than those in Rawat et al. (Reference Rawat, Cossu, Hwang and Rincon2015, Reference Rawat, Cossu and Rincon2016). This follows from the higher values of
$C_{s}$
and, in particular, the much higher Reynolds numbers considered here.
Table 2 summarises the propagating speed of the computed invariant solutions with their friction velocity and centreline velocity at
$Re_{m}=20\,133$
and
$Re_{m}=38\,133$
. The propagating speeds of the upper- and the lower-branch solutions are respectively found to be
$c\simeq 0.66U_{c}$
and
$c\simeq 0.62U_{c}$
, and they scale more closely with their centreline velocity
$U_{c}$
than with their friction velocity
$u_{\unicode[STIX]{x1D70F}}$
. Although the centreline velocity
$U_{c}$
of each of the invariant solutions is not the same as that from full turbulent statistics given in figure 1, this finding supports the observation by del Álamo & Jiménez (Reference del Álamo and Jiménez2009) and Song et al. (Reference Song, Barkley, Hof and Avila2016), who showed that the advection velocity of the outer structures scales with the centreline velocity. This behaviour is quite intriguing, because the single turnover time scale of self-sustaining process of outer coherent structures has been found to scale well with the friction velocity
$u_{\unicode[STIX]{x1D70F}}$
(Hwang & Bengana Reference Hwang and Bengana2016). It is currently too early to draw any firm conclusion, but this indicates that the velocity scale of propagation of an outer structure might not be the same as the turnover time scale of the self-sustaining process.
Scaling of the speed of the invariant solutions. Here,
$c$
and
$U_{c}$
are the propagating speed and the centreline velocity of each travelling wave solution, respectively. In the first column,
$U$
and
$L$
respectively indicate the upper- and lower-branch solutions.
The first- and second-order statistics of the two computed invariant solutions are compared with those of simulation
$S950$
in figure 4. Here, the data of
$S950$
with the shift-reflect symmetry (3.1) are plotted together. The good agreement of the statistics of simulations
$S950$
and
$S950$
with (3.1) suggests that the present minimal unit for the large-scale structures in the shift-reflect subspace does not greatly limit the generation of good statistics. The statistics of the upper-branch solution
$U950$
are found to show reasonably good agreement with those of
$S950$
roughly below
$y\simeq 0.4h$
–
$0.5h$
, despite the fact that the solution itself is a relative equilibrium, not able to fully describe the chaotic, quasiperiodic dynamics of simulation
$S950$
(i.e. bursting (Flores & Jiménez Reference Flores and Jiménez2010; Hwang & Bengana Reference Hwang and Bengana2016)). We also note that the difference in statistics between an upper-branch solution and full simulation is expected. This level of difference also appears in other works that make the analogous comparison at low Reynolds numbers (e.g. Kerswell & Tutty Reference Kerswell and Tutty2007; Schneider, Eckhardt & Yorke Reference Schneider, Eckhardt and Yorke2007; Park & Graham Reference Park and Graham2016). As these authors indicate, this observation suggests that computation of unstable time-periodic orbits would be a more promising way to represent the dynamics of coherent structures in simulation
$S950$
(Kawahara & Kida Reference Kawahara and Kida2001; Gibson et al.
Reference Gibson, Halcrow and Cvitanovic2008; Willis et al.
Reference Willis, Cvitanovic and Avila2013). Also as expected, the statistics of the lower-branch solution
$L950$
do not show such a level of agreement – the lower-branch solution sits on the edge state of simulation
$S950$
. These results indicate that only the upper-branch solution is statistically similar to simulation
$S950$
.
Finally, the invariant solutions and the solution trajectories of simulations
$S950$
and the symmetry-constrained
$S950$
are projected onto the planes
$E_{streak}$
–
$E_{vor}$
and
$I$
–
$D$
, reported in figure 5. Here,
$E_{streak}$
and
$E_{vor}$
respectively represent energy of the streak and that of the streamwise vortices, and are defined by
$E_{streak}=[1/(2V)]\int _{V}(u^{\prime }/U_{m})^{2}\,\text{d}V$
and
$E_{vor}=[1/(2V)]\int _{V}(v/U_{m})^{2}+(w/U_{m})^{2}\,\text{d}y$
, where
$u^{\prime }$
is the streamwise velocity fluctuation.
$I$
and
$D$
are respectively the energy input and dissipation of the system, defined by
$I=-(1/V)\int _{V}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p\,\text{d}V$
and
$D=-(1/V)\int _{V}\boldsymbol{u}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }((\unicode[STIX]{x1D708}_{T}/2)(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})))\,\text{d}V$
with
$\boldsymbol{u}=(u,v,w)$
and
$\unicode[STIX]{x1D708}_{T}=\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{t}$
. We note that if
$E\equiv [1/(2V)]\int _{V}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}V$
, then
$\text{d}E/\text{d}t=I-D$
. In the plane
$E_{streak}$
–
$E_{vor}$
(figure 5
a), solution trajectories for both simulations
$S950$
and the symmetry-constrained
$S950$
reveal that their
$E_{streak}$
and
$E_{vor}$
are found to be slightly negatively correlated, indicating the presence of the self-sustaining process given by the interaction between the streak and the streamwise vortices. The solution trajectories of these simulations are maintained roughly around
$20I_{l}$
and
$20D_{l}$
in the
$I$
–
$D$
plane (figure 5
b); these values are an order of magnitude larger than those at low Reynolds numbers (e.g. Kawahara & Kida Reference Kawahara and Kida2001; Gibson et al.
Reference Gibson, Halcrow and Cvitanovic2008; Willis et al.
Reference Willis, Cvitanovic and Avila2013, Reference Willis, Cvitanovic and Avila2016). This is essentially due to the high Reynolds number considered in the present study, which leads to the energy input to the system being substantially larger than that at low Reynolds numbers. This indicates that the role of the eddy viscosity introduced here differs from that of the molecular one, because it enables us to maintain the large energy input to the system occurring in a high-Reynolds-number flow unlike the molecular viscosity. In both the planes
$E_{streak}$
–
$E_{vor}$
and
$I$
–
$D$
, the upper-branch solution
$U950$
is found to be placed in the middle of the turbulent trajectories of simulations
$S950$
and the symmetry-constrained
$S950$
. The lower-branch solution
$L950$
is almost completely separated from the trajectories, consistent with the observation on the first- and second-order statistics made with figure 4.
Phase portrait in the planes (a)
$E_{streak}$
–
$E_{vor}$
and (b)
$I$
–
$D$
at
$Re_{m}=38\,133$
, where
$I_{l}$
and
$D_{l}$
indicate the corresponding production and dissipation of laminar flow at the same Reynolds number. Colours: red,
$S950$
; dashed green, the symmetry-constrained
$S950$
; blue,
$U950$
; orange,
$L950$
.
3.3 Bifurcation with the Reynolds number
Finally, to understand the connection between the invariant states of the Navier–Stokes equation at transitional Reynolds number and those of this study at fairly high Reynolds numbers, numerical continuation is further performed by gradually lowering the Reynolds number using the computed invariant solutions with high resolution at
$Re_{\unicode[STIX]{x1D70F}}\simeq 950$
(i.e.
$U950$
and
$L950$
). Figure 6 reports bifurcation of the invariant solutions with the Reynolds number. The solutions for
$C_{s}=0.30$
reveal a saddle-node bifurcation at a fairly low critical Reynolds number around
$Re_{m,c}\simeq 1693$
–
$1707$
(
$Re_{\unicode[STIX]{x1D70F},c}\simeq 71$
). The Smagorinsky constant
$C_{s}$
is subsequently lowered around this low
$Re_{m}$
, such that invariant solutions of the Navier–Stokes equation are retrieved. The exact solutions, directly linked to the large-scale structures in the minimal unit, are found to emerge approximately at
$Re_{m,c}\simeq 1120$
(
$Re_{\unicode[STIX]{x1D70F},c}\simeq 52$
) via saddle-node bifurcation. We note that this
$Re_{m,c}$
is much lower than
$Re_{m,c}\simeq 1867$
(
$Re_{\unicode[STIX]{x1D70F},c}\simeq 68.2$
) of the invariant solutions found by Park & Graham (Reference Park and Graham2016) with a similar box size, due to the different symmetry imposed in the present study (their solution
$P4$
with
$L_{x}\times L_{z}=\unicode[STIX]{x03C0}h\times \unicode[STIX]{x03C0}/2h$
).
Bifurcation diagram of the invariant solutions in the plane
$Re_{m}{-}Re_{\unicode[STIX]{x1D70F}}$
: ○,
$C_{s}=0.0$
; ●,
$C_{s}=0.30$
.
The exact solutions of the Navier–Stokes equations (
$C_{s}=0.0$
) are finally continued by increasing the Reynolds number. While the upper-branch solution could not be continued for
$Re_{m}\gtrsim 2223$
, the lower-branch solution is obtained up to
$Re_{m}=8889$
. Interestingly, the lower-branch solution of the Navier–Stokes equation is found to generate smaller skin friction than those solutions with elevated eddy viscosity (
$C_{s}=0.30$
) on increasing the Reynolds number. However, it should be noted that the eddy viscosity used here is designed to model all the surrounding smaller-scale structures, including the near-wall motions generating an appreciable amount of turbulent skin friction. The difference between the solutions with
$C_{s}=0.0$
and
$C_{s}=0.30$
probably originates from this nature of the eddy viscosity.
4 Concluding remarks
Invariant solutions corresponding to large-scale turbulent motions at high Reynolds numbers have been obtained following the approach of Rawat (Reference Rawat2014) and Rawat et al. (Reference Rawat, Cossu, Hwang and Rincon2015, Reference Rawat, Cossu and Rincon2016), where surrounding small-scale structures are modelled with an eddy viscosity (i.e. Smagorinsky model with the artificially elevated
$C_{s}$
as a continuation parameter). Here, we show that these solutions, and in particular the upper-branch solution, can be obtained at much higher Reynolds numbers (up to
$Re_{\unicode[STIX]{x1D70F}}\simeq 1000$
), a different regime where the large-scale structures in the form of long streaks (VLSMs) and quasistreamwise vortical structures (LSMs) emerge energetically in turbulent channel flow. This finding suggests that the large-scale structures are probably organised around these invariant solutions, providing direct evidence of their significance at sufficiently high Reynolds numbers. This finding also tightly establishes their relation to the self-sustaining process of the large-scale structures (Hwang & Cossu Reference Hwang and Cossu2010b
; Rawat et al.
Reference Rawat, Cossu, Hwang and Rincon2015; Hwang & Bengana Reference Hwang and Bengana2016), while implying that the dynamical systems approach, plus visualisation of the relevant phase space, could be used to enlighten at least some aspects of the given coherent structures at fairly high Reynolds numbers. However, the invariant solutions found here do not fully represent the dynamics occurring in real flows (i.e. bursting (Flores & Jiménez Reference Flores and Jiménez2010; Hwang & Bengana Reference Hwang and Bengana2016)). In this respect, extending the present approach to computation of the relative periodic orbits would be a fruitful path to follow towards more accurate modelling of the large-scale dynamics (Kawahara & Kida Reference Kawahara and Kida2001; Willis et al.
Reference Willis, Cvitanovic and Avila2013, Reference Willis, Cvitanovic and Avila2016). Finally, it should be mentioned that the self-sustaining energy-containing motions in wall-bounded turbulent flows at high Reynolds numbers appear in a self-similar form throughout the entire logarithmic region (Hwang & Cossu Reference Hwang and Cossu2011; Hwang Reference Hwang2015), as originally hypothesised by Townsend (Reference Townsend1976) (i.e. attached eddy hypothesis). Each of the energy-containing motions is typically characterised by its spanwise length scale (Hwang & Cossu Reference Hwang and Cossu2011; Hwang Reference Hwang2015; Hwang & Bengana Reference Hwang and Bengana2016), suggesting that the invariant solutions with different spanwise length scales may be linked to these self-similar motions. This also implies that a thorough investigation of the invariant solutions with different sets of the length scales needs to be carried out in order to clarify the relation between the ‘real’ and ‘exact’ coherent structures using a numerical experiment given in the present study. Exploring such a link between the invariant solutions and the concept of the attached eddies might unveil the nature of the self-sustaining coherent structures. This would be an important step towards a consistent theoretical description of wall-bounded shear flows in a wide range of Reynolds numbers, because the attached eddy hypothesis has provided an important theoretical basis for a consistent statistical description of high-Reynolds-number wall turbulence (e.g. Townsend Reference Townsend1976; Perry & Chong Reference Perry and Chong1982).
Acknowledgement
Y.H. was supported by the Engineering and Physical Sciences Research Council (EPSRC) in the UK (EP/N019342/1).
