1 Introduction

The effects of vegetation on water flow, turbulence, and sediment transport have attracted considerable interest from researchers in a variety of fields [1,2,3,4]. For example, it has been demonstrated that the characteristics of canopy flows should be clearly determined to properly understand the behaviors of wind flows involving forest canopies [5, 6] and groups of buildings in urban areas [7], as well as the behaviors of water flows involving coral reefs [8] and in gravel-bed rivers with boulders [9,10,11]. It is also important to know the resistance characteristics of a group of cylinders in fluid flow because of the occurrence of this structure in a number of applications of fluid dynamics, such as heat exchangers, where arrays of cylinders are a common configuration [12].

Generally, in flow calculations, multiscale obstacles are divided based on two characteristic scales: large-scale topography, which can be resolved using a computational mesh explicitly as a boundary shape; and small-scale topography, which is evaluated using subscale models [9, 11]. In a number of studies, numerical simulations of flows through vegetation have been conducted using fine meshes to resolve the effects of obstacles [7, 13]; in others, the effects of vegetation have been modeled using temporal-spatial averaging of the governing equations [3, 5, 6, 8, 14,15,16,17,18]. Based on the complexity of vegetation conditions in rivers, the second aforementioned approach is considered to be suitable for practical applications [18,19,20]. When averaged over space in Reynolds-averaged Navier–Stokes equations for local flows, in addition to the dispersion term due to sub-average-scale flows, the fluid force terms appear in conjunction with Gauss’ divergence theorem in the momentum equation integrating the local pressure distribution and shear stress acting on the objects [6, 15]. The problem, as with the closure problem for turbulence, is how to represent the hydrodynamic force terms induced by sub-average-scale flows around the boundaries and distinguish them from the average-scale flow [21]. The drag-force terms in the momentum equations are evaluated using drag coefficients, which determine the accuracy of the models [22]. The drag coefficients also appear in equations for evaluating the turbulence intensity and dispersion coefficients of passive particles in water flows [3, 14, 17, 23, 24].

Meanwhile, several studies have contributed to determining the drag induced by vegetation [25,26,27,28,29]. However, although actual vegetation tends to be irregularly shaped and flexible [18], circular cylinders have been employed to model rigid vegetation in calculations regarding the characteristics of canopy flow. The drag coefficient CD of a circular cylinder in a free-stream flow (defined as a uniform flow that lacks a free surface, and typically distinguished from a uniform open-channel flow) is known to be a function of the Reynolds number \( \mathcal{R}\) = Vd/v (V, free-stream velocity; d, cylinder diameter; v, kinematic viscosity); the same is true for a sphere among several representative immersed bodies [30,31,32]. The variation in the drag coefficient of an array of cylinders depends on the space between the cylinders and the density of the array. Several researchers have contributed to correlating drag coefficients of cylinder arrays with well-established relationships between the drag-force coefficient and Reynolds number for a single cylinder. For example, Koch and Ladd [25] and Tanino and Nepf [26] applied the equation proposed by Ergun [33], which was devised for pressure drops in packed columns, to the drag coefficient and Reynolds number function. Kothyari et al. [34] proposed a relationship between the drag coefficient and Reynolds number for subcritical flows as a positive logarithmic function of the density, showing that the drag coefficient increases rapidly for low density values. Cheng [35] proposed a formula for the drag coefficient of a cylinder array by expanding the relationship for an isolated cylinder based on the concept of a pseudo-fluid model [36] and validated the formula against several drag-force coefficient datasets.

On the other hand, in contrast to most canopy flows, open-channel flows, such as in rivers, involve the Froude number \( \mathcal{F}\) = V/(gh)0.5 (g, gravity; h, water depth), an important dimensionless parameter that influences the drag coefficient, especially for emergent obstacles in shallow-water flows. In fact, the drag coefficients of hydraulic structures are related to the Froude number [31]. However, only a few studies have investigated the effect of the Froude number on the drag coefficients of emergent cylinders in open-channel flows. Furthermore, Kothyari et al. [34] indicated that, whereas the drag coefficient is not substantially influenced by the Froude number in subcritical flows, the coefficient decreases with respect to the Froude number in supercritical flows. The results of a study by Huai et al. [37], in which a large-eddy numerical simulation model was used in conjunction with three identical-density conditions, showed that CD largely decreases as \( \mathcal{F}\) increases. It is significant that CD is smaller in supercritical flows than in subcritical flows; however, the underlying mechanism for this distinction has not yet been elucidated.

Most previous studies on the evaluation of drag forces and flow resistance were conducted based on uniform open-channel flows. However, to fulfill the objective of applying drag-force coefficients to numerical simulations of flows in natural rivers [19, 20], a method for evaluating the effects of non-equilibrium flow in average scales on the resistance must be developed. For example, Busari and Li [38] investigated the interference effects of finite-length vegetation in gradually varied flows, whereas Uchida et al. [11] proposed a dynamic rough-wall law to account for the effect of non-equilibrium flow near a bed with submerged boulders on the flow resistance. However, despite these efforts, there has been insufficient research on methods for evaluating the direct effects of non-equilibrium open-channel flows on the drag force.

To take into account the resistance of vegetation in rivers, it is necessary to elucidate the effects of variations in the water surfaces around obstacles. For non-equilibrium flows with pressure gradients in an average scale, because the separation regions behind the obstacles exhibit low flow velocities [12, 39, 40], pressures acting on the backs of the obstacles are inferred to be affected by the pressure downstream in the scale of the separation length.

