Flow and coherent structures generated by a circular array of rigid, emerged cylinders in a shallow channel

Abstract

The study investigates the flow structure, dynamics of the large-scale coherent structures, drag forces and sediment entrainment mechanisms generated by a circular array of diameter D containing rigid emerged circular cylinders of diameter d placed in a smooth-bed channel of depth h under strong shallow flow conditions (D/h = 20, d/D = 0.0125 and 0.025). Eddy resolving simulations are conducted with different values of the solid volume fraction 0.025 ≤ SVF ≤ 0.1 and non-dimensional frontal area per unit volume (2.56 ≤ aD ≤ 10.2, where for the present configuration, a = SVF(1/d)4/${\rm \pi}$) and a fixed channel Reynolds number (Re_h = 10 000). The flow conditions are such that a vortex street (VS)-type of shallow wake is expected to form for a solid cylinder of diameter D. Findings are compared with previous results obtained for cases with moderately shallow conditions 2.5 ≤ D/h ≤ 3.5 and with the limiting case of flow past a solid cylinder with D/h = 20. For moderately shallow conditions, the core of the main horseshoe vortex (HV) forming around the array occupies a small fraction of the water column and its coherence is the largest in front of the array. By contrast, for very shallow conditions, multiple HVs form around the array and their cores occupy a large fraction of the water column. Moreover, the coherence of most of these HVs peaks close to the sides of the array. Only for relatively large aD values, the main HV extends over the whole upstream face of the array, similar to the limiting case of a solid cylinder. For sufficiently high aD, a secondary instability is present inside the near wake that leads to the formation of parallel horizontal vortices in the vicinity of the wake roller vortices. The changes in the wake structure with decreasing SVF and aD are qualitatively similar to those observed for cases with moderately shallow flow conditions, with the antisymmetric wake shedding mode being suppressed for aD ≤ 2.5. However, the Strouhal number associated with the shedding of wake rollers, which is still close to 0.2 for a solid cylinder with D/h = 20, can be as high as 0.4 for aD = 2.5. The paper also discusses how the steady wake and total wake lengths and the strength of the bleeding flow vary with the SVF. Simulation results show that the capacity of the flow to entrain sediment inside and around the array peaks for SVF = 0.05 and aD = 5.2, with sediment entrainment monotonically increasing with the SVF outside the array and monotonically decreasing inside the array. Despite the differences in the flow structure next to and inside the array, the variation of the mean, time-averaged streamwise drag coefficient of the solid cylinders ${\bar{C}_d}$ with aD is close to that observed for arrays with moderately shallow flow conditions. The combined drag coefficient for the array decreases with increasing flow shallowness, with the decay being stronger for relatively large values of aD.

1. Introduction

The flow past isolated arrays of solid vertical cylinders placed on a bottom surface has been made the object of numerous experimental (e.g. see White & Nepf 2007; Liu et al. 2008; Rominger, Lightbody & Nepf 2010; Zong & Nepf 2010, 2011, 2012; Rominger & Nepf 2011; Chen et al. 2012; Kim, Kimura & Shimizu 2015; Lei & Nepf 2021) and numerical investigations (e.g. see Huai, Xue & Qian 2015; Chang, Constantinescu & Tsai 2017, 2020; Koken & Constantinescu 2020, 2023; Liu et al. 2021). These studies were mainly motivated by the need to understand flow and transport processes past patches of aquatic canopies present in the main channel of rivers and/or over their floodplains. Aquatic canopies in fluvial systems alter the local flow dynamics and create coherent structures in their wake thus influencing transport and mixing processes (Nepf 2012).

Aquatic canopies in the form of patches of emerged and submerged vegetation occupying various parts of river channels and floodplains (Carr, Duthie & Taylor 1997; Sukhodolov & Sukhodolova 2010) play several important ecological roles (Jones, Lawton & Schachak 1994; Costanza et al. 1997; Corenblit et al. 2007; Gurnell 2014). They create regions with enriched nutrients, enhanced light penetration and oxygen levels, and increased plant productivity (Wilcock et al. 1999; Horppila & Nurminen 2003; Schultz et al. 2003; Schoelynck et al. 2012; Camporeale et al. 2013). Vegetation canopies reduce the transport capacity of rivers due to increased plant-induced drag (Tal & Paola 2007; Plew, Cooper & Callaghan 2008). Canopies also create flow acceleration zones near the patch boundaries which causes bed erosion (Follett & Nepf 2012; Liu & Nepf 2016) and nutrient limitation (negative feedback for patch growth). Finally, canopies generate large-scale coherent structures and increased turbulence in the river reach (Righetti & Armanini 2002; Liu et al. 2008; Zong & Nepf 2012; Pan et al. 2014; Chang et al. 2017; Koken & Constantinescu 2020, 2023). Nutrient transport within canopies and plant nutrient uptake depend on the vertical momentum, secondary flow and turbulence intensity within the canopy. Turbulence, especially near the canopy front and its sides, inhibits deposition and resuspension, leading to bed scouring (Souliotis & Prinos 2010; Pan et al. 2014; Sukhodolov & Sukhodolova 2014; Kim et al. 2015; Yang, Chung & Nepf 2016). Inside canopies, sediment deposition and nutrient accumulation are enhanced, while sediment erosion is reduced (Gurnell et al. 2005; Bouma et al. 2007; van Wesenbeck et al. 2008; Schoelynck et al. 2012; Liu & Nepf 2016; Liu et al. 2017). This benefits ecosystem productivity by promoting canopy growth (Sukhodolov & Sukhodolova 2014; Liu et al. 2017).

In the case of patches of emergent vegetation, the width of the patch is generally much larger than the flow depth, especially for vegetation growing over the floodplain at low flooding conditions. In the case of braided rivers, the ratio between the (maximum) width of these patches, D, and the mean flow depth, h, is often larger than 10. For such conditions, flow shallowness effects associated with bed friction and large values of D/h are expected to strongly affect flow, transport and sediment entrainment mechanisms around, inside and downstream of the patch. For very shallow conditions, the dynamics of the large-scale eddies generated in the wake of porous obstacles (e.g. plant stems) may be strongly affected by their interaction with the channel bed (Prinz, Brevis & Socolofsky 2012). Moreover, the vertical development of the large-scale structures forming in the wake and around the patch (e.g. horseshoe vortices) is constrained by both the free surface and the channel bed. These phenomena are expected to be qualitatively similar to those observed for shallow wakes developing behind solid obstacles (e.g. solid circular cylinders with $D/h \gg 1$) which made the object of several experimental, numerical and stability theory studies (e.g. see Chen & Jirka 1995; Grubisic, Smith & Schar 1995; Balachandar, Ramachandran & Tachie 2000; Socolofsky, Carmer & Jirka 2004; Negretti et al. 2005; Rockwell 2008; Carmer, Rummel & Jirka 2009; Zeng & Constantinescu 2017). One of the most striking differences between a deep wake and a shallow wake is that the rate of growth of the latter decays with increasing distance from the obstacle. This is accompanied by a gradual loss of coherence of the wake billows, or wake roller vortices, once their average diameter becomes much larger than the flow depth and they start behaving as quasi-two-dimensional coherent structures. As these structures are advected over the (no slip) bed surface, their energy decays mainly due to dissipation at the bed surface. At sufficiently large distances from the obstacle, the wake region contains only small-scale turbulence and its width becomes close to constant (Balachandar, Chu & Zhang 1997). The dynamics is similar to what is also observed in shallow mixing layers (Cheng & Constantinescu 2020).

Several experimental investigations and stability theory studies have shown that three main types of flow patterns are observed inside the wake with increasing flow shallowness (Chen & Jirka 1995; Socolofsky et al. 2004; Socolofsky & Jirka 2004; Carmer et al. 2009). The change from one type to the other is mainly a function of the value of a non-dimensional number called the stability number S_s = c_fD/h, where c_f is the bed friction coefficient. This number characterizes the stabilizing influence of the channel bed on the large-scale eddies developing inside the near wake (Ingram & Chu 1987; Chen & Jirka 1995). For S_s < 0.2, quasi-two-dimensional wake rollers are shed in the wake. This vortex street (VS) regime resembles the von Kármán vortex street observed for deep wakes but wake stabilization occurs with increasing distance from the obstacle. The unsteady bubble (UB) regime is observed for 0.2 < S_s < 0.5. In this regime, a large unsteady bubble forms at the back of the obstacle and a turbulent wake of sinusoidal shape is present behind the bubble, but no wake rollers form. In the steady bubble (SB) regime observed for S_s > 0.5, the recirculation bubble is quasi steady and the wake contains no growing instabilities.

