3.1. Tracking the edge

Using the edge-tracking procedure, we attempt to identify an invariant solution (Kawahara et al. Reference Kawahara, Uhlmann and van Veen2012) to square-duct flow and the associated coherent structure. The edge tracking is started with the bisection (Itano & Toh Reference Itano and Toh2001; Skufca et al. Reference Skufca, Yorke and Eckhardt2006) followed by the custom Newton–Krylov solver (Viswanath Reference Viswanath2007, Reference Viswanath2009) developed within the Nektar (Moxey et al. Reference Moxey2020) framework. The initial condition for the bisection is formed as a superposition of a snapshot of the turbulent flow at ${Re}=4000$ and the laminar solution. We initially attempted bisection using the full-space configuration (i.e. no assumption on the solution symmetry has been made), but the process failed to arrive at a regular solution despite repeated attempts and using different initial conditions. Instead, each attempt resulted in shadowing of an irregular trajectory, with both the perturbation energy and friction factor decreased from respective turbulent values but changing in a recurrent manner. Similar behaviour of the full-space bisection process has been reported for the case of pipe flow (Duguet et al. Reference Duguet, Willis and Kerswell2008, Reference Duguet, Willis and Kerswell2010; Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013). We attribute the problem with the bisection in the full space to the fact that the edge state itself might be chaotic.

We expect that similarly to the pipe flow (Duguet et al. Reference Duguet, Willis and Kerswell2008, Reference Duguet, Willis and Kerswell2010; Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013), restricting the turbulent dynamics to one of the symmetric sub-spaces would allow us to identify a simple invariant state laying on the edge. Consequently, in the next attempt, we limit the bisection to consider only the $\pi$ -rotationally symmetric sub-space of the full state space. This is achieved by slicing the computational domain across the wall bisector $x=0$ and imposing ‘periodic’ boundary conditions on the cut $x=0$ , of the form:

(3.1) $ \begin{equation} \begin{cases} u(0,-y,z)=-u(0,y,z), \\ v(0,-y,z)=-v(0,y,z), \\ w(0,-y,z)=+w(0,y,z), \end{cases} \end{equation}$

which correspond to the $\pi$ -rotational symmetry with respect to the duct centre $(x,y)=(0,0)$ . In this symmetric sub-space, at ${Re}=4000$ , we observe no relaminarisation within observation time, similar to the full-space configuration. Initial conditions for the bisection are the same as in the full-space attempt and, this time, consecutive bisections quickly settle to shadow what seems to be a regular trajectory. The convergence of the bisection is illustrated in figure 2 via variations of $E_{3D}$ during consecutive bisection iterations. At the initial stage, when the $\pi$ -rotationally symmetric sub-space is considered (depicted by dash-dotted lines), bisection settles to shadow a state that features a significant decrease in perturbation energy and friction factor (not shown) compared with the turbulent flow. With the progress of the bisection, the solution autonomously exhibits mirror symmetry of the flow velocity across the wall bisector $y=0$ , leading to the mirror symmetry with respect to the other wall bisector $x=0$ . At this stage, we reinforce the symmetry requirement by imposing double mirror symmetries across both wall bisectors. This is done by casting the current solution onto a quarter of the domain and imposing a mix of Dirichlet–Neumann boundary conditions onto the velocity vector field, of the form:

(3.2) $ \begin{equation} \begin{cases} u(0,y,z)=0, \\ \dfrac {\partial v}{\partial x}(0,y,z)=\dfrac {\partial w}{\partial x}(0,y,z)=0 \end{cases}\quad \text{ and }\quad \begin{cases} v(x,0,z)=0, \\ \dfrac {\partial u}{\partial y}(x,0,z)=\dfrac {\partial w}{\partial y}(x,0,z)=0. \end{cases} \end{equation}$

Variations of $E_{3D}$ from this stage are depicted using solid grey and later black (the change corresponds to the restart of the bisection done to speed up the process) lines in figure 2. The reader might note that in figure 2 for some time, dash-dotted (initial, $\pi$ -rotational symmetry (3.1)) and solid grey (double mirror symmetry (3.2)) shadow the same trajectory, indicating that both types of symmetry restrictions cause the bisection to converge onto the same solution. Edge tracking is finalised with the custom Newton–Krylov iteration procedure (Viswanath Reference Viswanath2007, Reference Viswanath2009), using the initial condition corresponding to the bisection solution around $t=4000$ , which converges rapidly onto the solution marked by a thick dashed, green line in figure 2.

Figure 2. Variation of the perturbation energy $E_{3D}$ during bisection (curves) and the converged steady travelling wave (dashed, green horizontal line) through Newton–Krylov iterations with the initial guess obtained from the final edge tracking in the double mirror-symmetric sub-space. The edge-tracking process is started using the laminar solution and the turbulent snapshot cast onto the $\pi$ -rotationally symmetric configuration. Dash-dotted lines depict the initial bisection stage with only the $\pi$ -rotational symmetry enforced with respect to the duct centreline; solid grey lines correspond to the imposition of mirror symmetries across both wall bisectors, which are observed to appear autonomously in the initial bisection stage; and solid black lines correspond to consecutive restarts of the process. The Newton–Krylov iteration is performed in the double mirror-symmetric sub-space using the initial guess taken from the final edge tracking at $t\approx 4000$ time units. The dashed, green horizontal line represents the converged steady travelling wave. Variation of $E_{3D}$ from the turbulent simulation is provided for reference using a thin curve.

