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 Reference Spalart2009). 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 Reference Chang, Constantinescu and Park2007; Rodi, Constantinescu & Stoesser Reference Rodi, Constantinescu and Stoesser2013). The one-equation Spalart–Allmaras (SA) model (Spalart & Allmaras Reference Spalart and Allmaras1992) 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 Reference Zeng, Constantinescu and Weber2008). 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. (Reference Chang, Constantinescu and Park2007, Reference Chang, Constantinescu and Tsai2017).

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 (Reference Cheng and Constantinescu2020) and Constantinescu (Reference Constantinescu2014) 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 (Reference Kirkil and Constantinescu2012, Reference Kirkil and Constantinescu2015) 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. (Reference Chang, Constantinescu and Tsai2017) 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 (Reference Zong and Nepf2012) and Chen et al. (Reference Chen, Ortiz, Zong and Nepf2012).

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. (Reference Chang, Constantinescu and Tsai2017, Reference Chang, Constantinescu and Tsai2020) 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. (Reference Chang, Constantinescu and Tsai2017), 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. Reference Chen, Ortiz, Zong and Nepf2012).

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 (Reference Cheng and Constantinescu2020) and Zeng & Constantinescu (Reference Zeng and Constantinescu2017) 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. Reference Rodi, Constantinescu and Stoesser2013). 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. (Reference Chang, Constantinescu and Tsai2017) 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 Reference Kirkil and Constantinescu2015).

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 Reference Kirkil and Constantinescu2015). 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. (Reference Chang, Constantinescu and Tsai2017) 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. (Reference Chang, Constantinescu and Tsai2017).

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. Reference Chang, Constantinescu and Tsai2017). 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 Reference Kirkil and Constantinescu2012). 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 Reference Zeng and Constantinescu2017, 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 Reference Kirkil and Constantinescu2015). 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. (Reference Chang, Constantinescu and Tsai2017) 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 (Reference Zong and Nepf2011) and Chen et al. (Reference Chen, Ortiz, Zong and Nepf2012). 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. (Reference Chang, Constantinescu and Tsai2017) 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 (Reference Chen and Jirka1995) 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 (Reference Kirkil and Constantinescu2015). 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 (Reference Zeng and Constantinescu2017) 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 Reference Zong and Nepf2011). 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. Reference Chang, Constantinescu and Tsai2017), 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 Reference Chang and Constantinescu2015; Chang et al. Reference Chang, Constantinescu and Tsai2017), 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. (Reference Chang, Constantinescu and Tsai2017) 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. (Reference Chang, Constantinescu and Tsai2017), 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. (Reference Chang, Constantinescu and Tsai2017) 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. Reference Chang, Constantinescu and Tsai2017). 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 Reference Zong and Nepf2012). 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 (Reference Zong and Nepf2012) 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 (Reference Chang and Constantinescu2015), 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. (Reference Chang, Constantinescu and Tsai2017) 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. (Reference Chang, Constantinescu and Tsai2017).

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. Reference Yang, Chung and Nepf2016). 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. Reference Sumer, Chua, Cheng and Fredsøe2003). 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 Reference Cheng, Koken and Constantinescu2018). Following Chang et al. (Reference Chang, Constantinescu and Tsai2017), 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. (Reference Chang, Constantinescu and Tsai2017). 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. (Reference Chang, Constantinescu and Tsai2017) 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 (Reference van Rijn1984) formula which is often used to estimate the volumetric flux in eddy-resolving simulations with a movable bed (e.g. see Chou & Fringer Reference Chou and Fringer2010). 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. (Reference Chang, Constantinescu and Tsai2017) 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 (Reference Chang and Constantinescu2015) conducted for an array of long porous cylinders (solid blue circles) and the 3-D DES of Chang et al. (Reference Chang, Constantinescu and Tsai2017) 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. (Reference Chang, Constantinescu and Tsai2017) and Chang & Constantinescu (Reference Chang and Constantinescu2015) 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 Reference Chang and Constantinescu2015). 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 Reference Kirkil and Constantinescu2009). 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 Reference Zeng and Constantinescu2017).