In addition to the wake region, flow shallowness also affects the physics of the horseshoe vortices (HVs) forming around the base of solid obstacles. While in the case of a deep channel flow the average diameter of the main HVs is much smaller than the flow depth, the two length scales become comparable for $D/h \gg 1$. Moreover, the structure of the flow inside the HV region also changes with increasing D/h. Zeng & Constantinescu (2017) have shown that the dynamics of the HVs forming around a solid cylinder with D/h ≥ 10 resembles the one associated with the laminar breakaway regime (Kirkil & Constantinescu 2012). For such very shallow cases, the primary HV is not subject to aperiodic bimodal oscillations as generally observed for deep and moderately shallow flow conditions (e.g. see Kirkil & Constantinescu 2015). The other consequence of the cores of the HVs extending until close to the free surface is that the interactions between their legs and the separated shear layers forming on the two sides of the obstacle are not limited to the near-bed region. For very shallow conditions, the legs of these HVs can strongly affect the interactions between the two separated shear layers over the whole flow depth. Zeng & Constantinescu (2017) have shown that one consequence of these interactions is the generation of a secondary instability in the form of co-rotating, parallel, horizontal vortices forming in the vicinity of the wake rollers for cases where the VS regime is observed in the wake. Akilli & Rockwell (2002) and Carmer et al. (2009) also found that three-dimensional (3-D) effects in the wake become more important with increasing flow shallowness. Their study found that horizontal (braid) vortices form in between successively shed roller vortices and that the flow moves upwards inside the core of the roller vortices.

Despite their relevance for fluvial hydraulics, the physics of shallow wakes past porous, isolated obstacles in a channel was not yet investigated using modern experimental and numerical techniques that allow a quantitative characterization of the dynamics of the coherent structures forming around and inside the wake of the obstacle, and of the flow interactions with the individual elements of the porous obstacle. A main reason is that most experimental investigations of flow past isolated patches containing emerged plant stems or rigid cylinders were caried out in flumes that are not wide enough to accommodate patches with $D/h \gg 5$ while maintaining a fairly low blockage ratio (e.g. width to diameter ratio larger than 5), a sufficiently high channel Reynolds number and being able to perform accurate measurements over the flow depth. Eddy resolving simulations that resolve the flow past the individual elements (e.g. plant stems, cylinders) of the porous obstacle, though computationally very expensive, are not subject to the aforementioned limitations (Chang & Constantinescu 2015; Monti, Omidyeganeh & Pinelli 2019; Monti et al. 2020). Such an approach allows obtaining detailed quantitative information on the drag forces acting on individual cylinders, the flow and turbulence inside the array, and the bed shear stresses. Such simulations also allow characterizing the coherent structures, their dynamics and their interactions with the channel bed, something that is very complicated to achieve via experimental investigations (Chang et al. 2020).

The present study uses detached eddy simulation (DES) to investigate flow past a circular array of diameter D containing N rigid, emerged circular cylinders of diameter d placed in a straight and wide open channel (figure 1) for very shallow flow conditions (D/h = 20). The incoming flow is fully developed and turbulent. The channel Reynolds number is Re_h = 10 000. The circular shape of the array eliminates the angle of attack as a major variable that needs to be considered for arrays of other shapes. Moreover, close-to-circular shapes are very common for vegetation patches forming away from the river banks especially during their early growth stages (Sand-Jensen & Pedersen 2008). Though vegetated canopies are often characterized by volume heterogeneity (Hamed et al. 2017) and flexibility (Hong et al. 2022, Houseago et al. 2022), the study of a constant porosity or, equivalently, of a constant solid volume fraction (SVF) array serves as a base case with respect to which the effects of such complexities can be better evaluated for shallow flow conditions. In addition to d/D, the other relevant non-dimensional length scale is aD, where a is the frontal area per unit volume (e.g. a = Nd/(${\rm \pi}$D ²/4) for circular cylinders). In the case of an array containing circular cylinders of constant diameter, a = SVF(1/d)4/${\rm \pi}$. So, varying a, or aD, is equivalent to varying the SVF if d and D are kept constant. Given that the total force acting on the array scales with a rather than the SVF, comparison of results for circular arrays with close values of aD is generally more relevant for cases where the diameter of the solid cylinders forming the array is different (Zong & Nepf 2012). For the geometrical and flow conditions considered, the value of the stability number is S_s = 0.13. So, the VS regime is expected to be observed in the wake of the array, at least for sufficiently high SVF.

Figure 1. Sketch of the computational domain containing a circular array of solid cylinders. The emerged array of diameter D is placed in an open channel with incoming bulk velocity U and flow depth h. The polar angle ϕ is measured with respect to the symmetry plane (y/D = 0) upstream of the centre of the array.

The main goal of the present paper is to investigate the flow structure around and inside arrays of emerged circular cylinders with D/h = 20 for which the flow conditions can be classified as very shallow. A second goal is to investigate differences with respect to wakes past arrays of cylinders generated for deep and moderately shallow conditions. In particular, the present study tries to answer the following questions.

1. (i) What is the structure of the horseshoe vortex system for circular arrays with a large D/h?

2. (ii) Does the suppression of the antisymmetric shedding mode for shallow wakes in the VS regime occur at approximately the same values of SVF and aD as those observed for wakes past circular arrays with D/h < 5 (e.g. for wakes in the deep or moderately shallow regime)?

3. (iii) Is the dependence of some of the main non-dimensional variables characterizing the wake region (e.g. lengths of steady wake and total wake regions non-dimensionalized using the array diameter, D, Strouhal number associated with formation of wake rollers) and of the total streamwise drag force (e.g. mean drag coefficient acting on the cylinders forming the array, combined total drag parameter) with aD similar to that observed for arrays with D/h<5?

4. (iv) What are the main sediment erosion mechanisms around and inside circular arrays with a large D/h? How does the erosive capacity of the flow vary with increasing SVF for such arrays?

The paper is organized as follows. Section 2 briefly describes the flow solver and the DES module, the test cases, the boundary conditions and the relevant validation studies conducted using the same DES flow solver. Section 3 analyses how the SVF and d/D affect the formation of the horseshoe vortices, their position, size, circulation and dynamics. Section 4 investigates how aD affects the wake structure behind the array, 3-D effects, as well as some of the variables used to characterize the wake flow behind porous obstacles (e.g. the lengths of the steady wake and total wake regions). Section 5 investigates the bed friction velocity distributions in the mean and instantaneous flow fields and the role of the large-scale coherent structures in explaining sediment erosion inside and around the array. Section 6 discusses how the total non-dimensional streamwise drag force, the mean drag coefficient for the cylinders forming the array and the combined drag parameter for the array vary with aD. Results are also compared with those obtained for arrays with moderately shallow flow conditions. Section 7 summarizes the main findings and presents some conclusions.

2. Numerical approach

2.1. Numerical method and governing equations

The eddy-resolving simulations are performed using DES. DES combines Reynolds-averaged Navier–Stokes (RANS) modelling near solid boundaries with large eddy simulation (LES) away from these boundaries (Spalart 2009). The turbulence length scale in DES is redefined away from solid boundaries such that it becomes proportional to the cell size, similar to LES (Chang, Constantinescu & Park 2007; Rodi, Constantinescu & Stoesser 2013). The one-equation Spalart–Allmaras (SA) model (Spalart & Allmaras 1992) is used as the base turbulence model in the present DES simulations.

The incompressible Navier–Stokes equations are solved together with the transport equation for the modified eddy viscosity, $\tilde{\nu }$. The continuity and momentum equations are

(2.1)\begin{gather}\frac{{\partial {u_i}}}{{\partial {x_i}}} = 0,\end{gather}

(2.2)\begin{gather}\frac{{\partial {u_i}}}{{\partial t}} + \frac{{\partial ({u_i}{u_k})}}{{\partial {x_k}}} ={-} \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_k}}}\left[ {(\nu + {\nu_t})\left( {\frac{{\partial {u_i}}}{{\partial {x_k}}} + \frac{{\partial {u_k}}}{{\partial {x_i}}}} \right)} \right],\end{gather}

where ρ is the constant density, u_i is the resolved velocity component in the ith direction, p is the resolved hydrostatic pressure, t is the time, x_i is the coordinate in the i direction, ν_t is the subgrid-scale viscosity provided by the SA-based DES model and ν is the kinematic viscosity. The transport equation for $\tilde{\nu }$ is

(2.3)\begin{equation}\frac{{\partial \tilde{\nu }}}{{\partial t}} + {u_j}\frac{{\partial \tilde{\nu }}}{{\partial {x_j}}} = {C_{b1}}\tilde{S}\tilde{\nu } + \frac{1}{\sigma }[\boldsymbol{\nabla }\boldsymbol{\cdot }((\nu + \tilde{\nu })\boldsymbol{\nabla }\tilde{\nu }) + {C_{b2}}{(\boldsymbol{\nabla }\tilde{\nu })^2}] - {C_{w1}}\,{f_w}{\left[ {\frac{{\tilde{\nu }}}{l}} \right]^2},\end{equation}