We verify the accuracy of our edge tracking by comparing it with an invariant state computed by one of the authors (Okino Reference Okino2014) using a global spectral, Chebyshev–Fourier Galerkin method resulting in the same state, and by varying both the polynomial order and the number of Fourier modes used for the approximation in our current attempt while repeating the Newton–Krylov iteration step. The polynomial order $p$ is varied from the selected $p=24$ down to $20$ and up to $28$ , while either $M=128$ or $144$ complex Fourier modes are applied in the streamwise direction. Each time, the Newton–Krylov process converges onto the same solution. However, as the number of degrees of freedom is increased, the limitation of the Newton iteration correction step needs to be taken as it seems that the radii of convergence decreases.

3.2. Characterisation of the identified state

The identified invariant state features a slight increase in the hydraulic resistance, compared with the laminar solution. In this state, the bulk velocity $w_b$ is decreased to $w_{b}\approx 0.474$ from $w_b\approx 0.477$ of the laminar solution. The state is mirror-symmetric across both wall bisectors and has a form of a wave with a distinct, streamwise localised four-vortex structure travelling downstream at the constant speed of $0.675$ units ( ${\approx}1.424w_{b}$ , where $w_b$ corresponds to the identified invariant state). Figure 3 illustrates the contours of the second invariant of the velocity gradient tensor $Q=(1/2) ( \|\Omega \|^2 - \|S\|^2 )$ (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988), where $S = ({1}/{2}) ( \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^T )$ is the symmetric part of the velocity gradient tensor and $\Omega = ({1}/{2}) ( \boldsymbol{\nabla} \boldsymbol{u} - \boldsymbol{\nabla} \boldsymbol{u}^T )$ the antisymmetric part of the velocity gradient tensor. The contours are taken at $Q/w^2_{b}=1.5\times 10^{-1}$ and show two pairs of streamwise vortices around the wall bisector $y=0$ , next to the opposite wall pair $x=\pm 1$ and captured around the middle of the duct ( $z\approx 4\pi$ ). The corresponding streamwise velocity perturbation with relation to the laminar velocity distribution $w_L$ is illustrated in figure 4 and shows streamwise slices at different positions away from the duct centreline and at selected duct cross-sections. We conclude that while the identified invariant state remains streamwise localised, the resulting change of the streamwise velocity from the laminar flow distribution is global.

Figure 3. Identified invariant state visualised by contours of the second invariant of the velocity gradient tensor normalised by $w_{b}$ , taken at $Q/w_{b}^2=1.5\times 10^{-1}$ .

Figure 4. Contours of the difference of the streamwise velocity component $w$ of the identified invariant state with respect to the laminar solution $w_L$ . Slices are taken (a) along the streamwise direction $z$ at $x=0.5$ , $x=0.8$ , $y=0.5$ and $y=0.8$ and (b) on duct cross-sections $(x,y)$ placed at $z=0, 2\pi , 4\pi , 6\pi$ . Figure 15(a) shows the mean streamwise velocity distribution of the identified state.

Streamwise localisation property of the identified invariant state can be extracted from the density of the cross-flow and the streamwise velocity perturbation energies

(3.3) $ \begin{equation} E_{\perp }(z) = \frac {1}{4}\int ^{1}_{-1}\!\int ^{1}_{-1} \frac {1}{2}\big(u^2+v^2\big)\,\textrm {d}x\,\textrm {d}y \quad\text{ and }\quad E_{||}(z) = \frac {1}{4}\int ^{1}_{-1}\!\int ^{1}_{-1} \frac {1}{2}(w-w_L)^2\,\textrm {d}x\,\textrm {d}y. \end{equation}$

$E_{\perp }(z)$ provides a norm-like quantification of the in-plane motions and $E_{||}(z)$ quantifies changes to the streamwise velocity. Variations of those quantities with the streamwise coordinate $z$ normalised by the streamwise length $L_z$ of the domain are depicted in figure 5 and show significant localisation around $z\approx 4\pi$ ( $z/L_z \approx 0.5$ ). The plots in figure 5 show that localisation is maintained as the Reynolds number is decreased. Continuing the solution with $L_z$ indicates that with the increase of the domain length, the solution’s streamwise length remains nearly constant so that it occupies relatively smaller portions of the domain. We note that our results suggest that $L_z=8\pi$ is close to the saddle-node bifurcation length, below which the identified solution does not exist. Also, we have not observed a connection to any of the streamwise extended, periodic solutions neither by continuing with ${Re}$ nor with $L_z$ .

