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 ($C_D$) 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 $C_D$ 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 $C_D$ are compared to a range of experimental data from previous studies. All functions reasonably estimate $C_D$ at low φ and high cylinder Reynolds numbers (${\mathsf{R}}_{\mathsf{d}}$). At higher φ they typically underestimate $C_D$. Estimates of $C_D$ utilizing the re-parametrization presented here match the experimental data better than estimates of $C_D$ made using the other functions evaluated, particularly at low φ and low ${\mathsf{R}}_{\mathsf{d}}$.

Keywords:

• Cylinder arrays
• drag coefficient
• estimating drag
• porous media
• vegetated flows

1 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 (Sonnenwald, Guymer, & Stovin, 2017). Determining vegetation drag is therefore of interest for a range of hydraulic engineering applications.

1.1 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.

Table 1. Summary of experimental data describing drag in arrays of emergent cylinders

Study	Distribution	φ (–)	d (mm)	${\mathsf{R}}_{\mathsf{d}}$ (–)	Method
Ferreira et al. (2009)	Random	0.022–0.04	10	1185–1452	Equation (1(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) )
Kim and Stoesser (2011)	Staggered	0.022–0.087	10	761–1738	Equation (1(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) )
Tanino and Nepf (2008a)	Random	0.091–0.35	6.4	24–507	Equation (1(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) )
Tinoco and Cowen (2013)	Random	0.005–0.08	3–25	55–3838	Direct force measurement
Meftah and Mossa (2013)	Square	0.003	3	234–607	Equation (2(2) $\frac{\sqrt{k}}{U_p}=\gamma {\left(\frac{1}{(1-\phi )}{C}_{D}ad\right)}^{1/3}$(2) )
Koch and Ladd (1997)	Random	0.05–0.4	-${}^{\mathrm{a}}$	12–99	2D CFD (LBM${}^{\mathrm{b}}$)
Stoesser et al. (2010)	Staggered	0.015–0.251	6.4	133–1789	3D CFD (LES${}^{\mathrm{c}}$)

${}^{\mathrm{a}}$No diameters available.

${}^{\mathrm{b}}$Lattice Boltzmann method.

${}^{\mathrm{c}}$Large eddy simulation.

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: (1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) 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² m⁻³), and $U_{\mathrm{p}}$ is mean interstitial velocity (Stone & Shen, 2002; Tanino, 2012). For cylinders $\phi =ad\pi /4$ where d is cylinder diameter. In low velocities or low cylinder densities Eq. (1(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(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. (1(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) ) (Tanino & Nepf, 2008a). Instead, drag may be measured directly using a force sensor (Dittrich, Aberle, Schoneboom, Rodi, & Uhlmann, 2012; James, Goldbeck, Patini, & Jordanova, 2008; Tinoco & Cowen, 2013). Measured force is then equated directly with the the right hand side of Eq. (1(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) ).

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: (2) $\frac{\sqrt{k}}{U_p}=\gamma {\left(\frac{1}{(1-\phi )}{C}_{D}ad\right)}^{1/3}$(2) where k is turbulent kinetic energy and $\gamma \approx 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(1) $\rho gS(1-\phi )=\frac{1}{2}{C}_{D}a\rho U_{\mathrm{p}}^2$(1) ) from the streamwise pressure gradient or by extracting the force on a cylinder by integrating the pressure acting on the cylinder wall.

1.2 $C_D$ estimation functions

When no physical measurements are available and CFD-based 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 ${\mathsf{R}}_{\mathsf{d}}$, where ${\mathsf{R}}_{\mathsf{d}}=U_{p}d{\nu }^{-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.

Table 2. Equations of functions that estimate $C_D$ for arrays of emergent cylinders

Reference	Description
Cheng (2012)	Drag coefficient based on pseudofluid model modification of a single-cylinder function${}^{\mathrm{a}}$
	$C_D=\left(11{\left({\displaystyle \frac{{\mathsf{R}}_{\mathsf{d}}}{1+80\phi }}\right)}^{-0.75}+0.9\left[1-\exp \left(-{\displaystyle \frac{1000(1+80\phi )}{{\mathsf{R}}_{\mathsf{d}}}}\right)\right]\right.$
	$\left.+1.2\left[1-\exp \left(-{\left({\displaystyle \frac{{\mathsf{R}}_{\mathsf{d}}}{4500(1+80\phi )}}\right)}^{0.7}\right)\right]\right)(1-\phi )$
Ghisalberti and Nepf (2004)	Wake shading numerical model${}^{\mathrm{b}}$
	$C_D={\displaystyle \frac{C_{D\mathrm{White}}}{1.16}}\{1.16-9.31(ad)+38.6(ad{)}^2-59.8\left(ad{)}^3\right\}$
Tanino and Nepf (2008a)	Fit to cylinder array experimental data using Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) )
	$C_D=2(0.46+3.8\phi )$
Tinoco and Cowen (2013)	Fit to cylinder array experimental data using Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) )
	$C_D=2(0.58+6.49\phi )$
