Estimating drag coefficient for arrays of rigid cylinders representing emergent vegetation

ABSTRACT Flow resistance due to vegetation is of interest for a wide variety of hydraulic engineering applications. This note evaluates several practical engineering functions for estimating bulk drag coefficient () for arrays of rigid cylinders, which are commonly used to represent emergent vegetation. Many of the evaluated functions are based on an Ergun-derived expression that relates to two coefficients, describing viscous and inertial effects. A re-parametrization of the Ergun coefficients based on cylinder diameter (d) and solid volume fraction (φ) is presented. Estimates of are compared to a range of experimental data from previous studies. All functions reasonably estimate at low φ and high cylinder Reynolds numbers (). At higher φ they typically underestimate . Estimates of utilizing the re-parametrization presented here match the experimental data better than estimates of made using the other functions evaluated, particularly at low φ and low .


Introduction
Vegetation occurs in many natural and engineered water systems (O'Hare, 2015). In rivers the additional drag caused by vegetation acts to increase flow depths, potentially increasing the risk of flooding (Darby, 1999). In stormwater ponds, the resistance of vegetation has a dominant impact on the flow field and therefore affects treatment potential . Determining vegetation drag is therefore of interest for a range of hydraulic engineering applications.

Existing measurements of C D
Arrays of rigid cylinders are often used to represent emergent vegetation, e.g. Bennett, Pirim, and Barkdoll (2002), Nepf (1999), Rameshwaran and Shiono (2007), Rowiński and Kubrak (2002), Serra, Fernando, and Rodríguez (2004), Tanino and Nepf (2008b) and Tinoco and Cowen (2013). Table 1 presents seven datasets where the bulk drag coefficient, C D , for emergent cylinder arrays has been experimentally or numerically derived. The experimental and numerical methods used for determining C D are described below.
Traditionally, C D is obtained from experimental results for emergent cylinder arrays by equating driving forces with resistance caused by cylinders (Ferreira, Ricardo, & Franca, 2009;Kim & Stoesser, 2011;Tanino & Nepf, 2008a). Assuming wall and bed stresses are negligible, for emergent cylinders this equates to the balance of gravity and drag forces: where ρ is density, g is acceleration due to gravity, S is channel or energy slope, φ is solid volume fraction, a is frontal facing area (the cylinder area perpendicular to the direction of flow per unit volume, m 2 m −3 ), and U p is mean interstitial velocity (Stone & Shen, 2002;Tanino, 2012). For cylinders φ = adπ/4 where d is cylinder diameter. In low velocities or low cylinder densities Eq. (1) is impractical to apply, as it becomes difficult to measure surface slope. Bed and free surface stresses also become more important, eventually invalidating Eq.
As an alternative to direct measurement, Nepf (1999) assumed that turbulence production in vegetation (arrays of cylinders) is equal to dissipation, that drag dominates energy dissipation, and therefore that turbulence intensity can be equated to drag force as: where k is turbulent kinetic energy and γ ≈ 1 (Tanino & Nepf, 2008b). Thus, instantaneous velocity measurements, e.g. from acoustic Doppler velocimetry, may be used to determine k and hence C D (Meftah & Mossa, 2013). For simple geometries, such as a single cylinder or periodic arrays of cylinders, C D may be evaluated using computational fluid dynamics (CFD) tools (Kim & Stoesser, 2011;Koch & Ladd, 1997;Marjoribanks, Hardy, Lane, & Parsons, 2014;Rahman, Karim, & Alim, 2007;Stoesser, Kim, & Diplas, 2010), either by determining S in Eq. (1) from the streamwise pressure gradient or by extracting the force on a cylinder by integrating the pressure acting on the cylinder wall.

C D estimation functions
When no physical measurements are available and CFDbased approaches are infeasible (e.g. a complex geometry) C D must be estimated. It is well established that C D for a single-cylinder is dependent on cylinder Reynolds number R d , where R d = U p dν −1 and ν is kinematic viscosity (Schlichting, Gersten, Krause, Oertel, & Mayes, 1960;White, 1991). For cylinder arrays, C D is also dependent on array characteristics (Nepf, 1999). Table 2 lists several functions that estimate C D depending on array (or vegetation) characteristics. These functions all have a basis in experimental observations and it is of interest to evaluate how successfully they estimate C D . The White (1991) function is included as a base comparison.
The Tanino and Nepf (2008a) and Tinoco and Cowen (2013) functions share a common derivation. Koch and Ladd (1997) showed the Ergun (1952) expression for pressure drop in packed columns to successfully predict drag force. Tanino and Nepf (2008a) related this expression to drag coefficient giving: where α 0 and α 1 are coefficients describing viscous and inertial drag effects respectively. Tanino and Nepf (2008a) and Tinoco and Cowen (2013) used their experimental C D data to estimate α 0 and α 1 . Linking their values of α 0 and α 1 to the physical characteristics of their cylinder arrays led both to propose linear relationships for predicting α 1 as a function of φ. Tanino and Nepf (2008a) noted that α 0 appeared to be independent of cylinder array characteristics and omitted the viscous term from Eq. (3) in their function estimating C D . Tinoco and Cowen (2013) also excluded the viscous component from their function estimating C D and suggest it is most suitable at R d > 1000. Therefore, in both functions C D is solely a function of φ.
The similarity of the methods used in the Tanino and Nepf (2008a) and Tinoco and Cowen (2013) studies presents an opportunity to combine their results and create enhanced estimates of C D from Eq. (3). Sonnenwald, Hart, West, Stovin, and Guymer (2017) re-parametrized α 0 and α 1 in terms of φ and d. The objectives of this note are (i) to improve the re-parametrizations of  by including additional experimental data; (ii) to demonstrate the validity of these re-parametrizations by comparing estimates of C D made using Eq. (3) to experimental data; and (iii) to compare alternative estimates of C D with the re-parametrized Eq. (3) and with experimental data.
(3) C D = 2(0.58 + 6.49φ) White (1991) Fit to single-cylinder experimental data The pseudofluid model has an S term, but this is omitted here. For S ≤ 1/100, differences to estimated C D value are within ± 1%. b C DWhite refers to the White (1991) function.  Figure 1a does not suggest any systematic relationship between α 0 and φ, which is consistent with the conclusions of Tanino and Nepf (2008a). Figure 1b shows a positive correlation between α 0 and d. This is mainly due to the results of Tinoco and Cowen (2013), who varied both φ and d. Tanino and Nepf (2008a), who varied only φ, did not find a relationship with α 0 . It is therefore reasonable to conclude that the variation in α 0 observed by Tinoco and Cowen (2013) is due to d. The results of Koch and Ladd (1997) show a similar trend. Together they suggest a linear relationship between α 0 and d and as a result the data (excluding that of Koch & Ladd, 1997) presented in Fig. 1b have been fit to a linear function, Eq. (4a), shown in Fig. 1b. Tanino (2012) suggested that viscous drag, the component described by α 0 , is proportional to d/s, where s is cylinder spacing. A linear relationship with d is consistent with this. No relationship between α 0 and s (either on its own or with d) was found. Figure 1c shows a positive correlation between α 1 and φ, which is consistent with both Tanino and Nepf (2008a) and Tinoco and Cowen (2013) who both suggested a linear relationship between α 1 and φ. Tanino (2012) suggested that the inertial drag (described by α 1 ) is strongly linked to flow-field heterogeneity, and that φ provides a reasonable estimate of this. Figure 1d also shows a positive correlation between α 1 and d. A linear relationship between α 1 and d is suggested by the results of Tinoco and Cowen (2013) and Koch and Ladd (1997), similar to α 0 . If α 1 also depends on d, then d may serve to indicate flow-field heterogeneity.

A re-parametrization of the Ergun (1952) coefficients
Together, Figs 1c and 1d suggest that α 1 is a function of both φ and d and all data shown in these two figures (excluding that of Koch & Ladd, 1997) have been used to fit a single function (not shown in Fig. 1). Combining these two parameters gives a  Sonnenwald et al. Journal of Hydraulic Research Vol. 57, No. 4 (2019) variation in values of α 1 for the same d. Least-squares curvefitting was undertaken assuming α 0 = f (d) and α 1 = f (d, φ) are linear functions giving: where Eq. (4a) Figure 2 shows estimates of C D from the functions in Table 2 and Eq. (5) plotted for a range of R d , a representative selection of φ and d, and with the experimental data from Table 1. Each sub-figure shows increasing φ. Most functions show the expected dependency of C D on R d except the Tanino and Nepf (2008a) and Tinoco and Cowen (2013) functions, which exclude a viscous term and are therefore poor estimators Cheng (2012) Ghisalberti & Nepf (2004) Equation (5) Ghisalberti and Nepf (2004), (c) Tanino and Nepf (2008a), (d) Tinoco and Cowen (2013), (e) White (1991), and (f) Equation (5); -is a line of equality of C D at low R d . At φ 0.02, Fig. 2a-d, most of the functions provide good estimates of C D for R d > 200, also suggesting that the standard C D ≈ 1 is not unreasonable in this range. At φ 0.01 and R d < 200, Fig. 2a and 2b, the White (1991), Ghisalberti and Nepf (2004), and Cheng (2012)  As φ increases, φ 0.04 in Fig. 2e-h, the differences between the C D values estimated by each function become greater. The White (1991) and Ghisalberti and Nepf (2004) functions consistently underestimate C D . The Ghisalberti and Nepf (2004) function predicts decreasing C D with increasing φ, which is unique among the functions presented here. The Cheng (2012) function, in contrast, fits the data reasonably well at higher φ.

A comparison of estimates of C D against experimental results
The differences between the Tanino and Nepf (2008a) and Tinoco and Cowen (2013) functions become more apparent at higher φ, with the latter estimating greater values of C D . Compared to the experimental results, the Tinoco and Cowen (2013) function performs better at lower values of R d (Fig. 2g) while the Tanino and Nepf (2008a) function performs better at higher values (Fig. 2h). Despite their suggestion otherwise, the Tinoco and Cowen (2013) function performs well at R d < 1000. The estimates of C D made with Eq. (5) fit the data well at higher values of φ.
There are several instances where experimental configurations from the studies in Table 1 overlap such that measurements of C D were taken at the same φ but at different d. Figure 2a, 2b, and 2f show that across multiple R d , C D increases with d. Equation (5) reproduces this trend, justifying the dependence of Eq. (5) on d. It is the only function that consistently fits the experimental data in Fig. 2. Figure 3 provides a direct comparison between measured and estimated C D for each function. The Ghisalberti and Nepf (2004) and White (1991) functions ( Fig. 3b and 3e) consistently underestimate C D with RMSE values of 3.09 and 2.40. The Tanino and Nepf (2008a) function (Fig. 3c) performs better with an RMSE value of 1.66. The Tinoco and Cowen (2013) function (Fig. 3d) appears to perform well with an RMSE value of 1.16, but shows significant scatter. Horizontal bands in Fig. 3c and 3d indicate that the same C D value is estimated at the same φ despite different values of d and R d . The Cheng (2012) function (Fig. 3a) also appears to estimate C D reasonably well, with an RMSE value of 1.28 as it performs less well at higher C D and φ.
Equation (5) (Fig. 3f) has the tightest clustering around the line of equality with an RMSE value of 0.52, showing that of the six functions evaluated it estimates values of C D closest to experimental measurements. Therefore, the dependence of α 0 and α 1 on φ and d suggested by  is reasonable. Note that these functions have only been tested over the range of 0.003 ≤ φ ≤ 0.4, 0.003 ≤ d ≤ 0.025, and 12 ≤ R d ≤ 3838 and care must be taken applying them outside of this range.

Conclusions
A re-parametrization of the Ergun-derived coefficients α 0 and α 1 has been presented. This resulted in a function for estimating drag coefficient (C D ), which has been compared to experimental data alongside several other functions estimating C D in arrays of rigid cylinders representing emergent vegetation. All functions perform well for low solid volume fractions (φ) and high cylinder Reynolds number (R d ), and generally the standard C D ≈ 1 is not unreasonable here. As R d decreases, only those functions that include viscous drag effects provide reasonable results. As φ increases, many of the functions underestimate C D . The function detailed in this study, which includes viscous effects, φ, and also a dependency on cylinder diameter (d), provides improved estimates of C D .