Figure 5. Streamwise profile of the cross-flow energy $E_{\perp }$ and of the streamwise velocity perturbation energy $E_{||}$ of the identified invariant state illustrating streamwise localisation of the invariant solution. The plots in panel (a) correspond to cases which differ in the streamwise length and Reynolds number. Solid lines depict variation for $L_z=8\pi$ for a range of Reynolds numbers, while dashed lines illustrate the influence of varying the computational domain length at a fixed ${Re}=3050$ . Plots in panel (b) show $E_{\perp }$ and $E_{||}$ for different streamwise lengths and show that the streamwise length of the velocity perturbation does not vary much with the streamwise period $L_z$ .

At this stage, we would like to point out that the identified invariant solution has been found under the restriction of reflectional symmetry across wall bisectors and that solutions identified within the symmetric sub-space are necessarily solutions in the full space (Avila et al. Reference Avila, Mellibovsky, Roland and Hof2013) corresponding to physical (symmetric) states of the flow. We examine the consequences of the imposition of symmetry in the following sections.

3.3. Symmetric sub-space edge

Following the identification of an invariant solution, we examine the stability of this travelling wave by means of the Arnoldi iterations (Viswanath Reference Viswanath2009). We test the stability of the state, both in the full as well as in the symmetric sub-space. Examined symmetries include the $\pi$ -rotational symmetry (3.1) with respect to the duct centre and the double mirror symmetry (3.2) across both wall bisectors. This analysis shows that at ${Re}=4000$ , when confined to $\pi$ -rotationally symmetric sub-space, the state has one purely real unstable eigenvalue, while in the full space (no restriction on symmetry), the state already has multiple real, unstable eigenvalues. In addition, at ${Re}=4000$ under the $\pi$ -rotational symmetry (3.1), the initial condition formed by the identified LB solution perturbed in the unstable direction leads either to an uneventful laminarisation or towards a persistent non-laminar state (we do not observe relaminarisation even after a long time), depending on the sense of the perturbation. In the non-laminar state, the perturbation energy $E_{3D}$ , zero mode energy $E_0$ and the friction factor $f$ that we monitor have values in the range that we have observed in a long-time turbulent flow simulation, suggesting the state of the flow lands on the turbulent attractor. We observe this type of behaviour down to approximately ${Re}=3600$ , where we observe possible relaminarisation to happen within ten-thousand time units, indicating that as the Reynolds number is decreased, turbulence becomes a transient phenomenon within the observation time. Consequently, in the symmetric sub-space (both (3.1) or (3.2)), the identified state is an edge state, and its stable manifold of co-dimension one splits (at least locally) the state space into two parts and forms an edge. This edge separates the laminar attraction basin from the turbulent attractor under the $\pi$ -rotational symmetry (3.1) and at sufficiently high Reynolds numbers. As the Reynolds number is reduced, it seems that only the laminar attractor remains, while the turbulent one is replaced by a chaotic transient with the edge that separates them only locally.

Limiting considerations to the investigation of the edge in the symmetric sub-space, it is easy to track this state in ${Re}$ with an arc-length continuation method (Keller Reference Keller1977; Dijkstra et al. Reference Dijkstra2014). The determined state exists down to slightly below ${Re}\approx 3030$ , where it is created. It turns out that it is also possible to identify the upper branch (UB) of this travelling-wave solution. The bifurcation diagram is shown in figure 6 and also outlines the change of stability properties with ${Re}$ . We have tested the stability of the identified state while traversing both of the solution branches. For most cases, we applied double mirror-symmetric sub-space (3.2), only occasionally extending the solution space to include $\pi$ -rotationally symmetric state (3.1). On the UB, the invariant solution is also streamwise localised, but has multiple unstable eigenvalues. While on the lower branch (LB), the travelling wave remains to be an edge state in the symmetric sub-space, with a single, purely real unstable eigenvalue. Exactly at the turning point at approximately ${Re}=3030$ and upon turning onto the UB, there appears a second unstable, also real eigenvalue. This behaviour suggests that the bifurcation does not result in one of the solution branches becoming stable, as, for example, reported by Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013) for the case of circular-pipe flow. Rather, both the LB and UB solution branches are created unstable. Somewhere between ${Re}=3040$ and ${Re}=3046$ , on the UB, an additional pair of unstable complex-conjugated eigenvalues appear. Eventually, around ${Re}=3058$ , the two purely real, unstable eigenvalues disappear and are replaced by an additional complex-conjugated pair. We have tested that this character of stability holds at least up to ${Re}\approx 4000$ . Flow topology of the UB solution at ${Re}=4000$ is shown in figures 7 and 8. Comparison of the characteristic flow quantities at ${Re}=4000$ for the laminar, turbulent and both the UB and LB invariant solutions is given in table 1.

Figure 6. Bifurcation diagram for the identified travelling wave with Reynolds number ${Re}$ showing (a) the $E_{3D}$ and (b) friction factor. The thick dashed line in panel (b) represents the laminar flow friction factor. Stability in the symmetric sub-space is indicated with different line styles: the solid line identifies the lower branch (LB) with one purely real unstable eigenvalue up to the bifurcation point, the dashed fragment corresponds to two purely real unstable eigenvalues on the upper branch (UB) starting exactly at the bifurcation point, the thick fragment to two purely real unstable and an additional complex-conjugate unstable pair on the upper branch, and the solid, thin line to two complex-conjugate unstable pairs on the upper branch. Most of the presented stability results have been obtained within the double mirror-symmetric sub-space (3.2). Selected cases have also been examined in the $\pi$ -rotational symmetric sub-space (3.1).

