2.1 Large-eddy simulation

The LES code has been used for many years to model high- $Re_{\unicode[STIX]{x1D70F}}$ turbulent channel flow (resembling conditions within planetary atmospheric boundary layers, although for the present application we consider high Rossby number and neutral thermal stratification, i.e. turbulence is sustained exclusively by mechanical shear). Appendix A contains a comprehensive description of the LES code, and specific details of the computational mesh for the simulations performed for this study. Thus, we refer the interested reader to appendix A for further clarification on the LES code. Note that quantities resolved during LES are denoted by tilde, i.e. $\tilde{\boldsymbol{u}}=G_{\unicode[STIX]{x1D6E5}}\star \boldsymbol{u}$ , where $G_{\unicode[STIX]{x1D6E5}}$ is a filtering kernel at grid-filter scale, $\unicode[STIX]{x1D6E5}$ , and $\star$ denotes the convolution operator. As per (A 3) and accompanying text, the total (Reynolds) stress tensor, $\unicode[STIX]{x1D64F}$ , is attained via the addition of its constituent resolved and subgrid-scale components. The LES code is used for 24 simulations, which enables parametric testing on the role of element height, $h$ , and spanwise spacing, $s_{2}/\unicode[STIX]{x1D6FF}$ , where the attributes of different cases are summarized in table 1 and figure 1.

Table 1.

Summary of LES case attributes.

Figure 1 shows the geometric arrangements considered here, where the surfaces are composed of streamwise-aligned, vertically truncated, square-based pyramids, as highlighted in figure 1(h). The element heights, $h/\unicode[STIX]{x1D6FF}$ , are summarized in table 1, where the values vary in the range $5\leqslant \unicode[STIX]{x1D6FF}/h\leqslant 20$ ; these values indicate outer-layer similarity may only be exhibited by Cases A5 and A6 (i.e. those most closely resembling a roughness; see table 1 for details) (Flack et al. Reference Flack, Schultz and Connelly2007). Beyond the element ‘steepness’, denoted by $\unicode[STIX]{x1D703}$ in table 1 and figure 1(h), all other element attributes are summarized in § A.1 (i.e. $w_{b}$ , $l_{b}$ , $l_{t}$ , $w_{t}$ and  $s_{1}$ ). As per appendix A, results of computational mesh resolution sensitivity have not been included in this study. This topic has been thoroughly addressed in recent studies, for example Yang & Anderson (Reference Yang and Anderson2017) or, more recently, Awasthi & Anderson (Reference Awasthi and Anderson2018).

In fact, Yang & Anderson (Reference Yang and Anderson2017) used LES to model flow over geometric arrangements similar to those studied for this work, and the flow statistics (including the spatial attributes and coherence of turbulence above the element crests and troughs) exhibited no discernible response to changing computational mesh resolution. Moreover, the intensity and spatial nature of topography-induced, Reynolds-averaged turbulent secondary flows were shown to be insensitive to resolution (at least, over the range of resolutions considered in Yang & Anderson Reference Yang and Anderson2017, which encompass the mesh resolution considered for the present article). Note, also, that Awasthi & Anderson (Reference Awasthi and Anderson2018) present results of LES of flow over surfaces that exhibit geometric resemblance to those considered here – at least, in the sense that they considered large-scale elements with large spanwise spacing – and they presented statistics including spectrograms and correlation profiles between small- and large-scale content via modulation of small-scale frequency and amplitude. Resolution sensitivity was of foremost importance in Awasthi & Anderson (Reference Awasthi and Anderson2018), yet their results exhibited virtually no response to changing resolution. For these reasons, and in the interest of brevity, we consider the role of resolution to have been sufficiently addressed in preceding studies, and we have thus not included such results as part of this article.

Figure 1.

Visualization of LES cases. (a–f) Show visualization in the spanwise–wall-normal plane, while (g) shows visualization in the streamwise–wall-normal plane. (h) Shows perspective visualization of the roughness elements with annotations for geometric attributes, while (i) shows a planform visualization of Case A2, for perspective. In (a–g), line colours correspond to cases A $i$ (black), B $i$ (dark grey), C $i$ (grey) and D $i$ (light grey), where $i=1$ , 2, 3, 4, 5 and 6 (see table 1 for summary of case attributes). (a–f) Show cases with spacing, $s_{2}/\unicode[STIX]{x1D6FF}=6$ , 3, 2, 1, 0.5 and 0.25, respectively. Vertical red line in (a) denotes the location at which the (g) streamwise–wall-normal transect is shown; vertical blue line in (g) denotes the location at which the (a) spanwise–wall-normal transect is shown. In (a–f), the discrete location, $\boldsymbol{x}_{HMP}$ , denotes the location at which time-series data were recorded. $\boldsymbol{x}_{HMP}$ was selected to be precisely above the crest of the row of elements.

Figure 1(a–f) presents visualization of topographic transects in the spanwise–wall-normal plane, where the cases are denoted by different colours, as per the figure caption. As is apparent from the figures, this selection of surfaces provides a means to measure the influence of changing surface element height, $h/\unicode[STIX]{x1D6FF}$ , and spanwise spacing, $s_{2}/\unicode[STIX]{x1D6FF}$ , where both values are summarized in table 1. Figure 1(a–f) provides visual context for the surfaces. As outlined in appendix A, the LES code uses spectral discretization in the horizontal direction, thus rendering the horizontal boundary conditions periodic. The perception of a lateral domain termination is incorrect, and instead one can envision infinite repetition of the domain in the spanwise and streamwise direction. For this reason, the spanwise spacing for Cases A1, B1, C1 and D1, for example, is $s_{2}/\unicode[STIX]{x1D6FF}=6$ , since periodic boundary conditions impose that the nearest adjacent rows are located at $x_{2}/\unicode[STIX]{x1D6FF}=-3$ and $x_{2}/\unicode[STIX]{x1D6FF}=9$ .

Note that Cases A1, B1, C1 and D1 are designed to capture a limiting spacing, $s_{2}/\unicode[STIX]{x1D6FF}=6$ , at which turbulent secondary flows are able to occupy the depth of the flow without any interference from adjacent cells (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Yang & Anderson Reference Yang and Anderson2017). Surfaces in this regime have previously been described as a topography, owing to the global response they impart in the above flow. In contrast, cases with the smallest spacing (A5, A6, B5, B6, C5, C6, D5 and D6; see also figure 1 e,f and table 1) are designed to behave as a roughness, in the sense that any element-induced Reynolds-averaged turbulent secondary flows will be located in the immediate vicinity of the elements, but with spatial extent constrained by the lateral spacing between adjacent rows. Turbulence statistics in the inertial layer, or outer layer, for these cases should exhibit no Reynolds-averaged heterogeneity (beyond those associated with the wall-normal mechanical shear production). For these cases, then, the Reynolds-averaged secondary cells are best understood as those typical of the roughness sublayer – the region of the flow occupied by element-scale vortices (Goldstein & Tuan Reference Goldstein and Tuan1998; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2018).

It has been systematically shown that secondary cells within the roughness regime remain attached to the elements such that upwelling of low-momentum fluid occurs above the individual rows of elements (Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Yang & Anderson Reference Yang and Anderson2017; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2018; Vanderwel et al. Reference Vanderwel, Stroh, Kriegseis, Frohnapfel and Ganapathisubramani2019). In contrast, secondary cells within the topography regime are located in the outer layer such that downwelling of high-momentum fluid occurs above the rows of elements. When downwelling is directed towards relatively more rough regions of the surface, it is consistent with prior interpretation of high- and low-momentum pathways (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015), although we acknowledge that others have found opposite direction sign in complementary research. For this reason, we have added the discrete spatial marker, $\boldsymbol{x}_{HMP}$ , to figure 1(a–f), fully cognizant of the fact that in the roughness regime, the marker may not correspond to a HMP. Nevertheless, for Cases A1, B1, C1 and D1, $\boldsymbol{x}_{HMP}$ indeed corresponds to the location of a HMP, while for the transitional regime cases, we will demonstrate the probability for both to occur.

The selection of asymptotic roughness and topography regime cases ensures that the ensuing parameter variation is guaranteed to encompass the transitional regime, where both element-scale cells and HMP/LMP-type cells are possible (in §§ 3 and 4, we will demonstrate that a truncated Galerkin projection of transitional regime cells yields a system of ordinary differential equations, with orbits that undergo chaotic reversals about a phase-space saddle node). Results from LES of flow over the transitional cases selected for this study – Cases A2–A4, B2–B4, C2–C4 and D2–D4 – provide strong support for chaos as an underlying justification for the ostensible reversal in rotational sign from roughness to topography. The experimental results further this argument (§ 2.2). But, in both cases, the natural complexity of the problem and the presence of broadband mixing confounds efforts to isolate such chaotic mechanisms with clarity.

In figure 1(a), the vertical red line denotes the spanwise location at which the streamwise–wall-normal transect is shown in figure 1(g); likewise, the vertical blue line on figure 1(g) denotes the streamwise location at which the figure 1(a) transect is shown. Figure 1(h) provides further graphical context and perspective on the nature of the roughness elements; figure 1(i) shows a planform visualization of Case A2, for visual perspective.

Figure 2.

Instantaneous (a,c,e) and Reynolds-averaged (b,d,f) flow structure in the spanwise–wall-normal plane. Colour flood contour corresponds to streamwise velocity, $\tilde{u} _{1}/u_{\unicode[STIX]{x1D70F}}$ , while line contours represent positive (red) and negative (blue) components of vertical velocity, $\tilde{u} _{3}/u_{\unicode[STIX]{x1D70F}}$ . (a,b) Case C6; (c,d) Case A3; (e,f) Case D1. In all panels, location of discrete sampling point used for conditional averaging denoted by red circle with black fill. In a subsequent result, probability density functions of velocity at $\boldsymbol{x}_{HMP}$ will be used to demonstrate that $-2\lesssim \tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)\lesssim 4$ (see figure 6 and accompanying text), while figure 3 shows a wall-normal profile for Reynolds-averaged wall-normal velocity for Case D1 (f of the present figure). Note that (a,c,e) are shown over region, $2\leqslant x_{2}/\unicode[STIX]{x1D6FF}\leqslant 4$ , for clarity.

Figure 2 shows instantaneous (panels a,c,e) and Reynolds-averaged (panels b,d,f) contours of streamwise and wall-normal velocity in the spanwise–wall-normal plane. The specific cases are noted in the figure caption. The vertical velocity topographic contours were prepared such that regions of upwelling (positive wall-normal velocity; red) and downwelling (negative wall-normal velocity; blue) are emphasized. This approach aids interpretation of the streamwise velocity contours, since downwelling and upwelling are associated with streamwise momentum excesses and deficits, respectively. As per the caption, figure 3(b) shows wall-normal profiles of Reynolds-averaged vertical velocity for Case D1 (figure 2 f), while figure 6 shows p.d.f.s at $\boldsymbol{x}_{HMP}$ for a range of cases, which provide additional information on the mean and variability of $\tilde{u} _{3}^{\prime }(\boldsymbol{x},t)$ . For now, we note that the limiting negative and positive values of vertical velocity are $\tilde{u} _{3}^{\prime }(\boldsymbol{x},t)\approx -4$ and $\tilde{u} _{3}^{\prime }(\boldsymbol{x},t)\approx 4$ , respectively.

