*Article* **A Study on the Friction Factor and Reynolds Number Relationship for Flow in Smooth and Rough Channels**

**Yeon-Moon Choo , Jong-Gu Kim and Sang-Ho Park \***

Department of Civil and Environmental Engineering, Pusan National University, Busandaehak-ro, 63beon-gil, Geumjeong-gu, Busan 46241, Korea; chooyean@naver.com (Y.-M.C.); koreaws@empas.com (J.-G.K.) **\*** Correspondence: sangogo@nate.com; Tel.: +82-051-510-7654

**Abstract:** The shear velocity and friction coefficient for representing the resistance of flow are key factors to determine the flow characteristics of the open-channel flow. Various studies have been conducted in the open-channel flow, but many controversies remain over the form of equation and estimation methods. This is because the equations developed based on theory have not fully interpreted the friction characteristics in an open-channel flow. In this paper, a friction coefficient equation is proposed by using the entropy concept. The proposed equation is determined under the rectangular, the trapezoid, the parabolic round-bottomed triangle, and the parabolic-bottomed triangle open-channel flow conditions. To evaluate the proposed equation, the estimated results are compared with measured data in both the smooth and rough flow conditions. The evaluation results showed that R (correlation coefficient) is found to be above 0.96 in most cases, and the discrepancy ratio analysis results are very close to zero. The advantage of the developed equation is that the energy slope terms are not included, because the determination of the exact value is the most difficult in the open-channel flow. The developed equation uses only the mean velocity and entropy M to estimate the friction loss coefficient, which can be used for maximizing the design efficiency.

**Keywords:** friction coefficient; open-channel flow; entropy; Reynolds number

## **1. Introduction**

Head loss, *hf*, is a very important physical parameter for both the experimental and the theoretical analyses of fluid phenomena. The mechanism of the head loss in the openchannel flow is very complex and is not clearly explained yet. Usually, friction losses can be calculated whenever the friction coefficient, f, is clearly defined by the Darcy-Weisbach equation [1].

The friction coefficient is a key factor to determine the flow velocity in channel flows, which is also important to ensure the optimum hydraulic design. Furthermore, most studies on the friction coefficient f are not performed in an open channel but in a circular pipe flow. However, theoretical research for the open-channel flow was performed in the case of a relatively slow flow with no secondary current and small distribution. The theoretical analysis of the pipe flow was relatively easy compared to that of the open-channel flow [2].

The experimental data from Bazin [3,4] in the late 1800s and Varwick [5] in the early 1900s showed relationships between the friction coefficient and Moody curves in a pipe flow. For the open-channel flow, similar results by Bazin and Varwick [5] were presented. These results showed that the flow characteristics in the laminar, the transition, and the turbulence flows were similar to those in the pipe flow.

The Bazin, Manning [6], and Ganguillet-Kutter [7] equations were developed by using experimental data performed by Bazin and Varwick [5]. Many researchers, including Chezy [8] and Manning [6], developed empirical equations by using observed data in the open channel. However, these equations do not adequately represent the flow characteristics in the open channel. In order to improve the accuracy for representing the flow

**Citation:** Choo, Y.-M.; Kim, J.-G.; Park, S.-H. A Study on the Friction Factor and Reynolds Number Relationship for Flow in Smooth and Rough Channels. *Water* **2021**, *13*, 1714. https://doi.org/10.3390/w13121714

Academic Editor: Mouldi Ben Meftah

Received: 1 June 2021 Accepted: 16 June 2021 Published: 21 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

characteristics in an open channel, a sufficient study on the adequate and reliable analyses is required.

The Manning's equation [6] was less accurate, even in well-flowing controlled manmade waterways. Chow [9] also recommended the adjustment of Manning's roughness coefficient as a function of depth. Studies on the calculation of the friction coefficient and the head losses in the past are still in use or are not actively underway because of difficulties in the physical solution. This is based on studies done about 100 years ago. Additionally, Choo [10] used the entropy concept to derive the mean velocity using Chiu's velocity formula [11,12], which was also used in this paper to derive the friction coefficient.

For the safe design, operation, and management of an increasingly developed and complex water resource facility, research on the calculations of a more accurate friction coefficient should be based on this study.

Therefore, this paper proposes a new theoretical equation to reflect the probabilistic entropy concepts and hydraulic properties. This article proposes a theory using the twodimensional velocity formula and the probabilistic entropy to get the equation of the friction coefficient calculations. Equations based on these theories can be expected to be much more reliable than the empirical equations. The relationship between smooth and rough waterways is shown by means of comparison of the measured and estimated values to verify the accuracy of the determined values.

The results from the developed equation are based on a theoretical background. The friction coefficient was calculated directly based on the guidance equation combined with the physical factors. We found that it can be used to calculate the friction coefficients with very high accuracy without using energy slopes.

## **2. Theoretical Background**

#### *2.1. Existing Friction Loss Coefficient*

The Hagen-Poiseule equation is shown as Equation (1), which calculates the friction head loss and can be written as:

$$h\_L = \frac{8\mu\overline{u}}{\omega\_0 R\_H^4} l \tag{1}$$

where *hL* is the friction head loss, *μ* is the fluid viscosity coefficient, *ω*<sup>0</sup> is the water unit weight, *RH* is the hydraulic radius of a pipe, and *l* is the pipe length.

In this case, the friction head loss coefficient takes form by rearranging Equation (1) with *RH = d*/4, generating Equation (2) (Darcy-Weisbach [1]):

$$h\_L = f \frac{l}{d} \frac{\overline{u}^2}{2g} \tag{2}$$

The representative equation for the friction loss coefficient in uniform flow can be written as Equation (3):

$$f = \frac{8gR\_hS}{V^2} \tag{3}$$

where *f* is the friction coefficient, *V* is the mean velocity, *g* is the acceleration of gravity, *Rh* is the hydraulic radius, *S* = *hf* /*L* is an energy slope, *hf* is the friction head losses, and *L* is the length of a given section.

This equation can be applied in streams close to the uniform flow, because it is difficult to calculate the value accurately in a nonuniform or unsteady flow, because *S* is an energy slope.

In a smooth pipe flow, the relationship can be expressed by using the Blasius [13] equation, written as Equation (4):

$$f = \frac{0.223}{\text{Re}^{0.25}}\tag{4}$$

where *Re* is the Reynold's number.

The equation is limited to valid values of *Re* between 750 and 25,000. For higher values, von Karman [14] developed a general expression modified by Prandtl [15], which matches the data measured by Nikuradse more accurately [16]. The resulting Prandtl-von Karman formula [5,14,15] can be written as Equation (5):

$$\frac{1}{\sqrt{f}} = 2\log\left(\text{Re}\,\sqrt{f}\right) + 0.4\tag{5}$$

Therefore, it is possible to establish a relationship between *f* and *Re* by an experiment in an open-channel flow using the above equations. However, the relation factors to *f* are different between the open-channel flow and the pipe flow, because it is affected by multiple factors, such as free water surfaces in the waterways, hydraulic radius, and water surface slope. Figure 1 shows the relationship between smooth and rough channel flows by analyzing various overseas experimental data [17].

**Figure 1.** Relationship between f and Re for smooth channels (Flow conditions: Laminar, Transitional, and Turbulent [14]).

The shortcomings of the above studies are that it is difficult to calculate the energy slopes correctly in an open-channel flow. In addition, equations should be applied differently according to the scope of the *Re*. For the example, the Prandtl-von Karman's equation has to consider an uncertainty when using Equation (3). This is because it is hard to obtain an accurate flow velocity at the bottom of the channel.

#### *2.2. New Friction Coefficient Using Entropy*

Shannon [18] first defined entropy by function *H*(*x*) and can be written as Equation (6):

$$H(\mathbf{x}) = -\int\_{-\infty}^{+\infty} p(\mathbf{x}) \ln p(\mathbf{x}) d\mathbf{x} \tag{6}$$

where *p*(*x*) is the probability density function, and ln *p*(*x*) is dimensionless but *dx* has dimension.

Equation (6) means maximizing the entropy that represents the uncertainty for *x*, given *p*(*x*) for the continuous state variable *x*. Applying this concept to the water velocity can be written as Equation (7):

$$H(u) = -\int\_0^{u\_{\max}} p(u) \ln p(u) du\tag{7}$$

The following constraints are used, such as the average value and probability, which are available information about *u*, which is instantaneous (point) velocity and can be written as Equations (8) and (9):

$$\int\_{0}^{u\_{\text{max}}} p(u) du = 1 \tag{8}$$

$$\int\_{0}^{u\_{\text{max}}} \mu p(u) du = \overline{u} \tag{9}$$

Arranging the independent constraint conditions can be given as Equation (10):

$$\int\_{a}^{b} \Phi\_{i}(u, p) du \quad i = 1, \ 2 \tag{10}$$

where *a* is the minimum value of *u*, *b* is the maximum value of *u*, *i* is the constraint number (*i* = 1 is Equation (8) and *i* = 2 is Equation (9)).

Therefore, *p*(*u*), which maximizes the entropy, can be obtained using the method of Lagrange as Equations (11)–(13):

$$\frac{\partial I(\boldsymbol{\mu}, p)}{\partial p} + \sum\_{i=1}^{2} \lambda\_i \frac{\partial \phi\_i(\boldsymbol{\mu}, p)}{\partial p} = 0 \tag{11}$$

$$I(u, p) = -p(u) \ln p(u) \tag{12}$$

where *φ*1(*u*, *p*) = *p*(*u*), *φ*2(*u*, *p*) = *u* × *p*(*u*).

$$\frac{\partial \phi\_1(u, p)}{\partial p} = 1, \quad \frac{\partial \phi\_2(u, p)}{\partial p} = u \tag{13}$$

where *λ*<sup>1</sup> and *λ*<sup>2</sup> are Lagrange multipliers.

Substituting Equations (12) and (13) into Equation (11) can be constructed as the following Equation (14):

$$-1 - \ln p(\mu) + \lambda\_1 + \lambda\_2 \mu = 0\tag{14}$$

where *λ*<sup>1</sup> − 1 = *a*<sup>1</sup> and *λ*<sup>2</sup> = *a*<sup>2</sup> are the Lagrange multipliers.

Differentiating Equation (14) with respect to *p*(*u*) results in the velocity as Equation (15):

$$p(\mu) = e^{a\_1 + a\_2 \mu} \ (0 \le \mu \le \mu\_{\max}) \tag{15}$$

Equation (15) and *M* = *a*2*umax* (entropy coefficient) are substituted into Equation (8) to obtain Equation (16):

$$
\mu\_{\max} e^{a\_1} = \frac{M}{(e^M - 1)}\tag{16}
$$

Then Equation (15) and *M* = *a*2*umax* are substituted into Equation (9) to obtain Equation (17) (This is the two-dimensional average velocity equation, which is Chiu's velocity equation [11,12]):

$$\frac{\overline{u}}{u\_{\max}} = \phi(M) = \left(\frac{e^M}{(e^M - 1)} - \frac{1}{M}\right) \tag{17}$$

Equation (17) can be used to restructure Equation (16) to obtain Equation (18):

$$e^{a\_1} = \frac{\phi(M)M}{(e^M - 1)\overline{\mu}}\tag{18}$$

The shear stress is the product of the dynamic viscosity (kinematic viscosity is the dynamic viscosity divided by density) and the velocity gradient, which can be expressed as Equation (19):

$$
\pi = \rho \upsilon \left(\frac{du}{\overline{h}\_{\xi} d\xi}\right) \tag{19}
$$

where *τ* is the shear stress, *ρ* is the density of the fluid, *ν* is the kinematic viscosity of the fluid, *h<sup>ξ</sup>* is the mean value of *h<sup>ξ</sup>* , and *h<sup>ξ</sup>* is the scale factor, which has the length dimensions.

The shear stress at the channel boundary (bottom) is the shear stress when *ξ* is *ξ*0, as in Equation (20):

$$
\sigma\_0 = \rho \upsilon \left(\frac{du}{\overline{h}\_\xi d\xi}\right)\_{\xi = \xi\_0} = \rho g R\_h S\_f \tag{20}
$$

where *τ*<sup>0</sup> is the waterway boundary shear stress, *g* is the gravitational acceleration, and *Sf* is the energy gradient.

The velocity cumulative probability of *u* is suggested by Chiu [11,12] as Equation (21):

$$P(u) = \int\_0^u p(u) du = \frac{\mathfrak{f} - \mathfrak{f}\_0}{\mathfrak{f}\_{\max} - \mathfrak{f}\_0} \tag{21}$$

where *ξ* is the spatial coordinates (0 ≤ *ξ* ≤ 1), *u* is the velocity at *ξ*, *ξ*<sup>0</sup> is the minimum value of *ξ* (occurring at the channel boundary where *u* = 0), and *ξmax* is the maximum value of *ξ* (where *u* is at its maximum (i.e., *umax*)) (see Figure 2).

**Figure 2.** *ξ*-*η* of an open-channel flow (Chiu [11,12]).

The *ξ*-*η* coordinates are the isovel system, which was first developed by Chiu [11,12] to explain two-dimensional velocity distribution in the cross-section of an open channel.

In Equation (21), if *ξ* is *ξ*<sup>0</sup> at the bottom of the channel, *u* is 0 and *ξmax* − *ξ*<sup>0</sup> is 1, and by differentiating the velocity gradient in the channel bed, where *p*(*u*) = *ea*<sup>1</sup> , Equation (22) can be obtained:

$$\left(\frac{du}{\overline{h}\_{\xi}d\xi}\right)\_{\xi=\xi\_{0}} = \frac{1}{\overline{h}\_{\xi}e^{a\_{1}}}\tag{22}$$

For a wide channel, *ξ*/*ξ*max can be simplified to *y*/*D* and, hence, *h<sup>ξ</sup>* = *D*. For the hydraulic radius, *Rh* ∼= *D*. Therefore, substituting Equation (22) into Equation (20), Equation (23) is obtained:

$$\mathcal{e}^{a\_1} = \frac{\nu}{\mathcal{g}\mathcal{R}\_h^2 \mathcal{S}\_f} \tag{23}$$

Equating Equations (18) and Equation (22) expresses the average water velocity in the open-channel flow as Equation (24):

$$
\overline{\mu} = \frac{\text{g} \mathbb{R}\_h^2 \mathbb{S}\_f}{\nu F(M)} \tag{24}
$$

where *F*(*M*) = *<sup>e</sup><sup>M</sup>* − <sup>1</sup> /(*φ*(*M*)*M*).

For the friction velocity (*u*∗), the relationship between the average water velocity and the friction velocity is shown as Equation (25):

$$
\mu\_\* = \sqrt{\frac{\tau\_0}{\rho}}\tag{25}
$$

To calculate the friction term in Equation (25), Choo's mean velocity distribution [10] is used for Equation (26):

$$\mu = \frac{\overline{\pi}}{K(M)} \ln \left[ 1 + \left( e^M - 1 \right) \left( \frac{\widetilde{\xi} - \widetilde{\xi}\_0}{\widetilde{\xi}\_{\max} - \widetilde{\xi}\_0} \right) \right] \tag{26}$$

where *K*(*M*) = *φ*(*M*)*M*.

Choo's mean velocity was used earlier for calculating the discharge. However, in this paper, it will be used for converting friction velocity, since it has already been modified for the average water velocity in the open-channel flow [10].

The water velocity slope *du*/*hξdξ* is differentiated from Equation (26), and *ξ*<sup>0</sup> 0 and *ξmax* − *ξ*<sup>0</sup> = 1 are applied as Equation (27):

