Drag on a Square-Cylinder Array Placed in the Mixing Layer of a Compound Channel

Abstract

There are no studies specifically aimed at characterizing and quantifying drag forces on finite cylinder arrays in the mixing layer of compound channel flows. Addressing this research gap, the current study is aimed at characterizing experimentally drag forces and drag coefficients on a square-cylinder array placed near the main-channel/floodplain interface, where a mixing layer develops. Testing conditions comprise two values of relative submergence of the floodplain and similar ranges of Froude and bulk Reynolds numbers. Time-averaged hydrodynamic drag forces are calculated from an integral analysis: the Reynolds-averaged integral momentum (RAIM) conservation equations are applied to a control volume to compute the drag force, with all other terms in the RAIM equations directly estimated from velocity or depth measurements. This investigation revealed that, for both tested conditions, the values of the array-averaged drag coefficient are smaller than those of cylinders in tandem or side by side. It is argued that momentum exchanges between the flow in the main channel and the flow in front of the array contributes to reduce the pressure difference on cylinders closer to the interface. The observed drag reduction does not scale with the normalized shear rate or the relative submersion. It is proposed that the value of the drag coefficient is inversely proportional to a Reynolds number based on the velocity difference between the main-channel and the array and on cylinder spacing.

1. Introduction

1.1. Motivation

Floodplains are formed when the river becomes able to store its excess sediment load outside the main channel [1]. They extend from the riverbanks to the edge of the valley, forming, in each cross-section, a two-stage geometry. Floodplains are periodically flooded, which occurs when the conveyance of the main channel is exceeded. Floodplains are generally composed of fertile soil where a variety of natural vegetation, from meadows to shrubs and forests, can develop [2]. The combination of fertile land and waterways encouraged early humans settle in floodplains. In spite of the flooding risk, many of todays’ towns and cities are located in floodplains. The present study addresses the problem of determining drag forces and drag coefficients on finite arrays of emergent cylinders located at the interface between floodplain and the main-channel, in flood (overbank flow) events. Through finite arrays, it is intended that they do not cover the full width of the floodplain and that its longitudinal extension does not allow for the establishment of equilibrium conditions or even gradually varied flow conditions. To the best of the authors’ knowledge, there are no studies specifically aimed at characterizing and quantifying drag forces and drag coefficients on finite cylinder arrays when these are placed in the mixing layer of compound channel flows. This study was carried out keeping in mind, primarily, the problem of determining hydrodynamic forces on buildings during urban floods—the array of square-cylinders is a simple model of an urban mesh on a river waterfront. The outcome of this study should be useful to modelers of river flooding in the built environment, as it provides conceptual tools and actual values of drag coefficients able to be used in preliminary shallow-water analysis.

1.2. State of the Art

The drag force on isolated cylinders and the characteristics of their wakes have been extensively investigated by [3,4,5,6], among others. This knowledge, however, is not immediately useful to determine hydrodynamic forces on arrays of cylinders, first of all because of the proximity of cylinder elements in the array. Addressing this issue, the hydrodynamic drag on two or more cylinders, arranged in different ways, has been studied by [7,8,9,10] or [11], among others. The latter investigated how the drag-coefficient of consecutive circular cylinders varies with the number of cylinders, their geometry and surface roughness. They measured directly and separately the force exerted on each cylinder. They reported 50% reductions of the average drag-coefficient of pairs of cylinders, in comparison to that of an isolated cylinder, when the between distance was less than three times their diameter. Yen et al. (2008) [12] observed similar behavior at low Reynolds numbers for two tandem square cylinders. Yen et al. (2011) [13] studied how drag varies in side-by-side cylinder-pair arrangements. They showed that the average drag-coefficient decreases as the spacing between the cylinders becomes smaller. However, this decrease is not greater than 20% (always in comparison with the coefficient of an isolated cylinder), demonstrating that a tandem arrangement appears to result in a larger drag reduction. In the case of square cylinders, Robertson [14] performed experimental tests in order to assess the average drag coefficient of pairs of cylinders arrays. The author found that the average drag coefficient of tandem pairs is reduced, featuring negative drag on the downstream cylinder for distances smaller than three times the cylinder width. However, the drag coefficient for cylinders side-by-side could, on the contrary, increase.

These results indicate that proximity contributes to reduce the pressure differences across each cylinder in tandem configurations. This effect was also investigated in larger arrays of cylinders. The limiting case is that of infinite arrays—those occupying the entire flow width and for sufficient length so that the flow loses the memory of any pre-array conditions. Under these conditions, Tanino et al. [15] studying flows with low cylinder Reynolds number (25 < Re_d < 685), for different solid-volume fractions, found that the array-averaged drag coefficient increases with solid fraction and is inversely proportional to the cylinder Reynolds number. Earlier, Nepf [16] suggested that for larger Re_d the array-averaged drag coefficient might decrease with the increase of the solid fraction. Ferreira et al. (2009) [17] performed measurements for Re_d larger than 1000 and, by analyzing the terms of the double-averaged momentum conservation equation, showed that it is likely that the average drag coefficient increases with the solid fraction for Re_d > 1000 if the gradient of form-induced stresses also increases.

If the array is infinite but not uniform, e.g., formed by patches that alternate periodically in the longitudinal direction its areal stem density, Ricardo et al. [18] experimental work uncovered variations of form-induced and Reynolds stresses associated with, but out-of-phase, the modulation of the solid fraction. These could explain the results of [19], where the array-averaged drag coefficient seems to be strongly influenced by the spatial memory of the flow, as defined by [20]: essentially, form-induced stresses preserve and carry downstream the features of the upstream flow (for which the solid fraction was different) thus influencing non-locally the mean drag coefficient. Consequently, one expects hysteresis, as flows with the same Re_d and solid fraction have different values of the patch-averaged drag coefficient depending on the upstream areal density of stems and solid fraction. Hence, in infinite arrays, the equivalent array-averaged drag coefficient may be smaller or larger than that of an isolated cylinder, depending on patch uniformity, areal cylinder density and Reynolds number.