Figure 7. Invariant state on UB visualised by contours of the second invariant of the velocity gradient tensor normalised by $w_b$ at ${Re}=4000$ , taken at $Q/w^2_{b}=1.5\times 10^{-1}$ .

Figure 8. Contours of the difference of the streamwise velocity component $w$ of the invariant state on the UB at ${Re}=4000$ with respect to the laminar solution $w_L$ . The position of slices is the same as in figure 4. The mean streamwise velocity distribution of the identified state is shown in figure 15(c).

The change of the solution stability from one unstable direction for the LB solution to two unstable directions of the UB solution precisely at the bifurcation needs to be commented on. The reader will find a more detailed discussion of our approach to this issue in Appendix A, with just a brief characterisation provided here. We test the evolution of the flow state in the symmetric sub-space at ${Re}=3040$ (close to the bifurcation point) from the LB and UB solution perturbed in the direction determined by eigenvectors associated with the unstable eigenvalues of those solutions. The unstable direction is a line for the LB and a plane for the UB. We recall that at this value of the Reynolds number, we found no persistent non-laminar solution and only a chaotic transient seems to be available, irrespective of the applied symmetry constraint.

Our results show that in the case of both the perturbed LB and UB solutions, provided the perturbation amplitude is sufficiently small, the flow stays on the respective equilibrium solution for some time and either follows with an uneventful laminarisation or goes onto a short-lived, non-laminar transient that seems to follow the remains of the turbulent attractor and possibly reaches towards other invariant solutions. It seems reasonable to assume that this behaviour, i.e. uneventful laminarisation or a chaotic transient detour available from both solution branches, is because the bifurcation that gives rise to both solution branches happens on the edge sub-space (co-dimension one sub-space) and that both the LB and specifically the UB solutions remain in the edge sub-space, the stable manifold of the LB solution. We think this situation persists for a range of ${Re}$ close to the bifurcation point. That is, in a sub-space confined to the stable manifold of the LB solution (where the LB solution is a relative attractor), the LB solution is a node and the UB, which also exists in this sub-space, is a saddle.

Following this reasoning, we make a conjecture that the bifurcation in which both branches are formed also takes place on the edge, where it is a saddle-node bifurcation (the LB solution is a node, the UB a saddle in the sub-space confined to the edge). This character of the stability further implies that other invariant solutions should exist and possibly form at still lower values of the Reynolds number than the current solutions. Eventually, there should exist a solution pair that forms in a saddle-node bifurcation similar to that reported by Avila et al. (Reference Avila, Mellibovsky, Roland and Hof2013) for the pipe flow. This solution remains to be found, which might be made difficult by the transient nature of the non-laminar solution. Finally, with the possible existence of other invariant solutions at lower values of the Reynolds number, we conjecture that similarly to Duguet et al. (Reference Duguet, Willis and Kerswell2008), there can be many separate solutions confined to the edge (multiple local edge states) bifurcating in respective saddle-node (limited to co-dimension one manifold) bifurcations. Each such state should then be a local edge state in the sense that its stable manifold would split the phase space only locally and the state itself would remain a relative attractor only on a portion of the edge.

4.1. Flow symmetrisation

The laminar solution and the time-averaged velocity fields possess a number of symmetries related to the duct geometry, and so does the invariant state identified in § 3. We consider a heuristic indicator to examine changes to the flow over time and quantify the deviation of the flow state from the $\pi$ -rotational symmetry. The measure, normalised with the energy of the mean mode $E_0$ is of the following form:

(4.1) $ \begin{align} f_\pi &=\frac {1}{\Omega E_0}\iiint _{\Omega } \big[(u(x,y, z)+u(-x,-y, z))^2+(v(x,y, z)+v(-x,-y, z))^2\notag\\&\quad + (w(x,y, z) - w(-x,-y, z))^2\big]^{1/2}\textrm {d}\Omega \end{align}$