The objective of this study is to develop a method for evaluating the drag coefficient of emergent cylinders in an array in a non-uniform open-channel flow. The proposed method is to be designed for numerical simulations. In addition to the conventional drag-force expression, which includes the drag coefficient, this paper proposes two new terms for evaluating the effects of water surface variation in a local scale and pressure gradient in an average scale on the drag force for emergent cylinders. For the first term, a drag-coefficient equation is derived theoretically as a function of the Froude number. This equation is compared with drag coefficients from previous experiments and our laboratory experiment on uniform flows without the introduction of any calibration parameters. For the second term, the representative length of the separation zone is introduced to evaluate the effect of the pressure gradient in the average scale on the drag force in non-uniform open channel flows. The calculation of gradually varied flows with emergent cylinders reveals the necessity of the term to account for the effect of the pressure gradient on the drag forces in accelerating flows. The validity of the term is supported by the results of supplementary experiments that focus on the distribution of pressure on a cylinder in a group. Based on this pressure distribution, the negative pressure gradient for accelerating flow decreases the pressure in the separation zone behind the cylinder. It is then demonstrated that several non-uniform water-surface profiles, including accelerating and decelerating flows, can be reproduced well using a simple unified expression incorporating these two terms.

2 Experimental setup

Two sets of experiments were conducted in this study. The first set of experiments investigated longitudinal water surface profiles using an experimental channel, in which several small-scale emergent cylinders were installed, and focused on developing a method for evaluating the resistance. To clarify the effects of water surface variation and pressure gradients on drag force, three types of flow behavior, i.e., uniform, accelerating, and decelerating flows, were examined. The second set of experiments involved a supplemental study that elucidated the mechanism of the hydrodynamic force and pressure distribution on emergent cylinders. This experiment was conducted to measure the local water surface variations around large-scale emergent cylinders and the distribution of pressure acting on a cylinder in the group.

The first set of experiments was conducted in a water-recirculating flume (length, 3.0 m; width B, 0.3 m) with varying channel slopes. The flows were adjusted to either uniform or non-uniform flow using a weir installed downstream. For this study, we regard uniform flow as that when the water surface gradients coincide with the channel gradient. The emergent cylinders were circular wooden dowels (diameter d, 10 mm; height, 20 cm). The Styrofoam floor of the flume was drilled, and the dowels were installed individually in longitudinal and lateral intervals Δx and Δy, respectively, of 0.0344 m in a staggered array from the downstream end of the channel (Fig. 1) to maximize the flow resistance. The density of cylinders per unit area, λ, is defined as λ = Nπd2/4, where N is the number of cylinders per unit area (present experimental conditions: N = 1/(ΔxΔy)), and d is the cylinder diameter. The density λ = 0.067 in this experiment is in the range of actual vegetation densities in the flood plains of Japan [41].

Fig. 1
figure 1

Experimental setup of flume and cylinders: a side view (cylinders omitted); b, c plan views; d magnified plan view

Full size image

The experiments consisted of Case L and S series: In the Case L series (Fig. 1b), the cylinders cover the last 2.0 m of the flume, whereas in the Case S series (Fig. 1c), the cylinders cover only 1.0 m of the flume between the upstream and downstream ends. Flow depth was measured using a point gauge with a 0.1-mm vernier scale mounted on a carriage that could be moved along the flume rails. The water depth h was measured along the longitudinal axis on the centerline of the flume at intervals of 0.2 m. The resistance evaluation in this study is based on the results of this water surface profile measurement, i.e., the focus of this study is the modeling of phenomena that can affect the larger-scale analyses of one-dimensional measurement approaches, such as water surface profiles. On the other hand, the effects of the upstream and downstream ends of a cylinder group on the internal structure, such as on the velocity distribution and turbulence in two- or three-dimensional flow structures, are not included in this study. Investigations regarding these aspects of flow behavior will require further advancement of the results of this research. The volumetric discharge was measured using a triangular weir with a tank 0.8 m in length, 0.4 m in width, and 0.4 m in height. The average cross-sectional velocity (bulk velocity) U was defined using the mean depth h as U = Q/(hB), where Q and B are the discharge and width of the flume, respectively.

There are several possible definitions of representative velocity for evaluating canopy drag, including bulk velocity U [16,17,18, 42, 43], pore velocity Up = U/(1 − λ) [18, 26, 34, 41], and constructed velocity Uc = U/(1 − d/Δy) [44]. The undisturbed velocity, which is defined as a local velocity removing the target obstacle, is considered appropriate to evaluate a local pressure change by installing the target obstacle [45,46,47]. Because the local pressure change around the obstacle induces the hydrodynamic force, the representative velocities include the bulk, pore, and constructed velocities, which can be regarded as attempts to estimate the undisturbed velocity. However, the undisturbed velocity is affected by other cylinders and hydraulic conditions, and thus a proper method for the evaluation of the undisturbed velocity has not yet been established [47]. The definition of representative velocity is also relevant to discussions regarding the representative flow velocity for the Froude number, which is relevant to deriving the depth profile equation (see Eq. 8). In this study, the vegetation density was relatively low, and its effect was small [18]. Furthermore, as in many numerical calculations on flows involving vegetation [16, 17, 20], the bulk velocity was employed to focus on the effects of the water surface and pressure gradient on the drag coefficient. Nevertheless, there remain questions regarding the representative velocity for future studies to resolve.

Table 1 summarizes the bed slope and discharge conditions for the experiments on uniform, accelerating, and decelerating open-channel flows. For Cases L1–L6, the differences in the gradients between the water surface and bed, calculated using measurements obtained at 0.2-m intervals over 1.4 m, were less than 0.4 × 10−4. Large values of \( \mathcal{R}\) were used to decrease the dependence of the drag coefficient on the Reynolds number [26,27,28,29, 34, 35, 48]. To distinguish the effects of the Froude number from those of the Reynolds number, a more detailed investigation will be necessary because the effect of the Reynolds number on the drag coefficient varies with cylinder conditions [25,26,27,28,29]. Thus, instead of analyzing a large amount of drag coefficient data, we examine the effects of the Froude number and water surface theoretically in the following section.

Table 1 Experimental conditions for staggered array of vertical, rigid circular cylinders
Full size table