Beginning with the topography cases (figure 2 e,f), the instantaneous visualization shows a large plume of relatively low-momentum fluid immediately above the element, while the Reynolds-averaged contour clearly reveals a large ( $\unicode[STIX]{x1D6FF}$ -scale) zone of streamwise momentum excess above the element (figure 2 f, $2\lesssim x_{2}/\unicode[STIX]{x1D6FF}\lesssim 4$ ), which is concurrent with a large area of downwelling. This flow pattern is a clear manifestation of a HMP – as has been documented in previous studies (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015); note, however, a small zone of upwelling at the top of the element, owing to the presence of element-scale vortices (i.e. those that prevail exclusively for roughness regime surfaces). The large downwelling region (and momentum excess) is flanked by two large-scale zones of upwelling, which are concurrent with streamwise momentum excesses; this result is entirely consistent with HMP/LMP formation observed in prior studies.

Figure 2(a,b) shows instantaneous and Reynolds-averaged flow quantities for a roughness regime case. The instantaneous panel exhibits flow processes resembling those observed in flow over a homogeneous roughness, while the Reynolds-averaged panels show minimal inertial layer disturbances, with the exception of an undulation at $x_{2}/\unicode[STIX]{x1D6FF}\approx 4.8$ and $x_{3}/\unicode[STIX]{x1D6FF}\approx 0.8$ . This is a consequence of some modest transience within the flow that does not vanish with Reynolds averaging; for this study, every effort was made to attain truly Reynolds-averaged quantities. The number of large-eddy turnover times used to compute Reynolds averages is summarized in table 1 in the column with heading, $TU_{0}\unicode[STIX]{x1D6FF}^{-1}$ . We emphasize, also, that for all cases the averaging period was preceded by far longer development time. Nevertheless, such modest undulations are inherent to these simulations, as has been widely reported in prior studies (Reynolds et al. Reference Reynolds, Hayden, Castro and Robins2007; Fishpool et al. Reference Fishpool, Lardeau and Leschziner2009).

Figure 3.

Vertical profiles of Reynolds-averaged streamwise (a) and vertical (b) velocity within a HMP (i.e. vertical transect coincident with $\boldsymbol{x}_{HMP}$ in figure 1 a). Profiles correspond to Case A1 (black), B1 (dark grey), C1 (grey) and D1 (light grey), where vertical shift to beginning of profiles is due to region occupied by solid topographic elements. For context, the canonical channel flow logarithmic streamwise velocity profile, $U_{1}(x_{3})/u_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D705}^{-1}\log (x_{3}/l)$ , and vertical velocity profile, $U_{3}(x_{3})/u_{\unicode[STIX]{x1D70F}}=0$ , are shown with a dashed black line on panels (a) and (b), respectively. Direction of increasing element height added for perspective.

Finally, figure 2(c,d) shows instantaneous and Reynolds-averaged flow patterns for a transitional case. The instantaneous visualization, again, shows large-scale mixing over the depth of the flow, while the Reynolds-averaged streamwise velocity colour flood shows large-scale heterogeneities that are driven by the spanwise surface heterogeneity. Although the structures exhibit a lateral inclination, the vertical velocity contours demonstrate that HMPs are collocated with the rows of elements – that is, regions of downwelling occupy the zone above the element, and this overlaps precisely with streamwise momentum excesses in the outer layer, which is evidence of a HMP. Likewise, regions of net upwelling are present and overlap regions of streamwise velocity deficit, which indicates manifestation of an LMP. This lateral inclination is commonly observed, and is thought to be associated with insufficient time averaging (although, long time-averaging periods were used, as per table 1), or with numerical attributes of the LES code (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013). The pattern does not exhibit the same level of organization as attained for the topography case (i.e. see figure 2 f), but the Reynolds-averaged visualization is consistent with prior results on the position of HMPs and LMPs. Conditional sampling will show, however, that large-scale reversals occur with high probability.

In order to highlight flow structure within a HMP, figure 3(a,b) shows Reynolds-averaged streamwise and wall-normal velocity profiles at a spanwise location intersecting the HMP sampling point in figure 1(a). The profiles are for Cases A1, B1, C1 and D1 – i.e. those for which a distinct HMP resides above the elements – and for context we have included profiles for a canonical channel flow (logarithmic profile in figure 3 a and zero wall-normal velocity in figure 3 b). The origin at which each profile begins is vertically displaced by virtue of the presence of the elements. The immersed boundary method used for the simulations sets the velocity at internal, ‘immersed’ grid points to zero (all details of the LES code and immersed boundary method are provided in appendix A). It is clear, firstly, that the streamwise profiles are not logarithmic – instead, they exhibit a wall-normal gradient ostensibly greater than that for a logarithmic profile, $u_{\unicode[STIX]{x1D70F}}(\unicode[STIX]{x1D705}x_{3})^{-1}$ .

Moreover, the profiles attain a zero-stress state far lower in the domain, for example at $x_{3}/\unicode[STIX]{x1D6FF}\approx 0.5$ for Case A1. This larger wall-normal gradient is entirely consistent with prior work by Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015), who demonstrate that a production–dissipation imbalance in turbulent kinetic energy, $tke$ , above the relatively more rough regions (i.e. at spanwise locations corresponding to HMPs) was the underlying mechanism responsible for sustenance of the secondary flows (Hinze Reference Hinze1967; Bradshaw Reference Bradshaw1987). As would be expected, increasing the element height imposes a greater momentum penalty upon the flow, thereby reducing the centreline velocity: the figure demonstrates a monotonic reduction in outer velocity.

For all cases, figure 3(b) clearly demonstrates a region of downwelling over the depth of the flow (for Case D1, this is also apparent from inspection of figure 2f). However, when the profiles are presented in this format, an apparent reduction in the magnitude of $\langle \tilde{u} _{3}\rangle _{xt}(x_{2,HMP},x_{3})<0$ is apparent, while the magnitude of $\langle \tilde{u} _{3}\rangle _{xt}(x_{2,HMP},x_{3})>0$ in the near-wall region increases. This is a product of local, element-scale roughness sublayer mixing, and is also visible in figure 2(e,f). When these element-scale vortices are produced in close proximity to the surface, ambient background shear associated with the wall serves to homogenize their influence (i.e. Case A1 profiles in figure 3). However, when produced due to the presence of larger elements, the element-scale vortices are naturally larger, thereby inducing higher wall-normal velocities above the crest and weakening the subsequent downwelling.

Figure 4.

Data points for post-normalization shear velocity, $u_{\unicode[STIX]{x1D70F}}=(\text{max}|\unicode[STIX]{x1D61B}_{13}|)^{1/2}$ . Grey panel added to denote transitional regime. Symbols denote group of cases: ‘plus’ sign (Cases A1–A6), asterisk (Cases B1–B6), circle (Cases C1–C6) and square (Cases D1–D6).

For completeness, figure 4 shows the shear velocity attained via $u_{\unicode[STIX]{x1D70F}}=(\text{max}|\unicode[STIX]{x1D61B}_{13}|)^{1/2}$ . For flow over a homogeneous sandpaper distribution, the LES shear velocity is, by definition, $u_{\unicode[STIX]{x1D70F}}=1$ , and these values actually occur in the roughness regime (in an approximate sense). However, with increased spacing, $u_{\unicode[STIX]{x1D70F}}$ increases as large-scale flow heterogeneities emerge in the flow, and local peaks in $\unicode[STIX]{x1D61B}_{13}$ are concentrated above the elements (Awasthi & Anderson Reference Awasthi and Anderson2018). By post-normalizing the statistics with this  $u_{\unicode[STIX]{x1D70F}}$ , the results can be properly compared. In a subsequent result, probability density functions from $\boldsymbol{x}_{HMP}$ are presented, but before doing so some additional results on the underlying driving mechanism are presented.

Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) showed that the turbulent secondary flows associated with HMPs and LMPs – i.e. for which a streamwise momentum excess and downwelling occurs above the relative peaks in surface roughness – are a manifestation of Prandtl’s secondary flow of the second kind (Bradshaw Reference Bradshaw1987). In doing so, Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) used the Reynolds-averaged transport equation for streamwise vorticity and $tke$ . In a later section of this article, the $tke$ transport equation is reviewed again; for this stage, however, the Reynolds-averaged streamwise vorticity transport equation is apropos:

(2.1) $$\begin{eqnarray}\langle \tilde{u} _{2}\rangle _{T}\unicode[STIX]{x2202}_{2}\langle \tilde{\unicode[STIX]{x1D714}}_{1}\rangle _{T}+\langle \tilde{u} _{3}\rangle _{T}\unicode[STIX]{x2202}_{3}\langle \tilde{\unicode[STIX]{x1D714}}_{1}\rangle _{T}=(\unicode[STIX]{x2202}_{3}^{2}-\unicode[STIX]{x2202}_{2}^{2})\langle \unicode[STIX]{x1D61B}_{23}\rangle _{T}+\unicode[STIX]{x2202}_{23}(\langle \unicode[STIX]{x1D61B}_{22}\rangle _{T}-\langle \unicode[STIX]{x1D61B}_{33}\rangle _{T}),\end{eqnarray}$$

where the tensor components, $\langle \unicode[STIX]{x1D61B}_{23}\rangle _{T}$ , $\langle \unicode[STIX]{x1D61B}_{22}\rangle _{T}$ and $\langle \unicode[STIX]{x1D61B}_{33}\rangle _{T}$ , are retrieved a posteriori from $\langle \unicode[STIX]{x1D64F}\rangle _{T}$ via the sum of contributions from the resolved and subgrid-scale components (see also (A 3) and accompanying text). Upon inspection of (2.1), it is immediately apparent that in the absence of Reynolds stress spatial heterogeneity, no secondary flow is possible. Figure 5 shows colour flood contours of Reynolds stress terms appearing in (2.1), where underlying spatial heterogeneity is responsible for the non-zero left-hand side terms.

Figure 5.

Reynolds (turbulent) stresses responsible for Prandtl’s secondary flow of the second kind (see also (2.1) and accompanying text). (a,c,e) and (b,d,f) show results for Cases A5 (roughness regime) and A1 (topography regime), respectively. Colour flood contours correspond to $\unicode[STIX]{x1D61B}_{23}/u_{\unicode[STIX]{x1D70F}}^{2}$ (a,b), $\unicode[STIX]{x1D61B}_{22}/u_{\unicode[STIX]{x1D70F}}^{2}$ (c,d), and $\unicode[STIX]{x1D61B}_{33}/u_{\unicode[STIX]{x1D70F}}^{2}$ (e,f).

Contours are shown for cases in the roughness (figure 5 a,c,e) and topography (figure 5 b,d,f) regime, where case details are summarized in the figure caption. Since $\langle \tilde{\unicode[STIX]{x1D714}}_{1}\rangle _{T}\rightarrow 0$ in the absence of spatial heterogeneity of $\langle \unicode[STIX]{x1D61B}_{23}\rangle _{T}$ , $\langle \unicode[STIX]{x1D61B}_{22}\rangle _{T}$ , and $\langle \unicode[STIX]{x1D61B}_{33}\rangle _{T}$ , the emergence of a channel-like inertial layer in the roughness regime is understood from inspection of figure 5(a,c,e), where all components approach zero for $x_{3}/\unicode[STIX]{x1D6FF}\gtrsim 0.5$ . For the topography regime, regions of high-intensity turbulent stresses are still present in the immediate vicinity of the elements (note that figures 5 a,c,e and 5 b,d,f correspond to the same element height). In addition, however, are the zones of elevated stresses over the depth of the flow, which are the drivers of HMPs and downwelling onto the elements.