There are no published systematic studies of the variation of the drag coefficient with the length of the array. It is thus not possible to unambiguously determine if an array can be considered infinite or finite. For relatively short arrays, considered finite by the authors, the results of [21], using a fall velocity proxy argument and regular stem configurations, seem to show that the array-averaged drag coefficient decreases with the increase of solid fraction. These results should be confirmed for emergent cylinders but they are in agreement with the findings of the wind tunnel experiments and computational fluid dynamics (CFD) simulations of [22]: upstream elements are subjected to higher drag and change the flow in such way that elements in the back rows are subjected to lower drag. Overall, the effect is that of decreasing array-averaged drag coefficient.

In the case of prismatic compound channels, the flow in the floodplain is slower (because the work of friction forces per unit flow weight is higher), generating a mixing layer in the main-channel/floodplain interface, with strong lateral gradients of the streamwise velocity [23,24,25,26]. The nature of the mixing layer depends on the relative submersion of the floodplain, on the lateral shear rate, and on floodplain confinement [27].

For a highly sheared flow, the mass and momentum mixing processes are essentially two-dimensional and the mixing layer is akin to an unstable planar shear flow [28]. In this case, it is useful to recall research on the influence of planar spanwise shear flows on cylinder wakes. This has been detailed by [5,29,30,31], among others. Drag and lift induced by these flows has been shown to depend on Re_d and on the shear parameter K, defined as the shear rate normalized by the cylinder diameter and the flow velocity in the centre-line of the cylinder [30]. The experimental work of [32] has shown that, for Re_d lower than 1600, the drag coefficient C_d decreases with Re_d for non-zero values of the shear parameter. For a given Re_d, C_d was also observed to decrease with the increase of the shear parameter. These results are strictly true for planar and linear shear flows and may not be applicable to cylinders in mixing layers of compound channel flows.

The shear flow in compound channel mixing layers is not linear, even when is essentially planar and dominated by two-dimensional processes. Furthermore, and unlike the case of linear shear flows, it is subjected to instabilities akin to the inviscid Kelvin-Helmholtz [33]. Under these conditions, Gymnopoulos et al. [34] investigated the mean drag-force behavior of an isolated circular cylinder placed at the main-channel/floodplain interface. Using the same experimental facilities and employing the same methodology of the current study, they reported a 10% reduction of the mean drag coefficient relatively to the case where the cylinder is found in symmetrical flow conditions. This reduction is similar to that found by [30] but under different values of the Reynolds number (five times smaller) and of the shear parameter (three times larger). This indicates that the fluid–cylinder interaction in a compound channel mixing layer should be more complex and not assimilable to that occurring in a linear planar shear flow. If the normalized velocity difference between the main-channel and the floodplain is smaller than 0.3, it is expected that the mixing processes are essentially three-dimensional [27], including the effect of secondary currents [35,36]. This underlines the point that results for planar linear shear flows may not be adequate for a quantitative discussion of flow-cylinder interaction in compound channel mixing layers.

1.3. Objectives

According to the above review, it is clear that the flow within cylinder arrays exhibits, in general, a superimposition of local effects and spatial memory effects. If the array is sufficiently long and uniform, the difference between local and distal collapses and the mean drag coefficient increases with solid fraction. On the other hand, if the array is small, longitudinally, it has been seen that drag is unevenly distributed—the front cylinders bear the heavier load while those in the back (downstream) rows experience a lower drag. It is also clear that single cylinders in simple planar shear flows are subjected to lower drag forces. Single cylinders in compound channel mixing layers are also subjected to lower drag forces, but the drag reduction does not show the same trends of the linear shear flow.

The above review also reveals that the characterization of hydrodynamic actions on finite arrays of cylinders subjected to shear flow has not received sufficient attention. In particular, the drag force on cylinder arrays in the mixing layer of a compound channel flow has not been quantified.

The current study is thus aimed at: (i) characterizing the drag force that a finite array of square-cylinders sustains, in overbank-flow conditions, when placed in the mixing layer of a compound channel; and (ii) at discussing the reduction of the corresponding array-averaged drag coefficient. Two overbank conditions are tested, featuring relative floodplain submergences of 0.41 and of 0.31. In both cases, the bulk normalized velocity difference between main-channel and floodplain is $\lambda =0.35$. It is expected that, for this value of $\lambda$, three-dimensional mixing processes should not be negligible, although the main mixing processes should remain two-dimensional.

To fulfill the objective, experiments were conducted in a laboratory prismatic compound-channel. Nine square cylinders were placed in the floodplain, next to the main-channel/floodplain interface. At the downstream location of the array the width of the mixing layer does not vary, which means that there are not net mass and momentum exchange between main-channel and floodplain.

The drag force on the array, at a certain elevation from the floodplain bed, was assessed experimentally by applying the integral equation of conservation of time-averaged momentum in a fluid control volume, as described in section “Theory”. All terms of the momentum-balance equation in the streamwise direction, except the fluid–solid interaction term, were computed from acoustic Doppler velocimetry (ADV) and water depth measurements. The description of the experimental procedure can be found in the third section. Results are shown in the fourth section: the array-averaged drag coefficient, defined as the bulk array drag force normalized by the product of fluid density, the frontal solid area and the square of the velocity that characterizes the flow upstream the array within its frontal area, are calculated for two values of the normalized velocity difference, λ. The corresponding values are discussed and compared with those found in the literature for isolated cylinders, infinite arrays and cylinder pairs.

2. Theory

2.1. Computing the Bulk Drag Force from an Integral Momentum Balance

A control volume analysis is employed to compute the bulk drag force—the drag force on the totality of the cylinders in the array. The fluid control volume, shown in Figure 1, is limited by outer, open sections (sections S₁ to S₆), and inner, solid, sections (S₀). There are momentum fluxes across the open sections. The forces applied on S₀ constitute the reaction of the solid elements on the fluid in the control volume.

For a statistically stationary turbulent flow, the Reynolds-averaged integral equations of conservation of linear momentum (derived in the Appendix A and herein referred to as RAIM equations) are applied to the fluid control volume of Figure 1. Employing the Einstein-summation convention, they can be written:

$$\begin{array}{l}\underset{S_c\backslash S_0}{{\displaystyle \int }}\rho U_i(U_k{n}_k)\mathrm{d}S\\ =\underset{V_c}{{\displaystyle \int }}\rho g_i\mathrm{d}V+\underset{S_c\backslash S_0}{{\displaystyle \int }}\left\{-P{n}_i\right\}\mathrm{d}S+\underset{S_c\backslash S_0}{{\displaystyle \int }}{T}_{ik}{n}_k\mathrm{d}S+\underset{S_c\backslash S_0}{{\displaystyle \int }}\left\{-\rho \overline{u_i^{\prime }{u}_k^{\prime }}{n}_k\right\}\mathrm{d}S\\ +\left(\underset{S_0}{{\displaystyle \int }}\left\{-P{n}_i\right\}\mathrm{d}S+\underset{S_0}{{\displaystyle \int }}{T}_{ik}{n}_k\mathrm{d}S\right)\end{array}$$

(1)

where the free index i = 1, 2, 3 stands for the Cartesian spatial directions, U_i is the ith component of the time-averaged velocity vector, ${u^{\prime }}_i$ is the ith component of the velocity fluctuations, g_i is the ith components of the acceleration due to gravity, P is the time-averaged pressure, ρ is the density of the fluid, T_ik is the time-averaged viscous stress tensor and $-\overline{{u^{\prime }}_i{u^{\prime }}_k}$ is the Reynolds stress tensor. The control volume is represented by V_c while S_c is the total control surface (the union of S₁ to S₆ and S₀), S₀ is the solid part of S_c, that coincides with the surface of the cylinders, and S_c\S₀ represents the open control surface (the union of S₁, …, S₆), through which there may be mass and momentum fluxes. The components of the outward pointing normal unit-vectors applied to each boundary are n_i. For the sake of conciseness, Equation (1) is written in Cartesian tensor notation, where the dummy index k is used in the summations inherent to the dot products in the first, fourth, fifth and seventh terms. The last two terms on the right-hand side of Equation (1) represent the force exerted on the flow by the cylinders:

$$R_i=-\left(\underset{S_0}{{\displaystyle \int }}\left\{-P{n}_i\right\}\mathrm{d}S+\underset{S_0}{{\displaystyle \int }}{T}_{ik}{n}_k\mathrm{d}S\right)$$

(2)

where ${\displaystyle \underset{S_0}{\int }\left\{-P{n}_i\right\}}\mathrm{d}S$ is the net pressure force and $\underset{S_0}{{\displaystyle \int }}{T}_{ik}{n}_{k}dS$ is the net viscous force.

If Equation (1) is applied in the streamwise direction, force R_x is the reciprocal of the bulk drag force applied on the cylinder array (F_D), i.e., R_x = −F_D and R_x = |R_x| = |F_D|.

Point-wise time-averaged variables were weight-averaged in the corresponding control sections (accounting for the area of influence of each point) to calculate the integrals in Equation (1), i.e., ${\displaystyle {\int }_{S_m}a}\mathrm{d}S\approx \left[a^{\left(m\right)}\right]S^{(m)}$, where S_m is an open control section, a^(m) is a generic time-averaged variable measured in that control section, [ ] is the section-average operator, [a^(m)] is the mean value of a in that section and S^(m) is the area of the control section. Viscous stresses at the open boundaries of the control volume were considered negligible in comparison to Reynolds stresses. Resolving the integral terms, Equation (1) may be written as:

$$\begin{array}{l}{R}_x=\rho g\sin \left(\theta \right)V_c\\ +\rho S^{(1)}\left(\left[U_x{}^{(1)}{U}_x{}^{(1)}\right]+\left[{\overline{u_x^{\prime }{u}_x^{\prime }}}^{(1)}\right]+\left[P^{(1)}\right]/\rho \right)-\rho S^{(3)}\left(\left[U_x{}^{(3)}{U}_x{}^{(3)}\right]+\left[{\overline{u_x^{\prime }{u}_x^{\prime }}}^{(3)}\right]+\left[P^{(3)}\right]/\rho \right)\\ -\rho S^{(2)}\left(\left[U_x{}^{(2)}{U}_y{}^{(2)}\right]+\left[{\overline{u_x^{\prime }{u}_y^{\prime }}}^{(2)}\right]\right)+\rho S^{(4)}\left(\left[U_x{}^{(4)}{U}_y{}^{(4)}\right]+\left[{\overline{u_x^{\prime }{u}_y^{\prime }}}^{(4)}\right]\right)\\ -\rho S^{(5)}\left(\left[U_x{}^{(5)}{U}_z{}^{(5)}\right]+\left[{\overline{u_x^{\prime }{u}_z^{\prime }}}^{(5)}\right]\right)+\rho S^{(6)}\left(\left[U_x{}^{(6)}{U}_z{}^{(6)}\right]+\left[{\overline{u_x^{\prime }{u}_z^{\prime }}}^{(6)}\right]\right)\end{array}$$

(3)

where θ is the angle between the channel bottom and a horizontal plane. If the thickness of the control volume h (Figure 1) is small compared to the other dimensions of the control volume, implying that the net mean momentum exchange in the vertical direction is negligible, then the last two terms on the right-hand side of Equation (3) may be omitted from the conservation equation:

$$\begin{array}{l}{R}_x=\rho g\sin \left(\theta \right)V_c\\ +\rho S^{(1)}\left(\left[U_x{}^{(1)}{U}_x{}^{(1)}\right]+\left[{\overline{u_x^{\prime }{u}_x^{\prime }}}^{(1)}\right]+\left[P^{(1)}\right]/\rho \right)-\rho S^{(3)}\left(\left[U_x{}^{(3)}{U}_x{}^{(3)}\right]+\left[{\overline{u_x^{\prime }{u}_x^{\prime }}}^{(3)}\right]+\left[P^{(3)}\right]/\rho \right)\\ -\rho S^{(2)}\left(\left[U_x{}^{(2)}{U}_y{}^{(2)}\right]+\left[{\overline{u_x^{\prime }{u}_y^{\prime }}}^{(2)}\right]\right)+\rho S^{(4)}\left(\left[U_x{}^{(4)}{U}_y{}^{(4)}\right]+\left[{\overline{u_x^{\prime }{u}_y^{\prime }}}^{(4)}\right]\right)\end{array}$$

(4)