The second set of experiments, a supplementary study using larger-scale cylinders (diameter 0.1 m), was conducted in a channel of width 0.8 m and length 24 m to measure the local water surface around an emergent cylinder and the distribution of pressures acting on it. The channel slope was 1/200. The cylinders, 3.3 m in length, were installed in a staggered array with Δx = 0.25 m and Δy = 0.225 m, as shown in Fig. 2. The experiments included uniform-flow and acceleration conditions with a constant experimental discharge Q of 0.035 m3s−1. The experimental conditions are summarized in Table 2. In this set of experiments, the representative velocity for calculating the pressure coefficient and drag coefficient was measured using a MicroADV (SonTek) probe based on an undisturbed flow, i.e., without the cylinder of interest, to eliminate the effect of high velocities along the sidewalls, which were due to the small number of cylinders in the transverse direction (each row comprised only two or three cylinders, as shown in Fig. 2). The vertical distribution of pressures acting on the center cylinder of the group was also obtained. Figure 3 shows the pressure-measurement cylinder installed in the center of the group (x = 0 in Fig. 2). This 10-cm-diameter cylinder, along the side of which six equidistant holes were drilled (3 cm apart), was used to measure the pressure at different water depths. The cylinder was equipped with 5-mm-diameter tubes, which were connected to manometer glass tubes with a diameter of 0.6 mm. The pressure p on the cylinder was measured at an arbitrary angle θ with respect to the flow using the manometers as the pressure-measurement cylinder was rotated. The measurements were performed at intervals of π/18. The accuracy of the pressure measurement using the manometers was estimated to be 0.2–0.3 mm based on averages over the vertical pressure points for each angle.

Fig. 2
figure 2

Experimental setup with large-cylinder array for measuring pressure distribution around cylinders and local water surface profiles

Full size image
Table 2 Experimental conditions for pressure distribution on emergent cylinder
Full size table
Fig. 3
figure 3

Pressure-measurement cylinder with manometers

Full size image

3 Equation of water surface profile for flow with cylinders

To determine the relationship between the water surface profiles and drag forces on the vegetation, we first consider the momentum equation in the streamwise direction for shallow-water flows (hydrostatic pressure and no vertical velocity distribution) in a wide rectangular channel:

$$ \rho \left( {\frac{\partial U}{{\partial t}} + U\frac{\partial U}{{\partial x}}} \right) = \rho gS_{0} - \rho g\frac{\partial h}{{\partial x}} - \frac{{f_{t} }}{h} $$
(1)

where ft is the external force acting on the water column in a unit area, including the shear stress on the channel bed and drag forces on the vegetation. For a steady-flow condition, we have

$$ q = Uh = const. $$
(2)

where q is the unit flow discharge. By substituting Eq. (2) into Eq. (1), we obtain

$$ \left( {1 - \mathcal{F}^{2} } \right)\frac{dh}{{dx}} = S_{0} - \frac{{f_{t} }}{\rho gh} $$
(3)

where \( \mathcal{F}\)2 = q2/gh3. For the case involving no cylinders, ft = τ0 = ρgSf (τ0, bed shear stress; Sf, friction slope); therefore, the well-known governing equation for gradually varied flows [30, 49] can be obtained. For flows with emergent-cylinder vegetation models, the total force ft includes a term for the drag force acting on the vegetation:

$$ f_{t} = \tau_{0} + Nf_{d} $$
(4)

where Nfd is the drag force acting on the cylindrical obstacles in a control volume with a unit area, and fd is the drag force on a cylindrical obstacle. Substituting the Manning roughness coefficient n in SI units and the drag coefficient CD, we obtain

$$ f_{t} = \frac{{\rho gn^{2} U^{2} }}{{h^{1/3} }} + NC_{D} \frac{{\rho U^{2} }}{2}hd $$
(5)

For evaluations involving bed shear stress, we did not consider deformations of the vertical velocity profile, as in many two-dimensional analytical studies on rivers [16,17,18,19,20]. It is known that the vertical velocity profile varies between accelerating and decelerating flows [50] and with the resistance by the cylinder group [44]. In addition, because the presence of the cylinder group results in horseshoe vortices, a secondary resistance due to the cylinder group, apart from the drag forces acting on them, is also generated [11]. To accurately evaluate the total resistance, an advanced model capable of calculating the non-equilibrium velocity profile [11] will be necessary. However, as indicated afterward, because the contribution of bottom shear stress to the total resistance is assumed to be much smaller than the drag force on the cylinders, the first term of Eq. (5) was simplified based on an assumption of a uniform velocity profile. When Eq. (5) is substituted into Eq. (3) for non-uniform flow, the water depth profile for a gradually varying flow with emergent vegetation can be obtained [38, 49]:

$$ \frac{dh}{{dx}} = \frac{{S_{0} - \frac{{n^{2} U^{2} }}{{h^{4/3} }} - \frac{{C_{D} NdU^{2} }}{2g}}}{{1 - \mathcal{F}^{2} }} $$
(6)

The results for uniform flow (dh/dx = 0) can then be deduced:

$$ S_{0} - \left( {2\frac{{n^{2} }}{{h^{1/3} }} + C_{D} Ndh} \right)\frac{{U^{2} }}{2gh} = 0 $$
(7)

For a flow that involves dense emergent vegetation, the drag coefficient term is known to be dominant over the bed friction term [26]. The bed friction term with n = 0.010 estimated for the channel bed material (Styrofoam) used in this study is much smaller than the drag force terms (< 1% for Cases L1–L6 outlined in Table 1). In this study, the drag force acting on the emergent cylinders in the uniform open-channel flows (Cases L1–L6) was evaluated using Eq. (7), using averaged gradients of the total head U2/2 g + h + z, instead of S0, to minimize the effects of non-equilibrium flow.

4 Drag in uniform open-channel flows

In this study, the drag fd, defined in Eq. (4), on an emergent obstacle in open-channel flows comprises three components: the base component f0, which excludes the effects of water surface variation and pressure gradients or drag forces acting on an infinitely long cylinder for deep uniform open-channel flows (\( \mathcal{F}\) = 0); and two additional components fs and fp, associated with variations in the water surface and the pressure gradient attributed to the drag force, respectively.

$$ f_{d} = f_{0} + f_{s} + f_{p} $$
(8)

This section discusses fs; the other component fp is discussed further in the paper in the section on the calculation of the water surface profile.