For example, for the topography case, the $\unicode[STIX]{x1D61B}_{23}$ contour (figure 5 b) has zones of high magnitude positive and negative $\unicode[STIX]{x1D61B}_{23}$ in the area, $2.8\lesssim x_{2}/\unicode[STIX]{x1D6FF}\leqslant 3$ and $3\leqslant x_{2}/\unicode[STIX]{x1D6FF}\lesssim 3.2$ , respectively, which is virtually indistinguishable from the roughness case (figure 5 a). However, these two zones are, themselves, encapsulated within two $\unicode[STIX]{x1D6FF}$ -scale zones of opposite sign in the region immediately above. Similarly, the contours for $\unicode[STIX]{x1D61B}_{22}$ and $\unicode[STIX]{x1D61B}_{33}$ exhibit near-element distributions identical to those present in the roughness regimes – but, the near-element regions of element stresses are laterally flanked by $\unicode[STIX]{x1D6FF}$ -scale undulations. Again, when viewed in the context of (2.1), the significance of this outer-layer heterogeneity is readily apparent (note that computation of right-hand side terms in (2.1) was the topic of Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015, and was thus not repeated for this study). For brevity, Reynolds stress tensor components are not shown for intermediate cases corresponding to transitional spacing; in preliminary testing performed, it has been found that the secondary zone of Reynolds stresses emerges monotonically with increased spacing. These results suggest that the element-scale vortices are responsible for setting the rotational sign of the outer-layer circulations; however, since the outer-layer circulations exhibit relatively weaker vortical intensity, they are susceptible to perturbations from, for example, meandering LSMs. We will show, and argue, that these disruptions are non-periodic, and a subsequent theoretical exercise will demonstrate that the sign reversal is dynamically chaotic.

Figure 6.

Probability density functions of fluctuating streamwise (a,c,e) and vertical (b,d,f) velocities at discrete location, $\boldsymbol{x}_{HMP}$ , shown in figure 1(a–f) and 2. (a,b) Show results for Cases A1, B1, C1 and D1; (c,d) show results for Cases A3, B3, C3 and D3; and (e,g) show results for Cases A6, B6, C6 and D6. In all panels, line colouring corresponds to A $i$ (black), B $i$ (dark grey), C $i$ (grey) and D $i$ (light grey), where $i=1$  (a,b), 3 (c,d) and 6 (e,f). Vertical lines in (b,d,f) denote thresholds, $\unicode[STIX]{x1D6FC}_{-}$ and $\unicode[STIX]{x1D6FC}_{+}$ , used during conditional sampling of flow. Annotation for direction of increasing element height, $h$ , added to (a), for perspective.

Figure 6 shows probability density functions (p.d.f.) of fluctuating streamwise (a,c,e) and wall-normal (b,d,f) velocity from the discrete sampling point, $\boldsymbol{x}_{HMP}$ , shown in figures 1(a–f) and 2. The p.d.f.s are shown for cases in the topography (a,b), transitional (c,d) and roughness (e,f) regimes, with the direction of increasing element height annotated in (a). For the canonical arrangement of turbulent channel flow over a homogeneous roughness, streamwise and wall-normal velocity skewness at the outer-layer sampling point, $\boldsymbol{x}_{HMP}$ , should be negative and positive, respectively (Raupach et al. Reference Raupach, Antonia and Rajagopalan1991). This pattern is, indeed, observed for the roughness regime p.d.f.s, where the differences between cases are modest. Differences emerge at the relatively low-probability extremes of the p.d.f.s (i.e. the ‘tails’), which should be expected due to perturbations associated with the large elements. Although the p.d.f. tails are, by definition, associated with relatively lower-probability fluctuations, we direct the interested reader to table 1; the right-hand column summarizes the simulation duration via number of ‘large-eddy turnovers’, $TU_{0}\unicode[STIX]{x1D6FF}^{-1}\sim O(10^{3})$ . This integration time is sufficient for well resolved p.d.f.s, and stationary conditional averages, even when considering the relatively smaller subset of instants used to capture conditional averages.

In the transition and topography regime, the p.d.f.s are strikingly different: while the streamwise and wall-normal velocity fluctuations remain negatively and positive skewed, respectively, the elements now impart a dramatic influence on the velocity fluctuations.

The figure 6 p.d.f.s demonstrate the potential for system reversals – that is, events with opposite directional sign, which in this case implies the interchange between upwelling of low-momentum fluid or downwelling of high-momentum fluid. With this, and since downwelling and upwelling is a defining marker on the presence of a HMP or LMP, respectively (see figure 2 b,d,f and accompanying text), we have used the figure 6 p.d.f.s of vertical velocity to introduce thresholds for use during conditional sampling of the flow. The $\boldsymbol{x}_{HMP}$ marker noted on figures 1 and 2 was selected since, at this location, secondary flows associated with the roughness sublayer will have attenuated; instead, the only secondary flows present will be those related to outer-layer secondary flows. This article is concerned, specifically, with non-periodic reversals of the outer-layer secondary flows, which are driven by spatial heterogeneity of the Reynolds stresses (equation (2.1), figure 5, and accompanying text), but prone to sudden reversals. This approach will, thus, provide an aggregate view of the flow structure present during times of pronounced upwelling ( $\tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)>\unicode[STIX]{x1D6FC}_{-}$ ) or downwelling ( $\tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)<\unicode[STIX]{x1D6FC}_{-}$ ).

Conditional averaging is a popular technique for characterizing the spatial attributes of turbulence (Antonia Reference Antonia1981; Kim, Moin & Moser Reference Kim, Moin and Moser1987; Adrian Reference Adrian1988; Finnigan, Shaw & Patton Reference Finnigan, Shaw and Patton2009; Sheng, Malkiel & Katz Reference Sheng, Malkiel and Katz2009; Hutchins et al. Reference Hutchins, Monty, Ganapathisubramani, Ng and Marusic2011; Richter & Sullivan Reference Richter and Sullivan2014; Jacob & Anderson Reference Jacob and Anderson2016). Unlike a Reynolds average, wherein the resultant quantity is based on an average over every state, a conditionally averaged quantity is based on a sub-set of all states and is predicated upon some pre-defined condition. The conditionally averaged velocity fields, $\widehat{\tilde{\boldsymbol{u}}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ , are attained with the operations:

(2.2) $$\begin{eqnarray}{\displaystyle \frac{\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{-}}(\boldsymbol{x})}{u_{\unicode[STIX]{x1D70F}}}}=\left\langle \left.{\displaystyle \frac{\tilde{\boldsymbol{u}}(\boldsymbol{x},t)}{u_{\unicode[STIX]{x1D70F}}}}\right|\tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)<\unicode[STIX]{x1D6FC}_{-}\right\rangle _{N_{\unicode[STIX]{x1D6FC}_{-}}},\end{eqnarray}$$

and

(2.3) $$\begin{eqnarray}{\displaystyle \frac{\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{+}}(\boldsymbol{x})}{u_{\unicode[STIX]{x1D70F}}}}=\left\langle \left.{\displaystyle \frac{\tilde{\boldsymbol{u}}(\boldsymbol{x},t)}{u_{\unicode[STIX]{x1D70F}}}}\right|\tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)>\unicode[STIX]{x1D6FC}_{+}\right\rangle _{N_{\unicode[STIX]{x1D6FC}_{+}}},\end{eqnarray}$$

where $\tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)$ is the sampling quantity, $\unicode[STIX]{x1D6FC}_{-}$ and $\unicode[STIX]{x1D6FC}_{+}$ are pre-defined thresholds and $N_{\unicode[STIX]{x1D6FC}_{-}}$ and $N_{\unicode[STIX]{x1D6FC}_{+}}$ are integer number of instances at which the threshold is exceeded. Sample p.d.f.s of $\tilde{u} _{3}^{\prime }(\boldsymbol{x}_{HMP},t)$ for table 1 cases, and sampling thresholds, are shown in figure 6(b,d,e). For turbulence characterization, conditional sampling is appealing since turbulence structure associated with high-magnitude, low-probability events can be revealed with relative clarity (Jacob & Anderson Reference Jacob and Anderson2016).

Figure 7.

Conditionally averaged flow structure in the spanwise–wall-normal plane, based on figure 6 p.d.f.s. Colour flood contour corresponds to streamwise velocity, $\hat{\tilde{u} }_{1}/u_{\unicode[STIX]{x1D70F}}$ , while line contours represent positive (red) and negative (blue) components of vertical velocity, $\hat{\tilde{u} }_{3}/u_{\unicode[STIX]{x1D70F}}$ . (a,c,e) Correspond to flow sampling based upon $\tilde{u} _{3}(\boldsymbol{x}_{HMP},t)/u_{\unicode[STIX]{x1D70F}}\leqslant \unicode[STIX]{x1D6FC}_{-}$  (2.2), while (b,d,e) correspond to flow sampling based upon $\tilde{u} _{3}(\boldsymbol{x}_{HMP},t)/u_{\unicode[STIX]{x1D70F}}\geqslant \unicode[STIX]{x1D6FC}_{+}$  (2.3), where both thresholds are shown in figure 6(b,d,f). (a,b) Case C6; (c,d) Case A3; (d,e) Case D1.

Figure 7 shows conditionally averaged contours of streamwise velocity (colour flood), with line contours for wall-normal velocity superimposed; the wall-normal velocity contours are divided into positive (red) and negative (blue). The figure format is identical to that used for figure 2. The figures are arranged such that (a,b), (c,d) and (e,f) correspond to roughness, transitional and topography regime surfaces, respectively, while (a,c,e) and (b,d,f) show statistics for $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{-}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ and $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{+}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ , respectively.

In the roughness regime, figure 7(a,b) shows local zones of downwelling and upwelling, respectively, surrounding the sampling point, beyond which the distribution rapidly approaches a channel-like distribution above a homogeneous rough surface (i.e. sandpaper). In this sense, the upwelling and downwelling are not indications of an LMP or HMP, but rather due to the passage of coherent motions within the flow: since the flow is responding to a uniform distribution of roughness elements, roughness-driven secondary flows in the outer layer cannot form, which differs intrinsically to large-spacing cases that exhibit domain-scale Reynolds-averaged secondary cells in the outer layer (see also figure 2). As such, any domain-scale patterns of coherence for these roughness cases can only be attributed to the passage of coherent motions (Yang & Anderson Reference Yang and Anderson2017). For the large element roughness cases with $\unicode[STIX]{x1D6FF}/h=5$ , 6.7 or 10 (see table 1), one hesitates to attribute these to large-scale motions, since such outer-layer coherence attenuates in the presence of large elements. However, for the small element roughness cases (i.e. A5 and A6), outer layer similarity is possible and the structures could be due to the passage of large-scale motions encapsulated by hairpins. As spacing increases, however, the picture changes entirely.