and quantifies the overall volume-averaged deviation of the flow state from the $\pi$ -rotational symmetry. This means that this measure decreases as the flow field becomes more symmetric. Similar symmetry measures quantifying deviation of the flow velocity from either of the mirror symmetries across wall bisectors ( $x=0$ and $y=0$ ) are also considered and noted as $f_x$ and $f_y$ . Time variations of all three measures are shown in figure 10(a). The thick dashed line corresponds to the time average $\langle f_{\pi } \rangle$ and thinner lines represent consecutive displacements ( $k=1,2,3$ ) by the value of the standard deviation $\sigma$ of $f_{\pi }$ from $\langle f_{\pi } \rangle$ . Data show that usually, flow symmetry measure remains within the $\langle f_{\pi } \rangle \pm \sigma$ range, but on occasion, $f_\pi$ may be substantially decreased, indicating that the flow enters periods of increased symmetry. Interludes of decreased $f_{\pi }$ measure, when it drops below the $\langle f_{\pi } \rangle -2\sigma$ level, correlate well with the decrease of the friction factor and perturbation energy shown in figure 9(a,b) and marked with B and C. Surprisingly, time variations of all three measures overlap almost completely, meaning that mirror symmetries appear autonomously at the same time as the flow state approaches the $\pi$ -rotationally symmetric form. The fact that both mirror-symmetries seem to follow closely suggests that the flow develops states that are mirror-symmetric against $x=0$ and $y=0$ bisectors, but does not seem to support states that are symmetric only against one of the wall bisectors (either $x=0$ or $y=0$ ). This further implies and is consistent with our observations that states with only two vortices (similar to the edge trajectory described by Biau & Bottaro (Reference Biau and Bottaro2009)) next to one of the duct walls (the two-vortex states) are unlikely at the Reynolds number $4000$ .

Figure 10. Temporal variation of (a) the symmetry indicator $f_{\pi }$ (solid line) (4.1) and the corresponding $f_x$ and $f_y$ (dashed lines) and of (b) the vortex placement indicator $I$ (4.4). In panel (a), the thick, horizontal and dashed line indicates the mean value $\langle f_{\pi } \rangle$ (the temporal average of $f_\pi$ ), while thinner dashed lines represent corresponding $\langle f_{\pi } \rangle +k\sigma$ ( $k=1,2,3$ ) levels with $\sigma$ the standard deviation of $f_\pi$ .

4.2. Eight- and four-vortex state

In addition to symmetries that might transpire in the flow state, one of the distinct characteristics of the turbulent duct flow is the onset of turbulence-driven, Prandtl’s (Prandtl Reference Prandtl1926; Pirozzoli et al. Reference Pirozzoli, Modesti, Orlandi and Grasso2018) secondary motions of the second kind that manifest by the onset of pairs of vortices placed symmetrically around corner or wall bisectors. In the time-averaged velocity field, those motions manifest as either eight- or four-vortex configurations (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007, Reference Uhlmann, Kawahara and Pinelli2010; Gavrilakis Reference Gavrilakis1992), depending on the selected averaging time window (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007). Following Uhlmann et al. (Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007), we examine the flow by computing the square of the mean streamwise component of vorticity:

(4.2) $ \begin{equation} S_i(t) = \iint _{\Omega _{\perp i}} \left [\frac {1}{L_z}\int _0^{L_z}\omega _z\textrm {d}z\right ]^2 \textrm {d}\Omega _{\perp i} \end{equation}$

over four triangular regions $\Omega _{\perp i}$ (see figure 11), defined by diagonals:

(4.3) $ \begin{equation} \begin{cases} \Omega _{{\perp 1}} : \{(x,y)|y\gt x \bigcap y\lt -x\}, \Omega _{{\perp 3}} : \{(x,y)|y\lt x\bigcap y\gt -x\}, \\ \Omega _{{\perp 2}} : \{(x,y)|y\lt x \bigcap y\lt -x\}, \Omega _{{\perp 4}} : \{(x,y)|y\gt x\bigcap y\gt -x\} \end{cases} \end{equation}$

and further normalised to produce the vortex placement indicator

(4.4) $ \begin{equation} I(t) = \frac {S_1+S_3-S_2-S_4}{S_1+S_2+S_3+S_4}. \end{equation}$

Values of $I$ that remain close to zero indicate states where, on average, streamwise vorticity is equally distributed between the $x=\pm 1$ and the $y=\pm 1$ wall pairs, suggesting an overall eight-vortex state. However, relatively large positive or negative values of $I$ suggest the concentration of vorticity along one of the two wall bisectors and in the proximity of either the $x=\pm 1$ ( $I\gt 0$ ) or the $y=\pm 1$ ( $I\lt 0$ ) wall pair and is an indication towards an overall four-vortex state of the flow. Temporal variation of $I$ is shown in figure 10(b). In general, the value of $I$ remains strictly positive (or negative) over extended times, suggesting that, on average, there is a dominating tendency of the flow to favour one of the two four-vortex configurations (either the $y=\pm 1$ or the $x=\pm 1$ wall pair). Change of the configuration is possible but does not happen often. The absolute value of $I$ usually remains below $0.5$ . This suggests that in the usual state of the flow, vortices remain more pronounced in the proximity of one of the wall pairs, but at the same time, vorticity magnitude over the other wall pair remains significant. The reader might note that occasionally, $I$ reaches larger absolute values, indicating the transitional concentration of vorticity around one of the wall bisectors and suggesting the onset of a relatively well-defined four-vortex state with vortices either in the $y=\pm 1$ or $x=\pm 1$ orientation with the relatively quiescent flow next to the remaining walls. Comparing temporal variations of $f_\pi$ (as well as $f_x$ and $f_y$ ) and $I$ , we observe that during periods over which $f_\pi$ (as well as $f_x$ and $f_y$ ) decreases below $\langle f_{\pi } \rangle -2\sigma$ and the flow field becomes increasingly symmetric, the absolute value of $I$ increases above $0.75$ . This indicates that during symmetrisation, vorticity concentrates around one of the two wall bisectors suggesting that the four-vortex state should be expected. Also, this state coincides with a decrease in perturbation energy and hydraulic losses.