The base component f0 has already been analyzed in several past studies. Typically, the drag force is composed of pressure and friction components, the former being dominant in flows with high Reynolds numbers [30, 32], in which the drag is attributed directly to the relative pressure Δp acting on an object in the flow. Herein, the relative pressure on a vertical cylinder is expressed as Δp = p − p0, where p is the actual pressure acting on the cylinder surface, and p0 is the pressure of the “free stream” or “undisturbed flow” [45,46,47] at that point without the objective cylinder. The pressure drag is often referred to as “form drag” because of its strong dependence on the shape or form of the object [30]. In general, pressure is described as a dimensionless pressure coefficient Cp based on the representative velocity U under free-stream or undisturbed-flow conditions as follows:

$$ C_{p} = \frac{\Delta p}{{\rho U^{2} /2}} $$
(9)

Figure 4 shows the differential force for the base component of the drag force. For the free-stream condition without a water surface, the drag is calculated as

$$ f_{0} = \int\limits_{0}^{2\pi } {\Delta p\cos \theta dA} $$
(10)
Fig. 4
figure 4

Differential force acting on cylinder in free-stream flow for base component of drag force

Full size image

where dA = h(d/2)dθ. By substituting Δp from Eq. (9) into Eq. (10), we obtain

$$ f_{0} = \frac{hd}{4}\rho U^{2} \int\limits_{0}^{2\pi } {C_{p} \cos \theta d\theta } $$
(11)

To obtain the drag coefficient for the base component, CD0, we divide both sides of Eq. (11) by 0.5ρU2hd, yielding

$$ C_{D0} = \frac{1}{2}\int\limits_{0}^{2\pi } {C_{p} \cos \theta d\theta } $$
(12)

It should be noted from the projected area hd that the drag coefficient CD0 does not account for the water surface variation effect (i.e., the drag coefficient for a very small Froude number \( \mathcal{F}\)).

As the flow in an open channel passes around a cylinder, local water surface variations occur. The water surface profile around the cylinder is assumed to be superposed on an undisturbed open-channel flow, i.e., without the objective cylinder, and its variation Δh, as illustrated in Fig. 5. In this figure, the water level of the undisturbed open-channel flow for the objective cylinder is represented by a horizontal dashed line. Here, there is no pressure gradient in the streamwise direction of the undisturbed flow based on an assumption of a uniform flow. The additional term associated with the water surface variation, fs in Eq. (8), is applied such that the other additional term is not considered (i.e., fp = 0). The water surface variations Δh are induced by the pressure variation Δp around the cylinder under the water surface. Given that negative gauge pressures are untenable, the water level is reduced to the height of the zero-gauge pressure.

Fig. 5
figure 5

Schematic of water surface profile around objective circular cylinder and definition of water surface variation Δh

Full size image

To derive the effect of water surface variation around an emergent obstacle on the drag coefficient, we considered a typical pressure distribution around a cylindrical obstacle at an arbitrary water depth h (Fig. 6a). When the hydrostatic pressure in the undisturbed open-channel flow ρg(h − z) (where z is the height of the section) is subtracted from the actual pressure p, the result is Δp, distributed as shown in Fig. 6b. Based on an assumption of a constant Δp in the vertical direction and in the hydrostatic pressure distribution, Δh(θ) is defined as

$$ \Delta h(\theta ) = \frac{\Delta p(\theta )}{{\rho g}} = \frac{{U^{2} }}{2g}C_{p} (\theta ) $$
(13)
Fig. 6
figure 6

Schematic figures of pressure distribution on objective cylinder at several heights defined in Fig. 5: a variation in pressure distribution, p, at depth h − z below water surface at its normal level;b relative pressure distribution Δp; c variation in pressure distribution at lowest water level of fluctuating part; d pressure distribution at a higher water level than in (c)

Full size image

where Cp(θ) is the pressure coefficient.

However, because of the range of variation in the water surface (z >  Δpmin/ρg, Δpmin: minimum value of the relative pressure Δp), the effect of the water surface needs to be considered in the relative pressure distribution. Figure 6c shows the pressure variation at z =  0 where Δp = p, whereas Fig. 6d shows an example between the highest water level (in front) and lowest water level. In this range, a zero-pressure value occurs above the water surface, and the distribution of the actual Δp values is different from that shown in Fig. 6b. This induces an additional force as an effect of the water surface variation in open-channel flow.

To formulate the problem, we focused on the water surface variation component of the drag fs in Eq. (8). Based on sections at heights below the undisturbed water depth h, an additional pressure force Δf is required to compensate for negative pressures at angles with negative Δh, at which ρg(h − z) + Δp is negative. The additional differential force due to the water surface variation component, dfs, is defined as

$$ df_{s} = \left( {\frac{d}{2}d\theta } \right)\frac{\rho g}{2}\Delta h^{2} $$
(14)

Conversely, based on sections above the undisturbed water depth h, the same additional pressure as that appearing in Eq. (14) was applied under the assumption of a hydrostatic pressure distribution. Therefore, Eq. (14) is considered to be a definition of the infinitesimal force required to induce the water surface variation component of drag fs for both negative and positive Δh, and the drag component fs is given by

$$ f_{s} = \int\limits_{0}^{2\pi } {\frac{\rho gd}{4}\Delta h^{2} \cos \theta d\theta } $$
(15)

Substituting Δh from Eq. (13) into Eq. (15), we obtain the water surface variation component of the drag, fs:

$$ f_{s} = \frac{{\rho U^{4} }}{16g}d\int\limits_{0}^{2\pi } {C_{p}^{2} \cos \theta d\theta } $$
(16)

Substituting Eqs. (11), (12), and (16) into Eq. (8) for uniform-flow conditions (fp = 0), we obtain the following, using the drag coefficients CD for fd:

$$ \frac{{\rho U^{2} }}{2}C_{D} hd = \frac{{\rho U^{2} }}{2}C_{D0} hd + \frac{{\rho U^{4} }}{16g}d\int\limits_{0}^{2\pi } {C_{p}^{2} \cos \theta d\theta } $$
(17)

Dividing both sides by ρU2CD0hd/2 and substituting Eq. (12) into Eq. (17) yields the drag coefficient for open-channel flows.