Figure 7(e,f) shows $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{-}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ and $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{+}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ , respectively, which confirm that large-scale system reversals are possible. Figure 7(e), in fact, resembles the Reynolds-averaged distribution (figure 2 f), with a zone of downwelling above the element overlapping a region of streamwise momentum excess (HMP), and flanked by adjacent zones of upwelling of low-momentum fluid (i.e. an LMP). However, when (2.3) is used to compute $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{+}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ , a distinctly different pattern is found, with upwelling overlapping a $\unicode[STIX]{x1D6FF}$ -scale streamwise momentum deficit, adjacent to two regions of downwelling. Note, however, that even for $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{+}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ , the two adjacent zones of upwelling (i.e. $0.2\lesssim x_{2}/\unicode[STIX]{x1D6FF}\lesssim 1.8$ and $4.2\lesssim x_{2}/\unicode[STIX]{x1D6FF}\lesssim 5.8$ ) remain – with structure very similar to that observed in the $\widehat{\tilde{\boldsymbol{u}}}_{\unicode[STIX]{x1D6FC}_{-}}(\boldsymbol{x})/u_{\unicode[STIX]{x1D70F}}$ contours. Figure 7(e,f) demonstrates that during vertical velocity excursions of relevant probability, the system can exhibit large-scale reversals from the ‘standard’ LMP/HMP distributions (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015) to element-driven upwelling (Goldstein & Tuan Reference Goldstein and Tuan1998; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2018). The propensity for this to occur is especially relevant in the transitional regime, where the presence of both mechanisms – element-driven upwelling, and outer-layer-driven downwelling – occurs in close proximity (figure 7 c,d).

In this section, we have demonstrated that although a Reynolds-averaged direction sign emerges within the flow, large-scale reversals are possible in the roughness-driven secondary flow. This has been demonstrated in the p.d.f.s (figure 6) and conditionally averaged turbulence statistics (figure 7). We argue that these system reversals occur in a non-periodic manner; that is, when the system is projected in its underlying phase space, it does not evolve to a fixed point or limit cycle, but instead evolves chaotically; in complementary prior work, others have seen preliminary evidence suggesting this also (Vanderwel et al. Reference Vanderwel, Stroh, Kriegseis, Frohnapfel and Ganapathisubramani2019). When the secondary cells are projected upon a spatially harmonic base flow and substituted into the governing equations, we recover a system resembling the Lorenz system (Lorenz Reference Lorenz1963). In order to provide further support for this argument, however, complementary results from experimental measurement of turbulent boundary layer flow over a complex roughness with a predominant streamwise elongation are shown. Although the experimental results are derived from flow over an intrinsically different rough surfaces, we report conceptual agreement with the LES results.

2.2 Wind-tunnel experiments

In this section, we show results of experimental measurement of turbulent boundary layer flow over a complex, multiscale roughness, and a canonical smooth wall. The roughness is derived from a severely damaged gas turbine blade digital elevation model (Bons et al. Reference Bons, Taylor, McClain and Rivir2001). The roughness has been considered in a series of articles, from efforts to establish the limits of outer-layer similarity (Wu & Christensen Reference Wu and Christensen2007, Reference Wu and Christensen2010) to, more recently, the nature of inner–outer interactions in rough-wall boundary layers (Pathikonda & Christensen Reference Pathikonda and Christensen2017). Table 2 contains salient details of the flow conditions at the measurement location, at which the boundary layer had attained self-similar conditions. In table 2, the ‘Rough-HMP’ and ‘Rough-LMP’ labels correspond to the spanwise locations denoted in figure 8(a,b) by thick black lines with the ‘HMP’ and ‘LMP’ annotation, respectively. In addition to the brief discussion of the measurement protocol described here, note that appendix B contains a comprehensive description of the flow facility. As per table 2, the high $Re_{\unicode[STIX]{x1D70F}}$ values for the smooth- and rough-wall boundary layers enable direct comparison with the LES datasets under dynamically similar conditions (Jimenez Reference Jimenez2004).

Figure 8.

Results from experimental measurement of turbulent boundary layer flow over complex roughness (Bons et al. Reference Bons, Taylor, McClain and Rivir2001) with predominant streamwise elongation (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2010). (a) Shows a colour flood contour of Reynolds-averaged streamwise velocity, with Reynolds-averaged spanwise–wall-normal velocity vectors superimposed; data retrieved from a streamwise location at which the boundary layer has attained self-similarity, and the profile beneath the contour denotes the roughness element height at the same location. (b) Shows statistics of low-pass filtered (at scale,  $2\unicode[STIX]{x1D6FF}$ , and denoted by $\overline{(\ldots )}$ ) streamwise average of topography, $\overline{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6FF}$ (solid black line), and spanwise gradient of topographic profile, $\unicode[STIX]{x2202}_{2}\overline{\unicode[STIX]{x1D702}}$ (red line). In all panels, locations of high- and low-momentum pathways are denoted, while spanwise locations corresponding to measurement locations are denoted by a thick black line. Locations at which time-series measurements retrieved within HMP and LMP denoted with white arrows. Results also presented in Barros & Christensen (Reference Barros and Christensen2014).

Table 2.

Summary of wind-tunnel experiment attributes.

The complex roughness is multiscale, but exhibits underlying features that enable comparison against the canonical LES cases (§ 2.1). Mejia-Alvarez & Christensen (Reference Mejia-Alvarez and Christensen2010) present p.d.f.s of the spatial gradient of the figure 8(a) surface, $p(\unicode[STIX]{x2202}_{1}\unicode[STIX]{x1D702}(x_{1},x_{2}))$ and $p(\unicode[STIX]{x2202}_{2}\unicode[STIX]{x1D702}(x_{1},x_{2}))$ , where $\unicode[STIX]{x1D702}(x_{2},x_{3})$ is the surface height distribution; the kurtosis of $p(\unicode[STIX]{x2202}_{2}\unicode[STIX]{x1D702}(x_{1},x_{2}))$ far exceeded $p(\unicode[STIX]{x2202}_{1}\unicode[STIX]{x1D702}(x_{1},x_{2}))$ , an indication of the relative concentration of streamwise gradients about the mean (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2010, their figure 2c,d). In fact, this can be further understood when interpreted in the context of spectral density of height (Parseval’s identity). Recall that the variance of the gradient (surface slope) is derived from integration of its spectral density, $(\unicode[STIX]{x1D70E}^{\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D702}})^{2}=\int _{2\unicode[STIX]{x03C0}/L}^{2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D702}_{h}}E_{s}(k)\,\text{d}k$ , where $\unicode[STIX]{x1D702}_{h}$ and $L$ denote the cutoff wavelength and largest energy-containing wavelength of $\unicode[STIX]{x1D702}(x_{1},x_{2})$ , respectively, $k=\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}$ and $E_{s}(k)$ is spectral density of $\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}(x_{1},x_{2})$ . Since $E_{s}(k)\sim k^{2}E_{h}(k)$ , where $E_{h}(k)$ is the spectral density of height, it follows that $\unicode[STIX]{x1D70E}^{\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D702}}\sim (k^{2}E_{h}(k))^{1/2}$ . If $\unicode[STIX]{x1D702}$ is a self-similar distribution, for which the spectral density of height follows a power law, $E_{h}(k)\sim k^{\unicode[STIX]{x1D6FD}}$ , where $\unicode[STIX]{x1D6FD}$ is the spectral exponent, $\unicode[STIX]{x1D70E}^{\unicode[STIX]{x1D6FB}\unicode[STIX]{x1D702}}\sim k^{(1/2)(2+\unicode[STIX]{x1D6FD}_{s})}$ . Thus, that the p.d.f.s of the streamwise and spanwise gradients of the complex roughness are narrower and wider, respectively, demonstrates a relatively steeper and flatter spectral distribution, respectively.

The figure 8 complex roughness has been compared against conceptually similar results from LES on a number of previous occasions, where in all cases favourable agreement has been found (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). Figure 8(a), bottom, shows a spanwise transect of the roughness height recorded at the flow measurement location. Figure 8(a), top, shows a contour of Reynolds-averaged streamwise velocity at a streamwise location at which the boundary layer has attained self-similarity. In order to provide additional detail, vectors for spanwise and wall-normal velocity have been superimposed upon the colour flood (figure originally presented in Barros & Christensen Reference Barros and Christensen2014, and included here with the permission of the authors). In the absence of spanwise surface gradients, the spanwise undulations in Reynolds-averaged velocity vanish, leaving only the wall-normal gradient (and, of course, the relatively smaller streamwise gradient). The Reynolds-averaged flow is clearly modulated by the surface characteristics, although establishing causality is not readily possible from figure 8(a) (or the height transect shown below the panel). In order to clearly demonstrate how the flow responds to salient characteristics of the surface, Barros & Christensen (Reference Barros and Christensen2014) developed the figure 8(b) surface statistics of $\unicode[STIX]{x1D702}$ following low-pass filtering at filter scale, $2\unicode[STIX]{x1D6FF}$ , which is denoted by $\overline{(\ldots )}$ .

The figure 8(b) profiles offer immediate insight into how the surface modulates the Reynolds-averaged flow. Considering firstly $\overline{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6FF}$ , where positive and negative values represent relatively ‘more rough’ and ‘less rough’ regions of the surface, respectively, it is clear that large-scale downwelling (and streamwise momentum excesses) occurs above $\overline{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6FF}>0$ . In contrast, large-scale upwelling (and streamwise momentum deficits) is positioned above regions of $\overline{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D6FF}<0$ – the ‘low drag’ regions. Locations of upwelling and downwelling are categorically LMPs and HMPs, respectively, and the self-organization of HMPs above the ‘high drag’ regions has been established now in a series of preceding articles (Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; see also figure 2 f and accompanying discussion). Moreover, in-depth analysis of the turbulent stress distribution for the complex roughness case demonstrated conceptual agreement with the canonical LES cases (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). Nevertheless, in this section we will show evidence that within the HMP and LMP – in this case, those denoted by heavy black line in figure 8(a) – large-scale system reversals from the established Reynolds-averaged flow patterns are possible.

Figure 9 shows p.d.f.s of streamwise velocity fluctuations about the Reynolds-averaged velocity within the LMP and HMP:

(2.4) $$\begin{eqnarray}u_{1}^{\prime }(\boldsymbol{x},t)=u_{1}(\boldsymbol{x},t)-{\textstyle \frac{1}{2}}(\langle u_{1}\rangle _{T}(x_{1},x_{2,LMP},x_{3})+\langle u_{1}\rangle _{T}(x_{1},x_{2,HMP},x_{3})).\end{eqnarray}$$

The spanwise locations at which the time-series measurements were performed (shown in figure 8) were selected for two reasons: (i) they capture a distinct HMP and LMP; and (ii) they were sufficiently far from the tunnel side walls. The long time-series datasets were limited to a single streamwise location, two fixed spanwise locations and a series of discrete wall-normal locations (denoted by white arrows in figure 8 a,b). For this reason, subtraction of the right-hand side average in (2.4) is appealing since within a HMP and LMP, $\langle u_{1}^{\prime }\rangle _{T}(x_{1},x_{2,HMP},x_{3})>0$ and $\langle u_{1}^{\prime }\rangle _{T}(x_{1},x_{2,LMP},x_{3},t)<0$ , respectively, and thus large-scale excursions from the Reynolds-averaged flow patterns can be more readily discerned.

Figure 9.

Probability density functions of fluctuating streamwise velocity at elevation $x_{3}/\unicode[STIX]{x1D6FF}=0.1$  (a), 0.2 (b) and 0.355 (c), where fluctuations are measured relative to the average velocity within the HMP and LMP (2.4) for the rough-wall cases, while fluctuations are based upon deviations about the Reynolds average for the smooth-wall case. Measurement location graphically denoted in figure 8(a). In the panels, black, grey and light grey lines correspond to measurement within turbulent boundary layer flow over a smooth wall, within a HMP and within an LMP, respectively (see table 2 for experiment details; see figure 8 a,b for visualization of measurement location).