Figure 11. Regions $\Omega _{\perp i}$ splitting the duct cross-section $\Omega _\perp$ used to compute $S_i$ ( $i=1$ , $2$ , $3$ , $4$ ).

Using symmetry quantification, perturbation energy variation and the vortex placement indicator $I$ , we select three time periods for further discussion and flow averaging. The first period marked A, spans the $2000$ time units from $t=10^4$ up to $t=1.2\times 10^4$ , and corresponds to the period over which the value of $I$ switches from positive to negative and leads to the usual eight-vortex state when averaged. At the same time, during this period, the flow, at least as quantified by symmetry, friction factor and perturbation energy, remains close to respective averaged values. The second period marked B and the third C span a much shorter duration of $200$ time units and focus on periods when the absolute value of the vortex placement indicator $I$ becomes large, while perturbation energy, the symmetry measure and hydraulic resistance decrease substantially. This suggests the onset of a relatively quiescent, regular state of the flow with the value of $I$ implying a four-vortex state. Period B ranges from $t=6.4\times 10^3$ to $t=6.6\times 10^3$ , and the positive value of $I$ corresponds to vortices oriented along the wall bisector $y=0$ and positioned next to the walls $x=\pm 1$ . Period C spans the time from $t=1.26\times 10^4$ to $t=1.28\times 10^4$ with vortices oriented around the bisector $x=0$ and positioned next to the $y=\pm 1$ wall pair.

The character of the flow state during selected time periods A, B and C is shown in figures 12 and 13, depicting instantaneous snapshots of the flow, representative of selected time periods and showing the iso-contours of the second invariant of the velocity gradient tensor (representing vortex structures) $Q$ , normalised by $\langle w_{b} \rangle$ for the turbulent flow in figure 12 and variation of the streamwise velocity $w$ in figure 13. Snapshots have been taken at $t=10^4$ (panel a), $t=6.4\times 10^3$ (panel b) and at $t=1.26\times 10^4$ (panel c). Comparing the three states depicted in figure 12, we note that during the relatively quiescent interludes B and C, the region indicated by the $Q$ iso-surface, occupies a much smaller portion of the domain in comparison to the more common flow character, represented by A. The distribution of the $Q$ iso-surfaces suggests streamwise localisation of the turbulent fluctuations in the cases B and C while, in the case of A, those seem to be rather uniformly distributed in the streamwise direction. We observe that, in the case of the snapshot shown in A, vortex structures occupy the domain uniformly, while in cases B and C, those structures seem to be organised along one of the wall bisectors of two opposing wall pairs and absent from the remaining pair. We suspect that this vorticity configuration stems from vortex pairs formed next to the respective wall pairs and influences the streamwise velocity distribution next to the duct walls, which we shall now discuss.

Figure 12. Turbulent states visualised by contours of the second invariant of the velocity gradient tensor normalised by $\langle w_{b} \rangle$ of the turbulent flow, taken at $Q/\langle w_{b} \rangle ^2=4\times 10^{-1}$ . A, B and C correspond to the time periods distinguished in figures 9 and 10. (a) Snapshot taken at $t=10^4$ and a quasi-uniform distribution of vortex structures along the streamwise and in-plane directions; (b,c) flow state at $t=6.4\times 10^3$ and $t=1.26\times 10^4$ , respectively, and the relatively quiescent structures with pronounced streamwise localisation and an apparent spanwise organisation of the structures, exemplified by ellipses in the cross-axial plane view on the right.

Figure 13. Snapshots of the streamwise velocity component $w$ taken along the streamwise direction $z$ at $x=\pm 0.8$ , $y=\pm 0.8$ and at duct cross-sections placed at $z=0, 2\pi , 4\pi , 6\pi$ . A, B and C correspond to the time periods distinguished in figures 9 and 10. (a) Snapshot at $t=10^4$ and noticeable low-velocity streaks in the proximity of the duct wall bisectors, more pronounced at $x=\pm 1$ and present but less distinct at $y=\pm 1$ . (b,c) Snapshot at $t=6.4\times 10^3$ and $t=1.26\times 10^4$ , respectively, which shows low-velocity streaks forming close to bisectors of opposing, $x=\pm 1$ (or $y=\pm 1$ ) walls of the duct i.e. at $x=\pm 0.8$ (or $y=\pm 0.8$ ), with the flow near the other two walls remaining relatively quiescent. For consistency, colour scale used for $w$ is the same as in figures 14 and 15.