$$\left(\frac{du}{\overline{h}\_{\tilde{\xi}}d\tilde{\xi}}\right) = \frac{\overline{u}(e^{M} - 1)}{R\_{h}K(M)(1 + \tilde{\xi}(e^{M} - 1))}\tag{27}$$

where, because the bottom boundary layer *ξ*<sup>0</sup> = 0, Equation (27) is equal to Equation (28):

$$\left(\frac{du}{\overline{h}\_\circ d\overline{\xi}}\right)\_{\overline{\xi}=\overline{\xi}\_0} = \frac{\overline{\pi}\left(e^M - 1\right)}{\overline{R}\_h \overline{\kappa}(M)}\tag{28}$$

Equation (28) is inserted into Equation (20), which is the shear stress at the channel boundary, to obtain Equation (29):

$$
\pi\_0 = \rho \nu \left( \frac{du}{\overline{h}\_\circ d\overline{\xi}} \right)\_{\overline{\xi} = \overline{\xi}\_0} = \rho \nu \frac{\pi \left( \epsilon^M - 1 \right)}{R\_h \mathcal{K}(M)} \tag{29}
$$

The relationship between the average water velocity and the friction velocity of the friction loss coefficient of the pipe flow is shown in Equation (30):

$$\frac{\overline{u}}{u\_\*} = \sqrt{\frac{8}{f}}\tag{30}$$

Equations (25) and (30) are substituted into Equation (29) to obtain Equation (31):

$$\tau\_0 = \frac{8\rho v^2}{f} \left( \frac{\left(e^M - 1\right)}{R\_h K(M)} \right)^2 \tag{31}$$