To provide more clarity for this argument, and to cast the rough-wall results within a broader context, figure 9 shows p.d.f.s of fluctuating streamwise velocity (via (2.4)) at three outer-scaled wall-normal positions (values shown in the figure caption; profile colour coding is explained in the caption also). Salient details of the smooth-wall experiment are summarized in table 2 and appendix B. Since (2.4) is used to compute fluctuations for the rough-wall cases, the averages are shifted laterally from the true average, as expected.

The p.d.f.s further emphasize that within the location of a HMP, there is a non-negligible probability of streamwise velocity values resembling those within an LMP, and vice versa. This is true at the three sampling locations considered, all of which reside within a clearly defined HMP and LMP (see figure 8 b). In fact, the HMP and LMP p.d.f.s are negatively and positively skewed, respectively, indicative of the potential for occurrence of these reversal values. The passage of LSMs through a HMP or LMP may be capable of inducing a system reversal – such that locations of upwelling and downwelling reverse – although, in truth, to establish such causality is beyond the scope of this study. For the present section, we only wish to demonstrate that for flow over the complex roughness seen in figure 8, the Reynolds-average HMP and LMP arrangement is not stable for all time. By virtue of (2.4), the time-series measurements are able to highlight times when the streamwise velocity within a HMP (in a Reynolds-averaged sense) is comparable to the value within an LMP (again, in a Reynolds-averaged sense). This, by itself, provides support for the existence of large-scale (non-periodic) reversals to the Reynolds-averaged flow pattern.

Inclusion of the smooth-wall p.d.f. helps to indicate how an LMP and HMP alters turbulence properties (at least, at the discrete locations shown in figure 8 b). In the near-wall region (figure 9 a), the p.d.f.s are relatively similar (again note that the lateral shifting of the rough-wall p.d.f.s is a product of (2.4)), though not equivalent. The sampling point is ‘close’ to the surface, such that roughness sublayer processes could modulate the tails of the p.d.f. However, at progressively higher wall-normal locations (figure 9 b,c), the rough-wall p.d.f.s change materially. The HMP p.d.f. exhibits strong negative skewness (figure 9 b), and then is almost Gaussian (figure 9 c), while the LMP p.d.f. is less sensitive to wall-normal elevation. In a previous article, Awasthi & Anderson (Reference Awasthi and Anderson2018) found that turbulence structure within a HMP was profoundly modified, while turbulence structure within an LMP is relatively unaffected. The p.d.f.s here, shown relative to a baseline outer-scaled smooth-wall turbulent boundary layer, offer further support for this.

Since the time-series measurements within the LMP and HMP were not simultaneously recorded, the variability in $u_{1}^{\prime }(x_{1},x_{2,HMP},x_{3},t)$ and $u_{1}^{\prime }(x_{1},x_{2,LMP},x_{3},t)$ cannot be correlated. Thus, the existence of ‘cross-overs’ – transitions in the system such that downwelling and upwelling locations reverse – cannot be directly demonstrated. We once again note the results from Kevin and company, who showed the instantaneous propensity for depth-scale reversals from the Reynolds-averaged secondary flow patterns (Kevin et al. Reference Monty, Bai, Pathikonda, Nugroho, Barros, Christensen and Hutchins2017, Reference Monty and Hutchins2019). They considered a herringbone surface for which the spanwise variability in surface drag is less obvious – at least relative to the cases considered herein – but nevertheless their results provide further support for the arguments derived from the present study. However, interpretation of these results within the broader context of the LES results certainly indicates that such reversals are possible. The experimental time series and p.d.f.s have shown how the streamwise velocity within a HMP and LMP vary in time; conditional sampling of the LES results – based on thresholds identified from p.d.f.s – shows the presence of reversals.

2.3 Chaos in roughness-driven secondary flows

The LES and experimental results provided in §§ 2.1 and 2.2, respectively, have demonstrated that the spatial patterns associated with HMPs and LMPs within the Reynolds-averaged flow can be disrupted and potentially reversed. The LES results were based on flow over a series of surfaces encompassing changes in spanwise spacing between adjacent rows of elements, and element height. The Reynolds-averaged results confirm downwelling (i.e. a HMP) above the elements, but conditional sampling of the flow showed that upwelling (and a streamwise momentum deficit) can occur too. Experimental results also showed conceptually similar processes, although owing to the inherent challenges of flow access, etc., a direct repetition of the LES post-processing analysis was not possible. Nevertheless, p.d.f.s show wide variability within the LMP and HMP, and flow states different to those exhibited within the logarithmic region of a smooth-wall turbulent boundary layer. Recent work by Kevin and company has shown that LSMs meandering through the domain can disrupt any outer-layer circulations due to spanwise surface heterogeneity (Kevin et al. Reference Monty, Bai, Pathikonda, Nugroho, Barros, Christensen and Hutchins2017, Reference Monty and Hutchins2019). In contrast, others have explored how LSMs are disrupted by secondary flows due to spanwise-surface heterogeneity (Pathikonda & Christensen Reference Pathikonda and Christensen2017; Awasthi & Anderson Reference Awasthi and Anderson2018). An interaction is undeniable.

In preceding articles by the author and company, it has been shown that the HMP and LMP are positioned such that there is a predominant downwelling onto the ‘more rough’ regions of the flow (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2010; Willingham et al. Reference Willingham, Anderson, Christensen and Barros2013; Barros & Christensen Reference Barros and Christensen2014; Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015; Yang & Anderson Reference Yang and Anderson2017). Others, however, have reported the precise opposite, where we note specifically results from Ganapathisubramani and company, who have found upwelling onto relative roughness peaks (Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Medjnoun et al. Reference Medjnoun, Vanderwel and Ganapathisubramani2018; Meyers, Ganapathisubramani & Cal Reference Meyers, Ganapathisubramani and Cal2019; Vanderwel et al. Reference Vanderwel, Stroh, Kriegseis, Frohnapfel and Ganapathisubramani2019). Results presented herein are not used to reconcile these ostensibly conflicting results; instead, we have leveraged the widely accepted result that the Reynolds-averaged flow patterns are not stable for all time. Numerical and experimental results from groups studying this problem have, indeed, all shown a propensity for system reversals. In the following section, a projection for the secondary flow – which is, itself, superimposed upon a background flow – is evaluated within the context of the governing transport equations. A series of algebraic developments yields a system of ordinary differential equations governing the evolution of time-dependent parameters used to define the secondary flow. The system resembles the Lorenz system (Saltzman Reference Saltzman1962; Lorenz Reference Lorenz1963), and time integration reveals chaotic, non-periodic orbits in phase space. This result, then, provides evidence that the reversals observed in simulations and experiments occur in a chaotic, non-periodic manner.

3.1 Reduction of unknowns and models

In order to further advance the development, the number of unknowns must be reduced. A streamfunction defining flow components in the spanwise–wall-normal plane, $\unicode[STIX]{x1D713}(x_{2},x_{3},t)$ , will be introduced, thereby simplifying (3.12) to (3.15), since $u_{2}=\unicode[STIX]{x2202}_{3}\unicode[STIX]{x1D713}$ , $u_{3}=-\unicode[STIX]{x2202}_{2}\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D714}_{1}=-\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}$ (where it is understood that the streamwise component of the Laplacian operator vanishes due to the streamwise-homogeneity condition). Substitution into (3.12) and (3.15) yields:

(3.16) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x2202}_{i}^{2}\unicode[STIX]{x1D713})+J(\unicode[STIX]{x2202}_{i}^{2}\unicode[STIX]{x1D713},\unicode[STIX]{x1D713})=(\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{t})\unicode[STIX]{x2202}_{i}^{2}(\unicode[STIX]{x2202}_{i}^{2}\unicode[STIX]{x1D713})+U^{\unicode[STIX]{x1D714}}J(u_{1},g),\end{eqnarray}$$

and

(3.17) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}u_{1}+J(u_{1},\unicode[STIX]{x1D713})=-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}\unicode[STIX]{x2202}_{1}p+(\unicode[STIX]{x1D708}+\unicode[STIX]{x1D708}_{t})\unicode[STIX]{x2202}_{i}^{2}u_{1}-\unicode[STIX]{x1D6F1}.\end{eqnarray}$$

The number of unknowns – $u_{1}$ , $\unicode[STIX]{x1D713}$ and $p$ – is substantially diminished. To further the development, it is firstly helpful to normalize the system of equations by salient parameters of the problem – shear velocity, $u_{\unicode[STIX]{x1D70F}}$ , flow depth, $\unicode[STIX]{x1D6FF}$ , and vortex forcing magnitude, $U^{\unicode[STIX]{x1D714}}$ . Herein, superscript  $+$ will denote a normalized variable (although, for convenience, this superscript will be omitted in all subsequent steps for brevity). To summarize, the normalization is attained with:

(3.18a-g ) $$\begin{eqnarray}t=Hu_{\unicode[STIX]{x1D70F}}^{-1}t^{+}\quad \!u_{1}=u_{\unicode[STIX]{x1D70F}}u_{1}^{+}\quad \!\unicode[STIX]{x1D713}=u_{\unicode[STIX]{x1D70F}}H\unicode[STIX]{x1D713}^{+}\quad \!\boldsymbol{x}=\unicode[STIX]{x1D6FF}\boldsymbol{x}^{+}\quad \!u_{1}^{\unicode[STIX]{x1D714}}=U^{\unicode[STIX]{x1D714}}g\quad \!p=\unicode[STIX]{x1D70C}u_{\unicode[STIX]{x1D70F}}^{2}p^{+}\quad \!\unicode[STIX]{x1D6F1}=u_{\unicode[STIX]{x1D70F}}^{2}\unicode[STIX]{x1D6FF}^{-1}\unicode[STIX]{x1D6F1}^{+}.\end{eqnarray}$$

Upon substitution, the resultant system yields:

(3.19) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}^{+}(\unicode[STIX]{x2202}_{i}^{+,2}\unicode[STIX]{x1D713}^{+})+J^{+}(\unicode[STIX]{x2202}_{i}^{+,2}\unicode[STIX]{x1D713}^{+},\unicode[STIX]{x1D713}^{+})=({\mathcal{R}}^{-1}+{\mathcal{R}}_{t}^{-1})\unicode[STIX]{x2202}_{i}^{+,2}(\unicode[STIX]{x2202}_{j}^{+,2}\unicode[STIX]{x1D713}^{+})+{\mathcal{L}}J^{+}(g,u_{1}^{+}),\end{eqnarray}$$

and

(3.20) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}^{+}u_{1}^{+}+J^{+}(u_{1}^{+},\unicode[STIX]{x1D713}^{+})=-\unicode[STIX]{x2202}_{1}^{+}p^{+}+({\mathcal{R}}^{-1}+{\mathcal{R}}_{t}^{-1})\unicode[STIX]{x2202}_{i}^{+,2}u_{1}^{+}-\unicode[STIX]{x1D6F1}^{+}.\end{eqnarray}$$