The component R_x is estimated after all other terms in Equation (4) are determined experimentally, through acquisition of the three-component instantaneous velocity of the fluid in all open control sections and the free-surface elevation along the outer rim of the control section. If the vertical distribution of pressure is approximately hydrostatic, the calculation of the mean pressure on the surfaces is trivial and the gradients of the free-surface elevation are sufficient for the calculation of the pressure forces. Selecting a very thin prismatic control volume with its main dimension oriented parallel with the channel bottom, placed at an adequate distance above it, so that the effect of the bottom boundary to vertical distribution of velocities is not relevant, allows also for the definition of a single Reynolds number of reference for the assessed drag force [34]. The details of the mathematical derivation of Equation (3) are shown in the Appendix A.

2.2. Drag Coefficient

The drag coefficient of an isolated square cylinder, C_d, is expressed as:

$$C_d=\frac{R_x}{\frac{1}{2}\rho U_0^{2}hd}=\frac{2R_x}{\rho U_0^{2}hd}$$

(5)

where $R_x=\left|{\mathit{R}}_x\right|=\left|{\mathit{F}}_{\mathit{D}}\right|$ is the absolute value of the bulk drag force, h is the thickness of the control volume, d is the width of the square cylinder and U₀ is the time-averaged longitudinal velocity that characterizes the undisturbed upstream flow, calculated as a spanwise average within the frontal area of the cylinder.

For the case of an array that consists of finite number of identical cylinders, the array-averaged C_d is defined as the normalised bulk drag force (R_x) averaged over the total number of cylinders (n_c):

$$C_d=\frac{2R_x}{\rho U_0^2{n}_{c}hd}$$

(6)

where $\frac{1}{2}\rho U_0^2$ is stagnation pressure, $U_0^{}$ is a reference velocity representing the undisturbed flow upstream the array and defined as the longitudinal velocity laterally-averaged over the span of the cylinder array, and where $hd$ is the exposed area of each cylinder.

3. Experimental Facility, Instrumentation and Experimental Procedure

3.1. Experimental Facility

The experimental tests were conducted in the 10 m long prismatic and symmetric compound channel at the National Laboratory for Civil Engineering in Portugal (LNEC). The cross-section of the main channel is trapezoidal, while that of the floodplain is rectangular. The height of the main channel is 0.10 m and the slope of the margins (V:H) is 1.0. The exact dimensions are shown in Figure 2.

The bottom of the main channel is coated with polished concrete (hydraulically smooth). The floodplain is covered by a layer of synthetic grass, as described in [26]. The longitudinal slope of the channel is 0.0011. The inflow rates in the main channel and in the floodplain are regulated independently in order to achieve uniform-flow conditions at a relatively small distances from the entrance of the flow [37].

To minimize vertical accelerations and large-scale fluctuating motions, the flow is directed to the compound channel through a wall of perforated bricks in a tranquilization chamber and barrier of 0.02 m diameter PVC tubes. Acetate sheets are fixed at the end of the tubes pile for further tranquilization of the free surface. The flow depth at the outlet is controlled with tailgates, allowing uniform-flow depth for a given total discharge to be established. The channel is fed through a closed-circuit system involving a lower reservoir that collects the outflow of the channel, a stable-level upper reservoir that discharges into the channel and a pump station that connects the two reservoirs. The circuit ensures steady water head and discharge. Flow depth is measured by point gauges at 1.00 m and 8.00 m downstream of the channel entrance. The inflow rates in the main channel and the floodplain are monitored by electromagnetic flowmeters. A detailed description of the facility can be found in [38].

3.2. Experimental Tests

The experimental array is composed of nine square cylinders with width d = 0.045 m. It was placed at the interface between the main-channel and the flood plain, forming a grid orthogonal to the channel axis, as depicted in Figure 3. The distance d₀ between two consecutive cylinders was kept equal to 0.10 m. Spacing and sizing were chosen so as to align with those of the experiments carried out in the context of FlowRes ANR, a project oriented to the investigation of compound-channel flows with transition in bed roughness, described in [39]. The tests were run under uniform-flow conditions in the compound channel. To establish uniform flow in the channel, the total discharge Q was initially distributed in the main channel (Q_mc) and in the floodplain (Q_fp), according to the weighted divided channel method proposed by [40]. The cylinders were emergent under the tested flow conditions.

Two tests were conducted, with different characteristics of the mixing layer. A summary of the flow parameters for both tests, named SA_03 and SA_04, can be consulted in

Table 1. Experimental uniform-flow parameters.

Test	hr (-)	hfp (m)	Q (ls−1)	Qmc (ls−1)	Qfp (ls−1)	Frmc (-)	Frfp (-)	Red (-)
SA_03	0.31	0.045	58.9	42.3	16.6	0.46	0.40	13,800
SA_04	0.41	0.070	95.4	63.2	32.2	0.53	0.40	18,100

. These include the Froude number (Fr) and the cylinder Reynolds number defined as Re_d = dU₀/v, where $v$ is the kinematic viscosity of the fluid. The Froude number characterizes separately the main channel and the floodplain, and is given by Fr_fp = (Q_fp/(1.40h_fp))/√gh_fp and Fr_mc = {Q_mc/((0.40 + 0.20)(h_mc − 0.10) + 0.05)}/√gh_mc. The relative flow depth h_r is defined as the ratio between the floodplain flow-depth h_fp and the main-channel flow depth h_mc. The normalized velocity difference is λ = (U_mc − U_fp)/(U_mc + U_fp) [27], where U_mc and U_fp are the time-averaged and depth-averaged longitudinal velocities of the flow in the main channel and in the floodplain, respectively. In both flows, λ = 0.35, which indicates that the mixing processes should be dominated by two-dimensional planar interactions, associated to large-scale structures akin to Kelvin-Helmholtz instabilities. However, given that the value of λ is not significantly higher than 0.3, it was expected that three-dimensional effects in the mixing layer would not be negligible.

3.3. Velocity Measurements

Time series of the three components of the instantaneous velocity were acquired simultaneously with a side-looking Nortek Vectrino ADV at a rate of 200 Hz. The sampling interval was 3 minutes and sampling locations are shown as grid points in Figure 3a (test SA_03) and Figure 3b (test SA_04). The grid is denser in the case of SA_04 to account for possible higher flow complexity due to more pronounced three-dimensional effects. The grid in both tests was applied at a distance z₀ = h_fp/3 from the floodplain bed and at z = h_fp/3 − 0.06h_fp.