Therefore, if Equation (31) is substituted for Equations (24) (*F*(*M*) = *<sup>e</sup><sup>M</sup>* − <sup>1</sup> /(*φ*(*M*)*M*)) and (26) (*K*(*M*) = *φ*(*M*)*M*) and use *Re* = (*ud*)/*υ*, Equation (32) can be obtained:

$$f = \frac{8dF(M)}{\text{Re}R\_h} \tag{32}$$

Equation (32) can be used to estimate the frictional loss coefficient (f) of the openchannel flow, which reflects its entropy. Equation (32) does not require the hydraulic factors used in the existing equations, such as shear velocity (*u*∗) or energy gradient (*Sf*). In addition, the friction loss coefficient (*f*) can be expressed with only the average water velocity and the entropy *M*, which are easy to obtain. In addition, there is also an advantage in that the energy gradient (*Sf*) can be estimated by using Equation (32).

Therefore, in this study, we proposed the friction loss coefficient of Equation (32) in the open-channel flow by using the concept of entropy, which has been used in many fields recently. The data used to demonstrate the utility of the equation were obtained by Yuen [19] and Babaeyan-Koopaei [20] for each stream of water. It is shown in the Figure 1 that the estimated friction loss coefficient was compared with the measured friction loss coefficient.

## **3. Experimental Data**

To evaluate the accuracy of the proposed equation, we calculated the friction coefficient based on the data measured at the rectangular channel. The estimated results were compared with the measured data, as shown in Figure 1. First, the data measured by Yuen [19] at the trapezoidal section were used. Then, the data were measured by Babaeyan-Koopaei [20] at the trapezoid, the parabolic round-bottomed triangle, and the parabolic-bottomed triangle trapezoidal channel.

Yuen obtained data in a fully developed turbulent flow of the smooth trapezoidal open-channel flow. The ranges of the data were 0.5 < *Fr* < 3.5, 1.9·104 < *Re* < 6.2·10<sup>5</sup> and 0.3 < 2*b*/*H* < 15 (where 2*b*/*H* was the aspect ratio). The subcritical flow was also studied for the compound trapezoidal channel, which ranged in depths of 0.05 < *Dr* < 0.5. Here, *Dr* is the relative depth ratio ((*H*−*h*) *<sup>H</sup>* , where *H* is the flow depth of the channel, and *h* is the depth of the lower main channel). Several series of experiments were undertaken by using the Preston tube technique. These experiments were performed in a 21.26-m-long tilting channel with a working cross-section of 0.615 m wide × 0.365 m deep. A total of three sets were measured under the equivalent conditions, varying the bed slope at 0.001, 0.004, 0.009, 0.015, and 0.023. In addition, the point velocities were measured across the whole cross-section for the selected flow depths. Particular attention was focused on understanding the Reynolds and Froude number effects on these distributions.

Babaeyan-Koopaei [20] measured the data in the trapezoidal, parabolic, round-bottomed triangle, and parabolic-bottomed triangle channel. For each section, the measured data were used with the changes in the flow velocity and water levels under three flow conditions: 1 m3/s, 10 m3/s, and 100 m3/s (see Babaeyan-Koopaei for more information).

The values of the measured Re data are shown in Table 1. It can be seen that the measured data in the rectangular section reflected the transition zone and the turbulence zone. In the trapezoidal, parabolic, round-bottomed triangle, and parabolic-bottomed triangle sections, the measured data reflected the full turbulence zone.


**Table 1.** The range of Reynolds numbers with the cross-section shape and the channel slope.

## **4. Estimation of the Entropy Parameter, M**

An estimate of the entropy parameter *M* is needed to use Equation (29). For the estimation of entropy parameter *M*, most researches used Equation (15) to calculate the entropy parameters in which the equation required, essentially, the maximum flow velocity in an open channel.

However, the maximum velocity occurred at the center of the pipe flow, but the location of the maximum velocity was unclear at the open-channel flow. Additionally, a lot of manpower, time, and effort were required to measure the maximum velocity in the openchannel flow. Moramarco [21] calculated the *M* values by using Equation (15). For that, he used data obtained from the average and maximum velocities at the upper river basin. Moramarco [22] proposed an equation for calculating *φ*(*M*) by substituting Chiu's theory and the Manning and Prandtl-von Karman equations. However, the disadvantage of these equations were that it was difficult to clearly identify the point where the maximum velocity occurred, *ymax*, and imaginary distance, *y*0, where the velocity was zero in the riverbed.

This study determined the entropy parameter *M* by using the expression developed by Choo [23]. The advantage of this method was that the entropy parameters in the stream could be obtained at any time without using the uncertain maximum velocity.

The entropy parameters *M* and the *Re* were calculated by using the same characteristics as those shown in Figures 3 and 4. As the entropy parameter, *M* was increased, and the *Re* was also increased. On the other hand, as the friction coefficient *f* increased, *Re* decreased. Based on the value of *Re*, the two flows were identified as turbulent flows.

**Figure 3.** Relationship between *M* and *Re* calculated using Yuen's data.

**Figure 4.** Relationship between *M* and *Re* calculated using Babaeyan-Koopaei's data.

## **5. Results Analysis**

The entropy parameter M, defined in Section 4, was used in Equation (32) to calculate the coefficient of friction *f* in an open-channel flow. The relationship between *f* and *Re* is shown in Figures 5 and 6.

**Figure 5.** The relationship between *f* and *Re* calculated using Yuen's data.

The coefficient of friction f in Figures 5 and 6 shows the same trend as in Figure 1, where f tends to decrease as Re increases. In addition, in Figures 7 and 8, the friction coefficient, *f* , shows a tendency to decrease as the discharge increases. The discrepancy ratio is the error ratio between the measurement and calculated values, separated by a range. The proportions on the y-axis are the ratio of the total comparison quantity to the range of the x-axis. Figures 9–12 show that the discrepancy ratio results of the proposed equation were all distributed near 0.

**Figure 6.** The relationship of *f* and *Re* calculated using Babaeyan-Koopaei's data.

**Figure 7.** The relationship of *f* and *Q* calculated using Yuen's data.

**Figure 8.** The relationship of *f* and *Q* calculated using Babaeyan-Koopaei's data.

**Figure 9.** The discrepancy ratio for *f* and *Re* calculated using Yuen's data.

**Figure 10.** The discrepancy ratio for *f* and *Re* calculated using Babaeyan-Koopaei's data.

**Figure 11.** The discrepancy ratio for *f* and *Q* calculated using Yuen's data.

The above results are summarized as follows: The entropy parameter *M* is a linear function of *Log*(*Re*), which is an increasing function of *Re*. The friction coefficient f is a linear decreasing function of *Log*(*Re*), which is a decreasing function of *Re*. Figure 13 shows the above results calculated using the proposed Equation (32), along with the relationship between the friction coefficients *f* and *Re*. Comparing the scale with a previous empirical study resulted in rough flows; the coefficient of determination was observed to be 0.8 within the range of the flow of the rectangular channel and 0.75 within the range of the flow of the compound trapezoidal channel.

**Figure 12.** The discrepancy ratio for *f* and *Q* calculated using Babaeyan-Koopaei's data.

**Figure 13.** The comparison between the empirical and calculated friction coefficients using Yuen and Babaeyan-Koopaei (rough channel).

Figure 10 shows the results determined with Yuen's data measured in the rectangular and trapezoidal channels. Particularly, the calculated *f* value was expressed as 1.7 × 105 < *Re* < 6.18 × <sup>10</sup>6, which exceeded the existing *Re* value. This range means that the proposed equation can represent the actual phenomenon in the natural stream.

The picture on the right shows the results determined with Babaeyan-Koopaei's data measured in the trapezoidal, parabolic, round-bottomed triangle, and parabolic-bottomed triangle channels. In this case, the calculated *<sup>f</sup>* value was expressed as 13.5 × 104 < *Re* < 46.84 × <sup>10</sup>5, which was close to the maximum range of the previous *Re* value. The important thing is that even the compound trapezoidal sections that were similar to the natural section, which were not available in the existing graph, as shown in Figure 1, can easily be calculated for *f* and *Re*. In Tables 2 and 3, the regression equations from the Yuen [19] and Babaeyan-Koopaei [20] data show the relationship between *M*, *f* , *Re*, and *Q*.


**Table 2.** The results for the relationship between *M*, *f* , *Re*, and *Q* using Yuen's data.

**Table 3.** The results for the relationship between *M*, *f* , *Re*, and *Q* using Babaeyan-Koopaei's data.


Comparing the scale with a previous empirical study results in rough flows; the coefficient of determination was observed to be 0.8 within the range of the flow of the rectangular channel and 0.75 within the range of the flow of the compound trapezoidal channel.

As can be seen in Figure 13, the values of the coefficient of friction determined from past experiences are properly correlated with the measured data. There is no bouncing value on the graph. For the rectangular channel, this seems expressed fairly well as Bazin no. 17. Other types of channels matched well with the extended lines of the Prandtl-von Karman equation. In other words, it is meaningful that the values of the friction coefficient determined from proposed equation can be accurately estimate based on a theoretical formula rather than on an empirical parameter under a number of conditions.

#### **6. Conclusions**

The results of a study conducted approximately 100 years ago are still used to estimate the friction coefficient in an open-channel flow. However, as with the pipe flow (perfusion), the friction coefficient must be correctly determined in order to interpret the correct flow. This paper proposes a new form of friction coefficient calculation by using the two-dimensional velocity formula of Chiu [11,12] and probabilistic entropy.

The advantage of this equation is that it eliminates the terms of energy slopes, which are difficult to measure or calculate in an open-channel flow, making their application simple and very accurate on a theoretical basis.

In uniform flow conditions, a channel bed gradient may be the same or almost the same as an energy slope or water surface gradient. The normal depth is maintained as long as the slope, cross-section, and the surface roughness of the channel remains unchanged; thus, the average flow velocity remains constant. However, in natural flow and human-made open channels, such as irrigation systems and sewer lines, there are mostly nonuniform or unsteady flows. Unlike uniform flow conditions, these varied flows do not share the same energy slope, bed gradient, and water surface gradient.

Based on the data measured in the rectangular section, the proposed equation was used to determine the entropy parameter *M* and the friction coefficient *f* . The induced entropy parameters were shown to be a linear function of *Log*(*Re*), and the friction coefficient was the decreasing function of *Log*(*Re*).

If this study is to be carried out continuously by hydraulic data measured in various channel shapes, laboratory channels, and natural streams, the friction coefficient value estimated from the proposed equation will be actively used in the flow analysis and the design of hydraulic structures.

**Author Contributions:** Y.-M.C. and S.-H.P. carried out the survey of the previous studies. Y.-M.C. wrote the manuscript. S.-H.P. conducted all the simulations. Y.-M.C., J.-G.K., and S.-H.P. conceived the original idea of the proposed method. All authors have read and agreed to the published manuscript.

**Funding:** This research was funded by the Institute of Industrial Technology, Pusan National University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** This research was supported by the Institute of Industrial Technology, Pusan National University.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Ying-Tien Lin 1,2,3 , Yu Yang 2, Yu-Jia Chiu 4,\* and Xiaoyan Ji <sup>5</sup>**


**Abstract:** This study experimentally and numerically investigated the hydrodynamic characteristics of a 180◦ curved open channel over rough bed under the condition of constant downstream water depth. Three different sizes of bed particles (the small, middle and big cases based upon the grain size diameter *D*50) were selected for flume tests. Three-dimensional instantaneous velocities obtained by the acoustic Doppler velocimeter (ADV) were used to analyze hydrodynamic characteristics. Additionally, the Renormalization-Group (RNG) turbulence model was employed for numerical simulations. Experimental results show that rough bed strengthens turbulence and increases turbulent kinetic energy along curved channels. The power spectra of the longitudinal velocity fluctuation satisfy the classic Kolmogorov −5/3 law in the inertial subrange, and the existence of rough bed shortens the inertial subrange and causes the flow reach the viscous dissipation range in advance. The contributions of sweeps and ejections are more important than those of the outward and inward interactions over a rough bed for the middle case. Flow-3D was adopted to simulate flow patterns on two rough bed settings with same surface roughness (skin drag) but different bed shapes (form drag): one is bed covered with thick bottom sediment layers along the curved part of the flume (the big case) as the experimental condition, and the other one is uniform bed along the entire flume (called the big case\_flat only for simulations). Numerical simulations reveal that the secondary flow is confined to the near-bed area and the intensity of secondary flow is improved for both rough bed cases, possibly causing more serious bed erosion along a curved channel. In addition, the thick bottom sediments (the big case), i.e., larger form drag, can enhance turbulence strength near bed regions, enlarge the transverse range of secondary flow, and delay the shifting of the core region of maximum longitudinal velocity towards the concave bank.

**Keywords:** acoustic Doppler velocimeter; rough bed; secondary flow; turbulent bursting; turbulence kinetic energy

## **1. Introduction**

Curved channels are commonly found in natural rivers around the world. Flow in curved open-channel follows a helicoidal path and induces a secondary flow, which can redistribute mass, momentum, boundary shear stress and thereby plays an important role in hydraulic engineering [1,2]. Considerable studies focused on the curved channel characteristics in different ways, including theoretical derivation, field observations and measurements, physical model experiments, and numerical simulations. For example, Zeng et al. [3] showed that the turbulent kinetic energy in the channel bend is significantly larger than that at the entrance and the exit according to laboratory experiments. Blanckaert [4] and Blanckaert and Vriend [5] found that water surface gradient and streamwise

**Citation:** Lin, Y.-T.; Yang, Y.; Chiu, Y.-J.; Ji, X. Hydrodynamic Characteristics of Flow in a Strongly Curved Channel with Gravel Beds. *Water* **2021**, *13*, 1519. https://doi.org/ 10.3390/w13111519

Academic Editor: Mouldi BEN MEFTAH

Received: 6 April 2021 Accepted: 26 May 2021 Published: 28 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

curvature variations are the major factors for velocity redistribution in sharp bend channels. Vaghefi et al. [6] indicated that there are two persistent clockwise vortexes from the streamlines drawn in cross sections and the maximum shear stress occurs from the entrance of the bend to the bend apex area near the inner wall. However, these studies only focus on flow over smooth curved channels, and the hydrodynamic characteristics over a curved and rough bed, more consistent with practical field conditions, has been paid less attention and is not well-understood [7].

Over the past few decades, most of studies focus on the flow patterns over a straight and rough bed open channel. Grass et al. [8] found that the typical scale of near-wall vortical structure is directly proportional to the bed roughness size under fully developed flow conditions. Ferro [9,10] revealed that the bed roughness can change the mean velocity profile near the wall and the friction coefficient. Nikora et al. [11,12] applied the doubleaveraging method to large and organized roughness element bed flows and showed that the intrinsic averaged streamwise velocity profile is linear inside the bed roughness with high submergence. Mignot et al. [13] concluded that the strongest turbulence activity occurs at gravel crest levels *zc*/*h* = 0.1, where *zc* is the distance from the bed and *h* is the flow depth, and the turbulent diffusion also reaches a maximum at this elevation. Dey et al. [14] showed that the friction factor decreases with bed-load transport substantiating the concept of reduction of flow resistance and the downward vertical fluxes of the *TKE* values increase in presence of bed-load transport. Qi et al. [15] used the Particle Image Velocimetry (PIV) techniques to measure the near-wall rough regions and performed that the *TKE* value firstly increases and then decreases as water depth increases. Most of these literatures aim to understand the effect of rough bed on straight open channels. In order to better know the flow patterns in natural rivers, the study of curved open channel flow over rough bed is highly required.

The flow patterns of curved open channel flow over rough bed are more complicated than those over smooth bed. Up to now, the related studies are relatively fewer in comparison with the straight and rough bed channels. The studies can be divided into two types. The first one focuses on the roughness of the side wall in curved channels. In this type, the longitudinal velocity decreases and the shear stress increases in the bank vicinity [16]. In addition, a reversed secondary flow occurs for all rough outer bank cases is found and the reversed secondary flow becomes is stronger as the roughness of the outer bank increases [17]. Hersberger et al. [18] suggested that the macroroughness placed on the outer side wall causes the maximum velocity cell move towards the center of the channel and also reduces the sediment transport capacity in the bend. The other one is that the roughness is due to bed sediment particles over curved channels. Jamieson et al. [19] showed that the magnitude and distribution of three-dimensional Reynolds stresses increase through a 135◦ channel bend with a mobile sand bed. Pradhan et al. [20] revealed that the resistance caused by the curvature is larger on smoother channels than that on fine gravel channels at the meander bends. These studies have provided some useful information on the flow patterns over a rough and curved channel. However, these results may not be able to fully reflect some practical situations. For example, on the purpose of navigation the water level at the downstream usually need to maintain constant, which is the main focus in this study.

In this paper, laboratory experiments and numerical simulations were conducted to investigate the flow field in a 180 degree U-shaped curved and rough bed flume under the condition of constant downstream water level. Experimental investigations were carried out to study the profile of longitudinal velocity, turbulent kinetic energy, and turbulent bursting. Gravel sediment particles with different grain sizes were used to present the level of bed roughness. ADV was used to obtain three dimensional instantaneous velocities at different cross sections along the curved channel. In addition, the Renormalization-Group (RNG) turbulence model in the FLOW-3D software was then utilized to discuss flow patterns on two rough bed settings with same surface roughness (skin drag) but different bed shapes (form drag): one is bed covered with a thick bottom sediment layer along the curved part of the flume as the experimental condition, and the other one is uniform bed

along the entire flume. The paper was organized as follows: in Section 2 we detailed the experimental set-up and the numerical modelling. In Section 3, the experimental results including longitudinal velocity distribution, turbulent kinetic energy and power spectral density, and turbulent bursting under different bed roughness were performed. In addition, numerical results including water depth, longitudinal velocity, *TKE* (*k*) and secondary flow were presented and used to discuss the effects of thick sediment layers. Finally, conclusion and future prospects were given in Section 4.

## **2. Methodology**

## *2.1. Experimental Setup*

The experiments were carried out in a recirculating Plexiglas flume of 0.40 m wide, 0.40 m deep and 33 m long in the Ocean Engineering Laboratory in Zhoushan Campus, Zhejiang University, China. The flume consists of a 12 m straight inflow reach, a 180 degree U-shaped curved reach with the radius of curvature *R* in the centre position of the channel equal to 1.4 m, and a 16 m straight outflow reach with a constant slope of 0.5 % of entire flume. A series of honeycomb grids were installed at the entrance of the channel to stabilize the flow and prevent the formation of large-scale flow disturbances. Three different discharges *Q* were set at 0.015 m3 s<sup>−</sup>1, 0.025 m3 s−<sup>1</sup> and 0.030 m3 s<sup>−</sup>1, respectively. The flume has an adjustable tailgate at the downstream end of the flume, where the water level was set to 35 cm, to regulate the flow depth. The test sections were located in the curved region of the flume, where the flow was in fully developed turbulent regimes.

Nortek Vectrino ADV (Nortek AS, Vangkroken 2, N-1351 Rud, Norway) was used to measure three-dimensional instantaneous velocities (*x*—longitudinal direction, *y* transverse direction, and *z*—vertical direction), i.e., longitudinal velocity (*u*), transverse velocity (*v*) and vertical velocity (*w*). In the experiment, the *x* and *y* direction was relative to the curved channel. The ADV sampling volume is 0.09 cm3 and the actual measurement point is 5 cm below the tip of the probe. For each cross-section, the ADV beam with red marking (*x*-axis indicator) was pointed to the streamwise direction of the curved channel. In order to ensure the quality of the measured data, SNR (signal to noise ratio) should be guaranteed to be above 20 and Correlation coefficient should be above 70 for the measured data [21]. Additionally, the measured velocity signals were despiked using the method of Goring and Nikora [22]. In the experiments, the sampling time was set to 60 s to ensure the accuracy and reliability of the data [6,23,24] with the sampling frequency of 25 Hz. The 60 s records were compared to a 10 min record both for the bare case and the rough bed case, and the differences for the time-averaged velocities and the *TKE* values were <0.9% and <1.1%, respectively. Thus, the 60 s ADV records were sufficient to obtain stable first (time-averaged) and second (variance) moments of turbulent statistics. Velocity measurements were conducted at five selected cross-sections (0◦, 45◦, 90◦, 135◦ and 180◦) with 5 or 7 vertical locations for each cross-section (5 cm, 10 cm, 15 cm, 20 cm, 25 cm, 30 cm and 35 cm from the inner wall) and 15 measuring points at each vertical line (measurements carried out at every 1 cm at 0 cm to 10 cm from the bed, and every 2 cm at 10 cm to 20 cm from the bed). The surface of the rough bed was uneven, so the first measuring point (1 cm from the bed) at the vertical line was confirmed by the average distance of three sections (*y/B* = 0.25, *y/B* = 0.5 and *y/B* = 0.75). The flow depth was measured using water level gauges at three vertical locations (*y/B* = 0.25, *y/B* = 0.5 and *y/B* = 0.75) for each cross-section along the flume. For the bare case (smooth bed), five vertical locations were enough to obtain the hydrodynamics characteristics of flow. For rough bed cases, velocities at seven vertical locations were measured to provide the detailed flow fields. The scheme of experiment setup is shown in Figure 1.

**Figure 1.** Scheme of experiment setup: (**a**) The shadow shows rough bed regions and the inset figure represents rough bed configurations; (**b**) ADV measurement setup.

The rough bed in the experiments was made by placing a thick layer of sediments throughout the curved region (Figure 1a). The thickness of sediment layers was around 2~3*D*50, where *D*<sup>50</sup> is the grain size diameter *D* at 50% in the cumulative distribution of particle size [25]. Three different sizes of bed particles were chosen: *D*<sup>50</sup> = 1.5 cm (refer to the small case hereafter), 3.0 cm (refer to the middle case hereafter) and 5.0 cm (refer to the big case hereafter) in the study. Sediment particles were placed in the regions of 1 m near the inlet and outlet of the curved channel to prevent flow instability as shown in Figure 1.

## *2.2. Fundamental Flow Conditions*

For determining the flow conditions in the experiments, Reynolds numbers *Re* and Froude numbers *Fr* are respectively calculated by

$$Re = 4\mathcal{U}\_{m0}\mathcal{R}\_h/v\tag{1}$$

and

$$Fr = \mathcal{U}\_{m0} / \sqrt{gh\_0} \tag{2}$$

where *Um*<sup>0</sup> is mean flow velocity of the 0◦ section for each case; *Rh*(=*Bh/*(*B+2h*0)) is the hydraulic radius; *B*(= 0.4 m) is the channel width; *h0* is the water depth of the 0◦ section; *v* = 10−<sup>6</sup> m<sup>2</sup> s−<sup>1</sup> is the kinematic viscosity of water; and *g* = 9.81 m<sup>2</sup> s−<sup>1</sup> is the gravitational acceleration. As listed in Table 1, the Reynolds numbers ranged from 64,516 to 148,148, indicating the flow conditions were turbulent. Additionally, the Froude numbers for all experimental conditions were less than 1, meaning that the flows belong to the subcritical flow regime.

For rough bed cases, viscous sublayer is absent [20], and thus the hydrodynamic characteristics of flow is solely affected by the bed roughness. Yalin [26], and Schlichting and Gersten [27] classified the rough bed by the means of friction velocity *u*<sup>∗</sup> (m s−1) as depicted below:

$$0 < \frac{u^\*k\_s}{v} < 5\tag{3}$$

$$5 < \frac{u^\* k\_s}{v} < 70\tag{4}$$

$$70 < \frac{u^\*k\_s}{v} \tag{5}$$

where *u*<sup>∗</sup> = *gRhS* is the friction velocity of 0◦ section for each case, *S* = 0.005 denotes the bed slope and *ks* refers to the roughness height in meters. Equation (3) refers to the condition of flow over a hydraulically smooth bed, while Equations (4) and (5) represent flow over transition region and hydraulically rough surfaces, respectively [20]. In these cases, *D*<sup>50</sup> is used to define the grain size. The roughness height for the Perspex sheet is determined as 0.0001 m [20]. For non-uniform bottom sediment, *ks* is expressed as a characteristic grain diameter, i.e., *ks* = *D*<sup>50</sup> [28]. The detailed experimental parameters are listed in Table 1. As shown in Table 1, 0 < (*u*∗*ks*)/*v* < 5 for the bare case, representing a hydraulically smooth bed. For the rough bed cases, (*u*∗*ks*)/*v* 70 represents hydraulically rough surfaces.


**Table 1.** Parameters of the experiment.

## *2.3. Numerical Analysis Using FLOW-3D*

FLOW-3D is a commercial CFD software that specializes in solving transient, freesurface problems. The finite volume method (FVM) in a Cartesian, staggered grid is employed in FLOW-3D to solve the Reynolds' average Navier–Stokes (RANS) equations. FLOW-3D contains a powerful meshing capability through the Fractional Area/Volume Obstacle Representation (FAVOR), which is used to illustrate the complex boundaries of the computational domain. In FLOW-3D, the Volume of Fluid (VOF) is used to simulate the fluid with the free surface [29–32]. In this study, the water depth, longitudinal velocity and secondary flow with and without the thick layer of sediments were numerically studied using FLOW-3D. The numerical domains were based on the experimental channel mentioned above. In order to make the numerical results more accurate, the model geometry was lengthened 20 m in the straight inflow and outflow reach, respectively. The total grid number reached up to 1.17 million, and the grid size was 2.5 cm in the curved reach and 5 cm in the straight inflow and outflow reach. The irrelevance of the number of grids was verified by means of reducing the grid size, i.e., increasing the grid number. The grid size was increased to 2 cm, and the difference of velocity profiles in comparisons with the 2.5 cm grid size was less than 3%.

The boundary conditions for the inflow and outflow reach were set as mean flow velocity and water depth measured in flume experiments. For free water surface, the atmospheric pressure was assigned. Two inter-block junctions of straight and curved reach were defined as the symmetrical condition [33]. In addition, a no-slip condition was applied to wall boundaries. The simulation ran for 1000 s to ensure the reach of steady state conditions. In FLOW-3D, the Renormalization-Group (RNG) turbulence model known to describe low intensity turbulence flows and flows having strong shear regions more accurately [34–36], was selected. The RNG model systematically removes all small scales of motion from the governing equations by considering their effects in terms of larger scale motion and a modified viscosity [35].

Three cases numerically simulated in the study were: (i) the bare case (the same as experiment); (ii) the big case (the same as experiment): the rough bed region was set as a 10 cm high and solid step (=2*D*50), simplified conditions, to represent the thick sediment

layer, which is the same as the experiment setup. Additionally, a 20 cm long slope before and after the step was used to connect the step and flat regions of the flume. The solid step should be set to porous media in accordance with the experimental settings. However, for the FLOW-3D software it is unable to set the solid step as porous media as well as assign a surface roughness on the surface of porous media. Instead, we calibrated the surface roughness *ks* and ensured the simulations close to flume measurements; (iii) the big case\_flat, in which no thick sediment layer was covered, i.e., no step was set up in numerical tests. The entire flume bottom is uniform, but bed roughness was the same as the big case. Cases (ii) and (iii) were used to compare the effect of thick sediment layers on secondary flows along a curved channel. The roughness height *ks* in the FLOW-3D model was a tuning factor for the three cases.

#### **3. Results and Discussion**

## *3.1. Water Depth*

The transverse gradient at the channel curve caused by the centrifugal force is small and negligible due to the small quantity of flow and small radius of curvature. Therefore, the average depth of three sections (*y/B* = 0.25, *y/B* = 0.5 and *y/B* = 0.75) is calculated and shown in Figure 2 for four cases.

**Figure 2.** Water depth variations along the curved reach for four cases: (**a**) *<sup>Q</sup>* = 0.015 m<sup>3</sup> <sup>s</sup><sup>−</sup>1, (**b**) *<sup>Q</sup>* = 0.025 m<sup>3</sup> <sup>s</sup><sup>−</sup>1, and (**c**) *<sup>Q</sup>* = 0.030 m3 <sup>s</sup>−<sup>1</sup> (• Bare case; - Small case; Middle case; Big case).

It can be seen from Figure 2 that the water depth changes less obviously along the curved reach for the bare case (smooth bed). As bed particles become larger, the sediment layer over the curved channel is thicker. At the conditions of the constant downstream water level at the tailgate, the thicker sediment layers lead to shallow water depth over the curved reach, i.e., *hbig case* > *hmiddle case* > *hsmall case* in general. Therefore, the bed with the largest particle sizes leads to fastest mean flow velocity when the water depth at the tailgate is constant. In addition, the smaller discharge as shown in Figure 2a leads to less water depth differences for rough bed cases. As discharge becomes larger, water depth differences for rough bed cases are more apparent as provided in Figure 2c.

## *3.2. Longitudinal Velocity Distribution*

Figure 3 shows the longitudinal velocities along the depth of each section for smooth beds (bare case) and rough beds (small case, middle case and big case) for *Q* = 0.030 m<sup>3</sup> s<sup>−</sup>1. For smooth bed (bare case, Figure 3a), the longitudinal velocity along the water depth at 0◦ and 45◦ sections follows logarithmic distribution, similar to the case in straight open channels. At the cross sections of 90◦ and 135◦, the flow is affected by the channel bend. The longitudinal velocity gradually increases in the bottom region (*z*/*h* < 0.4) and reduces in the top region (*z*/*h* > 0.4), thus breaking the logarithmic distribution law. Consequently, the longitudinal velocity as shown in Figure 3a exhibits an approximate "constant" distribution along the water depth as Blanckaert [37] and Barbhuiya and Talukdar [38] found.

For rough bed cases (Figure 3b–d), the longitudinal velocity profile is significantly different as the bare case. For the bare case, longitudinal velocity ranges from 0.2 m s−<sup>1</sup> to 0.35 m s−1, while the values vary from 0.18 m s−<sup>1</sup> to 0.45 m s−<sup>1</sup> for the rough bed cases.

Since the water depth at the tailgate is constant for all cases, the water depth with larger bed sediment particles becomes shallower, leading to greater longitudinal velocity. At 0◦ sections, the longitudinal velocity along the water depth also follows distinct logarithmic distribution, but the longitudinal velocity near the bed greatly decreases due to the bottom friction from rough bed, causing the steeper slope of the longitudinal velocity profile and larger induced bottom shear stress as other studies [15,23,39,40] mentioned. For the middle and big cases, the slope of the longitudinal velocity profile has little difference, which possibly means that the longitudinal velocity profile changes a little when the grain size reaches a certain level. At the 45◦ section, the longitudinal velocity profile near convex bank is affected by the channel bend firstly, implying that the longitudinal velocity gradually increases in the bottom and reduces near the water surface. However, the velocity profile near concave bank keeps following the logarithmic distribution. At the cross sections of 90◦ and 135◦, the flow in the whole section is affected by the channel bend. Under the both effects of channel bend and rough bed, the longitudinal velocity distribution exhibits a trend of increasing near the bed and then decreasing near the surface along the water depth. It also can be found that the maximum longitudinal velocity at the central line occurs between 0.2 to 0.4 dimensionless depth from the bed. The grain size of bed particle only changes the magnitude of velocity, making the phenomena more remarkable. The trends of longitudinal velocity profiles for different rough bed cases are similar.

**Figure 3.** *Cont*.

**Figure 3.** Longitudinal velocity profiles for four cases at five cross sections: (**a**) Bare case, (**b**) Small case, (**c**) Middle case, and (**d**) Big case (*Q* = 0.030 m<sup>3</sup> s<sup>−</sup>1).

The velocity gradient between convex bank and concave bank for the bare case is unobvious. On the other hand, for the rough bed cases the transverse velocity gradient cannot be neglected, especially for the big case. The longitudinal velocity of concave bank is greater than that of convex bank for all rough bed cases, which indicates that the rough bed can enlarge the centrifugal effect of channel bend. Longitudinal velocity profiles for different flow discharge are similar and not discussed in detail here.

### *3.3. Turbulent Kinetic Energy and Power Spectral Density Analysis*

In the study, *TKE (k*), represented as the turbulence characteristics, is defined as

$$k = \frac{1}{2} \left( \overline{u'^2} + \overline{v'^2} + \overline{w'^2} \right) \tag{6}$$

where *u* = *u* − *u*, *v* = *v* − *v* and *w* = *w* − *w* are the fluctuation velocities of *u*, *v* and *w*, respectively.

Taking *Q* = 0.030 m<sup>3</sup> s−<sup>1</sup> as an example, the *TKE* profiles, normalized by the overall mean velocity *U* = *Q*/*Bh*, under different tough bed conditions are shown in Figure 4. For the bare case (Figure 4a), the normalized *k* values show a constant distribution along the water depth for each cross section, but the normalized *k* values gradually increase from

inlet (0◦ section) to the apex of bend (90◦ section) and decrease from the apex to outlet (180◦ section). The intense turbulence mainly concentrates in the cross section of 90◦–135◦. The reason is that the friction of the wall and the momentum exchange between the water flow and the wall leads to the increase of *k* values [4].

For the rough bed cases (Figure 4b–d), the normalized *k* values are smaller near the free surface, similar to the bare case. However, the normalized *k* values are larger near the bed region at the 0◦ section, implying that rough bed will enhance fluid mixing and strengthen turbulence as the results found in straight open channels [15,41]. The effect of bed roughness is not only confined for near bed regions but can be gradually extended to the whole water depth. As the particle size increases, the normalized *k* values near the bed increase a little because of the greater friction from the rough bed. Yet, the normalized *k* values near the surface remain basically unchanged. The slope of the variation of the normalized *k* values, i.e., <sup>Δ</sup>(*z*/*h*) <sup>Δ</sup>(normalized *<sup>k</sup>* values), decreases, which means that at the same difference Δ *z h* , Δ(normalized *k* values) becomes larger. As a result, the rate in loss of *k* for larger particles becomes faster. After the 45◦ section, the normalized *k* values are also affected by the channel bend. Similar to the bare case, the normalized *k* values of the rough bed on the convex bank are larger than those for concave bank. Additionally, the vertical gradient of the variation of normalized *k* value increases, indicating that the rate in loss of *k* becomes slow. The *k* distributions for three rough bed cases are similar in curved channels. Therefore, the bed roughness and curved channel both improve turbulence along the water depth, especially near the bed regions. The effect of roughness is essential for the transport of the turbulent kinetic energy along the vertical direction.

Figure 5 presents the variations in the power spectral density (PSD) of longitudinal velocity fluctuation *u* with frequency *f* for different experimental conditions measured at the mid-width of the flume (20 cm from inner bank). Since the flow discharge has little effect on PSD, the discharge *Q* equal to 0.030 m<sup>3</sup> s−<sup>1</sup> is used as an example to reveal the variations of PSD. The spectral analysis used here is based on the Welch method with Hamming type windowing [42]. In Figure 5a,c,e,g, the feature point was measured at mid-depth (12 cm from the bed). In order to explore the effect of distance from the channel bed on PSD, three characteristic positions, i.e., 2 cm (bottom), 12 cm (middle) and 16 cm (top) from bed along the vertical direction, are provided at the 90◦ section in Figure 5b,d,f,h. Due to the limitation of ADV frequency, the energy input, inertial subrange and part of dissipation range can be observed. However, the higher frequency (part of dissipation range) is unable to present. Therefore, some of the small-scale turbulent structure can still be investigated.

**Figure 4.** *Cont*.

**Figure 4.** TKE (*k*) profiles for four cases at five cross sections: (**a**) Bare case, (**b**) Small case, (**c**) Middle case, and (**d**) Big case (*Q* = 0.030 m3 s<sup>−</sup>1).

**Figure 5.** Power spectral density for different cases: (**a**) bare case, (**b**) different depth (2 cm, 12 cm, 16 cm) from bed for the bare case of 90◦ section, (**c**) small case, (**d**) different depth (2 cm, 12 cm, 16 cm) from bed for the small case of 90◦ sectioI (**e**) middle case, (**f**) different depth (2 cm, 12 cm, 16 cm) from bed for the middle case of 90◦ section, (**g**) big case, (**h**) different depth (2 cm, 12 cm, 16 cm) from bed for the big case of 90◦ section.

For the bare case (Figure 5a), the PSD value for the high-frequency part at the 45◦ section is the greatest, while the PSD value at the 180◦ section is the smallest, indicating that the turbulence becomes strong, and the high-frequency (small-scale) vortex increases due to the influence of channel bend. For Figure 5b, there are similar trends for different positions for the bare case, which means the structures of turbulent vortex along vertical directions are approximately identical over smooth bed.

For the rough bed cases (Figure 5c,e,g), the PSD value for the high-frequency part at the 90◦ section is obviously larger than other sections. For Figure 5d,f,h, the PSD values at near-bed positions are larger than those at middle and top positions, and the trend is more prominent with the increasing bed particle size, implying that the turbulent mixing near the bed regions are more intense. The PSD of the longitudinal velocity fluctuations satisfy the classic Kolmogorov −5/3 law in the inertial subrange. Comparing to the bare case, the existence of bed roughness shortens the inertial subrange, causing the energy access to the viscous dissipation in advance. The peak frequency for each rough bed case is similar and about 0.8 Hz, irrespective of bed particle sizes.

## *3.4. Turbulent Bursting*

It is found that fluid motion near a bed is not completely chaotic in nature, but it is a clear "sequence of ordered motion" [43]. Such coherent motions are called the bursting process [43]. To quantify the intermittent instantaneous Reynolds stresses as well as identify turbulence structures within a turbulent bursting sequence, one of the widely used conditional sampling techniques is the quadrant analysis of the Reynolds shear stress [44–46]. Here, quadrant analysis was conducted to study the effect of rough bed on turbulent bursting. In the quadrant analysis, the Reynolds stress has four types of contributions according to the signs of the instantaneous velocity fluctuations [42]. Accordingly, turbulent events are defined by the four quadrants as outward interactions (*i* = *I*, *u* > 0, *w* > 0), ejections (*i* = *I I*, *u* < 0, *w* > 0), inward interactions (*i* = *III*, *u* < 0, *w* < 0), and sweeps (*i* = *IV*, *u* > 0, *w* < 0). Results in quadrants *II* and *IV* mean the positive downward momentum flux and are involved in turbulence near-bed bursting [43].

In order to describe the turbulent event accurately, the hole (instead of zero) concept is used to eliminate smaller Reynolds stresses. The hole is formed by four hyperbolas |*u* (*t*)*w* (*t*)| = *G*0|*u w* |, where *G*<sup>0</sup> is a threshold value. By using the threshold, small values can be ignored in the *i*-th quadrant [47]. The contribution of each quadrant can be represented as *Sk* (*k* = *I*, *II*, *III*, *IV*, indicate the four quadrants), where

$$S\_k = \begin{cases} -1, & |u'(t)w'(t)| > \mathcal{G}\_0|u'w'|, \ [u'(t), w'(t)] \text{ in the same quadrant} \\ & 0, \text{ otherwise} \end{cases} \tag{7}$$

In the study, the threshold *G*<sup>0</sup> was set to 1.0 [48]. Then, the revised occurrence frequency, *fk*, for the four turbulent events is given as:

$$f\_k = \frac{\sum\_{t=0}^{T} \mathcal{S}\_k}{\sum\_{t=0}^{T} \mathcal{S}\_I + \sum\_{t=0}^{T} \mathcal{S}\_{II} + \sum\_{t=0}^{T} \mathcal{S}\_{III} + \sum\_{t=0}^{T} \mathcal{S}\_{IV}} \tag{8}$$

where *T* is the length of measurement time.

Figure 6 shows characteristic distributions of the occurrence frequency *fk* in different cases for *Q* = 0.030 m3 s−<sup>1</sup> as an example. The feature point is measured at the distance of 0.01 m from the bed in the vertical direction and 0.20 m from the inner bank (center of the flume) in the transverse direction in each section to investigate the turbulent event near bed. For the bare case (Figure 6a), it can be observed that sweep and ejection events (event *II* and *IV*) occur with comparable frequency and are much higher than outward and inward interactions at the 0◦ section, indicating that upwards motion of low-speed fluid and downwards motion of high-speed fluid are dominant. After the 45◦ section, the frequency of occurrence of events of *I* and *III* quadrants become higher. At the 180◦

section, the frequency of occurrence for four events is similar. This reveals that high-speed fluid reflects by the bottom and low-speed fluid is pushed back as the influence of channel bend, showing a strong interaction of the turbulent structures with the main flow. For rough bed cases, it is noted that higher magnitude of sweep and ejections (event *II* and *IV*) contributions as shown in Figure 6b (small case) similar to the bare case at the 0◦ section. When flow enters the curved channel, the frequency of occurrence for event *II* and *IV* increases, indicating that sweep and ejection are the most dominant processes, which can potentially influence the sediment transport in the stream, causing the exchange of energy and momentum in the flow and the bed formation. The middle case (Figure 6c) performs the similar trend as the small case. However, for the big case (Figure 6d) the frequency of occurrence for event *II* and *IV* decreases first (0◦ section to 90◦ section) and then increases (90◦ section to 180◦ section). Therefore, the contributions of sweeps and ejections are more important than those of the outward and inward interactions over a rough bed, which are the same results obtained by other researches in straight channels [43,49]. Bed roughness is the leading factor to the turbulent bursting in curved channels when the grain size of bed particles is moderate (the small case and middle case), but for the big case the turbulent bursting is weaken before the apex of the bend. The possible explanation is the larger grain size of bed particles possible disturbs the bottom boundary layer, changing turbulent structure near the bed.

**Figure 6.** Frequency of occurrence of coherence turbulent events for different cases (*Q* = 0.030 m3 s<sup>−</sup>1).

#### *3.5. Comparisons with Previous Experimental Results*

All of the experiment cases in this paper were done in a 180 degree U-shaped curved and rough bed flume under the condition of constant downstream water level, close to the practical situations on the purpose of navigation. Thus, the experimental results in comparison with previous studies may have some differences. Firstly, since the downstream water level is fixed, as the bed particles become larger, the water depth along the curved channel decreases, leading to greater velocity compared with that over smooth curved channels. Without the constant downstream water depth, the water depth along the curved channel may keep similar [18]. Secondly, the longitudinal velocity distribution is relatively different. Pradhan et al. [20] found that the velocity remains higher towards the inner wall for rough bed. In our experimental results, longitudinal velocity near the concave bank is greater than that near the convex bank for all rough bed cases, which indicates that the rough bed can enlarge the centrifugal effect of channel bend. The reason may be the different shape of cross-section of channel (trapezoidal channel for Pradhan et al. [20]). Thirdly, sweeps and ejections are more important than those of the outward and inward

interactions over a rough bed when the grain size of bed particles is intermediate [43,49]. For the big case, the proportion of the outward and inward interactions is higher than that of the sweeps and ejections in 90◦ section, which was not shown in the results of Najafabadi et al. [43] because of the limits of the grain size in the study.

## *3.6. Numerical Results*

Numerical tests aim to understand the effect of thick sediment layers on secondary flow over a curved channel. The roughness height *ks* of the FLOW-3D simulation was adjusted until the numerical results close to experimental data. The roughness height *ks* = 0.005 m for the bare case and the roughness height *ks* = 0.05 m for the big case can match well with the experimental water depths and longitudinal velocity as shown in Figures 7 and 8. The differences between the simulation and measurements can be possibly attributed to three aspects: (1) the isotropic assumptions in the RNG turbulence model are not suitable for curved channel flows, (2) the uneven bottom surface from sediment particles and flows through pores between sediment particles are unable to reflect in the numerical tests, where the solid step is used to represent the sediment layers, and (3) the slope of the entire channel may not be constant. For big case\_flat, the water depth increases a little due to the bottom friction from rough bed, and the water depth of each section is approximately the same, i.e., uniform flow conditions, meaning that the friction exerted by the rough bed mainly balances gravity due to sloping bed.

**Figure 7.** Comparisons of simulated and measured water depth for *Q* =0.030 m<sup>3</sup> s<sup>−</sup>1.

Figure 8 performs the longitudinal velocity profiles at the mid-width of the flume (20 cm from inner bank) under different cases. The simulated velocity distribution at the 90◦ and 135◦ section does not fit very well with the experimental results, especially for the big case. For big case\_flat, the longitudinal velocity near bed decreases due to the bed friction when compared with the bare case. On the other hand, the longitudinal velocity for the big case\_flat is smaller than that for the big case owing to larger water depths.

The distributions of longitudinal velocity in three typical cross-sections (0◦, 90◦ and 180◦) are shown in Figure 9. For the bare case (Figure 9a), the core region of maximum longitudinal velocity has obviously shifted toward the concave bank at the 90◦ and 180◦ section in comparison with that at the 0◦ section. This is because of the advective transport of streamwise momentum by the cross-stream circulation in open-channel bends. At the bend entrance, the shorter distance in the inner bend than that in the outer bend would lead to longitudinal velocity *u* decreasing from the convex toward the concave bank. The flow acceleration/deceleration is induced by streamwise pressure gradients related to the sudden transverse tilting of the water surface in the bend. In addition, the sudden disappearance of the transverse tilting of the water surface leads to pronounced flow accelerations/deceleration in the concave/convex bend at the bend exit [3,50]. For the big case (Figure 9b), the core region of maximum longitudinal velocity stays around the center region of the cross-section at the 0◦ and 90◦ sections. At the 180◦ section, the core region shifts to the concave bank. For the big case\_flat (Figure 9c), numerical results are similar to that for the bare case at the 0◦ and 90◦ sections. However, at the 180◦ section, the maximum longitudinal velocity for the big case\_flat concentrates on the right middle region rather than the lower right region for the bare case. Based upon the results, it can be concluded that thick sediment layers delay the shifting of the core of maximum longitudinal velocity towards the concave bank. The bed roughness only reduces the longitudinal velocity magnitude.

**Figure 8.** Comparisons of simulated and measured longitudinal velocity profiles on 5 cross sections: (**a**) 0◦ section, (**b**) 45◦ section, (**c**) 90◦ section, (**d**) 135◦ sectionInd (**e**) 180◦ section.

**Figure 9.** The distributions of the longitudinal velocity at the 0◦, 90◦ and 180◦ section: (**a**) Bare case, (**b**) Big case, and (**c**) Big case\_flat.

The distributions of *k* at the mid-width of the flume (20 cm from inner bank) at three typical cross-sections (0◦, 90◦ and 180◦) are shown in Figure 10. For the 0◦ section (Figure 10a), the *k* values of three cases are smaller near the free surface but become larger near the rough bed regions, implying that rough bed will enhance fluid mixing and strengthen turbulence. Additionally, the *k* values of the big case\_flat are slightly larger than those of the big case. For the 90◦ and 180◦ sections (Figure 10b,c), the *k* values firstly increase (from 0◦ section to 90◦ section) and then decrease (from 90◦ section to 180◦ section), indicating that the channel bend can also improve turbulence along the water depth, especially near the bed regions. However, the *k* values of the big case (red line) increase more than the other two cases along the curved channel. The result suggests that the bed roughness and curved channel both improve turbulence along the water depth, especially near the bed regions, which is the same as the experimental results. Additionally, the thick sediment layer can promote *TKE* more than the curved channel does.

**Figure 10.** The profiles of *TKE* (*k*) at the 0◦, 90◦ and 180◦ sections.

In order to explore the effect of distance from the convex bank on the distributions of *k*, three characteristic positions, i.e., <sup>1</sup> <sup>4</sup> cm (=1/4 *<sup>B</sup>*), <sup>1</sup> <sup>2</sup> cm (=1/2 *<sup>B</sup>*) and <sup>3</sup> <sup>4</sup> cm (=3/4 *B*) from the convex bank along the transverse direction, are provided at the 90◦ section for three cases as shown in Figure 11. For the bare case (Figure 11a), the *k* values of three lines basically collapse together, meaning that the transverse positive has little effect to the *k* values for the smooth bed. For the big case (Figure 11b), the *k* values near convex bank is slightly larger when *z/h* = 0.1–0.6. The thick sediment layer with bed roughness strengthens the turbulence more near convex bank, which is the same as the experimental result. However, for the big case\_flat (Figure 11c), the *k* values near concave bank at *z/h* > 0.1 is larger. The thick sediment layer increases the *TKE* (*k*) values, especially near convex bank.

**Figure 11.** The distributions of *TKE* (*k*) at the 90◦ section for three cases.

The secondary flow that can interact dynamically with the primary flow is important for meandering rivers and useful for studies of diffusion and navigation in natural waterways, because it is partly responsible for the large-scale bed topography of natural alluvial channel bends [51]. The vector of the cross-stream motion (*v*, *w*) at the 90◦ section is shown in Figure 12. For the bare case (Figure 12a), there is a classical helical motion obviously, called the center-region cell. In addition, a weaker counterrotating cell of cross-stream circulation, called the outer-bank cell, exists in the corner formed by the water surface and the outer bank as Blanckaert and Graf [2] revealed. The center-region cell plays an important role in the advective momentum transport, and the outer-bank cell is believed to be crucial to bank erosion processes [52]. For the big case (Figure 12b), the center-region cell and outer-bank cell still exist. The center-region cell is confined to the near-bed area, and

the outer-bank cell is stronger than that for the bare case, possibly accelerating the bank erosion. This is because the rough bed disturbs the flow, mainly reducing the longitudinal flow velocity near bed, and also affecting the distribution of the transverse velocity along the water depth. For the big case\_flat (Figure 12c), the center-region cell is confined to a smaller area compared to the bare and big cases, which means that the thick sediment layer can influence more riverbed regions. Furthermore, the outer-bank cell is weaker than that for the big case.

**Figure 12.** Numerical results of velocity vector for cross-stream motions: (**a**) the bare case, (**b**) the big case, and (**c**) the big case\_flat. Note: the arrowed circle represents the secondary flow.

The transverse velocity distribution can be used to display the intensity of secondary flow. Figure 13 shows the secondary flow distributions of different cases at the mid-width of the flume (20 cm from inner bank) at the 90◦ section, where X-axis is *vr*/(*uh*), Y-axis is *z*/*h*, *r* is curvature radius of bend, and *u* is average longitudinal velocity. When the calculated values (*vr*/(*uh*)) distribute on the same side of *vr*/(*uh*) = 0, i.e., all *vr*/(*uh*) values are larger or smaller than zero, it means that the direction of transverse velocity is the same at this position. In other words, there is no secondary flow. Once the calculated values (*vr*/(*uh*)) fall on the two sides of *vr*/(*uh*) = 0, there is secondary flow existing. Additionally, the secondary flow is stronger as the slope of transverse velocity distribution increases.

**Figure 13.** The transverse velocity distribution for experimental results (symbols) and numerical simulations (solid lines).

From Figure 13, it can be seen that numerical simulations (solid line) are close to experimental measurements (individual symbols). The slope of transverse velocity distribution for the big case is steeper than that for the bare case, which means that the rough bed would strengthen the secondary flow and possibly causes more riverbed erosion. Besides, the position of calculated values (*vr*/(*uh*) = 0) for the big case is lower than that of the bare case, representing the secondary flow is confined to the near-bed area as the results discussed above. For the big case and big case\_flat, the slope of transverse velocity distribution is comparable, i.e., similar secondary flow intensity.

#### **4. Conclusions**

Experiments and numerical simulation of curved open channel flows over rough beds under the condition of constant downstream water depth were carried out. Based on the ADV data, hydrodynamic characteristics of flow over rough bed such as longitudinal velocity distribution, turbulent characteristics, secondary flow, turbulent bursting and so on was analyzed. Additionally, the numerical results were used to show different flow patterns for the rough bed with and without thick sediment layers. The following conclusions are drawn:


classic Kolmogorov -5/3 law in the inertial subrange, and the existence of rough bed shortens the inertial subrange and reaches to the viscous dissipation in advance.


The findings of the research could improve the understanding of the interactive effects between the rough bed and the strongly curved channel flow. Based on the reported results, future researches can focus on studying the effect of flume curvature on the hydrodynamic characteristics of curved channel flow with gravel beds.

**Author Contributions:** Conceptualization, Y.-J.C., Y.-T.L. and X.J.; Methodology, Y.-J.C. and Y.-T.L.; Laboratory measurement, Y.Y.; Data analysis, Y.-T.L. and Y.Y.; Writing—Original Draft Preparation, Y.-T.L. and Y.Y.; Writing—Review and Editing, Y.-J.C.; Software, X.J.; Funding Acquisition, Y.-J.C. and X.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Key Research and Development Project of China [grant number 2017YFC0405205]; Ministry of Science and Technology of Taiwan [grant number MOST 109-2625-M-002-016]; Zhejiang Provincial Natural Science Foundation of China [grant number LY20A020009); State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University [grant number HESS-2014]; Fundamental Research Funds for the Central Universities [grant number 2020QNA4038].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data generated or analyzed during this study are included in this article.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Mouldi Ben Meftah \* , Diana De Padova , Francesca De Serio and Michele Mossa**

Department of Civil, Environmental, Land, Building Engineering and Chemistry, Polytechnic University of Bari, Via E. Orabona 4, 70125 Bari, Italy; diana.depadova@poliba.it (D.D.P.); francesca.deserio@poliba.it (F.D.S.); michele.mossa@poliba.it (M.M.)

**\*** Correspondence: mouldi.benmeftah@poliba.it; Tel.: +39-080-5963508

**Abstract:** Most studies on local scouring at grade control structures have principally focused on the analysis of the primary flow field, predicting the equilibrium scour depth. Despite the numerous studies on scouring processes, secondary currents were not often considered. Based on comprehensive measurements of flow velocities in clear water scours downstream of a grade control structure in a channel with non-cohesive sediments, in this study, we attempted to investigate the generation and turbulence properties of secondary currents across a scour hole at equilibrium condition. The flow velocity distributions through the cross-sectional planes at the downstream location of the maximum equilibrium scour depth clearly show the development of secondary current cells. The secondary currents form a sort of helical-like motion, occurring in both halves of the cross-section in an axisymmetric fashion. A detailed analysis of the turbulence intensities and Reynolds shear stresses was carried out and compared with previous studies. The results highlight considerable spatial heterogeneities of flow turbulence. The anisotropy term of normal stresses dominates the secondary shear stress, giving the impression of its crucial role in generating secondary flow motion across the scour hole. The anisotropy term shows maximum values near both the scour mouth and the scour bed, caused, respectively, by the grade control structure and the sediment ridge formation, which play fundamental roles in maintaining and enhancing the secondary flow motion.

**Keywords:** scour; equilibrium condition; velocity field; secondary currents; turbulence

## **1. Introduction**

The presence of natural or man-made structures on riverbeds plays an important role in the evolution of river morphology and sediment entrainment. Flow turbulence properties and secondary currents play a crucial role in sediment transport, and, in turn, suspended particle motion influences turbulence, such as Reynolds shear stress and velocity [1]. In addition to the fluid–particle interactions, which definitely influence the flow velocity distribution, the fluid–structure interactions, i.e., with natural vegetations, riverbed debris, bridge piers and abutments, sills, sluice gates, spillways, weirs, spur dikes, off-shore platforms, wind turbines, energy converters, etc., cause additional complex effects on the flow hydrodynamic characteristics. Local scouring is produced due to these complex flow patterns occurring in the surroundings of such structures. The local scouring process has attracted the attention and interest of many scientists for decades [2–18].

Experimental studies on the scouring process at grade control structures (GCSs) in riverbeds [3,5,9,17,18] showed that the scour often developed downstream of the structure. The extension of the scour hole is strongly influenced by the properties of the incoming flow, which is usually a two-dimensional jet-like flow. Owing to the high velocity of the entering jet flow, a large amount of sediment erosion locally occurs downstream of the GCS, forming the scour hole. Because of the large velocity gradients among the jet flow and that in the scour pool, the jet diffuses near the bed, and is redirected at a reduced bed velocity. The equilibrium state occurs when the path of the impinging jet becomes long enough and its diffused velocity is reduced to values below the minimum value required to move

**Citation:** Ben Meftah, M.; De Padova, D.; De Serio, F.; Mossa, M. Secondary Currents with Scour Hole at Grade Control Structures. *Water* **2021**, *13*, 319. https://doi.org/10.3390/ w13030319

Academic Editor: Roberto Gaudio Received: 25 December 2020 Accepted: 23 January 2021 Published: 28 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the sediments [17,18]. The entering jet produces a very complex flow field in the scour hole. For structures around which the flow passes (bypassed structures), i.e., bridge piers, abutments, offshore platforms, and wind turbines, the scour profile is considerably affected by the flow turbulence structures, generated in the form of a horseshoe vortex, wake vortex, and surface roller [14].

Most studies on local scouring downstream of GCSs have essentially focused on the prediction of the maximum scour dimensions at equilibrium, especially the maximum scour depth and the maximum length. For the sake of simplicity, in these studies, flows have been analyzed in the plane of flow symmetry (along the channel axis), neglecting the effect of the wall-normal velocity component. However, flows within scour holes, regardless of the GCS's geometry, are three-dimensional, and therefore, secondary currents can occur in some cross-sections of the scour hole. Secondary currents represent circulation of fluids around the axis of the primary flow [19].

In open channel flows, secondary currents affect the distribution of bed shear stress, Reynolds stress, and turbulence intensities across the channel [20]. In large channels, the secondary currents consist of a series of counter-rotating cells through the sectional planes. Between the cells, upwelling and downwelling alternating movement zones may occur, which extend over the entire water column. Secondary current cells often generate a sort of undulating bottom shear stress distribution in the transverse direction, affecting the whole depth-flow field and the free-surface flow pattern [20,21]. According to Albayrak and Lemmin [20], the dynamics of secondary currents are strongly affected by the channel aspect ratio. The authors [20] observed that the number of secondary current cells changes proportionally with the aspect ratio. Within the range of their experimental conditions, Albayrak and Lemmin [20] also argued that, for a given aspect ratio, the number of secondary current cells is not affected by changes in the Reynolds number or the Froude number. Papanicolaou et al. [22] proved that the presence of secondary currents increases the sidewall shear stress and affects the turbulent production within a channel flow.

Secondary flow circulation and scouring processes were also observed in stream confluences [23–25]. The increase of flow velocity due to tributary junction generates a zone of maximum scour located near the center of the confluence. This zone is characterized by dominant flow convergence and a consistent pattern of secondary circulation.

Despite the numerous studies conducted on local scouring, a lack of information regarding the structures and generation of a secondary current in scour holes has been noted in literature. The main factors generating secondary currents in straight non-circular channel flows have been and remains a topic of much discussion and conjecture. In this study, we essentially focus on the hydrodynamic characteristics of flow across scour holes developed downstream of a GCS in sand riverbeds. Thanks to several measurements of the flow velocities performed through the scour cross-section at the position of maximum equilibrium scour depth, the secondary current patterns in the scour hole were achieved. This paper aims to: (i) check the development of secondary currents in scour hole downstream of a GCS, (ii) analyze the evolution of the turbulence structure in the scour hole at the equilibrium condition, and (iii) try to understand the physical origin of secondary current cells across the equilibrium scour hole.

## **2. Theoretical Considerations**

Sand ridges, also termed sediment strips, which are usually aligned parallel to the direction of the mean flow, have been widely observed in nature, i.e., in rivers with a bed composed of non-cohesive material, continental shelves, estuaries, and deserts [26]. These phenomena are strongly related to the development of secondary currents. In the literature (e.g., [20,21,26]), steady secondary currents have been classified into two categories: (i) secondary currents of the first kind (skew-induced stream-wise vorticity), taking origin from mean flow, but driven by curvature effect, and (ii) secondary currents that are generated by anisotropy of flow turbulence.

For steady flow and incompressible fluid, the time-averaged continuity and Navier– Stokes equations are (using index notation):

$$\frac{\partial \mathcal{U}\_i}{\partial x\_i} = 0 \tag{1}$$

$$\frac{\partial \mathcal{U}\_{\dot{j}} \mathcal{U}\_{\dot{i}}}{\partial \mathbf{x}\_{\dot{j}}} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_{\dot{i}}} + \nu \frac{\partial^2 \mathcal{U}\_{\dot{i}}}{\partial \mathbf{x}\_{\dot{j}}^2} + \frac{\partial \tau\_{\dot{i}\dot{j}}}{\partial \mathbf{x}\_{\dot{j}}} + g\_{\dot{i}} \tag{2}$$

where *xi* and *xj* = (*x*, *y*, *z*) is the direction tensor, *Ui* and *Uj* = (*U*, *V*, *W*) is the time-averaged velocity tensor, in which *U*, *V*, and *W* are the velocity along *x*, *y*, and *z*, respectively, *ρ* is the water density, *p* is the pressure force, *ν* is the water kinematic viscosity, *gi* is the gravitational acceleration tensor, *τij* = −*U iU <sup>j</sup>* is the shear stress tensor, and *U iU <sup>j</sup>* is the time-averaged shear stress of *U iU <sup>j</sup>*(*t*) over the length of the time series. The instantaneous velocity is defined as *ui*(*t*) = *Ui + U <sup>i</sup>*(*t*), where *U <sup>i</sup>* = (*U* , *V* , *W* ) is the velocity fluctuation tensor of *u*, *v*, and *w*, respectively.

From Equation (2), one can obtain the equation of the transport of stream-wise momentum as:

$$\underbrace{-\frac{1}{\rho}\frac{\partial p}{\partial \mathbf{x}}}\_{\mathbf{I}} - \underbrace{\left(\mathcal{U}\frac{\partial \mathcal{U}}{\partial \mathbf{x}} + \mathcal{V}\frac{\partial \mathcal{U}}{\partial \mathbf{y}} + \mathcal{W}\frac{\partial \mathcal{U}}{\partial \mathbf{z}}\right)}\_{\mathbf{II}} + \underbrace{\nu\left(\frac{\partial^{2}\mathcal{U}}{\partial \mathbf{x}^{2}} + \frac{\partial^{2}\mathcal{U}}{\partial \mathbf{y}^{2}} + \frac{\partial^{2}\mathcal{U}}{\partial \mathbf{z}^{2}}\right)}\_{\mathbf{III}} - \underbrace{\left(\frac{\partial \mathcal{U}^{\prime}\mathcal{U}^{\prime}}{\partial \mathbf{x}} + \frac{\partial \mathcal{U}^{\prime}\mathcal{V}^{\prime}}{\partial \mathbf{y}} + \frac{\partial \mathcal{U}^{\prime}\mathcal{W}^{\prime}}{\partial \mathbf{z}}\right)}\_{\mathbf{IV}} = 0\tag{3}$$

In Equation (3), *I* is the pressure gradient driving the flow, *II* indicates the convective transport of stream-wise momentum, *III* represents the transport due to viscous shear stresses, and *IV* is the transport due to turbulent stresses. The term *III* is generally negligible, except close to boundaries (bed and banks/walls). In the case of uniformly developed turbulent flow, *∂*/*∂x* = 0 in the terms *II*, *III*, and *IV*.

Since, in this study, the channel is straight, and the anisotropy of the flow turbulence could be the source of the secondary current formation. In this case, the secondary currents (steady and incompressible flow) can be described by the stream-wise vorticity equation as:

$$\underbrace{\mathcal{U}\frac{\partial\omega\_{x}}{\partial x} + V\frac{\partial\omega\_{x}}{\partial y} + W\frac{\partial\omega\_{x}}{\partial z}}\_{A1} = \underbrace{\omega\_{x}\frac{\partial\mathcal{U}}{\partial x} + \omega\_{y}\frac{\partial\mathcal{U}}{\partial y} + \omega\_{z}\frac{\partial\mathcal{U}}{\partial z}}\_{A2} + \underbrace{\upsilon\left(\frac{\partial^{2}\omega\_{x}}{\partial x^{2}} + \frac{\partial^{2}\omega\_{x}}{\partial y^{2}} + \frac{\partial^{2}\omega\_{x}}{\partial z^{2}}\right)}\_{A4}$$

$$+\underbrace{\frac{\partial}{\partial x}\left(\frac{\partial(\mathcal{U}\mathcal{V}\mathcal{V})}{\partial z} - \frac{\partial(\mathcal{U}\mathcal{V}\mathcal{W})}{\partial y}\right)}\_{A4} + \underbrace{\frac{\partial^{2}}{\partial y\partial z}(\mathcal{V}^{2} - \mathcal{W}^{2})}\_{A5} + \underbrace{\left(\frac{\partial^{2}}{\partial z^{2}} - \frac{\partial^{2}}{\partial y^{2}}\right)V\mathcal{V}\mathcal{W}}\_{A6}\tag{4}$$

where *ω<sup>x</sup>* = *<sup>∂</sup><sup>W</sup> <sup>∂</sup><sup>y</sup>* <sup>−</sup> *<sup>∂</sup><sup>V</sup> <sup>∂</sup><sup>z</sup>* , *<sup>ω</sup><sup>y</sup>* <sup>=</sup> *<sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>z</sup>* <sup>−</sup> *<sup>∂</sup><sup>W</sup> <sup>∂</sup><sup>x</sup>* , and *<sup>ω</sup><sup>z</sup>* <sup>=</sup> *<sup>∂</sup><sup>V</sup> <sup>∂</sup><sup>x</sup>* <sup>−</sup> *<sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>y</sup>* are the stream-wise, the spanwise, and the vertical vorticities, respectively. The term A1 in Equation (4) represents the rate of change of *ω<sup>x</sup>* due to the convection of fluid by the mean flow velocities: primary flow *U* and secondary motions *V* and *W*. The term A2 represents the tilting and stretching of the three vorticities by the gradients of the mean flow velocities. The term A3 represents the viscous diffusion of *ωx*. The terms A4, A5, and A6 represent the rates of change of the normal and shear stresses in the cross-section plane, which are theoretically responsible for developing and maintaining the second kind of secondary currents [27,28].

## **3. Experimental Set-Up**

The experiments on the scour processes were carried out in a rectangular flume of closed-circuit flow at the Hydraulic Laboratory of the Mediterranean Agronomic Institute of Bari (Bari, Italy). The flume had glass sidewalls and a Plexiglas floor, allowing a good side view of the flow. It was 7.72 m long, 0.30 m wide, and 0.40 m deep. A pump with a maximum discharge of 24 L/s was used to deliver water from the laboratory sump to an upstream tank equipped with a baffle and lateral weir, maintaining a constant head upstream of a movable slide gate constructed at the inlet of the flume. The slide gate regulates channel flow discharge. To create a smooth flow transition from the upstream reservoir to the flume, a rectangular wooden ramp, playing the role of a first GCS, was installed at the inlet of the flume; it was 1.55 m long and 0.15 m thick, and of the same channel width (Figure 1). At the outlet of the flume, water was intercepted by a stilling tank, equipped with three vertical grids to stabilize water, and a triangular weir (V-notch sharp crested weir) to measure the discharge with a relative uncertainty of ± 8%. At the downstream end of the flume, a movable gate made of Plexiglas and hinged at the channel bottom was used to regulate the flow depth.

To simulate grade control structures to protect riverbeds against erosion, in this study, we used a series of sills that consisted of PVC plates 0.30 m wide and 0.01 m thick. The sills were installed on an experimental area extended 6 m along the channel, downstream of the wooden ramp. The sill height decreased progressively as it went downstream from the wooden ramp, respecting a determined initial slope *S0* of 0.0086. Different configurations were investigated, the difference between them being the distance between GCSs (sills).

The flume bottom between the GCSs was covered with an erodible bed material layer, consisting of almost uniform sand particles with a median size, *d50*, of 1.8 mm and a density of 2650 kg/m3. Along the experimental area, the sand layer was leveled respecting the maximum GCS heights, forming the original bed of the channel with a slope *S*<sup>0</sup> (Figure 1).