Figure 13 shows snapshots of the streamwise velocity component $w$ at slices located at $x,y=\pm 0.8$ in the proximity of duct walls $x,y=\pm 1$ (panel a) and at evenly spaced cross-sectional planes at $z=0, 2\pi , 4\pi$ and $6\pi$ (panel b). The distinctive feature that can be immediately noticed on the streamwise oriented slices are regions of decreased streamwise velocity corresponding to low-velocity streaks that form close to wall bisectors and, generally, run through the entire length of the duct. The corresponding cross-sectional slices show the low-momentum flow being displaced from the wall proximity and towards the duct centre, which is a manifestation of the rotational flow component due to the vortex pair forming near wall bisectors. We note that in the snapshots corresponding to B and C, the low-velocity streaks are well defined only on the two opposing walls. The streaks appear on the slices $x=\pm 0.8$ in the proximity of the walls $x=\pm 1$ for B and on the slices $y=\pm 0.8$ in the proximity of the walls $y=\pm 1$ for C, while in the proximity of the walls $y=\pm 1$ for B and of the walls $x=\pm 1$ for C, the flow remains relatively quiescent. However, the flow depicted in A shows a more disturbed and not so well-defined structure. The low-velocity streaks remain present and well defined next to the wall pair $x=\pm 1$ and to a lesser degree in the proximity of the other wall pair $y=\pm 1$ .

Averaging over a selected time period and streamwise direction leads to an ensemble average of the form:

(4.5) $ \begin{equation} \langle \!\langle \boldsymbol{u}\rangle \!\rangle (x,y) = \frac {1}{L_z}\int ^{L_{z}}_0 \langle \boldsymbol{u}(x,y,z,t) \rangle \,\textrm {d}z, \end{equation}$

where $\langle \!\langle \boldsymbol{\cdot} \rangle \!\rangle$ indicates average over time and the $z$ coordinate. Choosing the averaging time to correspond to periods A, B or C, the resulting ensemble average $\langle \!\langle \boldsymbol{u}\rangle \!\rangle$ depicts the dominant configuration that the flow adopts over those periods. Figure 14 shows contours of the average streamwise vorticity $\langle \!\langle \omega _z\rangle \!\rangle$ in the upper row and components of the mean velocity in the lower row. The streamwise mean flow $\langle \!\langle w\rangle \!\rangle$ is shown using contour lines, while arrows depict in-plane mean velocity components $\langle \!\langle u\rangle \!\rangle$ and $\langle \!\langle v\rangle \!\rangle$ . Averaging over period A yields a well-defined eight-vortex topology, while averages over periods B and C lead to four-vortex configurations that form along respective wall bisectors.

Figure 14. Ensemble average $\langle \!\langle \boldsymbol{u} \rangle \!\rangle$ defined by (4.5) applied over respective time periods A, B and C. (a) Long-time average taken over the period designated as (a) A, (b) B and (c) C. The upper row shows contours of the streamwise vorticity component $\langle \!\langle \omega _z\rangle \!\rangle$ , with negative values indicated with dashed lines, and contours of the streamwise velocity component $\langle \!\langle w\rangle \!\rangle$ are shown in the lower row with arrows illustrating the in-plane velocity components $\langle \!\langle u\rangle \!\rangle$ and $\langle \!\langle v\rangle \!\rangle$ . For consistency, colour scale used for $w$ is the same as in figures 13 and 15.

We find the structure of the mean flow shown in figure 14(b) reminiscent of the LB invariant state that we characterised in § 3. Figure 15(a) shows the streamwise (contours) and in-plane (arrows) velocity components at the slice through the invariant state velocity vector field taken at the streamwise position corresponding to the streamwise maximum of the $E_{\perp }$ . The illustrated velocity field corresponds to the solution lying on the LB at ${Re}=4000$ . The figure provides an insight into the in-plane position of the two vortex pairs (centres marked by red-filled triangles) that have formed around the wall bisector $y=0$ . We note that the effect of those vortices is in the local modification of the distribution of the streamwise velocity component, i.e. vortex pairs rotate in such a way as to push the low-momentum fluid from the wall centre and into the core of the flow. Locally, this leads to the formation of low-velocity streaks positioned around the wall bisector $y=0$ and next to the walls $x=\pm 1$ . For comparison, figure 15(b) shows the in-plane (arrows) velocity component of the ensemble-averaged $\langle \!\langle \boldsymbol{u}\rangle \!\rangle (x,y)$ flow from the relatively quiescent period B of the turbulent flow (see figure 14 b) with vortices (centres marked by purple-filled squares) placed around the $y=0$ bisector. For reference, positions of the centres of vortices from averaged flow from period C are also marked in the figure (cyan-filled squares) and compared with the orange-filled triangles that identify vortex centres of the identified LB invariant state, but in a configuration where vortices form around the $x=0$ bisector ( $\pi /2$ -rotation). It turns out that the four vortices in the identified LB invariant state appear at the cross-sectional positions close to those in the quiescent states with the four-vortex configuration of significant symmetrisation, observed in turbulent square-duct flow. Finally, for comparison, figure 15(c) shows the slice at the same position as in figure 15(a) through the UB solution at ${Re}=4000$ (see also figure 8), implying that the positions of the vortices in the UB invariant state are not consistent with those in the quiescent states of the turbulent flow.