where the turbulence length scale l is defined as the minimum between the distance to the closest wall surface, l_min, and C_DESΔ. The cell size Δ is defined as the minimum length of the cell along the three directions. The sub-grid scale viscosity is calculated as

(2.4)\begin{equation}{\nu _t} = {f_{v1}}\tilde{\nu },\end{equation}

where ${f_{v1}}$ is a damping function. The other variables and model parameters are

(2.5a)\begin{gather}\tilde{S} = \bar{S} + (\tilde{\nu }/{\kappa ^2}{l^2}){f_{v2}},\end{gather}

(2.5b)\begin{gather}{f_{v2}} = 1 - (\tilde{\nu }/\nu )/(1 + \chi {f_{v1}}),\end{gather}

(2.5c)\begin{gather}{f_{v1}} = \frac{{{\chi ^3}}}{{{\chi ^3} + C_{\nu 1}^3}},\end{gather}

(2.5d)\begin{gather}\chi = \frac{{\tilde{\nu }}}{\nu },\end{gather}

(2.5e)\begin{gather}{f_w} = g^{\prime}{\left[ {\frac{{1 + C_{w3}^6}}{{{{g^{\prime}}^6} + C_{w31}^6}}} \right]^{1/6}},\end{gather}

(2.5f)\begin{gather}g^{\prime} = r + {C_{w2}}({r^6} - r)\end{gather}

(2.5g)\begin{gather}r = \frac{{\tilde{\nu }}}{{\tilde{S}{\kappa ^2}{l^2}}},\end{gather}

where $\bar{S}$ is the magnitude of the resolved vorticity vector and κ is the von Kármán constant. The model constants are C_DES = 0.65, C_{b 1} = 0.135, C_{b 2} = 0.622, σ = 0.67, κ = 0.41, C_{v 1} = 7.1, C_{w 2} = 0.3, C_{w 3} = 2.0 and C_{w 1} = C_b1/κ ² + (1 + C_{b 2})/σ.

The governing equations were transformed to generalized curvilinear coordinates (Zeng, Constantinescu & Weber 2008). The use of a surface conforming structured mesh allows specifying no-slip boundary conditions on the surface of the cylinders. A blanking procedure was used to account for the presence of cylinders inside the computational domain. The 3-D incompressible Navier–Stokes equations were integrated using a fully implicit fractional-step method. The convective terms in the momentum equations and the transport equation for $\tilde{\nu }$ were discretized using a blend between the fifth-order accurate upwind biased scheme and the second-order central scheme. All other terms in the momentum equations and in the Poisson equation solved for the pressure were approximated using second-order central differences. The discrete momentum (predictor step) and turbulence model equations were integrated in pseudo-time using the alternate-direction-implicit approximate factorization scheme. Time integration was performed using a dual time-stepping algorithm by adding a pseudo time-derivative term in the momentum and transport equation for $\tilde{\nu }$. This term reduces to zero as convergence is achieved at each physical time step. The time discretization was second-order accurate. The mesh size in the wall-normal direction was sufficiently small to avoid the use of wall functions. Detailed descriptions of the non-hydrostatic, viscous flow solver and DES model are given in Chang et al. (2007, 2017).

The DES code was validated for several classes of open channel flows that are directly relevant for the applications discussed in the present study. Cheng & Constantinescu (2020) and Constantinescu (2014) discuss detailed validation for spatially developing shallow mixing layers in open channels with an aspects ratio of O(100) and a channel length to flow depth ratio of several hundreds. The code successfully predicted the velocity profiles and the passage frequency of the Kelvin–Helmholtz rollers at different locations as well as transition towards the equilibrium regime where the growth of the mixing layer width is negligible. Kirkil & Constantinescu (2012, 2015) discuss validation for flow past a solid cylinder placed in an open channel with $D/h \cong 1$ and 1000 < Re_D < 10⁶. Chang et al. (2017) discuss validation for flow past a circular array of emerged solid cylinders with D/h = 2.5 and 3.3. The profiles of the mean velocity and lateral velocity fluctuations along the array's centreline and the shedding frequency of the wake billows were found to be in good agreement with the experimental data of Zong & Nepf (2012) and Chen et al. (2012).

2.2. Model set-up and boundary conditions

The horizontal set-up of the present numerical simulations (figure 1) is similar to that used in the related studies of Chang et al. (2017, 2020) conducted for circular arrays with D/h < 3.5. In the present simulations, the bed is located at z = −h. The top boundary is situated at z = 0. The origin of the system of coordinates is located at the centre of the circular array. The streamwise direction is x, while the spanwise direction is y. The channel width is L_y = 140h, while the distances between the centre of the array and the inlet and outlet sections are L_{x 1} = 75h and L_{x 2} = 250h, respectively (figure 1). As for the simulations conducted by Chang et al. (2017), the solid cylinders are placed using a close-to-regular staggered pattern in the x–y planes with a random displacement of 0.5d–2.0d with respect to a regular staggered pattern. This arrangement is similar to that used in laboratory experiments of flows past vegetated canopies with rigid cylinders (e.g. see Chen et al. 2012).

The horizontal mesh contained close to three million cells and 36 points were used to resolve the flow in the vertical direction. The mesh was refined inside the array and some distance downstream of it. In the vertical direction, the mesh was refined close to the bed surface. Away from the array, the non-dimensional sizes of the cells in the horizontal and vertical directions were close to those used by Cheng & Constantinescu (2020) and Zeng & Constantinescu (2017) to perform simulations of shallow mixing layers and flow past a solid cylinder with 20 ≤ D/h ≤ 50, respectively. These simulations were also conducted in wide and long open channels with L_y/h = O(100) and (L_{x 1} + L_{x 2})/h > 200, and with Re_h ≈ 10 000. The grid spacing at the bed in the wall-normal direction was 0.0025h, corresponding to approximately one wall unit.