To clarify, the pressure-gradient forcing term within (3.20) goes to unity following normalization via the linear momentum balance conditions for a channel, viz. $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{-1}u_{\unicode[STIX]{x1D70F}}^{-2}\,\text{d}P_{\infty }/\text{d}x=\unicode[STIX]{x1D70F}^{w}\unicode[STIX]{x1D70C}^{-1}u_{\unicode[STIX]{x1D70F}}^{-2}=\unicode[STIX]{x1D6F1}^{+}=1$ . Within (3.19) and (3.20), we emphasize the emergence of system parameters ${\mathcal{R}}_{\unicode[STIX]{x1D708}}=u_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}^{-1}$ , ${\mathcal{R}}_{t}=u_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}_{t}^{-1}$ and ${\mathcal{L}}=U^{\unicode[STIX]{x1D714}}/u_{\unicode[STIX]{x1D70F}}$ , corresponding to the viscous Reynolds number, turbulent Reynolds number and a parameter loosely resembling the turbulent Langmuir number (Leibovich Reference Leibovich1983). In many engineering and geophysical flows of interest – for which the turbulent flow physics are completely dominated by inertia – it is commonplace to set ${\mathcal{R}}_{t}^{-1}({\mathcal{R}}_{t}{\mathcal{R}}_{\unicode[STIX]{x1D708}}^{-1}+1)\rightarrow {\mathcal{R}}_{t}^{-1}$ , although in the interest of generality we retain all terms for now and further simplify by setting ${\mathcal{R}}^{-1}={\mathcal{R}}_{\unicode[STIX]{x1D708}}^{-1}+{\mathcal{R}}_{t}^{-1}$ . Since all variables are now non-dimensional, the superscript  $+$ will henceforth be omitted.

Ordinarily, unknown pressure would be assessed dynamically during direct or large-eddy simulation via computation of the divergence of (3.20) and elimination of terms with $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}$ . This procedure yields a Poisson equation for pressure, for which the gradients are added to (3.20) as a pressure correction that serves to preserve a solenoidal vector, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ . However, in this work we will consider streamwise-homogeneous models for $u_{1}$ and  $\unicode[STIX]{x1D713}$ , automatically guaranteeing conservation of mass since:

(3.21) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=\unicode[STIX]{x2202}_{1}u_{1}+\unicode[STIX]{x2202}_{23}\unicode[STIX]{x1D713}-\unicode[STIX]{x2202}_{23}\unicode[STIX]{x1D713}=0.\end{eqnarray}$$

Under this condition, the streamwise pressure correction within (3.20) vanishes, further simplifying the system of equations and reducing the number of unknowns. This is somewhat similar to Waleffe (Reference Waleffe2001, Reference Waleffe2003), who implemented a ‘roll–streak’ or poloidal–toroidal projection of velocity within the governing momentum transport equations, thereby eliminating the pressure correction term. The need for such a projection is obviated herein by virtue of the streamwise-homogeneity condition.

Following elimination of the pressure correction, equations (3.19) and (3.20) can be solved, but direct solution is complicated by the need to compute partial derivatives in  $t$ , $x_{2}$ and  $x_{3}$ . For the Reynolds-averaged version of (3.19) and (3.20), one could quickly envision possible solutions via introduction of an appropriate similarity variable (Anderson et al. Reference Anderson, Yang, Shrestha and Awasthi2018), but this would not capture the transient nature of domain-scale rolls in the transitional regime. As an alternative, we will instead proceed with a severely truncated projection for the large-scale secondary rolls, which will yield a prognostic description for evolution of the projection coefficients.

4.1 Boundary conditions

Within the domain, $\boldsymbol{x}\in \unicode[STIX]{x1D6FA}$ , limiting conditions at the domain surface, $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ , are automatically defined through prescription of the models for $u_{1}$ and  $\unicode[STIX]{x1D713}$ . At the wall, $x_{3}=0$ , the no-slip condition imposes $u_{1}(x_{1},x_{2},x_{3}=0)=u_{2}(x_{1},x_{2},x_{3}=0)=u_{3}(x_{1},x_{2},x_{3}=0)=0$ , and quick inspection of the input functions ((4.1) and (4.2)) demonstrates that the model functions comply with this condition. At the domain top, $x_{3}=1$ , we can deduce that $u_{1}(x_{2},x_{3}=1)\rightarrow \langle a_{3}\rangle _{T}$ , where $\langle a_{3}\rangle _{T}$ is a centreline velocity superimposed upon the ambient channel profile. This could be defined by various means – from direct a priori prescription (for modelling purposes) to prescription based upon the existence of a plane-averaged profile, in which $\langle a_{3}\rangle _{T}$ is presumed to be some factor of the centreline velocity, $\langle a_{3}\rangle _{T}=\unicode[STIX]{x1D701}U_{0}/u_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D701}\unicode[STIX]{x1D705}^{-1}\log (l^{-1})$ , where $\unicode[STIX]{x1D705}$ is the von Kármán constant and $l$ is a roughness length defined a priori. Note that this can be further advanced by enforcing the existence of equilibrium conditions, in which $u_{\unicode[STIX]{x1D70F}}$ can be defined in terms of the ambient channel forcing and $U_{0}/(\unicode[STIX]{x1D70F}^{w}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6F1})^{1/2}=\unicode[STIX]{x1D705}^{-1}\log (l^{-1})$ . In addition, the flow models proposed are compliant with the channel centreline conditions, $u_{2}(x_{2},x_{3}=1)=u_{3}(x_{2},x_{3}=1)=0$ . Note, finally, that the standard zero-stress boundary condition typical of a channel is also attained; consider $\unicode[STIX]{x2202}_{3}u_{1}|_{x_{3}=1}=a_{2}\unicode[STIX]{x2202}_{3}\unicode[STIX]{x1D6FD}|_{x_{3}=1}+a_{3}\unicode[STIX]{x2202}_{3}\unicode[STIX]{x1D6FE}|_{x_{3}=1}=0$ .

In §§ 4.3 and 4.4, comprehensive results will be shown, all of which show that the system makes predictions consistent with physical expectations.

4.2 Numerical integration

Equations (4.10), (4.12) and (4.13) form a dynamical system,

(4.14) $$\begin{eqnarray}{\displaystyle \frac{\text{d}\boldsymbol{a}(t)}{\text{d}t}}=\boldsymbol{f}(x_{1},x_{2},\ldots ,x_{N})\quad \text{where}~\boldsymbol{a}(t)=\{a_{1}(t),a_{2}(t),a_{3}(t)\}~\text{and}~t\in [0,T],\end{eqnarray}$$

where $x_{N}$ parameters are not to be confused with the Cartesian spatial coordinates used hitherto, and $T$ denotes the integration period. Time integration of (4.14) predicts the phase-space trajectory (orbit) of system, $\boldsymbol{a}(t)$ . In this work, explicit forward Euler time integration is used to advance the system:

(4.15a,b ) $$\begin{eqnarray}\boldsymbol{a}(t+\unicode[STIX]{x1D6FF}t)=\boldsymbol{a}(t)+\unicode[STIX]{x1D6FF}t\boldsymbol{f}(x_{1},x_{2},\ldots ,x_{N})\quad \text{and}\quad t\in [0,T],\end{eqnarray}$$

where low-amplitude white noise is used to define initial conditions for the system, $\boldsymbol{a}(t=0)$ , as summarized in table 3. Note that the addition of unity to $a_{2}(0)$ is desirable, since this helps to initiate the system from a location in close proximity to a saddle node (results to follow). Table 3 contains a summary of cases, where for the present study we will quantify the system response to variable  ${\mathcal{R}}$ , for constant  ${\mathcal{L}}$ . Of course, further parameter variation is possible, but for the present article this range captures the values needed to attain chaotic orbits in the system. Similarly, the non-dimensional time step and non-dimensional integration time are provided, for reference. We will show that the system trajectory exhibits acceleration, and further generalizations could be adopted with local time steps, etc. but again such efforts would be superfluous to the core topic of this paper.

In the following sections, a series of results from (4.15) will be shown, which will highlight results from the system. In § 4, we will show instantaneous and Reynolds-averaged system results. In §§ 4.3 and 4.4, a series of scaling arguments for the magnitude of the Reynolds-averaged system parameters, $\langle \boldsymbol{a}(t)\rangle _{T}$ , are proposed, which provide further insights on how the intensity of streamfunction and streamwise velocity associated with the secondary cells should vary with system parameters. These results are used to recover scaling arguments on the intensity of compensated circulation, which agree with preceding results. Section 6 demonstrates recovery of the $tke$ imbalance arguments sustaining the Reynolds-averaged secondary flows, following substitution of the (4.15) model inputs. Finally, in § 7, Poincaré projection is used to identify sub-critical (i.e. dissipative system states, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{a}<0$ ) and conservative system forcing conditions, and to identify the fractal dimension (geometry) of the resultant attractor.

Table 3.

Summary of case attributes for dynamical system composed of (4.10), (4.12), and (4.13).

4.3 Results: non-periodic phase-space trajectories in governing equations

In order to demonstrate the potential for chaotic orbits from the system of ODEs, this section presents phase-space orbits of target cases in table 3. The cases capture parameter combinations corresponding to dissipative and conservative orbits, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{a}<0$ and $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{a}\approx 0$ , respectively, where the system conditions dictate the rate at which dissipative systems tend toward the fixed point, $\boldsymbol{a}(t)=0$ . Figures 11 and 12 show a series of phase-space orbits for increasing inertial conditions. As per table 3, the orbits are all initiated from approximately the same phase-space location (the addition of low-amplitude noise, $\unicode[STIX]{x1D702}$ , prevents initiation from precisely the same location, although the differences are negligible).

Figure 11.

Evolution of dynamical system composed of (4.10), (4.12) and (4.13), for ${\mathcal{L}}=250$ and ${\mathcal{R}}=100$ (Case E1; table 3). In phase-space visualization (a), blue and red squares – if visible – denote initial and final phase-space position.

Figure 12.

Evolution of dynamical system composed of (4.10), (4.12) and (4.13), for ${\mathcal{L}}=250$ and ${\mathcal{R}}=4000$ (Case E7; table 3). In phase-space visualization (a), blue and red squares – if visible – denote initial and final phase-space position.

For the weakest inertial conditions – Case E1, figure 11 – it is clear that the system is quickly overwhelmed by dissipation, as the orbit (a) amplitude quickly decreases while the system oscillates toward the origin (fixed point). The origin is a trivial fixed point, since it corresponds to a case for which no secondary circulations exist within the domain. Physically, this case corresponds to flow over a homogeneous roughness, for which there is no driving mechanism present to sustain roughness-driven secondary circulations. Gradual attenuation of the system is more clearly visible from the figure 11(b,c,d) time-series visualization, which shows how over-damping of the system suppresses its orbit amplitude (bearing resemblance to a damped sine wave). In fact, during early times, there are signs of chaotic orbits – most clearly visible in figure 11(d) for $tU_{0}\unicode[STIX]{x1D6FF}^{-1}\lesssim 3$ , but these oscillations are quickly damped as the system tends toward the origin. In a later section, we will demonstrate that these conditions correspond to dissipative systems (i.e. $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{a}<=0$ ).

Increasing the inertial conditions by an order of magnitude helps to preserve oscillations within the system, as can be seen in figure 12 for Case E7 ( ${\mathcal{R}}=4000$ ). Unlike for ${\mathcal{R}}=100$ (figure 11), here the phase-space orbit now exhibits chaotic reversals for all time, which thus corresponds to a conservative system (i.e. $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{a}\approx =0$ ). The attractor is shown in figure 12(a), where the perspective angle highlights orbits passing in close proximity to the saddle, and the occurrence of system reversals. The time-series evolution of parameters (figure 12 b,c,d) confirms the conservative nature of the system for ${\mathcal{R}}=4000$ , where large-scale variability in the system occurs in a non-periodic manner. In other testing, we have explored the limits upon continuous increases in ${\mathcal{R}}=4000$ . For ${\mathcal{R}}\gtrsim O(10^{4})$ , the system is unstable, with amplitude magnitudes continuing to increase before diverging rapidly to unphysical values.

