The analytical model consists of (1) an adapted drag formulation for closely-packed cylinder arrays, including blockage and sheltering, and (2) a turbulent kinetic energy balance, necessary to quantify sheltering. The turbulence model builds on the formulation suggested by Nepf25 for vegetation canopies, and incorporates a turbulence production term by flow expansion, and an extended drag formulation in the wake production term. The steps to derive the equations are presented below.
Drag model
The drag forces experienced by an array of cylinders, per unit mass, can be calculated as:
$$begin{aligned} F_{d} = frac{1}{2}c_D a |U|U end{aligned}$$
(1)
where (c_D) is the drag coefficient of a single cylinder, which can be estimated using the empirical expression of White30, given by:
$$begin{aligned} c_D = 1 + 10Re^{-2/3} end{aligned}$$
(2)
where Re is the Reynolds number based on the cylinder diameter and the depth-averaged local flow velocity U. a is the projected plant area per unit volume, defined by Nepf25 as:
$$begin{aligned} a = frac{d h}{h s^2} = frac{d}{s^2} end{aligned}$$
(3)
with d being the cylinder diameter, s the spacing between cylinders, and h the water depth.
The main unknown in Eq. (1) is the local flow velocity U. If a cylinder array is sufficiently sparse, the local flow velocity could be assumed equal to the depth-averaged incoming flow velocity, (U_{infty }), either measured or calculated with a free surface flow model. For denser configurations, the velocity will change as the flow propagates through the array due to (1) flow acceleration between the elements (blockage), and (2) flow deceleration due to the sheltering effects of upstream rows of cylinders. Both effects are illustrated in Fig. 1c. The changes in flow velocity are included by multiplying (U_{infty }) by a blockage factor, (f_b), and a sheltering factor, (f_s):
$$begin{aligned} U = f_b f_s U_{infty } end{aligned}$$
(4)
Inserting both factors in the expression for the drag force results in Eq. (5):
$$begin{aligned} F_{d} = frac{1}{2}c_D a |U|U = frac{1}{2}c_D a f_b^2 f_s^2 |U_{infty }|U_{infty } = frac{1}{2} c_{D,b} a |U_{infty }|U_{infty } end{aligned}$$
(5)
where the changes in velocity have been incorporated in the bulk drag coefficient, (c_{D,b} = c_D f_b^2 f_s^2). This expression provides a direct relationship between the drag coefficient of a single cylinder, (c_D), and bulk drag coefficients (c_{D,b}) measured for cylinder arrays in laboratory experiments. Predicting the drag force thus depends on determining the values of (f_b) and (f_s).
The blockage factor (f_b) can be estimated based on mass conservation through a row of cylinders11, considering that the velocity will increase as the same flow discharge travels through the smaller section between the elements:
$$begin{aligned} U_{infty } A = U_c A_c = f_b U_{infty } A_c end{aligned}$$
(6)
where the total frontal area is (A = h s_y), and (s_y) is the distance between cylinders perpendicular to the flow, center-to-center (see Fig. 1). Subtracting the frontal area of the cylinders from the total area gives the available flow area, (A_c):
$$begin{aligned} A_c = h s_y – h D = h (s_y-d) end{aligned}$$
(7)
Here we are assuming that the water depth is the same just upstream and in between the cylinders. Solving for (f_b) in Eq. (6) results in Eq. (8), see also Etminan et al.11:
$$begin{aligned} f_b = frac{h s_y}{ h (s_y-d)} = frac{1}{1-d/s_y} end{aligned}$$
(8)
The sheltering factor (f_s) can be estimated from the wake flow model by Eames et al.26, which predicts the velocity deficit behind a cylinder as a function of the distance downstream of the cylinder, (s_x), the cylinder diameter, the local turbulent intensity (I_t), and the drag coefficient:
$$begin{aligned} frac{U_{infty }-U_{w}}{U_{infty }} = frac{c_D d}{2sqrt{2 pi } I_t s_x} end{aligned}$$
(9)
where (U_{w}) is the velocity in the cylinder wake, (U_{infty }) is the incoming flow velocity, and (I_t) is the meant turbulent intensity, defined as (I_t = sqrt{k}/U_{infty })21,25. k represents the turbulent kinetic energy per unit mass, with (k = 1/2(overline{u’^{2}} + overline{v’^{2}} + overline{w’^{2}})), where (u’), (v’), and (w’) are the instantaneous velocity fluctuations in the streamwise, lateral, and vertical direction respectively, and where the overbar denotes time averaging. The turbulent velocity fluctuations are defined as the difference between the instantaneous velocities and their mean value over a measurement period. Here we consider the depth-averaged value of the turbulent intensity, in view of the uniformity of the turbulent properties over the vertical observed inside emergent arrays25.
Equation (9) was developed assuming turbulent flow. Viscous effects decrease the velocity deficit26, with the reduction factor being given by:
$$begin{aligned} f_{Re} = sqrt{frac{Re}{Re_{t}}} end{aligned}$$
(10)
where (Re_{t}) is the lowest Reynolds number corresponding to fully turbulent wake flow. Laminar effects are included in the wake flow model by multiplying the velocity deficit of Eq. (9) by the reduction factor (f_{Re}) for (Re < Re_t), where the the turbulent Reynolds number is assumed equal to (Re_t = 1,000). This value is based on the observation that although a wake starts becoming turbulent at (Re_{t} sim 200), drag coefficient measurements usually become constant at Reynolds numbers beyond (Re_{t} sim 1000), e.g. as shown in Figure 2.7 of Sumer and Fredsoe13. The influence of varying (Re_{t}) on the model results is investigated in “Results and discussion” section.
Defining the sheltering factor as (f_s = frac{U_{w}}{U_{infty }}), and including (f_{Re}) and the bulk drag coefficient in the definition of the velocity deficit results in Eq. (11):
$$begin{aligned} f_s = frac{U_{w}}{U_{infty }} = 1-f_{Re}frac{c_{D,b} d}{2sqrt{2 pi } I_t s_x} = 1-f_{Re}frac{c_{D,b} d}{2sqrt{2 pi } (sqrt{k}/U_{infty }) s_x} end{aligned}$$
(11)
Equation (9) also assumes that the downstream cylinder is placed inside the ballistic spreading region of the wake. The ballistic regime occurs for a downstream distance (s_x < L/It), where L is the integral length-scale of turbulence, and it is characterized by a rapidly decaying velocity deficit, and by a linear increase of the wake width with downstream distance. Inside the cylinder arrays, the length scale development is limited by the downstream spacing, resulting in (L < s_x). Considering that turbulent intensity measurements of Jansen29 varied between (I_t) = 0 and 0.8 inside cylinder arrays with n = 0.64–0.9, this would result in (L < s_x/It). This is a reasonable general assumption for the bamboo structures, since their porosity varies in a similar range. If the poles were sparsely placed, there would be a transition from ballistic to diffusive spreading of the wake. Eames et al.26 also developed expressions for turbulent flow under the diffusive regime, which could be used in place of Eq. (9).
In the opposite case of very high pole densities, there may be a point where the elements are so closely-packed that vortex shedding is inhibited by the presence of the neighboring cylinders. Considering an analogy with a cylinder placed close to a solid boundary, vortex shedding would not take place for spanwise spacings smaller than (s_y/d < 1.3)13, causing a decrease of the drag coefficient that would not be reproduced by the expression of White30. The application of the present model is thus restricted to (s_y/d > 1.3).
Balance of turbulent kinetic energy
Application of Eq. (11) requires predicting the turbulent kinetic energy. This is calculated by expanding the model developed by Nepf25, based on a balance between turbulence production and dissipation:
$$begin{aligned} P_w sim epsilon end{aligned}$$
(12)
where (P_w) is the turbulent production rate and (epsilon) is the dissipation rate. For a dense cylinder array, k is produced by (1) generation in the wakes of the cylinders25, and (2) shear production by the jets formed between the elements28. The total turbulence production term, (P_w), consequently has two parts:
$$begin{aligned} P_w = P_{w1}+P_{w2} end{aligned}$$
(13)
We assume that for dense cylinder arrays these two terms are much higher than turbulence production by shear at the bed, based on observations by Nepf25 for sparse arrays. This assumption is further tested in “Results and discussion” section.
The first term in Eq. (13), (P_{w1}), represents turbulence production at the wakes, and can be estimated as the work done by the drag force times the local flow velocity:
$$begin{aligned} P_{w1} = frac{1}{2}c_D a |U|U^2 = frac{1}{2}c_D a f_b^3 f_s^3 |U_{infty }|U_{infty }^2 end{aligned}$$
(14)
The second term, (P_{w2}), represents turbulence generation due to flow expansion28, and can be estimated from the Reynolds shear stresses:
$$begin{aligned} P_{w2} = overline{ u’ v’} frac{partial u }{partial y} end{aligned}$$
(15)
where the overbar denotes time averaging. The loss in mean kinetic energy (E_c) due to flow expansion is equal to:
$$begin{aligned} Delta E_c = frac{1}{2} U_{infty }^2 left( left( frac{A}{A_c}right) ^{2}-1 right) = frac{1}{2} left( f_b^{2}-1 right) U_{infty }^2 end{aligned}$$
(16)
where the energy loss due to flow expansion, (Delta E_c), is modelled using the Carnot losses. Assuming that the mean kinetic energy is transformed into turbulent kinetic energy (E_t), and assuming isotropic turbulence, gives Eq. (17):
$$begin{aligned} frac{1}{2} left( f_b^{2}-1 right) U_{infty }^2 = frac{3}{2}overline{ u’ u’} end{aligned}$$
(17)
Equation (17) enables expressing the normal Reynolds stress as a function of the incoming flow velocities and the blockage factor:
$$begin{aligned} overline{ u’ u’} = frac{1}{3} left( f_b^{2}-1 right) U_{infty }^2 end{aligned}$$
(18)
The Reynolds shear stress is estimated as (overline{ u’ v’} = Roverline{ u’ u’}), where the correlation factor R was given a constant value of 0.4 based on observations of Nezu and Nakagawa31. This value was derived for open channel flow conditions and is assumed acceptable as a first approximation, but it could vary inside a cylinder array. This is explored further in “Results and discussion” section.
The velocity gradient is estimated from the velocity difference between the side of the cylinders (dominated by blockage) and the wake of the cylinders (dominated by sheltering) resulting in Eq. (19):
$$begin{aligned} frac{partial u }{partial y} approx frac{U_{infty }(f_b-f_s)}{frac{1}{2} s_y} end{aligned}$$
(19)
Substitution into Eq. (15) gives Eq. (20):
$$begin{aligned} P_{w2} = frac{2}{3} R (f_b-f_s)(f_b^{2}-1)frac{U_{infty }^3}{s_y} end{aligned}$$
(20)
The dissipation term, (epsilon), is estimated as:
$$begin{aligned} epsilon sim k^{3/2} l^{-1} end{aligned}$$
(21)
The characteristic turbulent length scale l is limited by the surface-to-surface separation between the elements in the flow direction, (l = min(|s_x-d|, d)). This differs from the expression developed by Nepf25, who used the diameter as representative of the size of the eddies. We assume that in closely-packed cylinder arrays the spacing between cylinders may be smaller than the diameter, (|s_x-d| < d), which would limit turbulence development. The maximum value of l is set equal to the cylinder diameter. Here we also assume that for the dense cylinder arrangements, the spacing between cylinders is considerably smaller than the water depth, hence turbulence generated by bed friction is negligible.
Balancing the production and dissipation of turbulent kinetic energy results in Eq. (22):
$$begin{aligned} frac{k^{3/2}}{l} sim |U_{infty }|U_{infty }^2left( c_D a f_b^3 f_s^3 + frac{ 4R}{3s_y}(f_b^{2}-1)(f_b-f_s)right) end{aligned}$$
(22)
Taking the cubic root at both sides and introducing the scale factor (alpha _1) gives Eq. (23):
$$begin{aligned} frac{sqrt{k}}{U_{infty }} = alpha _1left( c_D f_b^3 f_s^3 a l + frac{4}{3}R(f_b^{2}-1)(f_b-f_s)frac{ l}{s_y}right) ^{1/3} end{aligned}$$
(23)
Where (alpha _1) is a coefficient of ({mathcal {O}}(1)), which is given a default value of (alpha _1 = 1). The sensitivity of the model to different (alpha _1) and R values is explored in “Results and discussion” section.
k can be calculated by solving Eq. (23) iteratively, using the incoming upstream velocity (U_{infty }) and the geometric characteristics of the structure, (s_y, s_x, d) and a, as an input. This enables determining the sheltering factor, (f_s = U_{w}/U_{infty }) from Eq. (11). The blockage factor (f_b=(1-d/s_y)^{-1}) can also be calculated from the geometric properties of each configuration. Both coefficients can be then combined to predict the bulk drag coefficient, with (c_D,_{b} =c_D(f_s)^2(f_b)^2). Deriving (c_D,_{b}) with the present approach relies on the assumption that the changes in water depth through the structure are small. This is a reasonable assumption given the short length of the bamboo structures in the streamwise direction, which varies between 0.7 and 1.5 m (see Fig. 1b). Longer structures that experience non-negligible changes in flow depth and velocity should be discretized, and the bulk drag coefficient should be calculated separately for the different sections. The model assumptions are discussed further in the following section.
Source: Ecology - nature.com