Post-processing involved de-spiking with the phase-space filter proposed by [41]. Time-series that exhibited bad data in more than 60% of their initial number of samples were eliminated. Time-averaged results of the corresponding grid points were replaced by space-interpolated values. These points represent less than 7% of the total number of grid points and are confined to the region close to the cylinders’ walls.

Free-surface time-series were acquired with ultrasonic distance-measuring sensors at a rate of 10 Hz during 2 min. These measurements were taken in the points of the grid that correspond to the outer rim of the control volume. The free-surface recordings were transformed into hydrostatic pressure forces at sections S₁ to S₄.

4. Results

4.1. Time-Averaged Flow and Main Wake Processes

The streamlines and the spatial distribution of the absolute value of time-averaged streamwise velocity are shown in Figure 4 for the SA_04 test. The control volume encompasses the near-wake field [42] of the array, where intense momentum exchange and turbulence occur. The spatial distribution of the velocities is indicative of the main wake-flow processes in both tests since similar patterns were also observed in the SA_03 test.

The ratio (d₀−d)/d of the array is equal to 1.22, a value similar to that for which [13] described the relevant wake patterns as the “gap-flow mode”. This wake flow features an alternatively deflecting jet along the streamwise direction between the cylinders and a nested vortex-shedding process within the one formed by the outward edges of two side-by-side cylinders. The streamlines observed in Figure 4 are compatible with this pattern with one notable exception—the formation length of the cylinder adjacent to the interface is smaller than that of the cylinder placed towards the floodplain. This indicates that drag is unevenly distributed in the lateral direction. The observed high velocities of the flow in the near wake and the relative small formation length of all downstream cylinders may indicate limited downstream velocity-deficits and, in turn, decreased bulk drag compared to that on an array with (d₀−d)/d >> 1.5 [13].

The patterns of the streamlines denote the main changes in the flow direction as it passes through the array. Upstream of the array, flow deflection at the frontal area of the cylinders is obvious. The deflected streamlines are directed both towards the main channel and the inner floodplain. The latter direction though, seems to present higher angles with the upstream-flow direction and the effect is traceable within a greater region in the floodplain (Figure 4). This may be attributed to the significantly lower flow velocities in the floodplain than those in the main channel, allowing the occurrence of greater yaw angles [43]. Downstream of the array the flow is attracted by the low-pressure region formed at the shaded areas (the formation length). Again, streamline angles with the main-flow directions in the floodplain are greater than the respective ones in the main channel, illustrating the effects of the time-averaged shear flow.

4.2. Analysis of Terms in Equation (4)

The terms of Equation (3) were calculated for a control volume with b = 2.7d₀, l = 6.25d₀ (note that the plan view of the control volume also shown in Figure 4) and h = 0.06h_fp. The momentum fluxes and the Reynolds stresses computed from the velocity measurements at z₀ = h_fp/3 (section S₅) and z = h_fp/3−0.06h_fp (section S₆) confirmed that net momentum flux across sections S₅ and S₆ is negligible and that the vertical variation of the turbulent stresses is small, justifying the application of Equation (4). All values are normalized with powers of the reference velocity U₀ estimated for each test. The values of the terms of Equation (4) are shown in Figure 5, Figure 6 and Figure 7.

The mechanisms underlying the intense increasing of the curvature of the streamlines involve momentum transport. The distributions of the convective transport U_xU_k at the vertical control sections S₁, S₂, S₃, S₄ are presented in Figure 6a,b. The product U_xU_k is also called as mean flux, as it represents the momentum flux of the time-averaged flow per unit mass and area, transported in the streamwise direction. The signs of the terms are chosen so that they represent simple summations in Equation (4).

As shown in Figure 5a, the maximum values of mean fluxes, in both control sections S₁ and S₃, are observed near the main-channel/floodplain interface, as a consequence of the higher velocities in the main channel. For the chosen locations of S₁ and S₃, both tests present rather balanced fluxes and a near zero net momentum flux. This was sought as a compromise, since placing the S₁ and S₃ very near the cylinder array would increase the uncertainty in all measured variables. This is especially true for the free-surface at S₃, because of the stronger vertical fluctuations generated by the interaction of the cylinder wakes and the lateral fluxes between the floodplain and the main channel.

Momentum transport along the lateral surfaces of the control volume S₂ and S₄ for both tests is shown in Figure 5b. Flow deflection upstream the array (out of the control volume-negative values in Figure 5b), observed in Figure 4, is also discernible here. The maximum values of the convective transport through these sections occurs just upstream the cylinder array. Flow reattachment is observed downstream the array. In both tests it seems that mean momentum transport out of the control volume upstream the array is higher close to the main channel (S₄) than in the inner floodplain (S₂). The different characteristics of the mixing layers of tests SA_03 and SA_04 is visible in the different values of the convective transport downstream the first row of cylinders. The outward convective fluxes upstream the cylinder are not influenced by the nature of the mixing layer. However, beyond the first cylinder row, the convective flux near the interface $U_x{U}_y$, oscillates between positive (inward in this context) and negative (outward) values. This is especially visible in the higher submersion test SA_04. The negative fluxes are much less expressive in the second and third rows, relative to those upstream of the first row. Outward convective fluxes are not seen at boundary $S_2$, next to the floodplain. This may indicate a perturbation in the array induced by the flow in the main-channel.

Note that the array spacing ratio is d₀/d = 2.22, compatible with the range 1.5 < d₀/d < 4 proposed by [12] for partially reattaching flow between two consecutive cylinders. Yet, the flow pattern and the convective fluxes at the boundary $S_4$ seem to indicate that the boundary layer separated from the upstream cylinder does not reattach to the side of the downstream cylinder but is driven inwards in the space between these two consecutive cylinders. It can be hypothesized that the relevant variables are the momentum difference between the main channel and flow in front of the array of cylinders and the space between the cylinders. In that case, a Reynolds number can be formed, ${\mathrm{Re}}_I=\frac{\left(U_{mc}-U_0\right)\left(d_0-d\right)}{\nu }$. Test SA_03 (lower submergence) is characterized by ${\mathrm{Re}}_I=13\times {10}^4$ while SA_04 (higher submergence) has ${\mathrm{Re}}_I=16\times {10}^4$.