Figure 15. (a) Slice through the velocity field of the LB ( ${Re}=4000$ ) invariant solution at the streamwise position of the maximum of cross-flow energy ( $z=4\pi$ , see figure 5), (b) in-plane components of the ensemble-averaged velocity field from period B and (c) a slice at $z=4\pi$ through the UB ( ${Re}=4000$ ) solution shown in figures 7 and 8. Contours in panels (a) and (c) depict streamwise velocity. In-plane velocity components are illustrated using arrows in all the panels. Centres (elliptic stagnation points) of the four vortices of the LB invariant state and their $\pi /2$ -rotated positions around the duct centre are marked by triangles, and compared with positions of vortex centres resulting from averaging the flows over the entire streamwise length and through periods B (purple squares) and C (cyan squares), as shown in figure 14. Colour scale used for $w$ is the same as in figures 13 and 14.

4.3. Streamwise localisation

Focusing on the interludes throughout which the flow features increased symmetrisation, i.e. times when $f_\pi$ drops below $\langle f_{\pi } \rangle -2\sigma$ , we attempt to examine the possible streamwise localisation property of the flow. In what follows, we focus on one of the quiescent interludes, marked B in figures 9 and 10. We note, however, that results remain qualitatively similar at other times (e.g. around $t=2000$ or at the period marked C). We start by considering the following averaging operator:

(4.6) $ \begin{equation} \boldsymbol{u}_s(x,y,z;s) = \frac {1}{(t_1-t_0)} \int ^{t_1}_{t_0} \boldsymbol{u}(x,y,z-st,t)\, \textrm {d}t \end{equation}$

over the fixed time $(t_1-t_0)$ with a constant, streamwise shift per unit time $s$ , where $(t_0,t_1)=(6.4,6.6)\times 10^3$ correspond to the quiescent period B. The meaning of (4.6) corresponds to the time averaging of the flow while continuously moving the averaging domain in the streamwise direction at speed $s$ . In principle, the averaging operation (4.6) highlights flow details that persist under the streamwise shift at speed $s$ , and blurs and attenuates remaining features. This should, in principle, highlight the persisting, coherent structures present in the flow and attenuate remaining fluctuations.

Applying (4.6) over the selected period while varying the speed of the streamwise shift $s$ , we look for the value of $s$ that maximises the extreme of $E_{\perp }$ taken over the streamwise direction $z$ . Plots in figure 16 show variations of the cross-flow energy $E_{\perp }$ for $\boldsymbol{u}_s$ with the streamwise direction for three values of streamwise shift per unit time $s$ . The reader might note that there seems to be a value of $s$ which maximises the extreme value of the cross-flow energy across the streamwise direction and at the same time, and around this particular value of $s$ , i.e. $s=0.397$ as will be shown later, the extreme value of $E_{\perp }$ becomes distinctly larger from neighbouring turbulent backgrounds in $E_{\perp }$ , suggesting a streamwise localisation of the coherent structures propagating downstream. Variation of the extreme value with the shift $s$ is shown in figure 17 and illustrates a well-defined overall maximum, which we evaluate to be approximately $s=0.397$ by bisection. Figure 18 illustrates iso-surfaces of the second invariant of the velocity gradient tensor $Q$ computed for the averaged velocity field $\boldsymbol{u}_s$ obtained for $s=0.397$ and normalised with $\langle w_b \rangle$ of the turbulent state. For comparison purposes among figures 3, 7 and 18, we take the same contour threshold so we can confirm the similar streamwise localisation as well as the in-plane orientation of the vortex structures. In figure 12(b), we can also observe less significant but visible streamwise localisation of the intense vortices even in the flow snapshot taken from the corresponding quiescent period B. Overall, the state resulting from the application of the averaging operation (4.6) for the appropriate value of $s$ suggests that the instantaneous flow structure in the quiescent state of turbulence could be related to the invariant state identified in § 3.

Figure 16. Cross-flow energy density $E_{\perp }$ for $\boldsymbol{u}_s$ along the streamwise direction for three values of the streamwise shift, $s=0.37$ , $0.397$ and $0.43$ . The extreme of $E_\perp$ for $s=0.397$ (thick curve) gives a maximum in figure 16.

Figure 17. Variation of the streamwise maximum of the cross-flow energy $E_{\perp }$ for $\boldsymbol{u}_s$ with the streamwise shift $s$ .

Figure 18. Contours of the second invariant of the velocity gradient tensor $Q$ computed for the vector field resulting from the averaging operator (4.6) for $s=0.397$ taken at $Q/\langle w_b\rangle ^2=1.5\times 10^{-1}$ . Contour level is the same as in figures 3 and 7.

Finally, we would like to stress that the evaluated shift speed $s$ with respect to $\langle w_b \rangle \approx 0.284$ for the turbulent flow is $s\approx 1.397\langle w_{b}\rangle$ , which is strikingly similar to the relation of the propagation speed of the symmetric sub-space edge state on the LB to the corresponding $w_b$ , i.e. $s\approx 1.424w_{b}$ .