$$ \frac{{C_{D} }}{{C_{D0} }} = 1 + \frac{{\beta_{Cp} }}{2}C_{D0} \mathcal{F}^{2} $$
(18)

The relative drag coefficient CD/CD0 is proportional to the square of the Froude number, with the proportionality factor composed of CD0 and βCp. The coefficient βCp is defined as the second-order moment of the pressure coefficient around the cylinder, expressed as

$$ \beta_{Cp} = \frac{{\int\limits_{0}^{2\pi } {C_{p}^{2} \cos \theta d\theta } }}{{\left( {\int\limits_{0}^{2\pi } {C_{p} \cos } \,\theta d\theta } \right)^{2} }} $$
(19)

5 Experimental results for uniform flows

The drag coefficient for base component of the drag force, CD0, and the coefficient βCp, calculated based on the distribution of the measured pressures on the cylinder, are listed in Table 2 and are compared with values calculated based on free stream-flow conditions [30]. Note that the coefficients βCp for accelerating flow that are outlined in Table 2 include the incremental drag due to the pressure gradient, which is described further in the paper, and are defined differently from the CD0 based on Eq. (8). Figure 7 shows the variations in CD0 and βCp with respect to vegetation density. Whereas the drag coefficient CD0 increases as the vegetation density increases, the coefficient βCp tends to decrease. This is because the pressure coefficient Cp used to determine βCp using Eq. (19) depends on the Reynolds number and the conditions of the cylinder array. Moreover, the symmetry between the variations in the drag coefficient CD0 and the variations in the coefficient βCp reveals a high correlation relationship between them. However, the available data are limited, and further study will be necessary to accurately predict the effect of water surface changes on the drag coefficient. Nonetheless, because all values of βCp calculated herein, as shown in Fig. 7, are negative, it can be inferred that the drag coefficient in an open-channel flow decreases as the Froude number increases, which is consistent with the findings of previous studies [34, 37].

Fig. 7
figure 7

Variations in the drag coefficient for the base component of drag force, CD0, and coefficient βCp with respect to vegetation density λ for high Reynolds number, including data from Etmian et al. [44] for \( \mathcal{R}\) = 1340. Data for λ = 0 are for single cylinder

Full size image

Figure 8a shows the drag coefficient CD as a function of the Froude number \( \mathcal{F}\) for the uniform-flow conditions listed in Table 1. Figure 8a enables comparisons between the experimental values in different hydraulic conditions and theoretical values based on Eq. (18). The drag-force coefficient based on the absence of water surface effects, or \(\mathcal{F} \to 0\), i.e., CD0, is necessary for Fig. 8. The value CD0 = 1.22 is obtained to minimize the root-mean-square error (RMSE) between the measured drag coefficient and the drag coefficient calculated using Eq. (18), with the data subject to \( \mathcal{F}\) < 0.5. The value CD0 = 1.79 in Table 2 is a result of a three-cylinder row array with high velocity along the sidewall, whereas CD0 = 1.22 when a sufficient number of cylinders in a row are used in the experiment. However, the differences in CD0 can be explained by the blockage effect, because the CD0 values 1.22 and 1.79 become 0.76 and 0.81, respectively, after correction for the blockage effect, i.e., C*D/CD= (1 −D/Δy)1.35 (C*D: corrected CD for blockage effect) [51]. For comparison, Fig. 8 also includes previous experimental results from Kothyari et al. [34] and Ishikawa et al. [43] on emergent cylinders in uniform flows, in which the values of CD0 were obtained in the same manner as described previously, i.e., using Eq. (20). In line with the variations in CD0 and βCp shown in Fig. 7, the experimental data are compared with calculation results based on possible ranges of values of CD0 and βCp. Most of the experimental data are explained well within the range of the calculation results, except for a number of scattered values with low \( \mathcal{R}\). This result demonstrated how the influence of the Froude number on the drag coefficient came to be overlooked. However, according to the relationship between CD/CD0 and \( \mathcal{F}\) (Fig. 8a) for higher values of \( \mathcal{F}\), the drag coefficient has an apparent tendency to decrease as the Froude number is increased, even in subcritical-flow conditions. This decrease in CD/CD0 with respect to increasing \( \mathcal{F}\) was also confirmed when we compared the drag coefficients for Cases L3–L7 (as outlined in Table 1), in which \( \mathcal{F}\) varied within a small range of \( \mathcal{R}\). The dependence of CD on \( \mathcal{R}\) and \( \mathcal{F}\) is discussed as follows.

Fig. 8
figure 8

CD/CD0 as a function of Froude number: a effect of \( \mathcal{F}\) on experimental CD data; b comparison between experimental data and theoretical results

Full size image

In Eq. (18), the drag coefficient CD is used in conjunction with CD0, βCp, and \( \mathcal{F}\). The values of CD0 and βCp are defined by the pressure coefficient Cp, as expressed in Eqs. (12) and (19), respectively. Based on an assumption of a hydrostatic pressure distribution, as described earlier, Cp is supposed to be independent of \( \mathcal{F}\), whereas CD is a function of \( \mathcal{F}\) and Cp. On the other hand, Cp around a cylinder can be expressed as a function of \( \mathcal{R}\) based on the similarity principle with respect to the Reynolds number for the Navier–Stokes equations [32]. From this, CD can be considered a function of the two independent variables \( \mathcal{F}\) and R, as CD(\( \mathcal{F}, \mathcal{R}\)). The empirical Ergun formula has been applied to the calculation of drag coefficient on aligned cylinders [25, 26], as follows:

$$ C_{D} = 2\left( {\frac{{\alpha_{0} }}{\mathcal{R}} + \alpha_{1} } \right) $$
(20)

Equation (20) represents the relationship between the drag coefficient and Reynolds number of a single body immersed in a free-stream flow with a relatively high Reynolds number [26]. The coefficients α0 and α1 of the equation have been investigated using experimental data and are known to be a function of the array characteristics.

If the CD0 in Eq. (18) is obtained using Eq. (20), and a constant βCp is assumed, the value of CD can be expressed as a function of \( \mathcal{R}\) and \( \mathcal{F}\) as follows:

$$ C_{D} = 2\left( {\frac{{\alpha_{0} }}{\mathcal{R}} + \alpha_{1} } \right) + \frac{{\beta_{Cp} \mathcal{F}^{2} }}{2}\left\{ {2\left( {\frac{{\alpha_{0} }}{\mathcal{R}} + \alpha_{1} } \right)} \right\}^{2} $$
(21)

Calculating the total derivative of dCD of Eq. (22) produces Eq. (23):

$$ dC_{D} = \frac{{\partial C_{D} }}{{\partial \mathcal{F}}}d\mathcal{F} + \frac{{\partial C_{D} }}{{\partial \mathcal{R}}}d\mathcal{R} $$
(22)
$$ \frac{{dC_{D} }}{{C_{D0} }}\, = - 2\frac{{\alpha_{0} }}{{C_{D0} \mathcal{R}}}\left( {1 + \beta_{Cp} C_{D0} \mathcal{F}^{2} } \right)\frac{{d\mathcal{R}}}{\mathcal{R}} + \beta_{Cp} C_{D0} \mathcal{F}^{2} \frac{{d\mathcal{F}}}{\mathcal{F}} $$
(23)

Referring to Tanino and Nepf [26] for the order of α0, we investigated the dependence of CD on \( \mathcal{F}\;{\text{and}}\;\mathcal{R}\) for CD0 = 1 and α0 = 100 (Fig. 9). Even in subcritical-flow conditions (log10(\( \mathcal{F}\)) < 0), the dependence of CD on \( \mathcal{F}\) in flows with large Reynolds numbers is not negligible relative to that on \( \mathcal{R}\). For example, the effect of \( \mathcal{F}\) becomes apparent at approximately \( \mathcal{F}\) = 0.3 (log10(\( \mathcal{F}\)) = 0.5) for \( \mathcal{R}\) = 103.

Fig. 9
figure 9

Dependence of CD on \( \mathcal{F}\) and \( \mathcal{R}\) for CD0 = 1 and α0 = 100. Labels indicate values of \(b/\sqrt {a^{2} + b^{2} }\), where a, b are the coefficients of d\( \mathcal{F}/\mathcal{R}\), d\( \mathcal{F}/\mathcal{F}\), respectively, when Eq. (22) is rewritten as \(dC_{D} /C_{D0} \, = a(d\mathcal{R}/\mathcal{R)} + b(d\mathcal{F}/\mathcal{F})\)

Full size image

In Fig. 8b, the theoretical lines calculated using Eq. (18) for free-stream, uniform, and accelerating open-channel flow conditions with constant values of βCp (Table 2) explain the decreasing trends exhibited by the experimental results. On the other hand, the decreasing trends exhibited by the experimental values of CD/CD0 with respect to increasing CD0\( \mathcal{F}\)2 became milder compared to the theoretical decreasing trend based on Eq. (18) with constant values of βCp. Overall, the experimental results demonstrated a tendency to deviate from the theoretical lines at higher Froude numbers. In addition to the variation in βCp with hydraulic conditions, this result implies that, at high Froude numbers, Eq. (18) with constant βCp cannot be applied to very shallow-water-depth conditions behind the cylinder. The limitation of Eq. (18) is evident in the negative drag coefficients at high Froude numbers. Based on the uniform-flow condition βCp =  − 0.4, Eq. (18) cannot be applied to flows with CD0\( \mathcal{F}\)2 >  5 because the calculation produces a negative value for CD. More specifically, for the constraint condition that satisfies the positive-water-depth condition h + Δh > 0, the limitation on Eq. (18) is derived from Eq. (13) using Cpmin, which takes the minimum water depth and value of Cp as F2 <  − 2/Cpmin (Cpmin < 0), based on an assumption of free-stream Cp(θ). For the condition F2 >  − 2/Cpmin, the bed behind each cylinder was dried. Equation (13) for the Δh behind the cylinder must be modified to maintain the water depth at h + Δh > 0; thus, a penalty function is applied:

$$ \Delta h = C_{p}^{\prime} \frac{{U^{2} }}{2g},\quad C_{p}^{\prime} = \max \left( {C_{p} ,\; - \frac{2}{{\mathcal{F}^{2} }}} \right). $$
(24)

Equation (24) can be introduced to satisfy the constraint condition h + Δh > 0 for the calculation of βcp using Eq. (19), in conjunction with the measurement of the pressure coefficient distribution based on a uniform open-channel flow (Table 2). At high Froude numbers (F2 >  − 2/Cpmin), the coefficient βCp must increase when the Cp′ in Eq. (24) is used as the Cp in Eq. (19). For this uniform open-channel flow, the value of βCp increases from a constant value (βCp =  − 0.24, − 0.40, − 0.51) for the supercritical-flow condition, eventually becoming positive at high Froude numbers (Fig. 10), indicating that the drag force on the emergent object will increase at high Froude numbers. Based on the variation in the coefficient βCp, the theoretical values for the drag coefficients are modified, as indicated by the dashed lines in Fig. 8b. The relative drag coefficient CD/CD0 attains its minimum value and gradually increases as the Froude number is increased. However, the theoretical lines based on Eq. (24) (dashed lines in Fig. 8) still underestimate the measured drag forces, which deviate from the solid theoretical lines over CD0\( \mathcal{F}\)2 = 1, where there is still little discrepancy between the solid and dashed lines. This discrepancy between the experimental results and the modified theoretical values is attributed to the occurrence of a vertical distribution of pressure change Δp induced by the vertical velocity distribution. By contrast, a uniform velocity distribution and a hydrostatic pressure distribution were assumed in the derivation of Eqs. (18) and (24) even for high Froude numbers. Because the velocity near the water surface is higher than the depth-averaged velocity, the actual effect of the water surface on the drag force is considered to be greater than that assumed in the present theory. Thus, the assumption that Cp is independent of \( \mathcal{F}\), which was used to introduce Eqs. (18) and (21), is no longer satisfied at high \( \mathcal{F}\) conditions. However, although detailed investigations will be necessary for more accurate estimations of the drag coefficient, particularly for flows with high Froude numbers, Eq. (18) can be adopted as a first approximation of the drag coefficient in open-channel flows.