White (1991)	Fit to single-cylinder experimental data
	$C_D=1+10.0{\mathsf{R}}_{\mathsf{d}}^{-2/3}$

${}^{\mathrm{a}}$The pseudofluid model has an S term, but this is omitted here. For $S\le 1/100$, differences to estimated $C_D$ value are within ±1%.

${}^{\mathrm{b}}{C}_{D\mathrm{White}}$ refers to the White (1991) function.

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: (3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) where ${\alpha }_0$ and ${\alpha }_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 ${\alpha }_0$ and ${\alpha }_1$. Linking their values of ${\alpha }_0$ and ${\alpha }_1$ to the physical characteristics of their cylinder arrays led both to propose linear relationships for predicting ${\alpha }_1$ as a function of φ. Tanino and Nepf (2008a) noted that ${\alpha }_0$ appeared to be independent of cylinder array characteristics and omitted the viscous term from Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(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 ${\mathsf{R}}_{\mathsf{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(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) ). Sonnenwald, Hart, West, Stovin, and Guymer (2017) re-parametrized ${\alpha }_0$ and ${\alpha }_1$ in terms of φ and d. The objectives of this note are (i) to improve the re-parametrizations of Sonnenwald et al. (2017) by including additional experimental data; (ii) to demonstrate the validity of these re-parametrizations by comparing estimates of $C_D$ made using Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) ) to experimental data; and (iii) to compare alternative estimates of $C_D$ with the re-parametrized Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) ) and with experimental data.

2 A re-parametrization of the Ergun (1952) coefficients

Figure 1 provides a comparison between values of ${\alpha }_0$ and ${\alpha }_1$ and the corresponding values of φ and d for a range of data. Results from Koch and Ladd (1997) are plotted taking lattice units as mm for comparison purposes.

Figure 1. Plots of the Ergun (1952) coefficients with respect to cylinder array characteristics, (a) ${\alpha }_0$ vs φ, (b) ${\alpha }_0$ vs d, (c) ${\alpha }_1$ vs φ, (d) ${\alpha }_1$ vs d; $▹$ Meftah and Mossa (2013),

Stoesser et al. (2010),

Tanino and Nepf (2008a) ○ Tinoco and Cowen (2013), − best-fit Eq. (4a(4a) ${\alpha }_0=6475d+32$(4a) ), ⋄ Koch and Ladd (1997) taking lattice units as mm

Figure 1a does not suggest any systematic relationship between ${\alpha }_0$ and φ, which is consistent with the conclusions of Tanino and Nepf (2008a).

Figure 1b shows a positive correlation between ${\alpha }_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 ${\alpha }_0$. It is therefore reasonable to conclude that the variation in ${\alpha }_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 ${\alpha }_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(4a) ${\alpha }_0=6475d+32$(4a) ), shown in Fig. 1b.

Tanino (2012) suggested that viscous drag, the component described by ${\alpha }_0$, is proportional to d/s, where s is cylinder spacing. A linear relationship with d is consistent with this. No relationship between ${\alpha }_0$ and s (either on its own or with d) was found.

Figure 1c shows a positive correlation between ${\alpha }_1$ and φ, which is consistent with both Tanino and Nepf (2008a) and Tinoco and Cowen (2013) who both suggested a linear relationship between ${\alpha }_1$ and φ. Tanino (2012) suggested that the inertial drag (described by ${\alpha }_1$) is strongly linked to flow-field heterogeneity, and that φ provides a reasonable estimate of this. Figure 1d also shows a positive correlation between ${\alpha }_1$ and d. A linear relationship between ${\alpha }_1$ and d is suggested by the results of Tinoco and Cowen (2013) and Koch and Ladd (1997), similar to ${\alpha }_0$. If ${\alpha }_1$ also depends on d, then d may serve to indicate flow-field heterogeneity.