The turbulent flow at the channel inlet was fully developed, with a bulk channel velocity, U. A convective boundary condition was used at the outflow of the domain (Rodi et al. 2013). As in the previous simulations of flow past arrays of cylinders performed using the present DES solver, the top boundary was modelled as a shear-free horizontal rigid lid, which is justified given the relatively low channel Froude number. The solid surfaces were assumed to be smooth. The modified eddy viscosity was set equal to zero at the solid surfaces. Symmetry boundary conditions were used at the two lateral boundaries. The estimated value of the bed friction coefficient in the incoming flow was c_f = 0.0064. The time step was dt = 0.0125h/U. Its value was sufficiently low to resolve the range of energetic frequencies associated with the dynamically important coherent structures in the flow (e.g. horseshoe vortices, vortices shed in the wake of the solid cylinders, vortex tubes in the array's separated shear layers). In all critical regions of the flow containing large-scale eddies, the Courant number was much less than one. Mean flow and turbulence statistics were collected over a time interval of approximately 1000h/U.

The number of cylinders in the array, N, and the diameter of the solid cylinders, d, were the two main array parameters varied in the simulations. For a circular array containing circular cylinders, SVF = Nd ²/D ² and aD = SVF*(D/d)*4/${\rm \pi}$. Thus, when comparing simulations with the same value of D/d, studying the effects of varying aD and SVF is equivalent. The two (non-dimensional) length scales controlling the total drag force acting on the array are 1/(aD) and d/D. The main geometrical and flow variables for the five simulations are given in table 1. Three of the simulations were conducted with constant d (=0.0125D) and varying number of cylinders (N = 120, 320 and 641). The corresponding values of the SVF were 0.1, 0.05 and 0.025, respectively. A fourth simulation was conducted with larger cylinders (d = 0.025D) and N = 160 such that SVF = 0.1. The range of SVFs was comparable to that (0.024 ≤ SVF ≤ 0.23) in the simulations of Chang et al. (2017) conducted with D/h = 2.5 and 3.33. Also included are results from a fifth simulation for the limiting case of flow past a solid cylinder (SVF = 1.0).

Table 1. Main geometrical and flow variables of the test cases (N is the number of cylinders in the array; d is diameter of the solid cylinders; D is the diameter of the circular array; h is the flow depth, U is the channel bulk velocity, Re_h = Uh/ν, Re_D = UD/ν, Re_d = Ud/ν, ν is the molecular viscosity, aD is the non-dimensional frontal area per unit volume for the array; U ₁ is the centreline streamwise velocity at the upstream end of the steady wake region; U ₂ is the maximum streamwise velocity of the flow on the sides of the array; ${L_{sw}}$ and ${L^{\prime}_{sw}}$ are the steady wake lengths defined based on mean velocity and spanwise normal stresses, respectively; ${L^{\prime}_w}$ is the total wake length; St_D is the Strouhal number defined with U and D corresponding to the most energetic frequency in the wake of the array; C_D is the mean solid cylinder streamwise drag coefficient; ${\varGamma _D} = {C_D}aD/(1 - \textrm{SVF})$ is the combined drag parameter for the array.

3. Horseshoe vortex system

The Q criterion is used in figure 2 to visualize the HV system in the mean flow associated with the array. The variable Q is defined as $Q = 0.5({\parallel} \varOmega {\parallel ^2} -{\parallel} S{\parallel ^2})({h^2}/{U^2})$ , where ${\parallel} \varOmega\!\!\parallel$ is the magnitude of the rotation rate tensor and ${\parallel} S\!\parallel$ is the magnitude of the rate of strain tensor calculated with the resolved velocities. Though horseshoe vortices form in all four simulations conducted with an array of cylinders, there are several important differences between the SVF = 0.1 case that has the largest value of aD and the other cases for which aD ≤ 5.12. It is only in the SVF = 0.1 case (figure 2b) that the primary HV forming in the immediate vicinity of the array remains strongly coherent until the symmetry plane corresponding to a polar angle ϕ = 0°. The polar angle is defined with respect to the symmetry plane y/D = 0 upstream of the centre of the array (figure 1). Its core and legs extend from the bed up to approximately 0.2h from the top boundary. For all the other cases shown in figure 2(c–e), no HVs are present in the mean flow in between the symmetry plane and |ϕ| = 30°–40° except for the small HVs generated around the base of the cylinders directly exposed to the incoming flow. These vortices are clearly captured in figure 2(e,f) for the SVF = 0.025 case. Despite the absence of strongly coherent HVs in front of the array, the region situated close to the symmetry plane and the front face of the array is still characterized by high turbulent kinetic energy (TKE) levels (e.g. k/U ² > 0.01) for aD > 5. This amplification of the TKE is mainly due to the smaller vortices that are shed from the region where the boundary layer separates (see TKE distributions in figure 3b,c). The size of this region and the peak TKE values inside it decay monotonically with decreasing aD (figures 3a–d). The turbulence amplification in the ϕ = 0° plane becomes negligible in the SVF = 0.025 case (figure 3d) where the HVs forming around the base of the cylinders directly exposed to the incoming flow are only subject to small oscillations around their mean flow position.

Figure 2. Visualization of the horseshoe vortex system using a Q isosurface (Q = 0.1). (a) SVF = 1.0 mean flow; (b) SVF = 0.1 mean flow; (c) SVF = 0.1(2d) mean flow; (d) SVF = 0.05 mean flow; (e) SVF = 0.025 mean flow; ( f) SVF = 0.025, instantaneous flow. BV indicates a bottom-attached vortex. The colour contours represent the distance from the top boundary, z/D. Note the large-scale hairpins generated along the legs of the largest HV, V1 in panel (f).

Figure 3. Visualization of flow structure inside the horseshoe vortex system in selected polar planes. (a) SVF = 0.1; (b) SVF = 0.1(2d); (c) SVF = 0.05; (d) SVF = 0.025; (e) SVF = 1.0. Top panels show the mean out of plane vorticity, ω_nh/U, 2-D streamlines for ϕ = 70°. Middle panels show the TKE, k/U ², for ϕ = 70°. Bottom panels show the TKE for ϕ = 0°. The vertical dashed line indicates the boundary of the array. The interior of solid cylinders is represented in grey. The horizontal to vertical aspect ratio in all frames is 1:2.5. The radial distance measured from the origin of the array is r.

Results in figure 3 show that the overall level of TKE amplification inside the HV system region scales with aD rather than the SVF. This is not surprising as the pressure on the cylinders directly exposed to the incoming flow and the magnitude of the adverse pressure gradients generated away from the frontal face of the array also scale with aD. For example, the mean non-dimensional pressure distributions on the top boundary (figure 4) show that the regions of high adverse pressure gradients are small and localized next to the upstream faces of the cylinders exposed to the incoming flow in the SVF = 0.025 case. In contrast, in the SVF = 0.1 case, the region of highest pressure magnitude extends from ϕ = −30° to ϕ = 30°, with strong adverse pressure gradients generated upstream of it. So, for high aD, the region of increased pressure is associated with the array as a whole (e.g. with the porous obstacle) rather than with the individual cylinders forming the array. The coherence of the HV system mainly scales with the magnitude of the adverse pressure gradients generated away from the obstacle (Kirkil & Constantinescu 2015).

Figure 4. Mean pressure, p/ρU ², at the top boundary (z/h = 0). (a) SVF = 0.1; (b) SVF = 0.025.

For the SVF = 1.0 and 0.1 cases, where the main HV remain strongly coherent for 0° < |ϕ| < 90°, the TKE levels are smaller in the SVF = 0.1 case compared with the SVF = 1.0 case (e.g. compare TKE distributions in figure 3a,e). This is a direct result of the decay of the adverse pressure gradients generated away from the upstream face of the array with decreasing SVF and aD. One interesting result is that the circulation of main horseshoe vortices V1 to V4 in the SVF = 0.1 case in the region where the coherence of these vortices is high (30° < |ϕ| < 90°) is comparable to that of the corresponding vortices V1 and V2 in the SVF = 1.0 case. For both cases, the interaction of the main HVs with the bed generates the formation of bottom-attached vortices (see BV vortices in figure 3a,e).

The TKE distributions in the ϕ = 70° plane are representative for the region where the coherence of the HVs is the largest in the shallow flow simulations conducted with an array of cylinders. The number of vortices part of the HV system is not the highest for the case where coherence of these vortices, as measured by the size of their cores and their circulation, is the largest (see vorticity distributions in figure 3) and does not vary monotonically with aD or SVF. For D/h = 20, the number of HVs is equal to five in the SVF = 0.01(2d) case. This number decreases to four in the SVF = 0.1 case and to three in the SVF = 0.05 and 0.025 cases (figure 2). As discussed later, the main reason for the TKE amplification in this region is their advection towards the HV situated the closest to the array.

In the region of highest coherence of the HVs, the cores of the largest vortices extend approximately 0.6h in the vertical direction in the SVF = 0.1 case. This vertical length reduces to approximately 0.25h–0.3h in the SVF = 0.05 and SVF = 0.025 cases (see vorticity distributions in figure 3). By comparison, the core of the main turbulent HV forming around solid cylinders with D/h ≈ 1 was approximately 0.25h (Kirkil & Constantinescu 2015). Large differences are also observed for the width of the HV system region. In the ϕ = 70° plane, HVs are present up to 7.5h (0.375D) away from the face of the array in the SVF = 0.1 case. This distance reduces to 3h (0.15D) in the SVF = 0.025 case. The value predicted for the SVF = 0.1 case is not very different from that (9h or 0.45D) predicted in the SVF = 1.0 case. By comparison, for moderately shallow conditions, Chang et al. (2017) found that the HV system region extended 0.25D away from the frontal face of the obstacle for a solid cylinder with D/h = 2.5 and approximately 0.2D for an array with SVF = 0.23, aD = 9.6 and D/h = 2.5. No HVs associated with the array were present in the other cases with aD < 4.0 investigated by Chang et al. (2017).

The presence of multiple vortices of comparable coherence inside the HV system and the fact that their coherence peaks for relatively high polar angle magnitudes are important differences with the HV system forming around arrays with D/h < 4 but with similar SVF and aD values (Chang et al. 2017). The structure of the HV system is also different from that observed in the case of flow around a solid circular cylinder with D/h < 1 where the HV system contains only one primary vortex that remain strongly coherent around the whole upstream face of the array. These differences are all induced by flow shallowness.

In the D/h = 20 cases, the dynamics of the cores of the main HVs resembles that of a laminar HV system in one of the breakaway sub-regimes containing two or more vortices (Kirkil & Constantinescu 2012). New HVs are shed from the region where the boundary layer separates while the cores of the other vortices advance towards the array. The core of V2 eventually merges with the core of the vortex situated the closest to the array (V1). As opposed to the case of solid cylinders with D/h > 10 (Zeng & Constantinescu 2017, there is no well-defined period associated with such a cycle in the present simulations conducted with aD < 10. This dynamic is also different from that observed for the turbulent HV system forming around solid bluff bodies for moderately shallow and deep flow ($D/h \ll 1$) conditions where the primary HV is subject to aperiodic bimodal oscillations and its interactions with small HVs shed from the region where the incoming boundary layer separates are fairly minor (Kirkil & Constantinescu 2015). One similarity with the dynamics of turbulent HVs forming around solid obstacles with D/h < 1 is that large-scale hairpins develop along the legs of the most coherent HVs in the instantaneous flow fields. For example, these hairpins can be observed in figure 2( f) for the SVF = 0.025 case, where their heads can raise up to approximately 0.4h above the bed surface. Only in the SVF = 0.1 case, the legs of the vortex situated the closest to the array, V1, strongly interact with the eddies shed in the separated shear layers over most of the flow depth. In all the other D/h = 20 cases, the HVs are situated relatively far from the array and/or their legs loose rapidly their coherence. This explains the weak interactions between the HVs and the separated shear layers in these cases.

4. Array and near wake regions

4.1. Wake instabilities and coherent structures

The distributions of the instantaneous vertical vorticity in horizontal planes show that a VS-type wake is present in the simulations conducted with SVF ≥ 0.05 and that the distance from the back of the array where roller vortices form increases with decreasing SVF and aD (e.g. see figure 5a,b). By contrast, wake billows do not form in the SVF = 0.025 case (figure 5c) where a turbulent wake undergoing undulatory motions forms behind a region of low velocities at the back of the array. This is somewhat similar to a UB-type wake behind a solid obstacle except that there is no recirculation bubble at the back of the array. So, the transition between the VS mode, where wake rollers form, and a wake mode, where the separated shear layers cannot interact, occurs at SVF values slightly smaller than 0.05. This value is very close to the threshold value observed by Chang et al. (2017) for moderately shallow conditions (the antisymmetric wake shedding mode was absent in their simulation conducted with SVF = 0.034) and also to the experimental results of Zong & Nepf (2011) and Chen et al. (2012). Based on this result, one can conclude that changes in the wake structure are not significantly affected by flow shallowness.

Figure 5. Non-dimensional vertical vorticity, ω_zh/U, in the instantaneous flow close to the top boundary (z/h = −0.1). (a) SVF = 1.0; (b) SVF = 0.1; (c) SVF = 0.025.

However, this is not the case for the dominant shedding frequency of the roller vortices, f. For moderately shallow conditions and sufficiently large SVF, for the antisymmetric wake shedding mode to be dominant, Chang et al. (2017) predicted values of the Strouhal number St =  f D/U between 0.22 (SVF = 1) and 0.2 (SVF = 0.092), consistent with values measured in experiments conducted for porous cylinders with similar values of D/h. In the present simulations, St = 0.22 for SVF = 1.0, which is very close to the range of values (St = 0.21–0.23) measured experimentally by Chen & Jirka (1995) for very shallow flow conditions with a VS-type wake forming behind the solid cylinder. The St = 0.22 value is slightly smaller than the values (St = 0.25–0.27) predicted for solid cylinders with D/h ≈ 1 and 16 000 < Re_D < 500 000 by Kirkil & Constantinescu (2015). Rather than being constant or slightly decreasing with decreasing aD and SVF, the Strouhal number in the simulations conducted with D/h = 20 increases monotonically with decreasing aD such that St ≈ 0.25 for aD ≈ 10 (SVF = 0.1 case) and St ≈ 0.36–0.4 for aD ≈ 5. So, flow shallowness has a major effect on the non-dimensional period associated with the interactions of the separated shear layers forming on the two sides of the array that result in the generation of wake rollers.

In the SVF = 1.0 (figure 5a) and SVF = 0.1 (figure 5b) cases, the legs of the main horseshoe vortices can get very close to the separated shear layers and interact with them (see figures 2a,b and 5a,b). For solid cylinders with 10 ≤ D/h ≤ 25, Zeng & Constantinescu (2017) observed the formation of parallel corotating horizontal vortices close to the roller vortices. This 3-D near-wake instability is driven by the strong interaction of the horseshoe vortices with the separated shear layers. The instability is also present in the SVF = 1 simulation (figure 6a) where these vortices are situated slightly behind and on the outer side, with respect to the symmetry plane, of the roller vortices. The length of these vortices can attain 10h or 0.5D and their spacing in the spanwise direction is of the order of 1h to 1.5h. The streamwise velocity inside the cores of these vortices is smaller than that in the surrounding wake flow (figure 7a). Given that the cores of these vortices extend over approximately 80 % of the flow depth, they can play an important role in vertical mixing inside the near wake. Moreover, these vortices can have a strong effect on the transport of fine sediments and nutrients by advecting near-bed fluid towards the top boundary (e.g. free surface for a patch of vegetation in an open channel). Despite the common belief that shallow wakes are quasi-two-dimensional, strong 3-D instabilities can be generated because of flow shallowness. The present study also shows that this instability is also present in the case of arrays of cylinders provided that the non-dimensional frontal area per unit volume is sufficiently large such that strong interactions take place between the legs of the most coherent horseshoe vortices and the separated shear layers. Such streamwise vortices form in the vicinity of the roller vortices only in the SVF = 0.1 simulation for which aD = 10.2 (figures 6b and 7b) but not in the other simulations with aD ≤ 5.2. So, the threshold value of aD is between 5 and 10. As expected, the average circulation of the streamwise vortices in the SVF = 0.1 case is smaller compared with that of the corresponding vortices forming in the SVF = 1.0 simulation but the cores of these vortices still occupy a large fraction of the channel depth.

Figure 6. Visualization of coherent structures in the instantaneous flow using a Q isosurface (Q = 0.002). (a) SVF = 1.0; (b) SVF = 0.1. The black lines in panels (a) and (b) show the position of the x/h = constant planes shown in figure 7.

Figure 7. Streamwise vorticity, ω_xh/U, (top panel) and streamwise velocity, u/U, (bottom panel) in an instantaneous flow field. (a) SVF = 1.0, x/h = 90; (b) SVF = 0.1, x/h = 120. The red arrows point towards patches of high streamwise vorticity associated with some of the streamwise-oriented vortices visualized in figure 6.

4.2. Flow inside the array

Despite the 2-D geometry, the flow inside the array is not quasi-two-dimensional. In the case of patches of vegetation, the secondary flow and vertical velocity fluctuations inside the array have important implications for trapping heavy metals and suspended sediments as well as for the transport of nutrients (Zong & Nepf 2011). The main 3-D effects are due to the downflow forming parallel to the upstream face of the array and, for sufficiently high SVFs, the rotational motion of the core of the horseshoe vortex situated in the immediate vicinity of the array. In the SVF = 0.1 case (figure 8a), which is the only case where a strongly coherent horseshoe vortex forms along the whole frontal face of the array, 3-D mean-flow effects are fairly limited and associated with the fluid moving with the downflow that penetrates slightly inside the array. As aD decreases, so do the coherence and diameter of the horseshoe vortices. As for the case of circular arrays with moderately shallow flow conditions (Chang et al. 2017), the flow entering the array has a strong velocity component oriented towards the bed. Given the incompressibility constraint, this downflow is compensated by strong upwelling taking place in the wake of the cylinders situated close to the frontal face of the array (e.g. see figure 8b for SVF = 0.025 where the peak vertical velocity magnitudes are of the order of 0.25U). For similar values of aD, the upwelling and downwelling motions inside the array weaken with increasing flow shallowness (e.g. with increasing D/h).

Figure 8. Non-dimensional mean vertical velocity, w/U, at mid depth (z/h = −0.5). Regions with |w/U| < 0.04 are blanked. (a) SVF = 0.1; (b) SVF = 0.025.

The level of amplification of the TKE inside the array is directly related to the formation and shedding of vortices in the wake of the cylinders and the velocity magnitude in the flow approaching the cylinders. Close to the top boundary, the largest TKE is always observed near the front of the array (figure 9). However, for sufficiently low aD values, the largest TKE values close to the bed surface occur some distance behind the front face (e.g. see figure 9d). This effect is probably related to the downflow penetrating the array.

Figure 9. Non-dimensional TKE, k/U ², close to the top boundary (x/h = −0.1, top half) and the channel bed (z/h = −0.9, bottom half). (a) SVF = 1.0; (b) SVF = 0.1; (c) SVF = 0.05; (d) SVF = 0.025. The red circle shows the boundary of the array region.

Similar to arrays of cylinders with D/h = O(1) and to arrays of long cylinders (Chang & Constantinescu 2015; Chang et al. 2017), the horizontal flow inside the array is not oriented along the streamwise direction. Rather, the flow entering the array diverges away from the symmetry plane, as shown by the mean streamline patterns in figure 10. This effect weakens with decreasing aD. The flow inside the array decelerates due to the drag induced by the cylinders. Table 1 reports the values of the bleeding flow velocity at the back of the array near the symmetry plane, U ₁. The relative deceleration of the flow over the length of the array can be approximated as (U − U ₁)/U. This value decreases with increasing aD from approximately 0.9 for aD ≈ 10 to approximately 0.65 for aD ≈ 2.5. Comparison with the simulations reported by Chang et al. (2017) shows that, for comparable values of aD, increasing the flow shallowness results in a slight reduction (e.g. by approximately 10 % as D/h increases from approximately 3 to 20) of the flow deceleration inside the array. This is mainly due to a slight decrease of the mean drag force acting on the cylinders with increasing flow shallowness, as will be discussed in § 6.

Figure 10. 2-D mean flow streamline patterns close to the top boundary (z/h = −0.1). (a) SVF = 1.0; (b) SVF = 0.1; (c) SVF = 0.05; (d) SVF = 0.025. The red circle shows the boundary of the array region.

4.3. Mean flow and turbulent kinetic energy inside the wake

As for flow past emerged arrays with moderately shallow flow conditions, the near wake contains a recirculation bubble in the mean flow for cases where strong rollers are shed in the wake. The position and size of this bubble are mainly a function of aD. No bubble forms for cases when the antisymmetric vortex shedding mode is very weak or absent (e.g. see figure 10c,d). This is confirmed by comparison of the two SVF = 0.1 simulations. A fairly large bubble forms in the SVF = 0.1 case (aD = 10.2). No bubble forms in the SVF = 0.1(2d) case (aD = 5.1) and in the SVF = 0.05 case (aD = 5.1). Based also on the data of Chang et al. (2017), one can estimate the threshold value of aD to be close to 6 for a recirculating bubble to form in the wake of circular emerged arrays.

Comparison of the present simulation results with those of Chang et al. (2017) shows that the vertical non-uniformity of the bubble increases significantly with D/h. For D/h = 20 and SVF = 1.0, the length of the recirculation bubble increases from 0.5D near the top boundary to 0.7D near the bed. By contrast, the length of the same eddy was 0.7D at all flow depths for D/h = 2.5 (Chang et al. 2017). The differences increase when cases with aD ≈ 10 are compared. For D/h = 20 and SVF = 0.1, the detached bubble starts at x = 2.25D near the bed and at x = 2.5D, near the top boundary. The bubble length is close to 1.2D near the top boundary and 1.5D near the bed. Meanwhile, for D/h = 2.5 and aD = 9.6, the bubble length was 1.1D at all flow depths. So, increasing the flow shallowness decreases the two-dimensionality of the near-wake mean flow. The non-uniformity over the vertical direction of the flow and TKE decreases with decreasing aD especially once the roller vortices become very weak (SVF = 0.05 case) or do not form (SVF = 0.025 case). This can be seen by comparing the TKE distributions near the top boundary and near the bed in figure 9.

The steady wake is defined as the region where the mean streamwise velocity behind the array decreases (Zong & Nepf 2012). Its length is measured along the symmetry plane starting at the back of the array (x/D = 0.5). The end of the steady wake region approximatively corresponds to the length of the separated shear layers for the array. As the flow starts recovering inside the wake, the mean streamwise velocity in the symmetry plane starts increasing toward the free stream value. Zong & Nepf (2012) defined the length of the total wake region as the streamwise distance between the back of the array and the streamwise location where the rate of increase of the streamwise centreline velocity becomes very small (e.g. smaller than a threshold value). The longitudinal profiles of the streamwise velocity along the centreline (y/D = 0) in figure 11(a) show that a relative minimum associated with the end of the steady wake region is only observed in the SVF = 0.1 simulation. Following Chang & Constantinescu (2015), we use the longitudinal profiles of the spanwise normal stress $\overline {v^{\prime}v^{\prime}} $ along the centreline to identify the end of the steady wake and total wake regions. The end of the steady wake is defined as the location where $\overline {v^{\prime}v^{\prime}} $ starts increasing monotonically. The end of the total wake corresponds to the location where $\overline {v^{\prime}v^{\prime}} $ reaches its maximum (figure 11b). The lengths of the steady and total wake regions are measured from the back of the array and are denoted ${L^{\prime}_{sw}}$ and ${L^{\prime}_w}$, respectively. This definition of ${L^{\prime}_w}$ has the advantage that it avoids defining an arbitrary threshold value for the streamwise velocity gradient along the centreline to identify the end of the total wake region. The streamwise variations of the centreline velocity and $\overline {v^{\prime}v^{\prime}} $ are close to identical in the SVF = 0.1(2d) and SVF = 0.05 simulations for which aD ≈ 5.1, which shows that aD rather than the SVF is the main parameter. Both variables are monotonically increasing with aD. Compared with arrays with 2.5 ≤ D/h ≤ 3.3, the wake recovery is slower in the simulations conducted with D/h = 20. Moreover, the change from a streamwise velocity profile containing a relative minimum to one displaying a fairly monotonic increase over the whole wake occurs between aD = 5 and aD = 10 for D/h = 20 compared with aD ≈ 2 for 2.5 ≤ D/h ≤ 3.3.

Figure 11. Longitudinal profiles of: (a) mean streamwise velocity, $\bar{u}/U$; and (b) spanwise normal stress, $\overline {v^{\prime}v^{\prime}} /{U^2}$, at the top boundary. The profiles are shown along the centreline, y/D = 0. The grey region shows the extent of the array region (−10 < x/D < 10). The dashed arrows in panels (a) and (b) show the end of the steady wake region based on the predicted values of $\; {L_{sw}}/h$ and ${L^{\prime}_{sw}}/h$, respectively. The solid arrows in panel (b) show the end of the total wake region based on the predicted values of ${L^{\prime}_w}/h$.

The variations of ${L^{\prime}_{sw}}/D$ and ${L^{\prime}_w}/D$ with the ratio between the bleeding velocity at the back of the array, U ₁, and the maximum streamwise velocity in the regions of flow acceleration on the two sides of the array, U ₂, show a similar trend for the present simulations (D/h = 20) and for the simulations conducted by Chang et al. (2017) for arrays with 2.5 ≤ D/h ≤ 3.3 (figure 12). However, for comparable aD values, ${L^{\prime}_w}/D$ is at least 30 % larger in the simulations conducted with D/h = 20.

Figure 12. Variation of the length of the near-wake region with the velocity ratio U ₁/U ₂. (a) Steady wake length, ${L^{\prime}_{sw}}/D$; (b) total wake length, ${L^{\prime}_w}/D$. In addition to the present test cases (D/h = 20) with d = 0.0125D (open red squares) and d = 0.025D (solid red squares), results are also included for DES conducted for an array with D/h = 2.5–3.3 (solid green diamonds) by Chang et al. (2017).

5. Bed friction velocity and volumetric flux of entrained sediment

For flows past arrays of surface-mounted cylinders, large-scale turbulent eddies are a main contributor to erosion inside and around the array. Sediment can be entrained in regions where the mean bed friction velocity, ${\bar{u}_\tau }$, is less than the threshold value for sediment entrainment provided that the turbulence is strong near the bed (Yang et al. 2016). To characterize in an average way the capacity of the flow to entrain sediment at a certain location, one needs to consider not only the value of ${\bar{u}_\tau }$ but also the deviations of the instantaneous bed friction velocity with respect to the time-averaged value. The effect of these fluctuations on the capacity of the flow to entrain sediment can be quantified in an approximative way by the standard deviation of the bed friction velocity, $u_\tau ^{SD}$ (Sumer et al. 2003). This is also why many RANS codes used to perform simulations with a deformable loose bed employ an ‘enhanced’ value of the bed shear stress, or equivalently of the bed friction velocity (${u_{\tau e}} = {\bar{u}_\tau } + C^{\prime}u_\tau ^{SD}$), in the equation used to estimate the local flux of entrained sediment, where the value of the constant C′ is of the order of one (Cheng, Koken & Constantinescu 2018). Following Chang et al. (2017), the magnitude of the local bed friction velocity in the instantaneous flow is calculated as$\; {u_\tau }/U = {((1/Re)(({u_m}/U)/({n_1}/h)))^{1/2}}$, where u_m is the instantaneous velocity magnitude in a plane parallel to the bed situated at a normal distance n ₁ from the bed. The distance n ₁ is sufficiently small for the plane to cut through the viscous sub-layer. The same definition is used to calculate ${\bar{u}_\tau }$ and $u_\tau ^{SD}$ with u_m being the mean velocity magnitude and the square root of the turbulent kinetic energy, respectively.

Similar to the case of solid circular cylinders placed on a surface with D/h ≤ O(1), most of the entrainment in the SVF = 1.0 case occurs beneath the regions of strong flow acceleration on the sides of the cylinder and beneath the cores of the wake roller vortices (see bed friction velocity distributions in figures 13a and 14a). When non-dimensionalized with the diameter of the cylinder, the width of the region of high u_τ/U forming on the side of the cylinder and the mean diameter of the regions of high u_τ/U forming beneath the cores of the wake rollers are approximately 30 % smaller in the simulation conducted with D/h = 20 compared with the values in the corresponding D/h = 2.5 solid cylinder simulation discussed by Chang et al. (2017). The peak bed friction velocity values inside these regions are also smaller in the D/h = 20 simulation. This means that the flux of entrained sediment from the bed around the solid cylinder decreases with increasing flow shallowness when normalized by the area of the cylinder.

Figure 13. Non-dimensional bed friction velocity magnitude, u_τ/U, in the instantaneous flow. (a) SVF = 1.0; (b) SVF = 0.1; (c) SVF = 0.05; (d) SVF = 0.025.

Figure 14. Non-dimensional bed friction velocity magnitude, $\overline {{u_\tau }} /U$, in the mean flow. (a) SVF = 1.0; (b) SVF = 0.1; (c) SVF = 0.05; (d) SVF = 0.025.

Similar trends are observed in the variation of the bed friction velocity in the vicinity and in the wake of the array for cases with SVF < 1. The main difference between the simulations performed with D/h = 20 and those with 2.5 ≤ D/h ≤ 3.3 is the relatively larger amplification of the mean and instantaneous bed friction velocity on the outer side of the separated shear layers (e.g. see figures 13b and 14b for SVF = 0.1) where strongly coherent horseshoe vortices are present in some of the simulations (figure 2). In general, the cores of the HVs have a larger capacity to locally amplify the bed shear stress with increasing flow shallowness. Their increased size relative to the flow depth in the D/h = 20 simulations allows them to advect high streamwise velocity fluid from the immediate vicinity of the top boundary towards the bed. This is the main mechanics for the amplification of u_τ in this region (e.g. see figure 13b). Another difference with the simulations of Chang et al. (2017) conducted with 2.5 ≤ D/h ≤ 3.3 is that the wake rollers and their trajectories are more regular in the D/h = 20 simulations. This also explains the fairly regular staggered pattern of the regions of high u_τ/U inside the near wake in figure 13(a–c). For cases where the antisymmetric wake shedding mode is very weak (SVF = 0.05 case) or absent (SVF = 0.025 case), most of the sediment entrained in the wake of the array originates in the thin elongated regions starting beneath the separated shear layers (figures 13c,d and 14c,d). The passage of the wake rollers and the movement of the HVs are the main mechanisms for the amplification the bed friction velocities outside of the array. Overall, the D/h = 20 simulations are characterized by rather regular trajectories of the larger-scale eddies originating inside the separated shear layers.

The regions of high $u_\tau ^{SD}/U$ in figure 15 correlate well with the regions of high TKE in figure 9. As the SVF and aD decrease, $u_\tau ^{SD}/U$ decreases around the upstream face of the array and inside its wake (figure 15), which is consistent with the observed decrease in the coherence of the HVs, the increase in the length of the separated shear layers and the decrease in the strength of the antisymmetric wake shedding mode. The positions and sizes of the regions of high $u_\tau ^{SD}/U$ inside the array also vary with decreasing SVF and aD. These changes are driven by the positions and number of cylinders that generate strong wake vortex shedding. A decrease of the SVF increases the percentage of the cylinders in the array that generate strong vortex shedding.

Figure 15. Non-dimensional root mean square of the bed friction velocity magnitude, $u_\tau ^{SD}/U$. (a) SVF = 1.0; (b) SVF = 0.025.

One way to assess in a more quantitative way the effect of varying aD and SVF on the sediment entrainment capacity of the flow inside the array and over the whole region affected by the presence of the cylinders inside the channel is to calculate the time-averaged flux of entrained sediment, $\overline {{E_\tau }} $. This variable is defined as the time-averaged value of the flux evaluated from the instantaneous flow fields, ${E_\tau }(t)$. Sediment entrainment formulae for non-cohesive sediment generally assume that ${E_\tau }$ is proportional to ${(u_\tau ^2 - u_{\tau cr}^2)^\gamma }\; $with γ > 1, where ${u_{\tau cr}}$ is the critical value for sediment entrainment for a given mean particle diameter, d ₅₀. The value of γ is 1.5 in van Rijn's (1984) formula which is often used to estimate the volumetric flux in eddy-resolving simulations with a movable bed (e.g. see Chou & Fringer 2010). In non-dimensional form, the volumetric flux of sediment entrained from the bed can be written as

(5.1)\begin{equation}{E_\tau }(tU/D) = \frac{1}{{{D^2}}}\int_{A^{\prime}} {\frac{{{{(u_\tau ^2 - u_{\tau cr}^2)}^{1.5}}\,\textrm{d}A}}{{{U^3}}}} ,\end{equation}

where A is the area over which entrainment is calculated and A′(t) is the region of the bed where ${u_\tau }(t,x,y) > {u_{\tau cr}}$. Figure 16 includes three values of $\overline {{E_\tau }} $ calculated assuming the area A is the whole channel bed, the bed region situated inside the array and the bed region situated outside the array. Results are reported for u_τcr/U = 0.065 (critical bed shear stress τ_cr/ρU ² = 0.0042) such that no sediment is entrained at large distances from the array. The values of $\overline {{E_\tau }} $ inside and outside the array are not close in the two simulations conducted with aD ≈ 1. This means that the entrainment flux is a function of both SVF and the relative size of the cylinders forming the array. For constant aD, the fluxes entrained inside and outside the array are larger in the simulation conducted with cylinders of larger diameter (e.g. in the SVF = 0.1(2d) case compared with the SVF = 0.05 case). For constant SVF, the flux of sediment entrained outside the array decreases with increasing diameter of the solid cylinder, while the opposite trend is observed for the flux of sediment entrained inside the array. For constant d/D, the flux of sediment entrained outside the array increases monotonically with the SVF and aD. The main reasons are the increased size of the region of strong flow acceleration forming on the sides of the array for |ϕ| > 30° and the increased near-bed velocity magnitudes inside this region. However, the velocity of the flow entering the array increases with decreasing SVF and aD. The opposite is true for the streamwise rate of decrease of the velocity in between the cylinders. These varying trends explain the monotonic decay of $\overline {{E_\tau }} $ with increasing SVF in figure 16 for the simulations conducted with constant d. Given the opposite variations of $\overline {{E_\tau }} $ inside and outside the array, the total sediment flux entrained over the whole bed reaches a minimum around SVF = 0.05. This value of the SVF roughly corresponds to the transition between a regime where counter-rotating roller vortices form and a regime where such vortices do not form.

Figure 16. Non-dimensional flux of entrained sediment, $\overline {{E_\tau }} $, as a function of the solid volume fraction, SVF. The red line and red solid symbols show $\overline {{E_\tau }} $ calculated inside the array. The blue line and blue solid symbols show $\overline {{E_\tau }} $ calculated outside the array. The black line and black solid symbols show $\overline {{E_\tau }} $ calculated over the whole channel bed. The open symbols correspond to the SVF = 0.1(2d) case. Results are plotted for a critical bed friction velocity, u_τc/U = 0.065.

6. Drag forces

For moderately shallow conditions and a sufficiently large value of aD such that the antisymmetric vortex shedding mode is strong, Chang et al. (2017) found that most of the total drag is due to drag forces acting on the cylinders directly exposed to the incoming flow and those situated in their immediate vicinity. Present results show this is also the case for very shallow flow conditions if aD ≥ 10. As aD decreases to values where the antisymmetric shedding mode is very weak or absent, all cylinders are subject to significant drag forces.

Figure 17 shows a plot of the mean streamwise drag force on the solid cylinders normalized by the drag force acting on an isolated cylinder. This information is especially valuable in the case of patches of aquatic vegetation where the movement of the plant stems, and the capacity of the plant stems to avoid dislocation from the bed are a function of the magnitude of the drag force acting on them. The distribution of the normalized streamwise drag force is qualitatively different in the SVF = 0.1 case (aD = 10.2) with respect to the other cases. In the SVF = 0.1 case, more than 60 % of the cylinders in the array are subject to fairly negligible drag forces defined as forces that are less than 20 % of the force acting on an isolated cylinder of same diameter for the same value of Re_d. Interestingly, most of the cylinders directly exposed to the incoming flow are also subject to streamwise drag forces smaller by approximately 20–40 % compared with the force acting on an isolated cylinder. This is because the cylinders situated close to the frontal face of the array contribute to the deceleration of the flow as it approaches the array and this deceleration is stronger than that for an isolated cylinder, especially for relatively large values of aD. As aD decreases, the percentage of the cylinders for which the drag force is less than 20 % of the force acting on an isolated one decreases to 20 % for aD ≈ 5 and to 0 % for aD ≈ 2.5. For aD < 5, 30 % to 40 % of the cylinders are subject to drag forces that are approximately 30 % to 50 % of the drag force acting on an isolated cylinder, with a very small number of cylinders (less than 5 %) being subject to drag forces that are comparable or larger than that acting on an isolated cylinder. It is only for sufficiently small values of aD, when no rollers form in the wake, that the number of cylinders subject to forces that are comparable to, or larger than, the one acting on an isolated cylinder becomes important. This percentage is close to 20 % for aD = 0.05. Though drag forces on the individual cylinders in the SVF = 0.1(2d) case are much larger than those in the SVF = 0.05 case, the distributions of the normalized drag forces are close to identical in the two cases. This shows that the distribution of the normalized drag forces inside the array is mainly a function of aD.

Figure 17. Distribution of the normalized mean streamwise drag force acting on the solid cylinders forming the array. The force on each cylinder is normalized by the drag force acting on an isolated solid cylinder. The vertical axis shows the percentage of the cylinders forming the array for which the normalized drag force is in between (2i)/10 and (2i + 2)/10 with i = 0 to 5.

The time-averaged, total streamwise drag force acting on the array is also mainly a function of aD. The variation of this force with aD is plotted in non-dimensional form in figure 18(a). A total drag coefficient for the array is defined, ${\bar{C}_{DG}} = {\bar{F}_x}/(0.5\rho {U^2}Dh)$, where ${\bar{F}_x}$ is the total streamwise drag force acting on the N cylinders forming the array. For the range of aD values considered in the present study, ${\bar{C}_{DG}}$ increases monotonically with aD, though the rate of increase gradually reduces for the larger values of aD. In addition to the total streamwise drag force acting on the array, a common measure of the drag acting on the array is the average non-dimensional streamwise drag force acting on the solid cylinders, ${\bar{C}_d} = ({\bar{C}_{DG}}/N)(D/d) = (1/N)\sum_{n = 1}^N {\bar{C}_D^n} $, where $\bar{C}_D^n = \bar{F}_x^n/(0.5\rho {U^2}dh)$ is the non-dimensional streamwise drag force acting on the rank n cylinder in the array (n = 1 to N). Though the instantaneous values of the streamwise drag force acting on the cylinders directly exposed to the incoming flow and to those in their vicinity can be a couple of times larger than the mean value, the peak values of the instantaneous total streamwise drag force, or equivalently of the total drag coefficient, ${C_{DG}}$, are of the order of 10 % of the corresponding time-averaged value (e.g. ${\bar{C}_{DG}}$). The largest fluctuations of the total non-dimensional streamwise drag force acting on the array are predicted in the SVF = 0.025 case and in the SVF = 0.1(2d) case. For these two cases, the maximum values of the instantaneous total spanwise force acting on the array are approximately 10 % of the time-averaged value of the total streamwise drag force.

Figure 18. Non-dimensional time-averaged streamwise drag force. (a) Time-averaged total drag coefficient for the array, ${\bar{C}_{DG}}$ versus aD; (b) mean, time-averaged solid cylinder drag coefficient, ${\bar{C}_d} = ({\bar{C}_{DG}}/N)(D/d)$ versus aD; (c) ${\bar{C}_d}$ versus aD in log-linear scale; (d) combined drag parameter for the array, ${\Gamma _D} = {\bar{C}_d}(aD)/(1 - \textrm{SVF})$ versus aD; (e) $\varGamma_{D}$ versus 1 − U ₁/U ₂. In addition to the present test cases (D/h = 20) with d = 0.0125D (open squares) and d = 0.025D (solid squares), some of the frames show data from the 3-D LES predictions of Chang & Constantinescu (2015) conducted for an array of long porous cylinders (solid blue circles) and the 3-D DES of Chang et al. (2017) conducted for arrays with D/h = 2.5 and 3.3 (solid green diamonds). The dashed line in panel (c) shows a best fit to the numerical data, ${\bar{C}_d} = 0.75\hbox{--}0.248\ast \textrm{ln}(aD)$.

As shown in figure 18(b,c), ${\bar{C}_d}$ decreases with aD and the decrease is close to logarithmic. Figure 18(b) also includes the predictions of ${\bar{C}_d}$ from the numerical simulations of Chang et al. (2017) and Chang & Constantinescu (2015) performed for arrays with aD < 5 and for arrays of long cylinders, respectively. Overall, the values for the present simulation performed for D/h = 20 are slightly lower than those in the other two numerical studies. So, increasing flow shallowness slightly reduces the mean non-dimensional drag force acting on the cylinders forming the array. This effect is even stronger for the combined drag parameter for the array, ${\varGamma _D} = {\bar{C}_d}aD/(1 - \textrm{SVF})$, shown in figure 18(d). This variable characterizes the relative influence of the drag term relative to the advection term (Chang & Constantinescu 2015). It also represents the ratio of the array diameter, D, to the array's drag length scale, ${L_c} = 2(1 - \textrm{SVF})/({\bar{C}_d} \cdot a)$. Increasing D/h from 3 to 20 results in a decrease of up to 25 % of $\varGamma_{D}$ for aD ≈ 10. The decrease is less than 7 % for aD < 3. So, for this variable flow, shallowness effects are also a function of aD. The combined drag parameter characterizes the overall flow deceleration inside the array, which is of the order of the difference between the incoming velocity, U ₀ and the bleeding velocity at the back of the array, U ₁. Results in figure 18(e) confirm that $\varGamma_{D}$ increases close to monotonically with 1 − U ₁/U ₀ and that a unique curve can be used to describe this variation for both moderately shallow and very shallow flow conditions.

7. Summary and conclusions

The present study showed that several important qualitative and quantitative changes occur in the flow and turbulence structure generated by a circular array of emerged cylinders in an open channel in between cases with moderately shallow flow conditions (1 < D/h < 5) and cases with very shallow flow conditions (D/h > 10). The paper focused on cases where the stability number was low enough such that roller vortices can still form in the wake for D/h = 20. This is the case of most significance for practical applications related to flow in shallow rivers containing vegetated patches away from their boundaries.

The main two qualitative changes associated with the increase in flow shallowness were the formation of a turbulent horseshoe vortex system containing multiple, strongly coherent vortices (up to 5) and the formation of arrays of parallel, co-rotating horizontal vortices in the vicinity of the roller vortices for sufficiently large values of aD. The dynamics of the horseshoe vortices was fairly similar to that of the laminar breakaway sub-regime observed for solid cylinders The present simulations showed that for sufficiently low aD (e.g. aD < 10), the horseshoe vortices lose their coherence for small polar angles, with the peak coherence occurring for polar angle magnitudes between 40° and 90°. This result is not as surprising given that main necklace vortices of low coherence near the symmetry plane were also observed for flow over a high-aspect-ratio rectangular cylinder. In the latter case, the coherence of the turbulent horseshoe vortices was the highest around the lateral edges of the cylinder rather than close to its symmetry plane (Kirkil & Constantinescu 2009). Only for relatively large values of aD, the main necklace vortex remained strongly coherent over all the frontal side of the array with D/h = 20. By contrast, for cases with D/h < 5, the turbulent horseshoe vortex system was found to contain only one main vortex whose coherence is the largest in the symmetry (ϕ = 0°) plane. The other main difference with arrays with D/h < 5 was that the cores of the main horseshoe vortices extended over most of the flow depth. So, their interactions with the eddies shed in the separated shear layers occurred over the whole flow depth rather than being confined to the near-bed region. The present study showed that this is a main condition for the eventual formation of a secondary, near-wake instability in the form of array of horizontal, parallel vortices in the wake of a surface-mounted porous obstruction.

For comparable values of aD, increasing the flow shallowness resulted in a decrease of the strength of the upwelling and downwelling motions inside the array. This was accompanied by a decrease of the combined drag coefficient for the array, $\varGamma_{D}$, and thus of the overall flow deceleration inside the array, with the effects being stronger for relatively large values of aD. Meanwhile, the variation of the mean, time-averaged streamwise drag coefficient of the solid cylinders, ${\bar{C}_d}$, with aD was only mildly affected by varying D/h. Increasing the flow shallowness also triggered an increase in the width (e.g. up to 7.5h or 0.375D in the SVF = 0.1 case) of the region where the erosive capacity of the horseshoe vortices is large. The near wake was also affected by flow shallowness. Most notably, an increase in the length of the separation bubble and of the degree of vertical non-uniformity of the mean flow were observed for relatively large values of aD. Meanwhile, the minimum value of aD (≈6.0) required for a bubble to form behind the array was found to be independent of D/h. The increase in the length of the separation bubble with D/h was accompanied by an increase of the non-dimensional length of the total wake, ${L^{\prime}_w}/D$. A strong increase of the Strouhal number associated with the formation and passage of roller vortices with D/h was observed for relatively low values of aD (e.g. St ≈ 0.4 for aD ≈ 5). This result is quite interesting given that the Strouhal number for solid cylinders was found to be largely unaffected by the increase in D/h (Zeng & Constantinescu 2017).