Fig. 10
figure 10

Variation in coefficient βCp with Froude number \( \mathcal{F}\) using pressure coefficient distributions for free stream, uniform open-channel, and accelerating open-channel flows (outlined as conditions in Table 2) in conjunction with the constraint condition of Eq. (24)

Full size image

6 Calculation of water surface profile for flow with emergent cylinders

Before examining the validity of Eq. (18) for calculating water depth profiles in gradually varied flows in conjunction with Eq. (6), we examine the accuracy of the equations for the uniform subcritical flows (L1–L5), as visualized in Fig. 11.

Fig. 11
figure 11

Calculation results based on Eq. (6) and RMSEs of water depth for water surface profiles in uniform-flow conditions. Equation (18) with CD0 = 1.22 is applied to Run U1, whereas constant CD = 1.22 is applied to Run U2. Symbols are measurement results; solid and dashed lines are calculated for Run U1 and U2, respectively

Full size image

Although the coefficient βCp varies with respect to the vegetation density, as indicated in Fig. 7, the coefficient is set to βCp =  −0.4 in the upcoming analysis because the variation with respect to the density does not appear to be significant when the vegetation density is greater than 0.04. Furthermore, we want to focus on the effects of the non-equilibrium open-channel flows at the constant-density condition λ = 0.067. The water surface profiles are calculated using Eq. (6) from the downstream end water depth at x =  − 0.1 m. The drag-force coefficient CD is calculated with Eq. (18) for Run U1, whereas the constant CD = CD0 is applied for Run U2. The RMSEs between the results of the calculation and experiment for L3–L5 are larger than those for L1–L2, indicating that the RMSEs depend on the velocity head. For the calculations, CD0 = 1.22 is adopted, which is obtained earlier in the study for Fig. 8. The RMSEs are affected by measurement errors and the uniform-velocity assumption introduced for Eq. (6). The velocity head causes the local water level to rise and is considered to be the representative length of the water surface measurement error, measured manually using a point gauge. The RMSEs are also affected by errors in the calculation model based on Eq. (6) from the momentum correlation coefficient in the advection term in Eq. (1) induced by the spatiotemporal variation in the velocity. The RMSEs of Run U1 are smaller than those of Run U2 except for Case L3, in which the effect of Eq. (18) for fs on CD is manifested in cases with relatively large \( \mathcal{F}\). The experimental data and discussions in this paper regarding water depth are considered to include an order of error similar to the velocity head, i.e., a few millimeters at the maximum. The water surface profiles for L1–L5 are well reproduced within this scale.

Subsequently, the same calculation method with the same coefficients as those in Run U1 and U2 were applied to gradually varied flows, including accelerating and decelerating flows (Fig. 12). The calculation results from Eq. (6), in conjunction with Eq. (18) for CD and constant CD = CD0, are indicated by dashed lines and small square symbols, respectively. Whereas the water depth profiles for the decelerating flows (b) were properly reproduced by the calculation results, those for the accelerating flows (a) were underestimated. Although the calculation method was designed to consider an error of a few millimeters even for uniform-flow conditions (Fig. 11), the underestimations by the calculations (Fig. 12a) were still significant relative to the error. On the other hand, compared to the differences between the measurement and calculation results of Eq. (6), the differences in the calculated results between Eq. (18) for CD and constant CD = CD0 were negligible compared to the deference from the measured results because relatively high Froude number flow restricted to near the boundary. This indicates that accelerating flows involve a mechanism that increases the drag coefficient or resistance of the obstacles.

Fig. 12
figure 12

Variation in water surface profile for gradually varied flows. Solid symbols are measurement results; solid lines are calculated using Eqs. (27) and (18); dashed lines are calculated using Eqs. (6) and (18) (f = f0 + fs); no-fill, small square symbols are calculated using Eq. (6) with constant CD (f = f0)

Full size image

As shown in Table 2, the pressure measurement around a cylinder also revealed that the base drag coefficient CD0 for the accelerating flow (CD0 = 2.68) was larger than that for the uniform flow (CD0 = 1.79). The base drag coefficients are calculated using the pressure coefficient Cp, which excludes the effect of the water surface from the coefficient, even though the coefficients βCp for accelerating flow include the effect of the pressure gradient, as mentioned in a previous section. Figure 13 compares the distributions of the pressure coefficients for a solitary cylinder in a free stream [30], and for a single cylinder in a cylinder group in uniform and accelerating open-channel flows. The Cp values of these experiments were calculated based on an assumption of Cp = 1 at the stagnation point [44] to obtain an undisturbed water surface level, which is necessary for calculating the Δh distribution. Whereas the distributions of the positive pressure coefficients in front of the cylinder among the three cases are similar, the negative pressure coefficients in the separation zone differ. For the accelerating flow, the pressure in the separation zone decreases to values lower than those for the uniform flow, contributing to an increase in the drag coefficient for the accelerating flow. The mechanism whereby the pressure decreases in the separation zone is explained in Fig. 14, which compares the longitudinal water surface profiles for uniform and accelerating flows. It is assumed that, in the separation zone, the velocity is low [12, 39, 40], and the water depth is nearly constant; furthermore, there is exposure to the water depth just downstream from the separation zone. Therefore, the pressure in the separation zone is less for the accelerating flow than for the uniform flow because of the pressure gradient in the direction of flow.

Fig. 13
figure 13

Distributions of pressure coefficients around a solitary circular cylinder in free-stream flow (circles), and around a single cylinder in a group in uniform flow (squares) and accelerating flow (triangles)

Full size image
Fig. 14
figure 14

Water surface profiles adjacent to circular cylinders in uniform and accelerating flows

Full size image

The pressure reduction in the separation zone investigated in the experiment and visualized in Figs. 13 and 14 indicates that the additional drag force associated with the pressure gradient, fp, in the accelerating flow should be considered in evaluating the drag force fd of Eq. (8). The additional pressure due to the decrease in the separation area is expressed as

