Channel flow with large longitudinal ribs

Abstract We present data from direct numerical simulations of flow through channels containing large, longitudinal, surface-mounted, rectangular ribs at various spanwise spacings, which lead to secondary flows. It is shown that appropriate modifications to the classical log-law, predicated on a greater wetted surface area than in a plane channel, lead to a log-law-like region in the spanwise-averaged axial mean velocity profiles, even though local profiles may be very different. The secondary flows resulting from the presence of the ribs are examined and their effects discussed. Comparing our results with the literature we conclude that the sense of the secondary flows is largely independent of the particular rib spacing whether normalised by channel depth or rib width. The strength of the secondary flows, however, is shown to depend on the ratio of rib spacing to rib width and on Reynolds number. Topological features of the secondary flow structure are illustrated via a critical point analysis and shown to be characterised in all cases by a free stagnation point above the centre of the rib. Finally, we show that if the domain size is chosen as a ‘minimal channel’ size, rather than a size which allows adequate development of the usual outer layer flow structures, the secondary flows can be affected and this leads inevitably to differences in the near-rib flows so that for ribbed channels, unlike plain channels, it is unwise to use minimal domains to identify details of the near-wall flow.


Introduction
There have been a number of recent papers exploring the nature of the flow in either a boundary layer or a channel when the wall surface contains longitudinal ribs whose height h is not a very small fraction of the boundary layer thickness, δ (or, in the case of a channel, the half-height, H). Typical examples include Vanderwel & Ganapathisubramani (2015), Hwang & Lee (2018), Vanderwel et al. (2019) and, most recently, Zampiron, Cameron & Nikora (2020) for an open channel flow. The ribs are much larger (usually O(0.1H), say) and more widely spaced than the 'riblets' classically studied in the context of seeking to reduce wall drag, and they generate secondary flows which may stretch through a significant portion of the boundary layer height. The secondary flows are driven by spanwise Reynolds-stress gradients and are thus secondary flows of Prandtl's second kind (Bradshaw 1987;Anderson et al. 2015). An example of the kind of flow produced is shown in figure 1(a), which is from Hwang & Lee (2018), who used direct numerical simulation (DNS) to compute a developing boundary layer flow over longitudinal roughness. Note, incidentally, that in some ways 'roughness' is a misnomer for this kind of surface spanwise heterogeneity: all surfaces are smooth and the ribs continue for the entire fetch so there is no form drag generated by pressure differences. The latter is what the usual kinds of surface roughness produce and a fully rough surface is normally defined as one for which this (pressure drag) component of surface stress completely swamps any contribution there may be from viscous stresses. For smooth longitudinal ribs on a smooth base, in contrast, the surface stress is entirely a result of viscous forces. This is an important caveat whenever such ribs are referred to as 'longitudinal roughness'.
Note that the particular rib height in figure 1(a) is approximately 0.08δ with centre-to-centre spacing of S/W = 4 (with W the rib width) and S/δ ≈ 0.8. This leads to downflow at the centre of the span and upflow above the centre of each rib, a pattern also noted by Vanderwel et al. (2019). On the other hand, Stroh et al. (2016) used changes in the flat-surface boundary condition, rather than ribs, and found secondary flows whose circulations could lead to vertical flows (e.g. over the centre of the 'quasi-ribs') of either sign depending on the value of H/S (for S/W = 2). It has been suggested that different geometrical arrangements can lead to rather different flow structures (e.g. Medjnoun, Vanderwel & Ganapathisubramani 2018;Yang & Anderson 2018;Stroh et al. 2020). For flows with longitudinal ribs, apart from the boundary layer DNS of Hwang & Lee (2018) and the channel large eddy simulation (LES) of Yang & Anderson (2018), there appears to be only the study containing data from (channel) DNS and a laboratory boundary layer (Vanderwel et al. 2019) and, even more recently, wind tunnel boundary layer studies by Medjnoun, Vanderwel & Ganapathisubramani (2020) and a similar set of experiments, but in a water flume, by Zampiron et al. (2020).
By presenting numerical studies of nominally two-dimensional (2-D) smooth-wall channel flow with smooth-wall longitudinal ribs, this paper seeks mainly to complement these earlier works, including comparisons where appropriate and exploring inter alia the effect of the rib spacing with respect to rib width, S/W and domain (or boundary layer) height, H/S (see figure 1b), and also Reynolds number. In particular, we consider the secondary flows by examining the critical points in the cross-stream flow, guided by the necessary topological constraints. We also make limited comparisons with some of the studies available in the literature which consider, by contrast, inhomogeneous flat surfaces (i.e. having no physical ribs).
A feature of all the previously published computational or laboratory studies is that almost none show spanwise-averaged mean velocity profiles. The exceptions are the recent works of Medjnoun et al. (2018) and Zampiron et al. (2020). In all cases, including the present, it turns out that the profiles lie significantly below the usual log-law line expected in regular smooth-wall channels and, not unnaturally, have a 'kink' at the rib height (i.e. at z = h). This will also be explored here, but it is worth noting immediately that the profiles obtained at specific spanwise locations in laboratory experiments also lie below the regular log-law, as shown by Medjnoun et al. (2020) and Zampiron et al. (2020). It is much more difficult to obtain spanwise-averaged profiles below z = h in such experiments and that was not attempted in these latter works. Medjnoun et al. (2018Medjnoun et al. ( , 2020 denote the profile offset by ΔU + in viscous units, in common with genuine rough-wall boundary layers; they call this the 'roughness function' but we will avoid this terminology, given that the boundary layer is not a genuinely rough-wall flow, as noted above. The present paper also discusses the spanwise-averaged turbulence stresses.
The next section outlines the three quite different DNS codes used and is followed, in § § 3.1 and 3.2, by a presentation of the basic flow statistics -the mean axial flow and the corresponding turbulence field. The influences and the nature of the secondary flows are discussed next, in § § 3.3 and 3.4, respectively, with topological considerations in § 3.5 and final discussion and conclusions in § 4.

Methodologies
The first set of computations were undertaken using a modern version of the in-house DNS code CGLES, originally written and described by Thomas & Williams (1997). This is a finite-difference, parallel, multiblock Navier-Stokes solver written so that its efficiency generally increases with mesh count. Cartesian, uniform meshes were used and second-order central differencing was applied to all spatial derivatives, with a second-order Adams-Bashforth scheme employed for time advancement using the pressure projection method. Continuity at the next time step was enforced implicitly by solving a Poisson equation for pressure using a parallel multigrid method. The initial mesh had 192 × 512 × 320 nodes in the x, y, z directions, respectively (see figure 1b for the coordinate system), with a mesh size of h/32. The domain size was thus 6h × 16h × 10h (i.e. 0.6H × 1.6H × 1.0H) and it contained two ribs symmetrically placed with S/W = 4. It should be emphasised that the domain was both too short and too narrow, with respect to its height (H/h = 10), to allow capture of the larger-scale motions common in the outer layer of channel flows. This was expected to lead, in the outer layer, to a rise in the mean velocity profile above the log-law, as has been frequently shown in the context of 'minimal flow channel' explorations (e.g. Jiménez & Moin 1991;Lozano-Durán & Jiménez 2014;MacDonald et al. 2017). It was nonetheless sufficient to allow generation of the secondary motions, since these are driven by processes near the bottom wall. The non-negligible impact on outer layer profiles and the secondary flows themselves will be discussed in due course. Using the same mesh size (h/32), a further computation was undertaken using a much larger domain size, 24h × 32h × 10h (i.e. 2.4H × 3.2H × 1.0H), in order to reduce domain size effects and clarify the effects of a limited domain on the near-wall flow.
No-slip conditions were applied at the bottom surface and on the ribs, which were captured naturally by the body-conforming Cartesian mesh. Periodic conditions were applied in the axial and spanwise directions and the top of the domain was free-slip. The flow for the above computations was driven by an applied pressure gradient, chosen to yield a spatially averaged surface stress equivalent to a friction velocity of approximately u τ = 0.92 ms −1 in the converged flow. With the specified viscosity, this led to a channel Kármán number, Re τ = Hu τ /ν, of around 850, with a normalised mesh size of Δ + ≈ 2.7 (in all directions, since the mesh was uniform), sufficiently low to ensure that the flow was adequately resolved by the DNS whilst high enough to ensure full turbulence. This was confirmed by estimating the smallest Kolmogorov length scale, η, in the flow. Assuming that energy production roughly balances dissipation, one can write for the latter = −u w dU/dz, where −u w is the Reynolds shear stress. Around z = 2h (where the stress was a maximum) this was approximately 0.7u 2 τ so that ≈ 1.3u 3 τ /h. With the Reynolds number hu τ /ν of 85 this leads to η ≈ 0.0056h = 0.18Δ, or Δ ≈ 0.6η. Near the surface, where the shear stress was close to its maximum value at z = 2h, and one again expects production to roughly balance dissipation, a similar calculation yields Δ ≈ 1.0η. Moin & Mahesh (1998) pointed out that the smallest length scale that must be resolved can be typically significantly larger than the Kolmogorov scale (O(10η), say) so we are confident that the present simulations capture practically all of the dissipation spectrum. We identify the two computations with Re τ = 850 as LC4ml and LC4ms, for the large and small domain cases, respectively. Another case (on the small domain) having a lower Reynolds number, Re τ = 500, was also computed; this will be termed LC4ls (see table 1).
A second set of computations was undertaken using the same numerical method as that described in Vanderwel et al. (2019). Only the salient details will thus be summarised here. The Navier-Stokes solver is based on a spectral method for the velocity-vorticity equations. Convective and viscous terms are discretised using the third-order Runge-Kutta and Crank-Nicholson methods, respectively. An immersed boundary-type method (IBM) was used to enforce the bottom surface morphology. Four cases were simulated, having S/W = 2, 4, 7 and 14, all with h/H = 0.082 and on a fixed domain of size 8H × 4H × H. In these cases, therefore, the domain was more than adequate for capturing the outer layer flow structures. A mesh size of 768 × 384 × 301 nodes in the x, y, z directions, respectively, was employed yielding a grid size of Δ + x = Δ + y = 5.2 and, in the vertical (for which a Chebyshev polynomial grid distribution was used), a minimum Δ + z of approximately 0.014. Again, a constant axial pressure gradient was applied, designed in these cases to yield a channel Kármán number of approximately 500. We identify the four cases as VS2, VS4, VS7 and VS14. The VS7 case, having S/W = 7 and four ribs within the domain span, is similar to the case presented by Vanderwel et al. (2019); in that paper, however, the ribs used were Lego blocks (to match the experiments) so had small pimples on the top surface. Again, no-slip conditions were applied on the bottom surfaces, axial and spanwise periodic conditions were applied at the domain ends and sides, respectively, but in these cases symmetry (rather than free-slip) was applied at the upper boundary. In view of some of the results shown later, it should be mentioned that the ribs were introduced using an IBM, based on the technique proposed by Goldstein, Handler & Sirovich (1993). This method can make the precise boundary location difficult to determine, so the boundary is in some sense rather 'fuzzy', especially around the rib corners. The implications will become evident in due course.
A final computation was undertaken using a very different DNS code -CANARD (compressible aerodynamics and aeroacoustics research code) developed at the University of Southampton. As indicated by the name, CANARD is a compressible flow solver for which the Mach number should be specified. For the present work, this was set to 0.25 in order to keep the compressibility effects minimal and not to cause excessive computational cost. The solver is based on a fourth-order pentadiagonal compact finite-difference scheme on a seven-point stencil that has been optimised for the maximum wavenumber resolution attainable (Kim 2007). Fourth-order Runge-Kutta time stepping was carried out with a Courant-Friedrichs-Lewy number of 1.0. Numerical stability was maintained by implementing sixth-order pentadiagonal compact filters for which the cut-off wavenumber (normalised by the grid spacing) was set to 0.89π (Kim 2010). Characteristics-based wall boundary conditions (Kim & Lee 2004) were applied. Data for cases of S/W = 4 with H/S = 1.25, 1.6 and 1.75 were obtained. The cases are denoted by KC4a,b,c, respectively, in table 1. A domain size of 8H × 4H × H was used and the number of grid cells was 510 × 1000 × 240 where each of the four ribs was resolved by 50 and 40 cells in the spanwise and vertical directions, respectively. The first wall-normal grid spacing was maintained at Δ + y = Δ + z ≈ 1.1. Boundary conditions identical to those used in the first set of computations described above were applied (i.e. a slip wall at the domain top). Parallel computing based on the message passing interface (known as MPI) was implemented, for which a precise and efficient technique specially designed for the compact finite-difference schemes and filters was used (Kim 2013).
For all methodologies, a long integration time was used to reach statistical stationarity and all the results shown herein were then obtained using a sufficiently extensive averaging time to ensure convergence. Table 1 lists all the various cases considered here and includes the salient parameter values for each. The formalism used for spatial averaging is well 915 A92-5   known (Raupach & Shaw 1982;Finnigan 2000) and follows from decomposing prognostic variables, like u i , into three components. Thus, u i = U i + u i +ũ i , where U i = u i is the time and space averaged velocity (the mean velocity),ũ i = u i − U i is the spatial variation of the time mean flow and u i = u i − U i −ũ i is the turbulent fluctuation; the overbar denotes a time average and angle brackets denote a spatial average. Applying time and spatial averaging to the momentum equations then leads to an additional dispersive stress term ( ũ iũj ) in addition to the usual Reynolds stress (u i u j ).

Basic statistics -the axial mean flow field
Before considering the secondary flows, it is of interest to explore some of the mean and turbulence statistics. In all cases except where indicated, for the region below the rib height extrinsic averaging is used -i.e. spatial averages taken over the axial direction and the entire span (including both fluid and solid regions). In every computation, the domain span covers an integral number of rib wavelengths. To set the scene, we show first, in figure 2(a), the temporally and spatially averaged profile of mean velocity for the KC4a computation, having S/W = 4, H/S = 1.25 and Re τ = 550. A number of points should be noted. First, the (green) profile displays a reasonable region of log-law behaviour (and has a mild wake region, as expected for a channel flow) but it lies well below the normal smooth-wall log-law, U + = (1/κ) ln(z + ) + A, with κ = 0.384, A = 4.65. These values of the constants fit the Hoyas & Jiménez (2008) smooth channel data (at Re τ = 550) but differ somewhat from those at higher Reynolds numbers, as in the more recent evaluations of (e.g. Marusic et al. 2010). Although, interestingly, none of the published (computational) works on flows with spanwise surface heterogeneities actually show U + profiles, it turns out that all of them yield significant profile offsets from the regular log-law; some of them are shown in figure 2(b) (discussed below).
Apart from this offset the most obvious feature of the U + profile is the kink that occurs at a z + around the top of the ribs (z + = h + = 55). The velocity on the rib's top surface is zero, leading to the relatively rapid fall in velocity as z approaches h from above. This kink is inevitable and would only disappear as W/S → 0 or W/S → 1. It becomes even more obvious if, below z = d, the spanwise averaging is done over the fluid region only, which is usually termed intrinsic averaging; this gives the dashed green line in the figure. Profiles at a specific location either between or above the ribs naturally do not show such kinks and we consider these later. It is of interest to compare the U + data in figure 2(a) (KC4a) with those obtained for the other computations including those from other workers (although, as noted above, these are not presented in their own papers). Figure 2(b) shows U + profiles obtained for the S/W = 4 case for different Re τ . The DNS boundary layer data from Hwang & Lee (2018) are included (HL4), along with the laboratory boundary layer data of Vanderwel et al. (2019) (VS4lego, but note that these do not extend below approximately z = 1.5h). Naturally, as Re τ decreases, the kink in the profile, which is at z + = h + = (h/H)Re τ , moves to the left. More interestingly, the log-law offset also seems to vary from case to case. For the channel computations the offset increases monotonically with Re τ . It is tempting to label this offset by ΔU + , as classically done for genuine rough-wall surfaces. Fitting the data to a shifted log-law yields ΔU + values of around 1.9, 2.0, 2.3 and 2.7, for the cases having Re τ = 490, 500, 550 and 850, respectively. But as argued in § 1, these flows are not rough-wall flows -there is no pressure contribution to the surface drag -and the fits (not shown) are at best rather mediocre, not least because the log-linear slopes are not close to 1/κ. Note that the small domain LC4 profiles (LC4ms and LC4ls) show U + eventually rising above the log-law; this is entirely a result of the limited extent of the domain, and the locations where the profile peels off from the log-law are largely consistent with the findings of, for example, Lozano-Durán & Jiménez (2014) and Abe, Antonia & Toh (2018). In contrast, the LC4ml, VS4 and KC4 simulations all have much larger domains, which (if the channel surface were flat) should be sufficient to maintain the log-law virtually all the way to the upper boundary. The computation of Hwang & Lee (2018) also used a sufficiently large domain but is of a boundary layer so the U + profile exhibits an outer layer wake region.
Close inspection of the mean flow profiles in figure 2 shows that, in each case, the straight-line region is not closely parallel to the regular log-law (so a value for ΔU + is anyway somewhat arbitrary). This can be explained at least partly by the way the ribs influence the effective height of the surface. Vanderwel et al. (2019) suggested an effective zero plane displacement height equal to the geometric average of the surface height across the span. For S/W = 4 this is h/4. However, it is arguably more appropriate to use the height at which the surface drag appears to act -as explained by Jackson (1981). An exact value of d/h on this basis requires a computation of the moment generated by the surface frictional forces on the rib's horizontal and vertical surfaces and the Appendix A shows how this is done. For the particular cases in figure 2 (all having S/W = 4) the result is d/h = 0.34, significantly above the geometric average. The U + profile for KC4a is plotted against (z − d) + with this value of d in figure 2(a) and it clearly increases the extent of the straight-line region which, nonetheless, remains non-parallel to the regular log-law. No physically reasonable value of d/h (i.e. between zero and unity) would improve things. This, together with the large offset from the regular log-law, leads one to consider whether u τ is in fact an appropriate normalising friction velocity for log-laws in ribbed channels. The following analysis suggests that it is not because, crucially, the wetted area on the bottom surface of the channel is larger, and the cross-sectional area is smaller, than if it were a regular channel (i.e. with h = 0).
For an axially straight ribbed channel, the axial force balance is written as where τ w,rib = ρu 2 τ,rib is the effective wall stress, A x is the cross-sectional area of the channel in the axial direction and A w is the wetted surface area of the wall at the bottom. These areas are given by where Δx is the axial length of the channel and Δp the pressure difference across that length. Substituting these into the balance equation and rearranging gives For a plain channel (h = 0), the above equation is, as usual, So the friction velocity of a plain channel that yields the same pressure gradient to that of a ribbed channel is The classical log-law for the plane channel is If we assume that the same log-law can also be applied to a ribbed channel when normalised by u τ,rib , i.e.
the following rearrangement can be made in order to compare it directly with the plain channel case: This gives an equivalent log-law for a ribbed channel: Note that β is always below unity, so the results indicate that the appropriate log-law line for a ribbed channel will always have both a lower slope and a significant offset when compared with the plain channel log-law. The analysis must clearly fail as S/W approaches unity, for the sidewall contribution to the surface wetted area, the A w in (3.2a,b), is still included.
For the specific case shown in figure 2(a), for which h/H = 1/10, W/S = 1/4 and h/S = 1/8, the above expressions give u τ,rib = 0.883u τ , κ rib = 0.435 and A rib = 3.821. The modified log-law is shown as the dotted red line in the figure. The fit with the data remains imprecise, but it is significantly improved if account is taken of the zero plane displacement. Assuming that the flow beneath z = d can be ignored for the purpose of defining the flow's effective cross-sectional and surface-wetted areas, it is straightforward to repeat the analysis, which leads to a modified value of β, given by (3.11) With d/h = 0.34 this changes the values of the parameters to u τ,rib = 0.903u τ , κ rib = 0.425 and A rib = 3.958. These are minor adjustments but lead to a noticeably improved fit to the data, as shown by the solid red line in figure 2(a). A more extensive test of this modification is provided by the four VS cases, which each had a different value of S/W ranging from two to 14 (see table 1). Figure 3 shows that provided S/W is large enough, the data collapse well to the modified log-law. As noted above, one might anticipate that for small S/W the analysis would be less satisfactory, as is demonstrated by figure 3(d) (S/W = 2). Even quite large adjustments to the d/h used in this case do not improve the fit.
The boundary layer data of Hwang & Lee (2018) are included in figure 2(b) and the profile has an even larger offset then the channel flow cases (as does the boundary layer data of Vanderwel et al. (2019)). In any case, the analysis clearly needs adjustment for zero-pressure-gradient boundary layers. Consider the momentum integral equation for a 2-D boundary layer, which can be written as (3.12) in the usual notation. If, at a specific axial location in a zero-pressure gradient flow, a 'virtual' pressure gradient can be considered to model the growth of the momentum thickness and yield the same wall stress, we can apply the following equation locally: Comparing this with (3.4), a channel analogy can thus be achieved by replacing H in the latter by (2θ + δ * ). Using the values of momentum and displacement thickness given in Hwang & Lee (2018), for this particular flow (at x/θ in = 300 for their P12S3 case -i.e. S/W = 4 in our notation) we find that (2θ + δ * )/h = 5.4. Taking the same value of d/h (0.34) as in the corresponding channel case and using κ = 0.41 and A = 5 as appropriate for this boundary layer having Re τ = 292, and noting also that Avasarkivos et al. (2014) found that A did not vary within the range 125 < Re τ < 550, we obtain β = u τ,rib /u τ = 0.775, κ rib = 0.495 and A rib = 3.393. The resulting modified log-law is shown along with the data in figure 3(e). Given the uncertainty in the precise value of d/h and the possible implications of a non-zero d for the momentum and displacement thicknesses derived by Hwang & Lee (2018), the fit to the data is acceptable -not least in confirming that the offset is greater for this boundary layer than for the channel case at the same S/W. We conclude that in a ribbed channel, or indeed in a similarly ribbed boundary layer, a reasonable fit of the velocity profile to a log-law can be obtained provided one modifies the normalising friction velocity to account for the fact that the wetted area is larger and the cross-sectional area is smaller than for the corresponding plane channel. Accounting for  the zero plane displacement further improves the fit. Nonetheless, we point out that there is no really convincing reason why such log-laws should appear in spanwise-averaged profiles when there are significant secondary flows, especially if the strength of these is large. This is discussed further in § 3.3, when profiles at specific spanwise locations are presented.
3.2. Basic statistics -the turbulence field Stress profiles (normalised by u 2 τ ) are shown in figures 4 and 5. Recall that extrinsic averaging (i.e. with solid regions included) is used below z = h. This is the only kind of averaging that ensures continuity in the momentum fluxes across z = h, as fully explained by Xie & Fuka (2018), who explored the whole issue in some detail, recognising the different possible kinds of averaging ( Raupach & Shaw (1982), for example). Note first that, in figure 4(a) showing shear stress data for KC4a, the dispersive shear stress (caused by the spanwise inhomogeneity in the mean flow, generated by the secondary motions) remains non-zero all the way to z/H ≈ 0.6 (z/h = 6). There is thus a noticeable deviation up to around this height between the Reynolds stress profile and the usual straight-line total stress for a regular channel (shown by the dashed line). In common with the U +  profile there are also rapid variations in the stress gradients around the top of the ribs (z/H = 0.1), with the viscous stress in particular only being significant in this region and near z = 0, as expected. Below z = h the total stress profile includes the contribution to the stress at height z from the viscous sidewall stresses (integrated downwards from the top of the rib and including the viscous stress on top of the rib). If the ratio of the rib volume to the computational domain volume were zero, the total extrinsic stress would continue along the expected straight line all the way to z = 0 (Xie & Fuka 2018). However, because the forcing term applied at every cell in the computation is applied also within the solid volume of the rib, one expects the total stress at a z below z = h to differ from the expected straight line by a factor equal to the ratio of the fluid volume (above z) to the total volume of the domain (above z). At height z the factor is (1 − W(h − z)/S(H − z)) which, in this case, is 0.975 at z = 0, increasing to 1.0 at z = h, so the data are close to those expected. Small differences, especially around the top of the rib and the bottom surface, are probably the result of errors arising in the estimation of the friction at the walls. Figure 4(b) shows the corresponding axial normal stress component and its dispersive counterpart for KC4a, along with data from the other S/W = 4 simulations. The large spike (especially in the dispersive contribution) centred around z = 0.1H arises because there is a step change in axial velocity across the horizontal plane where the velocity is very small just above the rib (and zero inside it). The major point to note is that, despite the significantly non-zero dispersive stress up to around z/H = 0.6, the LC4ml, VS4 and KC4a data are, above the rib, similar to the smooth-wall channel data of Hoyas & Jiménez (2008); data from the small domain computation, LC4ms, are not. Similar behaviour is evident in the profiles of the other two normal stress components, shown in figure 5. For these stresses the LC4ms simulations (not shown) yield values lower than those from LC4ml, VS4 and KC4a (and also the smooth channel), rather than higher as in the axial stress case (figure 4b). As anticipated, the limited domain size has a major influence on the stress components as well as the structure of the secondary motions (see later). Indeed, the changes in the stress profiles caused by using a minimal domain closely follow the findings of Abe et al. (2018), who demonstrated that the spanwise and vertical stresses were lower than for a full domain, whereas the axial stress was higher (as in figure 4b). This strengthens the earlier, not unexpected, conclusion that minimal flow channel domains cannot be relied on for flows of this kind, for if the stress profiles above the ribs are incorrect (as they are) then the spanwise stress gradients which are the essential cause of the secondary motions will be incorrect. It is also of interest that the dispersive components of the vertical and spanwise stresses (the lower red lines in figure 5) are very small; they provide a much smaller contribution to the total stress than is the case for the axial normal stress and the shear stress (figure 4).

The influence of the secondary flows
We turn now to consider the influence of the secondary flows. These will obviously affect the dispersive stresses. Figure 6(a) shows the dispersive shear stress and figure 6(b) the vertical component of the normal dispersive stress, for the four VS cases. Very similar plots arise if these stresses are normalised by their corresponding Reynolds stresses and note that above the rib height (z/H = 0.1) ww + never exceeds 6 % of w 2 + . As S/W rises the dispersive stresses become relatively larger in the outer region, although the strength of the secondary flows is much weaker in this region than nearer the ribs (see later). That their strength depends on the geometrical parameter S/W is at least suggested by figure 3(a-d), which show that the mean velocity profiles for the four (VS) data sets having closely constant Re τ ≈ 500 sink below the standard log-law by an amount which increases with increasing W/S. The trend must eventually reverse, since when W/S = 1 there are no ribs and the regular smooth-wall log-law must be recovered, as it must also for W/S = 0.  (see  table 1) because the domain span is fixed, W/h is fixed, and S/W is varied by changing the number of ribs within the span. The S/W = 4 case has H/S = 1.75 and the case on the other side of the (speculated) location of the peak has H/S = 3.45, W/S = 0.5. Different values of H/S for fixed S/W, particularly if it is below H/S = 1, might be expected to change the secondary flow strength somewhat, as well as changing the proportion of the channel height over which the secondary flows are particularly noticeable. This will be the topic of a future study. We show in § 3.5 that the secondary flows are essentially the same for H/S between 1.25 and 1.75 (with S/W = 4). It is worth emphasising, therefore, that S/W is clearly more significant than H/S, at least for the current range of the latter. Yang & Anderson (2018) and Medjnoun et al. (2018) both suggest that the maximum swirl strength occurs for H/S = O(1) as H/S varies and there is nothing in the present data which would contradict that; we have not studied cases for which H/S < 0.5. A direct measure of the size of the log-law shift could conveniently be taken as the value of β d , the factor by which the friction velocity is reduced in the modified log-law and deduced using (3.11) as discussed earlier. This is plotted (as 1 − β d ) in figure 7(b). We have made the assumption that d/h varies smoothly through the values used for the four modified plots in figure 3 and with values of zero and unity at W/S = 0 and 1, respectively. As noted earlier, we cannot expect the model to be appropriate close to, and certainly at, W/S = 1, for then the situation reverts to a simple plain-walled channel but of lower depth (H − h, 915 A92-13 I.P. Castro, J.W. Kim, A. Stroh and H.C. Lim (c), except that z-scaling with α = 0.5, see (3.14), is used for u τ . In panels (b), (c) and (d) no attempt is made to identify separately the closely clustered profiles.
in fact), even though the expression for the wetted surface area, A w still contains the 2h contribution (see § 3.1). Also, β d depends critically on the ratio of rib and channel heights, h/H. Smaller values of the latter than that used for the VS computations (h/H = 0.085) will clearly lead to smaller 1 − β d and this is illustrated in figure 7(b) by the solid red line, for which h/H = 0.04. Very similar plots are obtained if one uses (A − A rib ) rather than (1 − β d ) as a measure of the change from the classical log-law. One might anticipate that the strength of the secondary motions will affect the size of the variations in the mean axial surface friction across the span. Figure 8(a) shows this variation for S/W = 4 cases as an example, with wall stress values generally estimated by taking the local surface stress to be ν(∂U/∂z) calculated at the first grid point above the surface and normalised so that, in each case, the spanwise average (not including contributions from the rib's sidewalls) is unity. In all cases the frictional stress on the top of the rib is higher than the average and generally not dissimilar to values near the edges of the span (i.e. at the centre of the gaps between ribs). This could perhaps suggest that the vertical flow is downwards towards the ribs (and upwards just outboard of the rib), which is opposite to what was apparently found by, for example, Vanderwel & Ganapathisubramani (2015) and Vanderwel et al. (2019), whose visualisations suggested upward flow all the way from the centre of the rib to the top of the domain (but see later) with, presumably, consequent spanwise flows near the top surface directed towards the rib centre from its corners. However, it is not really possible to deduce the direction of near-surface spanwise flow from the relative strength of the axial flow there and we return to this issue in § 3.5. Note that the computation at Re τ = 850 yields a rather higher axial friction above the rib than in the other two cases (figure 8a) with correspondingly smaller friction outboard. This might suggest a Reynolds number effect and we return to this in § 3.4, where a possible reason for the differences in axial friction seen in the VS4 and KC4a cases is also suggested.
The secondary motions lead naturally to differences across the span in the vertical profiles of mean velocity. Before looking more closely at the secondary motions it is of interest to inspect the degree of collapse, or otherwise, in the axial mean velocity profiles at different locations, using a wall-scaling based on the local wall stress, rather than the spanwise-averaged one (used for the span-averaged velocity profiles in figures 2 and 3). Figure 8(b) shows that using the latter for local profiles does not yield much collapse (except at the top of the domain), as would be expected. Note that for profiles above the rib (labelled (5) and (6) in the legends of figure 8b-d), z is measured from the top surface so that profiles in the viscous sublayer can be expected to collapse. Note too that intrinsic averaging (below z = h) is used for the spanwise-averaged profile in figure 8(b), since only this would be expected to lead to the usual viscous law collapse close to the wall. The two most obvious outliers are for y/S = 0.3 and 0.35, which (as seen in figure 8a) are locations close to the rib where the local u τ is much lower than the spanwise average. Very near the bottom wall one would expect classic viscous scaling to yield collapse if the local u τ were used and this is confirmed in figure 8(c); below approximately z + = 3 the profiles do indeed collapse onto the viscous wall law. It might be tempting to try other forms of scaling, e.g. using a local u τ satisfying (3.14) where u τ ( y) is the local friction velocity varying with spanwise location and u τ is the spanwise average normally represented herein by u τ . This expression is chosen so that near z = 0 the local u τ dominates whereas with α > 0 the spanwise-averaged u τ dominates, at least in the upper region, to an extent depending on the value of α. However, figure 8(d) shows that although one can choose a value for α which leads to reasonable collapse over all z for most of the profiles (α = 0.5 is about the best), this approach is unlikely ever to work for the more extreme outliers (e.g. at y/S = 0.35) where the wall stress is far from the spanwise average. There seems, in any case, no physical reason to expect an exact overall collapse; it has been known for a long time that secondary motions within a boundary layer lead to distortion of the mean velocity profile (e.g. Mehta & Bradshaw 1988) and one could argue that even if a spanwise-averaged velocity profile appears to yield a reasonable (although modified) log-law region, as in fact was shown in § 3.1, this must be somewhat fortuitous. It seems obvious that for any flow like these, local velocity profiles -at least those within the gap between ribs but close to them -cannot have the usual log-law form, so there is no expectation that a spanwise-averaged profile could display such a log-law. Furthermore, since Mehta & Bradshaw (1988) actually studied a smooth-wall case in which the secondary motions were generated by upstream vortex generators, we suspect that this conclusion is independent of the nature of the surface. But, of course, as the secondary motions become weaker and weaker one expects the classical log-law to re-emerge. 3.4. The nature of the secondary flows Examples of the secondary flows are shown in figure 9, which presents contours of swirling strengths overlaid on velocity vectors in the spanwise plane, for the two cases KC4c and VS4, both having S/W = 4, H/S = 1.75 and approximately the same Re h . Swirling strength is here defined in the usual way as λ Ci ω x /|ω x | (Zhou et al. 1999), as this is generally recognised to be a more satisfactory way of identifying swirling motions. The vorticity itself, ω x , is less appropriate as it cannot distinguish between genuine vortex motions and regions of strong shear. The sign of the vertical velocity just above the centre of the rib seems to be positive (upwards) for VS4 but negative (downwards) for KC4c, which was confirmed by close inspection of the velocity field. A closer look at the secondary velocity fields reveals further features of the flow structure above the rib. Figure 10 shows two S/W = 4 cases at identical H/S and rather different Re h , but (in their plotting) more highly resolved around the rib region than the plots in figure 9 (or figure 1a). The crucial point is that there is an elevated critical point (a saddle) above the centre of the rib. Its vertical location seems to depend on the Reynolds number. In figure 10(a) it lies around z/S = 0.23 whereas in figure 10(b) it is significantly above that point. Beneath the saddle, the flow is downwards towards the rib centre and outwards on the surface either side of there, but the strength of that flow depends (at least partly) on where the saddle is. Above the saddle, the flow is always upwards, all the way to the top of the domain. This elevated saddle point seems not to have been clearly identified previously (but see § 3.5). Beneath it there are two sets of nodes, whose approximate centres are shown on the figure. These nodes are in somewhat different positions in the KC4a case (figure 10a, Re h = 55) and clearly rather weaker than those seen in figure 10(b) (Re h = 85) -one might expect the secondary flows to become rather 'tighter' at higher Reynolds numbers. Likewise, the outboard recirculating regions are more diffuse and extend higher in the lower Reynolds number case (figure 10a).
Recall now that the VS4 flow of figure 9(b) seems not to contain the elevated saddle point. This is almost certainly because the IBM used to model the rib in the VS4 case gives a rather fuzzy rib boundary, particularly at the rib corners, allowing the cross-flow to move smoothly over the corners rather than having to separate. For that case, therefore, the flow on the top surface of the rib is inward towards the centre from the corners, leading to a half-saddle of separation at the centre and flow upwards from there to the top of the domain. In fact, a further (VS) computation using a more refined IBM which reduces the fuzziness of the rib corners, did yield genuine separation at those corners so that the overall flow is more like that shown in figure 9(a). Incidentally, this rather fuzzy boundary probably explains the small differences in axial surface friction on the rib seen in figure 8(a) between the KC4a and VS4 cases; it is not so easy to determine wall friction when precise boundary locations are uncertain. It is worth emphasising, incidentally, that the main features of the secondary flows seen in figure 10 are not changed by the increasing Re h . This was the finding of Vanderwel et al. (2019). In fact, in their paper, despite the lack of resolution in the flow near the ribs, one can just detect a rise in the critical point location above the rib centre as Re h rises (compare their figures 2c and 2d). We comment finally on the effect of using a minimal domain on these secondary flows. Figure 11 shows two computations having the same Reynolds number, S/W and H/S, one obtained using a small domain size (LC4ms) and the other using a much larger domain (LC4ml). The differences are clear. Although the general topology of the recirculation pattern is the same, the small domain significantly constrains the secondary flows, making those just above the rib rather more intense, thus forcing the saddle point to move higher above the centre of the rib and weakening the recirculating velocities further aloft. This can be seen by comparing the vector lengths at, say, z/S = 0.8vertical velocities are noticeably larger when a large domain is used. Even the flow near the rib is thus influenced if too small a domain is used, so it cannot be argued that a minimal channel can be used to obtain adequate near-wall flows -unlike the case for plane channels. 3.5. Topological considerations Next, we consider the topological constraints on the nature and number of critical points in the spanwise plane. This is a great help in interpreting the cross-plane visualisations and reduces the likelihood of incorrect conclusions. Denoting saddle points by S and nodes by N, with half-saddles and nodes (on boundaries) as S and N , it is known (e.g. Hunt et al. 1978) that for a singly connected domain like this, A good way of clarifying the presence and location of critical points is to inspect the cross-plane velocity field -in particular to consider (V 2 + W 2 ). This quantity will ideally be zero at all critical points. Figure 12(a) shows a contour plot of γ = log 10 (V 2 + W 2 ) from the KC4a computation at Re h = 55. Critical points show up as concentrated regions of white or near-white, where γ → −∞. As noted earlier, there is a surface half-saddle at the centre of the top of the rib, below the saddle just above, and a pair of half-saddles on each of the two side surfaces of the rib (including ones at the rib corners). The critical point structure is sketched in figure 12(b). Note that on the central y = 0 line there are matching half-saddles at the top and the bottom of the domain, with downward flow beneath them. A half-saddle also exists at the top of the domain (z = H) on y/S = 0.5, matching the saddle below. Features close to the bottom surface in the gap between ribs are rather less easy to identify, but close inspection of the vector field (not shown here) indicated that there are half-saddles near y/S = 0.17 and 0.83, with two nodes between these and the two rib sidewalls. In total, there are six nodes, one saddle and 10 half-saddles, satisfying the topological requirement given above. From the visualisations shown in the various published papers on this topic it is difficult, if not impossible, to identify all these various critical points, although there has been one previous attempt (Stroh et al. 2016) for a related flow (having no ribs but spanwise changes in surface conditions). Indeed, often the resolution in the published plots shown is insufficient to clarify, for example, the structure near the rib's top surface or along the bottom surface between the ribs. This is sometimes because (in a laboratory study) the particle image velocimetry (PIV) field is not sufficiently resolved and (in a numerical study) the vector density shown is insufficient -even though the computational mesh resolution may be quite sufficient to delineate the flow structure accurately. Incidentally, at higher Reynolds numbers one expects the eventual appearance of further critical points near the bottom corners of the rib; these are not really visible in figure 10(b) (Re h = 85), for example, but would become more so at higher Reynolds number. Any attempt to construct the overall critical point structure in such cases would need to satisfy the topological constraint. Figure 13 shows the topological structure for the three KC4 cases, having H/S = 1.25, 1.6 and 1.75. It is clear that the secondary flow below z/S ≈ 0.8 is essentially independent of H/S. Increasing H/S does lead, however, to a just perceptible rise in the location of the saddle point above centre of the rib. The corresponding swirling strengths for these three cases are shown in figure 14 and indicate, similarly, that H/S is not an influential parameter. Indeed, an average of the (modulus of) the swirling strength across the span and from z = 0 to z = 2h differs by no more than approximately ±2 % across these three cases. Above z/S ≈ 0.8 the figures emphasise the extremely small values of the cross-stream velocities. Although figure 14 might suggest that the flow is practically homogeneous across the span in the upper region, it is not exactly so, as is obvious from figure 13. Nonetheless, the variations in axial velocity are sufficiently small that the dispersive stresses are, as discussed in § 3.2 and seen in figures 4 and 5, very close to zero above z/S ≈ 0.6 (i.e. z/H ≈ 0.5 for H/S = 1.25).

The mean velocity profiles
The present work has shown that there are regions across the span (particularly those close to rib walls) where the vertical profile of axial velocity (U + ) is far from having any classical shape and there is no obvious scaling which allows local profiles at different points across the span to be closely collapsed. This is not very surprising. Nonetheless, what is somewhat more surprising is that the spanwise-averaged profile does contain a reasonable log-linear region. However, this lies below the classical log-law by an amount which depends on W/S -or, more strictly, on the difference between W/S and either zero or unity. This region also has a different slope, which is not consistent with the classical Kármán constant. These differences from the usual log-law are probably greatest when W/S ≈ 0.4 (figure 7b), which is when the strength of the secondary flows is at its greatest (figure 7a). A major conclusion of our work is that if, using an appropriate analysis ( § 3.1), account is taken of the increase in wetted surface area and reduction in cross-sectional area (because of the ribs), and also of an appropriate zero-plane displacement, d, then a modified log-law can be predicted. This has different values of κ and A, the constants in the classical log-law relation, and it fits the data quite well provided S/W is not too small. Furthermore, d can be directly calculated using the computed frictional forces on all the rib walls (see the Appendix A). The analysis can be extended to zero-pressure-gradient boundary layers and this again leads to a modified log-law which fits the data reasonably well. The analysis inevitably fails as W/S → 1.
The fit to any kind of log-law is perhaps surprising because, as argued earlier, even on a flat surface secondary flows distort the U + profile (e.g. Mehta & Bradshaw 1988). But in the present case the secondary flows, although very distinctive are, compared with the strength of the axial flow, very weak indeed. They represent a negligible contribution to the total energy in the flow -the pointwise maximum of W 2 /U 2 , for example, nowhere exceeds approximately 1.2 × 10 −3 , so it seems reasonable that these secondary motions are not large enough significantly to negate the usual assumptions on which the appearance of a log-law is based.
The use of a minimal domain size leads, as is well known, to artificial behaviour in both the mean velocity and turbulence stresses in the outer region of the flow. However, in these cases, where secondary flows are significant, we have shown that it also leads, perhaps not surprisingly, to modification of the secondary flow structures near the ribs, so adequate determination of the latter requires appropriately large domains, just as required for proper characterisation of the outer flow in plane channels.

The secondary flows
Overall, in terms of the secondary flows and, in particular, the direction of the large-scale swirling motions above and outboard of the rib, the present results are similar to those of other investigators who have considered elevated ribs rather than changes in surface condition. The boundary layer DNS of Hwang & Lee (2018), for example, with S/W = 4, shows upflow over the ribs and downflow between them at the downstream location where H/S ≈ 1.25, exactly as seen in our H/S = 1.25 (figures 9b and 14a) and 1.6 and 1.75 (figure 14b,c) cases. Likewise, the channel LES of Yang & Anderson (2018) showed that for an H/S = 1.56, S/W = 8.5 case there was upflow over the ribs and downflow in the spaces between them. Furthermore, the H/S = 1.1, S/W = 5.9 case studied by Vanderwel et al. (2019) has a large-scale upflow above the ribs but with a (just) discernible critical point just above the ribs with downflow below it. Taken with the present cases (S/W = 4), these results suggest that S/W is not significant in setting the direction of the large-scale secondary flows. On the other hand, Yang & Anderson (2018) also considered an H/S = 1, S/W = 13.2 case; the visualisation (of vorticity and cross-flow velocities) tends to suggest a large-scale downflow over the ribs and upflow between them, which led the authors to propose a switch in orientation when H/S is somewhere between their two cases of H/S = 1 and 1.56. It should be noted, however, that the ribs used by Yang & Anderson (2018) were small 'house-shaped' obstacles, having sloping roof sides up to a narrow top, which reduces the likelihood of separations. In addition, since the obstacles were three-dimensional, pressure forces would have contributed to the surface drag and this could perhaps significantly change the topography's generation of secondary flows. Compared with the present scenario of smooth 2-D rectangular obstacles, their study may well, therefore, represent an isolated case. Zampiron et al. (2020), who studied flow over ribs in a water flume, found large-scale upflow over the ribs for all their cases, certainly down to H/S = 0.64. Furthermore, no switch was apparent in the Vanderwel et al. (2019) cases either. They used sharp-topped triangular ribs not too dissimilar to the ribs of Yang & Anderson (2018) (except that they were 2-D) and these encouraged the large-scale rotating flow near each side of the ribs to sweep up, meeting at the top of the rib and continuing upwards. This cannot happen for flat-topped ribs with sharp corners, even if the large-scale secondary flow still leads to upflow above the ribs. In such cases there will always be a smaller scale region near the top of the rib, encompassing separations at the corners with the concomitant downflow at the rib centre and thus a critical point aloft, as shown in figure 12(b). This is exactly the situation in the H/S = 1.25 case of Hwang & Lee (2018), the H/S = 1.1 case of Vanderwel et al. (2019) and the H/S ≈ 0.8 cases of Medjnoun et al. (2020). It should also be noted that one of the cases studied by Hwang & Lee (2018) had H/S as low as approximately 0.4 and the present VS14 case had H/S = 0.5; in both cases the secondary flows were in the same direction as for the larger H/S cases. We conclude that a directional switch never occurs for 2-D ribs (of any shape) as H/S changes. Nonetheless, in most of the visualisations of those authors mentioned above (who used rectangular ribs) it is possible to discern, with more or less difficulty, the critical points just above the rib centres, with very local downflow below and the larger scale upflow above. Although the authors did not discuss these smaller-scale features, they are very clear in all the present cases (0.5 ≤ H/S ≤ 3.45).
It would seem that in cases when the large-scale recirculations lead to upflow above the rib (i.e. for all 2-D rib cases), whether or not there is an elevated critical point with local downflow below will depend crucially on the shape of the rib. Indeed, Medjnoun et al. (2020) have shown that the rib shape can be important in setting what happens to the flow in its vicinity. In all their cases, H/S was in the range 0.8-0.87 and they only detected a downflow over the rib centre when it was of rectangular shape and unusually wide, having S/W of only 1.79. However, their PIV data did not always extend downwards enough (i.e. closer to the ribs) to detect the small-scale recirculating regions which must have existed for the larger S/W cases with rectangular ribs (i.e. cases with narrower ribs). These regions were inevitably of significantly smaller scale because of the more limited spanwise extent between a rib's two corners. The data are not inconsistent with the presence of a critical point above the rib centre, albeit too near the rib to be visible. It seems that the recirculating regions associated with corner separations and of opposite sign to the larger-scale motions aloft are relatively small for large enough S/W, and the larger scale contrary circulations above and outboard become more dominant, whereas at smaller S/W there may be insufficient room between the ribs to allow full development of the latter. This is essentially the argument of Hwang & Lee (2018), who suggested that it is S − W (i.e. the valley width) that determines the sizes and strength of the secondary flows; strictly, it should presumably be a normalised parameter (e.g. (S − W)/H) which is the relevant quantity.
Whether or not W/S is an important parameter controlling the flow just above the rib must clearly depend on rib shape; in extreme cases, like the triangular rib cases of Zampiron et al. (2020) and (some of) the ribs of Medjnoun et al. (2020), W is essentially zero at the top of the rib, so small-scale separation-driven recirculations cannot occur. The latter must always be a feature of rectangular ribs and, as Medjnoun et al. (2020) show, S/W can then be important. Our results (figure 7a) show that the peak strength of the larger scale secondary flows occurs somewhere in the range 2.7 ≤ S/W ≤ 3.3 and, by comparison with the related literature, these flows always correspond to upflow over the ribs (the high momentum pathways, (HMP, commonly mentioned in the literature)) and downflow between them (the low momentum pathways, (known as LMP)) for H/S ≥ 0.6 at least. We emphasise that in the present KC4 cases (S/W = 4, 1.25 ≤ H/S ≤ 1.75) the details of the secondary flow are essentially independent of H/S (figures 13 and 14). As mentioned in § 3.5, the swirling strength in the lower half of the flow varies very little. For H/S lower than O(1) the situation remains uncertain.
As noted in § 1, some authors have shown that changes in surface condition without any change in surface height also leads to significant secondary flows. Anderson et al. (2015), for example, used LES to study channel flow containing longitudinal strips of roughness (width W) having a higher z o (roughness length) than in the regions between them. They had H/S = 0.32, S/W = 5.2 and varied the ratio of the high to low roughness. In all cases, there was downflow over the high roughness regions (in the high momentum pathway regions) and upflow between them (in the low momentum pathway regions). Willingham et al. (2014) studied similar cases with H/S = 0.32 and in terms of the orientation of the secondary flows the results are essentially the same, for cases with 3.1 ≤ S/W ≤ 15.7. These findings are, therefore, contrary to those found in the present work and in all others using physical ribs, which all show upflow above the ribs and downflow between them, quite independently of H/S or W/S, albeit with the dominance of each matching pair of vortices varying significantly with W/S. Vanderwel & Ganapathisubramani (2015) argued that this difference in secondary flow direction arises because of the different way the spanwise inhomogeneity is imposed. This could well be the case, although the issue merits further study.

Topology and a final comment
A significant feature of the present work has been the use of topological constraints to guide the interpretation of flow visualisations. In particular, these have helped to ensure that the critical points in the cross-stream flow are identified properly and are consistent with a kinematically valid flow field. It is suggested that this should always be considered for these (and no doubt other) kinds of flows, just as recommended by Hunt et al. (1978), whether visualisations are from laboratory experiments (typically now PIV) or from computational studies (DNS or LES). Finally, it is worth emphasising again that (i) the results shown in this paper only relate to the mean flow field and (ii) the cross-plane mean velocities (i.e. V and W) are at every point very small compared with the mean axial velocity (U) at the same point. Taking point (ii) first, the ratios V/U and W/U nowhere exceed approximately 3.5 % and are usually much smaller, particularly nearer the top of the domain. The local mean flow energy ratio, (V 2 + W 2 )/U 2 , is thus extremely small everywhere. Note also that it is possible that because the cross-flow velocities are relatively so low, corner separations at the ribs might disappear at low enough Reynolds numbers. We have not explored this; recall that the apparent lack of separation seen in figure 9(b) has been found to be a result of the 'fuzzy' IBM used, rather than any Reynolds number effect. Regarding the point (i), we emphasise that in common with all turbulent flows the mean flow never actually exists. Figure 15 illustrates just how different the instantaneous flow is compared with the mean. Readers may be interested to see a video from which this snapshot was taken and which follows the flow in time; this is available at https://soton. ac.uk/engineering/about/staff/jwk.page. The dynamics of these kinds of flows compared with those in regular channels have begun to be studied (e.g. Wangsawijaya et al. 2020;Zampiron et al. 2020) and this is clearly a topic that merits further work. the force moments about that line as an axis, produced by the surface stresses, ensuring a balance between the two integrated force moments below z = d and the two above that line. The moments provided by the forces from the bottom surface stress and that on top of the rib are given by respectively. These act in opposite directions about z = d. (Note that the relationships define I 1 and I 2 as the two integrals.) Assuming that there is no stress on the sidewall (i.e. τ s = 0) and the stress on the horizontal surfaces is everywhere the same, then equating these two expressions leads simply to d/h = W/S which is in fact just the average height of the horizontal surfaces; for S/W = 2, for example, d/h = 0.5, as expected. The above assumptions are in general untrue, of course, so one has to consider both the variation of stress along the horizontal surfaces and the moments generated by the stress on sidewalls. The latter are given by for the moment of the forces below and above z = d, respectively, and again, these act in opposite directions. The appropriate balance required to ensure that z = d is the line about which there is no resultant moment (and is thus the height at which the total drag force acts) then becomes Defining I 3 and I 4 by respectively, and with some rearrangement of the two integrals over the sidewall -the second and third terms in (A3) -we eventually obtain d h = I 2 + I 4 I 1 + I 2 + I 3 .
Given the distribution of frictional stress along all the walls from the DNS data (i.e. τ o ( y), τ h ( y) and τ s (z)) for any given case, it is straightforward to determine the values of the four integrals and thus d/h from (A5). In the more general case of a surface containing three-dimensional obstacles (at any orientation), there will be additional terms to include in (A3), expressing the moments (above and below z = d) generated by, firstly, the axial frictional components on, say, any additional obstacle sidewalls parallel to the flow direction and, secondly, the axial components of all the pressure forces acting on the surfaces of the obstacles. As noted above, in such cases these latter terms will normally dominate, as shown in detail for arrays of cubes by Leonardi & Castro (2010), who used exactly this method to determine d.
As an example, we take the S/W = 4, H/S = 1.25 case KC4a in the present paper. Figure 17(a) shows the variation of axial friction along the span, which has been expanded in order to show the two sidewall stress variations at y = (S − W)/2 and (S + W)/2. Computation of the four integrals in (A3) yield I 1 = 2.28, I 2 = 0.875, I 3 = 0.602 and I 4 = 0.39. Equation (A5) then gives d/h = 0.34, which is somewhat above the geometrical average surface height of 0.25 (for this S/W = 4 case). This is an expected result; it is common for d to lie somewhere above the geometrical average surface height (Leonardi & Castro 2010) and, for cases like the present one, the difference tends to increase with decreasing S/W. The resulting velocity profiles with and without the inclusion of a zero-plane displacement are shown in figure 17(b), where it is clear that using (z − d) with d/h = 0.34 yields a much better fit to the modified log-law -better than if d/h = 0.25 were used (not shown).