In separate testing for intermediate inertial conditions (i.e. $10^{2}\lesssim {\mathcal{R}}\lesssim \times 10^{3}$ ), we show monotonic increase in the duration over which the system exhibits chaotic evolution with increasing  ${\mathcal{R}}$ . However, in these intermediate cases (Cases E1 to E5 in table 3), the system is always overwhelmed by dissipation such that the evolution is ultimately dissipative. In these intermediate cases, the orbit resembles those shown for Cases E1 and E7 (figures 11 and 12, respectively) – i.e. two limit cycles, connected at a saddle node at which orbit reversals occur – but the amplitude always tends to zero. We consistently found that ${\mathcal{R}}\gtrsim 3.5\times 10^{3}$ is required for conservative system evolution. For brevity, visualization of system orbit for these intermediate cases are now shown, although we note that this topic is revisited in § 7.

In order to summarize the system dynamics more succinctly, figure 13 shows an idealized sketch of system orbits in close proximity to the saddle seen on figure 12(a). As is depicted in the figure, orbits adjacent to the saddle point are considered stable, in the sense that they remain upon a limit cycle-like trajectory, and idealized boundaries demarcating the stable and unstable trajectories are denoted (orange and blue). Orbits residing in the stable regime eventually approach the saddle node from within the unstable region of the phase space (top and bottom), and in doing so the systems undergo global reversals. In the context of the LES and experimental results, such reversals correspond to the propensity of the systems to instantaneously transition between both downwelling and upwelling (i.e. between a HMP and LMP, respectively). The figure also includes idealized sampling points denoting locations where the system orbit is coincident with a plane parallel to the  $a_{2}{-}a_{3}$ plane, and placed at $a_{3}=0$ . By sampling a large collection of these points as the system evolves, a Poincaré projection can be attained, which provides additional insight on the dissipative or conservative nature of the system and the system fractal dimension (attractor geometry).

Figure 13.

Sketch of saddle point and idealized stable (grey) and unstable (black) phase-space trajectories. Orange and blue lines demarcate idealized boundary between trajectories likely to remain stable, and those corresponding to a system reversal. Shown also are idealized Poincaré projection points based upon an $a_{1}{-}a_{2}$ projection plane.

As a final assessment of the system dynamics, figure 14 shows the phase-space orbit for Case E6. Unlike the preceding figures, however, here we have added colour-coded data points at each discrete position on the orbit. The colour code is based upon the ‘system velocity’:

(4.16) $$\begin{eqnarray}\boldsymbol{u}_{s}(t)={\displaystyle \frac{\text{d}\boldsymbol{a}(t)}{\text{d}t}}\approx \unicode[STIX]{x1D6FF}t^{-1}\{a_{1}(t)-a_{1}(t-\unicode[STIX]{x1D6FF}t),a_{2}(t)-a_{2}(t-\unicode[STIX]{x1D6FF}t),a_{3}(t)-a_{3}(t-\unicode[STIX]{x1D6FF}t)\},\end{eqnarray}$$

where the error terms are presumed small and not consequential to the scientific deductions to be drawn from the present results. We then compute the magnitude, $\boldsymbol{u}_{S}\boldsymbol{\cdot }\boldsymbol{u}_{S}$ , which is used to select the symbol colour at each orbit point (warmer colours denote faster evolution, light grey highlights slowest points in orbit). This approach is helpful in terms of understanding where the orbits are accelerating and decelerating. It is clear from figure 14 that the orbit advances fastest in the outermost (in  $a_{3}$ ) regions of the phase space. The orbit decelerates as it approaches the saddle node, at which point the trajectory and ambient stochastic noise dictate the transition from stable to unstable orbit. Following this, the system rapidly accelerates as it moves outward from the saddle node.

Figure 14.

Evolution of dynamical system composed of (4.10), (4.12) and (4.13), for ${\mathcal{L}}=250$ and ${\mathcal{R}}=3500$ (Case E6; table 3). Here, phase-space system evolution at discrete locations is colour-coded based upon ‘system velocity’ (4.16), which graphically reveals the rate of system evolution. Here, the inner product is computed, $\boldsymbol{u}_{S}\boldsymbol{\cdot }\boldsymbol{u}_{S}$ , and hotter colours correspond to the largest values of $\boldsymbol{u}_{S}\boldsymbol{\cdot }\boldsymbol{u}_{S}$ .

Inspection of the figure 12(a) orbit demonstrates some similarities and differences with the Lorenz system. Of course, as was shown in § 4, the system of ODEs differ intrinsically to the Lorenz system, guaranteeing different orbits. In the Lorenz system, a constant heat flux drives a constant rise in temperature, which explains why $Z(t)$ – the counterpart to $a_{3}(t)$ in the present system – is always positive. Here, the system is predicated upon the balance of linear momentum associated with the large-scale secondary rolls. Yet the system predicts non-periodic transitions from states resembling that observed in figure 10 – a central LMP, flanked by adjacent HMPs – to a system with a central HMP and flanked by adjacent LMPs. These transitions induce differing net momentum fluxes, such that $a_{3}(t)$ cannot be positive for all times. Nevertheless, we will show that when the base flow is superimposed upon a canonical logarithmic profile – thus attaining the entire Reynolds-averaged flow – realistic profiles are recovered.

4.4 Results: Reynolds-averaged statistics

This section presents a series of Reynolds-averaged turbulence statistics based upon a posteriori outputs from the (4.15) system. The results are selected to demonstrate that realistic predictions are attained from the model. Recall that the system captures the large-scale (i.e. supra- $\unicode[STIX]{x1D6FF}$ ) rolls, and not the entire flow. For this reason, an additional logarithmic profile is added to the streamwise velocity base flow:

(4.17) $$\begin{eqnarray}\langle u_{1}\rangle _{T}(x_{2},x_{3})=\langle a_{2}\rangle _{T}\cos (k\unicode[STIX]{x03C0}x_{2})\sin (\unicode[STIX]{x03C0}x_{3})+\langle a_{3}\rangle _{T}\sin (\unicode[STIX]{x03C0}x_{3}/2)+\unicode[STIX]{x1D705}^{-1}\log \left({\displaystyle \frac{x_{3}}{l}}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D705}=0.4$ is the von Kármán constant and $l/\unicode[STIX]{x1D6FF}=10^{-4}$ is a pre-defined aerodynamic roughness length. This logarithmic profile captures the wall-normal profile due to sub- $\unicode[STIX]{x1D6FF}$ attached eddies, which were low-pass filtered from both variables before the development began. The first and second right-hand side terms within (4.17) represent the streamwise velocity distribution due to the turbulent secondary flows, both of which are additively combined with the logarithmic profile to capture the entire distribution.

Figure 15 shows Reynolds-averaged colour flood contours of streamwise velocity (a), and streamwise vorticity, $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}=-\unicode[STIX]{x1D6FB}^{2}\langle \unicode[STIX]{x1D713}\rangle _{T}$ (b). Superimposed upon both panels is a vector map for the spanwise and wall-normal velocity components. Results are shown for Case E6 (table 3), for which ${\mathcal{R}}=3500$ . As established in the preceding section, this Reynolds number is sufficient for conservative evolution of the system and the preservation of chaotic oscillations.

The $\langle u_{1}\rangle _{T}$ contour confirms that the model predicts a central LMP – evidenced by a pronounced plume of relatively low-momentum fluid at $x_{2}/\unicode[STIX]{x1D6FF}=1$ , and flanked by adjacent counter-rotating cells. By superimposing the vector map for the spanwise and wall-normal velocity components, the locations of upwelling and downwelling can be immediately visualized – in this case, upwelling and downwelling occurs within the LMP and adjacent HMPs, respectively, in precisely the same fashion observed from the LES and experimental results (§§ 2.1 and 2.2, respectively).

Figure 15(b) shows a corresponding contour of streamwise vorticity. To some extent, the distribution is evidenced from inspection of the vector field, but nevertheless the figure confirms that the rotational sign of the secondary rolls is such that positive and negative streamwise vorticity flanks the central LMP, as seen in the figure. Since the $x_{1}$ dimension is aligned in the ‘out of page’ direction, the sign of $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}$ is correct.

Figure 15.

Reynolds-averaged flow from numerical solution of (4.10), (4.12) and (4.13) system. (a) Colour flood shows streamwise velocity with $\langle a_{2}\rangle _{T}$ and $\langle a_{3}\rangle _{T}$ in (4.2), while in (b) colour flood shows streamwise vorticity, $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}=-\unicode[STIX]{x1D6FB}^{2}\langle \unicode[STIX]{x1D713}\rangle _{T}$ , for $\langle a_{1}\rangle _{T}$ in (4.1). In (a) a logarithmic profile accounting for the background, ambient shear associated with the channel flow has been added (4.17). In both panels, the vector map is spanwise–wall-normal velocity, $\{\unicode[STIX]{x2202}_{3}\langle \unicode[STIX]{x1D713}\rangle _{T},-\unicode[STIX]{x2202}_{2}\langle \unicode[STIX]{x1D713}\rangle _{T}\}$ , for $\langle a_{1}\rangle _{T}$ in (4.1). Result is shown for Case E6.

Figure 16(a) shows vertical profiles of Reynolds-averaged streamwise velocity, retrieved from transects corresponding to an LMP and HMP in figure 15(a) (i.e. from $x_{2}/\unicode[STIX]{x1D6FF}\approx 1$ and $x_{2}/\unicode[STIX]{x1D6FF}\approx 0$ , respectively). The HMP and LMP are denoted by black and grey, respectively, while the additional logarithmic profile (third right-hand side term in (4.17)) is denoted by dashed black profile. For the same system conditions, figure 16(b) shows Reynolds-averaged vertical velocity at spanwise locations equivalent to the locations at which the figure 16(a) profiles were retrieved, with the canonical profile, $\langle u_{3}\rangle _{T}/u_{\unicode[STIX]{x1D70F}}=0$ , included for reference (dashed black). The profiles exhibit physically sensible trends. Owing to the vertical trends within the profiles, the wall-normal gradient within the HMP exceeds the corresponding value within an LMP. This can be understood by the added drag beneath a HMP (Anderson et al. Reference Anderson, Barros, Christensen and Awasthi2015). The vertical velocity profiles are consistent with the position of downwelling and upwelling within a HMP and LMP, respectively. The profiles also agree with the profiles from LES (figure 3), especially for the smallest element height cases (i.e. Case A1).

Figure 16.

Reynolds-averaged vertical profile of streamwise velocity computed with (4.17) (a), and vertical velocity (b), for $\langle a_{2}\rangle _{T}=2$ and $\langle a_{3}\rangle _{T}=0.02$ , at discrete spanwise locations corresponding to an LMP (grey) and HMP (black) in figure 15(a). In (a,b), dashed black profiles represent logarithmic base flow for $l/\unicode[STIX]{x1D6FF}=10^{-4}$ (i.e. third right-hand side term in (4.17)), and canonical channel profile, $\langle u_{3}\rangle _{T}/u_{\unicode[STIX]{x1D70F}}=0$ , respectively.