$$ f_{p} = - \rho ghdL\frac{dh}{{dx}} $$
(25)

where L is the representative length of the separation zone. Equation (25) is not attributed to the longitudinal buoyancy, which is calculated using dL = πd2/4, but indicates an increasing form drag associated with a decrease in pressure in the separation zone for the negative pressure gradient. Based on an assumption that the representative length is expressed by the cylinder dowel diameter d as L = kd (k, a coefficient), and a substitution of Eq. (25) into Eq. (8), the momentum equation for Eq. (1) can be re-written for steady flow conditions as follows:

$$ U\frac{dU}{{dx}} + g\frac{dh}{{dx}} = gS_{0} - g\frac{{n^{2} U^{2} }}{{h^{4/3} }} - N\left( {C_{D} \frac{{dU^{2} }}{2} - gkd^{2} \frac{dh}{{dx}}} \right) $$
(26)

Equation (26) with a constant unit-width discharge (Eq. 2) yields the following depth profile equation for a flow with emergent vegetation:

$$ \frac{dh}{{dx}} = \frac{{S_{0} - \frac{{n^{2} U^{2} }}{{h^{4/3} }} - N\frac{{C_{D} dU^{2} }}{2g}}}{{1 - \mathcal{F}^{2} - Nkd^{2} }} $$
(27)

To calculate accelerating flows using Eq. (27), the coefficient k should be determined. In this study, the value k was obtained to minimize the discrepancies of water depth between the measurements and calculations based on Eq. (27).

Figure 15 shows the RMSEs of the calculation results at all measurement points against k for all accelerating and decelerating flows. For all cases, the RMSE decreases from k = 0 and then increases after reaching a minimum, except for L10 in decelerating flow, which takes the minimum value at k = 0. The k value yielding the minimum RMSE varied by case, ranging from 0.0 to 2.1. Therefore, for each case, the k value that yielded the minimum RMSE is applied to the calculation using Eq. (27). According to Fig. 15, the calculations based on Eq. (27) with appropriate values k yielded low RMSEs in all cases, and the minimized RMSEs, i.e., 0.29–1.23 mm, were not significant relative to those in uniform flows (Fig. 11). The water surface profiles calculated using Eq. (27) with appropriate k values from Fig. 15 were then compared with those based on the measurement results, as shown in Fig. 12. The results obtained using this calculation method are in good agreement with the measurement results. Although introducing Eq. (25), for the calculation of fp, to the examination of decelerating flows had little effect on the validity of the analysis, the introduction of the additional drag force fp from Eq. (8) into the examination of accelerating flows was demonstrated to be a necessity. It must be recalled that the calculations with Eq. (6), which adjusts the drag coefficient for the uniform flow and accounts for water surface effects with Eq. (18), could not reproduce the water depth profiles for non-uniform flows. Because the practical use of the numerical calculation of flows with obstacles necessitates evaluation of the drag force for non-uniform flows, it is significant that several non-uniform flows can be reproduced using a unified calculation method. It is known that the relative separation length to cylinder diameter varies depending on the hydraulic conditions, including the Reynolds number and obstacle arrangement condition [12, 39, 40], whereas the representative length introduced to evaluate the additional drag force in Eq. (25) is indirectly related to the separation length. On the other hand, with regard to the proper evaluation of L with the appropriate value of k, that remains a question for future research efforts to resolve.

Fig. 15
figure 15

Root-mean-square error (RMSE) of calculated water depths based on Eq. (27) as a function of coefficient k for the representative separation length

Full size image

7 Conclusions

This study proposed additional terms for the drag-force expression with the drag coefficient to include the effects of the water surface and the pressure gradient in an average scale. The effects of adding these terms were investigated through two types of laboratory experiments on uniform and non-uniform open-channel flows with emergent cylinders, focusing on the entire longitudinal water depth profile and the distribution of local pressure around the cylinder. The main results of the study are as follows:

  1. (1)

    Based on the assumption of hydrostatic pressure distribution, it was determined analytically that the drag coefficient CD decreases in proportion to the square of the Froude number. The proportionality factor in the derived formula is composed of the drag coefficient of the base component of the drag force, CD0, and the second-order moment of the pressure coefficient, βCp. The derived formula can reproduce the characteristics of the drag force on emergent cylinders in a uniform flow, as obtained from our experiments and from others previously conducted.

  2. (2)

    For low Froude number flow, the water surface effects on the drag force for the depth profile calculation were negligible in this study. Although the drag coefficient formula for uniform flows can be applied to the calculation of water depth profiles in decelerating flows with emergent cylinders, the same formula led to an underestimation of resistance in accelerating flows, where the pressure in the separation zone decreased because of the negative gradient of the water depth in the longitudinal direction. Therefore, a modified method with additional drag-force terms accounting for the pressure drop in the separation zone was proposed. The accuracy of the acceleration flow reproduced therewith was improved considerably through the incorporation of this additional drag force into the calculation.

As reported in several past publications on the drag force in open-channel flows, the drag coefficient depends on hydraulic conditions and on several factors pertaining to the obstacles, such as their arrangements and shapes. Thus, estimating an appropriate value for the drag coefficient of the base component of the drag force is challenging. In fact, parameterization of the channel resistance coefficients is inevitably necessary in practical applications such as flood flow simulations of rivers [19, 20]. The most practical use of numerical calculations for flows with obstacles necessitates an evaluation of drag forces in non-uniform flows. The main contribution of this research is the development of an appropriate modification of the drag force from the base component for several non-uniform-flow conditions using a simple unified expression, which accounts for the open-channel effects due to water surface variation and pressure gradients in accelerating flows on the drag-force coefficient. It is expected to be particularly useful for resistance evaluation of obstacles in high Froude number flows with strong non-equilibrium and unsteadiness, such as dam-break flows and tsunamis, where both effects of the water surface and pressure gradient on the drag force are considered to be significant. In order to improve the validity of the analysis of the non-equilibrium flow field, it is necessary to clarify the effect of the pressure gradient on the drag force, i.e., the representative length L.