The perturbation imposed by the flow at the interface on the flow in the space between cylinders should be mainly a reduction of the pressure difference across the first cylinder, similar to what is found in arrays affected by cross-flow [20]. On the second and third cylinders, the effects of the perturbation are more difficult to envisage: for instance, while lateral intrusion of momentum may contribute to increase the pressure downstream the second cylinder, it may also increase the pressure upstream that cylinder, which, in the reattachment regime, could be lower than the downstream pressure. In the absence of data to clarify this issue, it is further hypothesized that the overall effect of the perturbation induced by the flow at the interface on the flow in the space between the cylinders is to lower bulk drag force. It follows that the higher the value of ${\mathrm{Re}}_I$, the larger the perturbation and the lower the normalized drag force. Test SA_04, with the higher submergence and the higher momentum at the interface, should be the most affected.

On the side of the floodplain, section $S_2$, the regularity of the convective fluxes is compatible with cylinders in the reattachment regime characterized by none-to-weak vortex formation and negative drag acting on the downstream cylinder. The occurrence of negative drag is highly probable in the last row, which was also observed by [14] for d₀/d < 3 and Re_d~5600–12,800.

The values of the turbulent transport $\overline{{u^{\prime }}_x{u^{\prime }}_k}$, normalized with the scale velocity U₀, are shown in Figure 6a,b. The variable $\overline{{u^{\prime }}_x{u^{\prime }}_k}$ is also referred as stress, since it expresses the Reynolds stress per unit fluid density. The values of the stresses, both normal and shear, are smaller than those of the convective fluxes. This observation is highlighted by the use of the same scale of the vertical axis in Figure 5a and Figure 6a. The $\overline{{u^{\prime }}_x{u^{\prime }}_x}$ stress values in S₁ and S₃ (Figure 6a) resemble the distribution of the streamwise convective momentum transport (Figure 5a) featuring higher values at the main-channel/floodplain. Stronger stresses near the interface (y/d₀ = 0) are typical of compound-channel flows [23], which suggests that these peaks are not due to the presence of the array. The corresponding peak for S₃ is larger in both tests, compared to that of S₁. This observation implies that there is an excess in the $\overline{{u^{\prime }}_x{u^{\prime }}_x}$ stress, generated by interaction of the wake flow in the near-interface cylinder and the shear flow. This effect seems more pronounced in the case of test SA_03–lower submersion.

The shear stresses $\overline{{u^{\prime }}_x{u^{\prime }}_k}$ on the lateral control sections S₂ and S₄ are seen in Figure 6b. It is noteworthy that the peak values are obtained at the second row of the array, following a turbulence suppression caused by the accelerated flow upstream the first row. It is also relevant to note that in the case of test SA_03 (lower submergence) turbulent transport is much higher on the side of the interface (S₄). In the case SA_04 (higher submergence), turbulent transport is higher in section S₄, relatively to section S₂, just upstream the second row and for $x/d_0=1$ onwards. Downstream the latter section the characteristics of the compound flow mixing layer become again prevalent, relative to those of the array wake flow. In the space between the first and the second rows, the separated boundary layer from the first cylinder may seem to cross the $S_4$ section, in the case of test SA_03. However, this is not the case of SA_04, which is line with the hypothesis that, for this test, boundary layer reattachment occurs inside the space between cylinders (see discussion of convective fluxes at $S_4$). This justifies the stronger difference between shear stresses at $S_2$ and at $S_4$, in case of SA_03.

The distribution of the terms $g\left(z_{fp}-z_0\right)/U_0^2$, proxies for the pressure forces in Equation (4), assumed hydrostatic, are shown in Figure 7. Pressures upstream the array are larger than pressures downstream. It is recalled that the convective momentum flux measurements shown in Figure 5a revealed a very small imbalance between incoming (at $S_1$) and outgoing ($S_3$) convective momentum fluxes.

Lateral net later flux of momentum mostly compensates for the longitudinal momentum lost due to deflection in front of the array. However, the irreversible loss of mechanical energy in the fluid–cylinder interactions results in the pressure drop visible in Figure 7.

4.3. Total Net Contribution of Terms in Equation (4)

The net contributions of the terms of Equation (4) for both tests referring to the fluid control-volume of Figure 4 are quantified and presented in Figure 8 as percentages of the absolute value of the drag force.

The term $\left\{\rho S^{(1)}\left[U_x{}^{(1)}{U}_x{}^{(1)}\right]-\rho S^{(3)}\left[U_x{}^{(3)}{U}_x{}^{(3)}\right]-\rho S^{(2)}\left[U_x{}^{(2)}{U}_y{}^{(2)}\right]+\rho S^{(4)}\left[U_x{}^{(4)}{U}_y{}^{(4)}\right]\right\}/R_x$ is referred to as net convective transport, the term $\left\{\rho S^{(1)}\left[{\overline{u_x^{\prime }{u}_x^{\prime }}}^{(1)}\right]-\rho S^{(3)}\left[{\overline{u_x^{\prime }{u}_x^{\prime }}}^{(3)}\right]-\rho S^{(2)}\left[{\overline{u_x^{\prime }{u}_y^{\prime }}}^{(2)}\right]+\rho S^{(4)}\left[{\overline{u_x^{\prime }{u}_y^{\prime }}}^{(4)}\right]\right\}/R_x$ is referred to as the net turbulent transport, the term $\left\{S^{(1)}\left[P^{(1)}\right]-S^{(3)}\left[P^{(3)}\right]\right\}/R_x$ as the pressure imbalance and the term $\rho g\sin (\theta )V_c/R_x$ as the gravity contribution. Figure 8 shows that the pressure imbalance between S₁ and S₃ is the single major contribution for the drag force. This raises concerns about the accuracy of the estimation of the drag coefficient, given the assumption of the hydrostatic hypothesis.

While the pressure drop is the dominant term, it is also clear that the effects of the mixing layer, expressed through the contributions of the net convective transport and the net turbulent transport, cannot be neglected. The different contributions of the net convective transport of SA_03 and SA_04 justify quantitative and qualitative differences in the value of the drag coefficient.