Since spatial heterogeneity of the Reynolds stresses is foremost in the sustenance of the present turbulent secondary flows – Prandtl’s secondary flow of the second kind, equation (2.1) – it is of interest here to consider the spatial distributions of turbulent stresses from the dynamical system. Fortunately, this can be readily accomplished using the Boussinesq assumption (equation (3.10), Pope Reference Pope2000), and the base flow projections based upon a posteriori Reynolds-averaged system parameters, $\langle \boldsymbol{a}\rangle _{T}$ . Only one component of $\unicode[STIX]{x1D64F}$ cannot be computed using the present projections, $\unicode[STIX]{x1D61B}_{11}$ , which is due to the streamwise-homogeneity conditions. Nevertheless, the other components can be found readily, viz.:

(4.18) $$\begin{eqnarray}\displaystyle & \langle \unicode[STIX]{x1D61B}_{12}\rangle _{T}=-{\mathcal{R}}^{-1}\unicode[STIX]{x2202}_{2}\langle u_{1}\rangle _{T}=k\unicode[STIX]{x03C0}{\mathcal{R}}^{-1}\sin (k\unicode[STIX]{x03C0}x_{2})\sin (\unicode[STIX]{x03C0}x_{3})\langle a_{2}\rangle _{T} & \displaystyle\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D61B}_{13}\rangle _{T} & = & \displaystyle -{\mathcal{R}}^{-1}\unicode[STIX]{x2202}_{3}\langle u_{1}\rangle _{T}=-{\mathcal{R}}^{-1} [\langle a_{2}\rangle _{T}\unicode[STIX]{x03C0}\cos (k\unicode[STIX]{x03C0}x_{2})\cos (\unicode[STIX]{x03C0}x_{3})\nonumber\\ \displaystyle & & \displaystyle +\,(\unicode[STIX]{x03C0}/2)\langle a_{3}\rangle _{T}\cos (\unicode[STIX]{x03C0}x_{3}/2)\!]\end{eqnarray}$$

(4.20) $$\begin{eqnarray}\displaystyle & \langle \unicode[STIX]{x1D61B}_{23}\rangle _{T}={\mathcal{R}}^{-1}(\unicode[STIX]{x2202}_{2}^{2}\langle \unicode[STIX]{x1D713}\rangle _{T}-\unicode[STIX]{x2202}_{3}^{2}\langle \unicode[STIX]{x1D713}\rangle _{T})={\mathcal{R}}^{-1}\unicode[STIX]{x03C0}^{2}\langle a_{1}\rangle _{T}(k^{2}-1)\sin (k\unicode[STIX]{x03C0}x_{2})\sin (\unicode[STIX]{x03C0}x_{3})\qquad & \displaystyle\end{eqnarray}$$

(4.21) $$\begin{eqnarray}\displaystyle & \langle \unicode[STIX]{x1D61B}_{22}\rangle _{T}-\frac{2}{3}tke=-2{\mathcal{R}}^{-1}\unicode[STIX]{x2202}_{23}\langle \unicode[STIX]{x1D713}\rangle _{T}=-2{\mathcal{R}}^{-1}k\unicode[STIX]{x03C0}^{2}\langle a_{1}\rangle _{T}\cos (k\unicode[STIX]{x03C0}x_{2})\cos (\unicode[STIX]{x03C0}x_{3}), & \displaystyle\end{eqnarray}$$

and

(4.22) $$\begin{eqnarray}\langle \unicode[STIX]{x1D61B}_{33}\rangle _{T}-{\textstyle \frac{2}{3}}tke=2{\mathcal{R}}^{-1}\unicode[STIX]{x2202}_{23}\langle \unicode[STIX]{x1D713}\rangle _{T}=2{\mathcal{R}}^{-1}k\unicode[STIX]{x03C0}^{2}\langle a_{1}\rangle _{T}\cos (k\unicode[STIX]{x03C0}x_{2})\cos (\unicode[STIX]{x03C0}x_{3}).\end{eqnarray}$$

Substitution of $\langle \boldsymbol{a}\rangle _{T}$ into (4.18) to (4.22) yields the turbulent stress distributions. Note that for $k=1$ , equation (4.20) would indicate that $\unicode[STIX]{x1D61B}_{23}=0$ . In order to circumvent this, however, we selected $k=0.98$ for this case, which enables non-zero stresses.

Figure 17.

Reynolds-averaged deviatoric turbulent stresses: $\langle \unicode[STIX]{x1D61B}_{12}\rangle _{T}$ (a, equation (4.18)), $\langle \unicode[STIX]{x1D61B}_{13}\rangle _{T}$ (b, equation (4.19)),  $\langle \unicode[STIX]{x1D61B}_{23}\rangle _{T}$ (c, equation (4.20)),  $\langle \unicode[STIX]{x1D61B}_{22}\rangle _{T}-(2/3)tke$ (d, equation (4.21)) and $\langle \unicode[STIX]{x1D61B}_{33}\rangle _{T}-(2/3)tke$ (e, equation (4.22)). Cases shown are for ${\mathcal{R}}=3500$ and $\langle \boldsymbol{a}\rangle _{T}=\{2,2,0.02\}$ .

The turbulent stresses are shown in figure 17, for ${\mathcal{R}}=3500$ and Reynolds-averaged system parameter values summarized in the figure caption. Note, again, that these are only the stresses associated with the base flow, and as such direct comparison with the figure 5 contours is not possible. Nevertheless, important conclusions can be drawn from inspection of the contours. Figure 17(a,b,c) corresponds to shear-stress components $\unicode[STIX]{x1D61B}_{12}$ , $\unicode[STIX]{x1D61B}_{13}$ and $\unicode[STIX]{x1D61B}_{23}$ , respectively, where the distributions can be readily interpreted in the context of the underlying flow. Since $\unicode[STIX]{x1D61B}_{12}$ is developed by the negative of spanwise gradients in $u_{1}$ ((4.18) and figure 15 a), the central LMP and adjacent HMPs necessitate the positive–negative distribution seen in figure 17(a).

Figure 17(b) shows the distribution of $\unicode[STIX]{x1D61B}_{13}$ . We again emphasize that this distribution corresponds to the base flow projection – not the total – but nonetheless $\unicode[STIX]{x1D61B}_{13}$ confirms a relative excess and deficit in turbulent stresses beneath the central LMP and lateral HMPs, respectively. This result is conceptually identical to prior results from this problem – indeed, Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) show shear-normalized streamwise–wall-normal surface stress in the spanwise direction, and results confirmed elevated stresses beneath a HMP. This can be understood, too, in the context of figure 16, since the wall-normal gradient of $\langle u_{1}\rangle _{T}$ is larger within a HMP, relative to an LMP.

As per (2.1), the distributions of $\unicode[STIX]{x1D61B}_{12}$ and $\unicode[STIX]{x1D61B}_{13}$ do not actually contribute to the balance of streamwise vorticity. However, the second-order spatial derivatives of $\unicode[STIX]{x1D61B}_{23}$ constitute the shear-stress contribution to sustenance of $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}$ , and figure 17(c) shows non-zero variability in $x_{2}$ and $x_{3}$ (in fact, as per (4.20), the stress is projected upon the basis function, $\unicode[STIX]{x1D6FC}(x_{2},x_{3})$ , equation (4.5), but scaled by different amplitudes). Similarly, figure 17(c,d) shows the distributions of $\unicode[STIX]{x1D61B}_{22}$ and $\unicode[STIX]{x1D61B}_{33}$ , respectively, and (2.1) indicates that the difference, $\unicode[STIX]{x1D61B}_{22}-\unicode[STIX]{x1D61B}_{23}$ , regulates the normal stress contribution to $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}$ . The normal stress contours are heterogeneous in $x_{2}$ and  $x_{3}$ , confirming a finite contribution to evolution of the dynamical system.

Figure 18.

Reynolds-averaged contributions to production and dissipation of streamwise vorticity, based upon truncated Galerkin projection (i.e. balance between right-hand side terms of (4.9)). (a) Shows production by vortex forcing term, ${\mathcal{P}}^{\unicode[STIX]{x1D714}}$ (4.23), while (b) shows dissipation, ${\mathcal{T}}^{\unicode[STIX]{x1D714}}$ (4.24). Cases shown are for ${\mathcal{R}}=3500$ , ${\mathcal{L}}=16$ and $\langle \boldsymbol{a}\rangle _{T}=\{2,2,0.02\}$ .

As a final result in this section, figure 18 shows spatial contours of production and dissipation of streamwise vorticity by the vortex forcing term and turbulent torque, respectively, but in a Reynolds-averaged sense:

(4.23) $$\begin{eqnarray}{\mathcal{P}}^{\unicode[STIX]{x1D714}}={\mathcal{L}}{\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\cos (\unicode[STIX]{x03C0}x_{3}/2)\unicode[STIX]{x1D6FC}(x_{2},x_{3})\langle a_{3}\rangle _{T},\end{eqnarray}$$

and

(4.24) $$\begin{eqnarray}{\mathcal{T}}^{\unicode[STIX]{x1D714}}=-{\mathcal{R}}^{-1}\unicode[STIX]{x03C0}^{4}(1+k^{2})^{2}\unicode[STIX]{x1D6FC}(x_{2},x_{3})\langle a_{1}\rangle _{T},\end{eqnarray}$$

where readers will quickly recognize (4.23) and (4.24) as the Reynolds-averaged right-hand side terms of (4.9). In this sense, the role of ${\mathcal{P}}^{\unicode[STIX]{x1D714}}$ and ${\mathcal{T}}^{\unicode[STIX]{x1D714}}$ in the balance of $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}$ can be more succinctly appreciated. For the system conditions summarized in the figure 18 caption, (a,b) show ${\mathcal{P}}^{\unicode[STIX]{x1D714}}$ and ${\mathcal{T}}^{\unicode[STIX]{x1D714}}$ , respectively. For the present case, in which a central LMP is flanked by adjacent HMPs (see also figure 15 a), once can deduce that a continuous production of $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}$ must occur in the region, $0\leqslant x_{2}/\unicode[STIX]{x1D6FF}\leqslant 1$ . Indeed, figure 18(a) shows ${\mathcal{P}}^{\unicode[STIX]{x1D714}}>0$ in this region, due to the vortex forcing term, and this production is balanced by an equal and opposite ${\mathcal{T}}^{\unicode[STIX]{x1D714}}$ . This aspect of the present development is a departure from the LES and experimental results, since (2.1) shows a non-zero left-hand side due to continuous advection of $\langle \unicode[STIX]{x1D714}_{1}\rangle _{T}$ . However, the projections for $u_{1}$ and $\unicode[STIX]{x1D713}$ result in a vanishing of the advective term, since $J(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713})\rightarrow 0$ ((4.9) and accompanying text). This modest shortcoming could potentially be remedied with the consideration of alternate basis functions, or with generalization for high-order modes beyond the current first-order projection. Both pursuits have been considered, but we ultimately anticipate that higher-fidelity results would not alter the primary scientific result presented here: in the transitional region, where inner- and outer-driven secondary flows can both occur in close proximity, the transition between each occurs in a chaotic manner.

Over the region, $1\leqslant x_{2}/\unicode[STIX]{x1D6FF}\leqslant 2$ , the turbulent torque term is responsible for production of vorticity (4.24), which is balanced here by the vortex forcing term. Over the domain, a balance is apparent in a local and integral sense, i.e. $\int _{\text{d}^{2}\boldsymbol{x}}({\mathcal{P}}^{\unicode[STIX]{x1D714}}-{\mathcal{T}}^{\unicode[STIX]{x1D714}})\,\text{d}^{2}\boldsymbol{x}\approx 0$ , which conforms to the underlying consideration of the Reynolds-averaged system.