Together, Figs 1c and 1d suggest that ${\alpha }_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 variation in values of ${\alpha }_1$ for the same d. Least-squares curve-fitting was undertaken assuming ${\alpha }_0=f\left(d\right)$ and ${\alpha }_1=f(d,\phi )$ are linear functions giving: (4a) ${\alpha }_0=6475d+32$(4a) (4b) ${\alpha }_1=17d+3.2\phi +0.50$(4b) where Eq. (4(4a) ${\alpha }_0=6475d+32$(4a) a) provides an estimate of the viscous effects and Eq. (4(4a) ${\alpha }_0=6475d+32$(4a) b) provides an estimate of the inertial effects of drag when used in Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) ). Note that the coefficients to the d terms must have units m⁻¹ to retain non-dimensionality. Root mean square error (RMSE) values of 38.0 and 0.131 were obtained respectively for ${\alpha }_0$ and ${\alpha }_1$. Substituting Eq. (4(4a) ${\alpha }_0=6475d+32$(4a) ) into Eq. (3(3) $C_D=2\left(\frac{{\alpha }_0}{{\mathsf{R}}_{\mathsf{d}}}+{\alpha }_1\right)$(3) ) gives a new function for estimating $C_D$: (5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(5)

3 A comparison of estimates of $C_D$ against experimental results

Figure 2 shows estimates of $C_D$ from the functions in Table 2 and Eq. (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(5) ) plotted for a range of ${\mathsf{R}}_{\mathsf{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 ${\mathsf{R}}_{\mathsf{d}}$ except the Tanino and Nepf (2008a) and Tinoco and Cowen (2013) functions, which exclude a viscous term and are therefore poor estimators of $C_D$ at low ${\mathsf{R}}_{\mathsf{d}}$. At $\phi \lesssim 0.02$, Fig. 2a–d, most of the functions provide good estimates of $C_D$ for ${\mathsf{R}}_{\mathsf{d}}>200$, also suggesting that the standard $C_D\approx 1$ is not unreasonable in this range. At $\phi \lesssim 0.01$ and ${\mathsf{R}}_{\mathsf{d}}<200$, Fig.  2a and 2b, the White (1991), Ghisalberti and Nepf (2004), and Cheng (2012) functions produce similar underestimates of $C_D$. All three are based on single-cylinder formulations of drag. Figure 2a–d show that only Eq. (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(5) ) estimates $C_D$ well at ${\mathsf{R}}_{\mathsf{d}}<200$.

Figure 2. Comparison of experimental values of $C_D$ (data) to estimates (functions) at a selection of different values of φ and d (shown in top right corner of plot)

As φ increases, $\phi ⪆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 φ.

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 ${\mathsf{R}}_{\mathsf{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 ${\mathsf{R}}_{\mathsf{d}}<1000$. The estimates of $C_D$ made with Eq. (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(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 ${\mathsf{R}}_{\mathsf{d}}$, $C_D$ increases with d. Equation (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(5) ) reproduces this trend, justifying the dependence of Eq. (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(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 ${\mathsf{R}}_{\mathsf{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 φ.

Figure 3. Measured $C_D$ compared with estimated $C_D$ using the functions of (a) Cheng (2012), (b) Ghisalberti and Nepf (2004), (c) Tanino and Nepf (2008a), (d) Tinoco and Cowen (2013), (e) White (1991), and (f) Equation (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(5) ); – is a line of equality

Equation (5(5) $C_D=2\left(\frac{6475d+32}{{\mathsf{R}}_{\mathsf{d}}}+17d+3.2\phi +0.50\right)$(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 ${\alpha }_0$ and ${\alpha }_1$ on φ and d suggested by Sonnenwald et al. (2017) is reasonable. Note that these functions have only been tested over the range of $0.003\le \phi \le 0.4$, $0.003\le d\le 0.025$, and $12\le {\mathsf{R}}_{\mathsf{d}}\le 3838$ and care must be taken applying them outside of this range.

4 Conclusions

A re-parametrization of the Ergun-derived coefficients ${\alpha }_0$ and ${\alpha }_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 (${\mathsf{R}}_{\mathsf{d}}$), and generally the standard $C_D\approx 1$ is not unreasonable here. As ${\mathsf{R}}_{\mathsf{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$.

Notations

a	=	frontal facing area (m2 m−3)
$C_D$	=	drag coefficient (–)
d	=	cylinder diameter (m)
g	=	gravity acceleration (m s−2)
k	=	turbulent kinetic energy (m2 s−2)
${\mathsf{R}}_{\mathsf{d}}$	=	cylinder Reynolds number (–)
S	=	channel slope (–)
s	=	cylinder spacing (m)
$U_{\mathrm{p}}$	=	mean interstitial velocity (m s−1)
${\alpha }_0$	=	coefficient describing viscous effects (m−1)
${\alpha }_1$	=	coefficient describing inertial effects (−)
γ	=	turbulence intensity scaling coefficient (–)
ν	=	kinematic viscosity (m2 s−1)
ρ	=	density (kg m−3)
φ	=	solid volume fraction (–)

Acknowledgments

The authors thank the anonymous reviewers for their feedback.

ORCID

Fred Sonnenwald http://orcid.org/0000-0002-2822-0406

Virginia Stovin http://orcid.org/0000-0001-9444-5251

Ian Guymer http://orcid.org/0000-0002-1425-5093

Additional information

Funding

This research was funded by the Engineering and Physical Sciences Research Council (EPSRC) [grants EP/K024442/1, EP/K025589/1 and via the Warwick Impact Acceleration Account EP/K503848/1].