The negative value of the net convective transport can be seen as an effect of the higher momentum imbalance between the flow in the main-channel adjacent to the interface and the flow in front of the array. The physical processes involved were presented in Section 4.1. It is hypothesized that the penetration of flow from the main channel into the space between cylinders contributes to reduce the pressure imbalance, particularly in the first cylinder row. The main attribute of this penetration is to reduce the value of the net convective momentum transport.

Finally, it should be underlined that the net turbulent momentum transport is not negligible, even if the fluxes associated to Reynolds stresses are one order of magnitude smaller than convective fluxes. Failure to include the turbulence contributions would result in errors of 10% to 15% in the value of the drag coefficient.

4.4. Sensitivity Analysis of the Values of the Drag Coefficient Relatively to the Position of Section S₃

Mean flow and turbulence at the wake of the array are seen to strongly vary longitudinally. To demonstrate the validity of the present integral method, it is necessary to assess the sensitivity of the values of the drag force and of drag coefficient to the changes in the wake flow. The drag coefficient was computed for different control volumes, differing only in the longitudinal positions of section $S_3$ and, hence, of length l (Figure 1). The bulk drag force on the array of cylinders was estimated according to Equation (4) and the drag coefficient was computed from Equation (6).

The different net contributions for the different control volumes are shown in Figure 9, for the case of test SA_04, for which the density of measurements was higher. Each contribution corresponds to its array-averaged value divided by the reference stagnation pressure $\frac{1}{2}\rho U_0^2$ and by the exposed cylinder area $hd$. As an example, the pressure contribution corresponds to $\left(\left\{S^{(1)}\left[P^{(1)}\right]-S^{(3)}\left[P^{(3)}\right]\right\}/n_c\right)/\left(\frac{1}{2}\rho U_0^{2}hd\right)$. The simple sum of all the contributions is the value of the drag coefficient. It is verified that the values of the drag coefficient are not significantly sensitive to the location of the downstream boundary. As examples, the values of the drag coefficient for $x/d_0=1.0$, $x/d_0=3.0$ and $x/d_0=5.2$ are 0.583, 0.578 and 0.580, respectively.

If the $S_3$ boundary is placed very near the cylinder array, the pressure imbalance is larger, but the negative contribution of the net convective transport is also larger. These negative values express the negative contribution of the outward momentum fluxes in the upstream reaches of the control volume, which, for short lengths, are not balanced by the positive inward contributions in the wake. As the length of the control volume increases, the positive momentum (inward) contributions increase, and the net momentum transport becomes less negative. Evidently, the further the $S_3$ section is placed, the smaller the pressure imbalance. The net turbulent transport increases with the length of the control volume as it picks up the lateral shear stress imbalance between the interface and the floodplain, downstream of the array.

4.5. Values of the Drag Coefficient and Discussion

The array-averaged values of C_d for tests SA_03 and SA_04 are presented in

Table 2. Drag coefficients estimated for various configurations of square cylinders.

Reference	Flow Type	Configuration		Red	Cd
SA_03	compound-channel flow	9-cylinder regular array (d0/d = 2.22)		13,800	0.74
SA_04	compound-channel flow	9-cylinder regular array (d0/d = 2.22)		18,100	0.58
Robertson (2016) [14]	open-channel flow	infinite regular array (x: d0/d = 5.26, y: d0/d = 2.63)		9670	2.04
Robertson (2016) [14]	open-channel flow	isolated cylinder in uniform flow		10,000–22,000	2.11
Lyn et al. (1995) [44]	closed water-channel flow	isolated cylinder in uniform flow		21,400	2.10
Yen and Liu (2011) [13]	air flow	isolated cylinder in uniform flow		21,000	2.06
Robertson (2016) [14]	open-channel flow	pairs side-by-side (2 < d0/d < 3)		5600–12,800	2.59–3.28
Yen and Liu (2011) [13]	air flow	pairs side-by-side (d0/d = 2.5)		21,000	1.90
Robertson (2016) [14]	open-channel flow	pairs tandem (2 < d0/d <3)		5600–12,800	0.92–1.15
Yen et al. (2008) [12]	air flow	pairs tandem (d0/d = 3)		900–1200	0.50
Yen et al. (2008) [12]	air flow	pairs tandem (d0/d = 1.5)		900–1200	1.15

. The same table summarizes the values of drag coefficients found in other studies, corresponding to different configurations of square cylinders and types of flow. Among these, some employ cylinders in open-channel flows [23], closed water-channel flows [44] and air flows [13]. One case concerns a regular array in open-channel flow [23] and the others regard cylinder pairs arranged in tandem or side-by-side configurations with relative distance between the cylinders d₀/d similar to that of the regular array examined herein. All C_d values in Table 2 are array-average values. As a general trend, the drag coefficient of square cylinders in tandem is smaller than that of isolated cylinders in uniform flow. The drag coefficient of cylinders in a side-by-side configuration may be increased relatively to that of isolated cylinders. The drag coefficients of the tests conducted here (SA_03 and SA_04) are smaller than those of isolated cylinders and also smaller than those of cylinders in tandem for the same cylinder spacing and range of Reynolds number. The larger cylinder spacing reported by [12] features lower values of the drag coefficient. This is attributed to the lower value of the Reynolds number and not discussed in this paper. The discussion below addresses the reason why the drag coefficient of the cylinders in the mixing layer is smaller than that of cylinders in tandem in uniform flow and clarifies the differences between the values obtained for the two tests SA_03 and SA_04.

Test SA_03, with the lower relative submergence, yielded the higher drag coefficient. The different values of the drag coefficient are mainly due to the differences in the net convective transport. Without this term, the summation of the normalized contributions for the drag coefficient are 0.72, for test SA_03, and 0.70 for test SA_04. The different values of the net convective transport express different hydrodynamics of the interaction between the compound channel mixing layers and the flow within the array.

It hypothesized that the key difference between test SA_03 and SA_04 is a result of the difference between the momentum in the main channel and the momentum of the flow in front of the array. If this difference is normalized by the total momentum transported compound channel, one obtains a parameter that is approximately the same as the normalized shear rate. This parameter is fundamental to understand the nature of the mixing processes between floodplain and main channel in the absence of localized perturbations. In the presence of the finite array, it is hypothesized that it is not appropriate to scale the momentum difference by any other momentum quantity, since the array locally disrupts the lateral exchange processes of the compound channel flow. Instead, it is proposed that the momentum difference between the flow in the main channel and the flow in front of the array should be taken as a measure of the ability of the flow in the main channel to penetrate in the space between consecutive cylinders. The velocity difference, a geometrical variable expressing the distance between consecutive cylinders and the fluid properties density and viscosity are grouped into a Reynolds number that may be more adequate to parameterize the differences in the net convective transport.

The differences between test SA_03 and test SA_04 are thus not a simple consequence of different relative submersions of the floodplain. The latter are a result of different discharges in the main-channel and in the floodplain. The higher momentum transported in the main-channel in test SA_04 (with the higher relative submersion) induces a change in the reattachment regime of the cylinders adjacent to the floodplain relatively to what would be expected in a simple array in uniform flow. This change, envisaged in the values of the momentum fluxes at section $S_4$ and in the pattern of streamlines shows, consists in a deflection of the separated boundary layers from the upstream cylinders. As per the effect of lateral momentum influxes, they do not re-attach at the outer walls of the downstream cylinders, as is expected in that range of Reynolds numbers and cylinder distances, but against the upstream faces of the downstream cylinders.

It is noteworthy that both tests produced significantly smaller values of C_d, comparatively to those of isolated, tandem or side-by-side cylinders in uniform flows (see Table 2). The array-averaged values of C_d of both SA_03 and SA_04 are closer to those of tandem cylinders than to those of the side-by-side reported by [14] or [13] for the same ranges of the Reynolds number and of the ratio d₀/d. It is thus hypothesized that downstream cylinders, especially those on the inner column or on the side of the floodplain, experience significantly less drag—or even negative drag—than that of the upstream cylinder. In short, it was expected that the array-averaged drag coefficient of the array under investigation would be smaller than that of an isolated cylinder, since each column is likely to behave as a system of cylinders in tandem, even if the side-by-side configuration may attenuate this effect. The additional drag reduction brought about by the mixing layer called for an explanation. It was hypothesized that it should be the effect of convective fluxes that equalize pressures across the second row of cylinders, especially in the column adjacent to the interface, an effect that is stronger in test SA_04.

It was earlier argued that the relative importance of the pressure forces to the total drag force poses the question of the accuracy of the results, given the hydrostatic hypothesis. Deviations from hydrostatic pressure distribution are expected only if streamlines exhibit strong vertical curvatures. It is unlikely that such strong curvatures are found in this flow. In the case of the upstream section, the values agree well with those obtained by application of the Bernoulli principle—conversion of kinetic head into pressure head. In the case of the downstream section, the results of Figure 9 show that the drag coefficient is essentially independent of the longitudinal position of that section. Given that vertical accelerations are more likely to occur just downstream of the array, it is concluded that streamline curvature in the vertical plane should be negligible everywhere.

As a last remark, concerning the precision of the results, if it is assumed that the precision of the free-surface measurements at the downstream section is 1 mm, the values of the drag coefficient would vary as follows: 0.742 ± 0.09, for test SA_03, and $0.577\pm 0.07$ for test SA_04.

5. Conclusions

This study was aimed at experimentally assessing the bulk drag force on arrays of cylinders in compound channel mixing layers and the corresponding array-averaged drag coefficients. The drag force was estimated by applying the Reynolds-averaged integral momentum conservation equations (RAIM) in a fluid control-volume around the array. The robustness of the methodology was demonstrated by carrying out a sensitivity analysis to the size of the control volume. It was shown that the values of the drag coefficients do not depend on the location of the downstream section.

Two values of the relative submergence were tested. It was found that the values of the array-averaged drag coefficients, of both tests, are smaller than those of cylinders in tandem or side-by-side, obtained for comparable values of the Reynolds number and cylinder spacing. Drag reduction in the downstream cylinders of a tandem arrangement has been verified in the reattachment regime. The drag reduction observed in the tested arrays should be partially explained by a strong drop of the pressure imbalance in the downstream rows of the array, leading, possibly, to negative drag.

Further drag reduction should be due to the effect of lateral convective fluxes of momentum from the main-channel to the space between consecutive cylinders. This perturbation is likely to strongly reduce the pressure imbalance on the first cylinder row, on the side of the main channel. It is argued that this effect should become more pronounced with the increase of the difference between the momentum in the main-channel flow and the momentum of the flow in front of the array. The value of the array-averaged drag coefficient should thus be inversely proportional to a power of the Reynolds number ${\mathrm{Re}}_I=\left(U_{mc}-U_0\right)\left(d_0-d\right)/\nu$. Test SA_04 featured the highest imbalance of momentum relatively to the kinematic scale $\nu /\left(d_0-d\right)$ and, hence, exhibited the stronger drag reduction. The normalized shear rate or the relative submergence, per se, are not adequate parameters to express drag reduction in cylinder arrays in the mixing layer of compound channel flows. Further attention should be given to finding a relation between the variation of the drag coefficient of finite arrays in mixing layers with ${\mathrm{Re}}_I$ for a wider range of normalized shear rates, Reynolds numbers and relative floodplain submersions.

This study is meant to provide researchers on urban flooding with a conceptual framework to derive drag forces from fundamental principles of hydrodynamics and to propose an actual quantification of drag coefficients susceptible to be used in simplified two-dimensional shallow-water models using, e.g., porosity approaches to determine momentum and energy losses due to urban meshes. The method used in this study provided values for the drag coefficient at a height of approximately 1/3 of the flow depth above the floodplain bed. In principle, the method can be applied to obtain depth-averaged estimates of the drag coefficient. For that purpose, it would be beneficial to obtain velocity data at several horizontal planes. The method is sensitive to the errors in the pressure forces. The assumption of hydrostatic vertical pressure distribution is an important limitation. The method can be improved if pressures are taken with pressure transducers. Another improvement would be the determination of all three velocity components in planar regions with stereo particle image velocimetry. This would allow for the calculation of instantaneous forces on the cylinders provided that the velocity measurements could be synchronized with pressure acquisitions. It would thus strengthen the argument that the additional drag reduction due to the mixing layer is mainly the effect of the lateral convective fluxes that equalize pressures across the second and following rows of cylinders.