**Numerical Analysis on Hydraulic Characteristics of U-shaped Channel of Various Trapezoidal Cross-Sections**

#### **Ruichang Hu and Jianmin Zhang \***

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; hrc@stu.scu.edu.cn

**\*** Correspondence: zhangjianmin@scu.edu.cn

Received: 9 November 2018; Accepted: 1 December 2018; Published: 5 December 2018

**Abstract:** Curved channel with trapezoidal cross-section is approximate to the common form in nature fluvial networks and its hydraulic characteristics are considerably complex and variable. Combined with volume of fluid (VOF) method, renormalization group (RNG) k-ε turbulence model was employed to numerically investigate the flow properties in the U-shaped channel with various trapezoidal cross-sections. Analyses were performed from the aspects of the water surface transverse slope in bend apex (WTS-BA), longitudinal velocity, secondary flow, shear stress and turbulent kinetic energy (TKE) under several scenarios, specifically, four types of radius-to-width ratio and seven types of slope coefficient with a constant aspect ratio. The calculated results suggested that the maximums of shear stress and TKE in the bend were observed in the convex bank and the maximal intensities of secondary flow were observed within the range of 60 to 75 degrees for various varieties. As the radius-to-width ratio increased, the maximums of shear stress, TKE and WTS-BA decreased; but increased with increasing slope coefficients. The intensity of secondary flow decreased as slope coefficients increased and the angle of maximum intensity of secondary flow moved to the upstream for the increasing radius-to-width ratios. In addition, a new equation concerning the vertical distribution of longitudinal velocity in trapezoidal cross-sectioned channel was presented.

**Keywords:** trapezoidal cross-section; U-shaped channel; radius-to-width ratio; slope coefficient; hydraulic characteristics

#### **1. Introduction**

The water surface transverse slope (WTS), violent turbulence, forceful secondary flow and shear stress aggravate the erosion and siltation of the channel [1]. Consequently, detailed study regarding the hydraulic characteristics of the channel is an essential aspect for solving problems related to river restoration and navigability.

The secondary flow was observed by the movement of impurities in straight and curved open channel [2] and then the investigations referring to the secondary flow and shear stress in the channel were extensively conducted [3–6]. The structure of secondary flow was changed for the varying shape of cross-section in straight channel, which was investigated using three slope coefficients of 0.58, 1 and 1.73 by Tominaga et al. [7]. And the spanwise distribution of boundary shear stress was extraordinarily affected by the secondary flow and the shape of the cross-section [8] in the straight trapezoidal channel. The value of shear stress increased in the bottom of the wall and a dramatical variation from the apex to the middle of the sidewall was observed with increasing slope coefficients [8,9] in straight trapezoidal channel.

The curved trapezoidal channel is approximate to the common form in nature. The phenomenon of multi-vortex was also observed in curved trapezoidal channel [2], however, the influence of slope coefficient had not been explained completely. Farhadi et al. [10] investigated the transformation law of multi-vortex in the cross-section for the varying Froude number with the slope coefficient of 2.375 and Termini and Piranino [11] suggested that counter rotating circulation was conspicuous only in the lower aspect ratio, while the effect of various slope coefficients and radius-to-width ratio on secondary flow had not been studied. Shukla et al. [12] predicted the free surface flow velocities and shear stress reasonably well, using the Reynolds-averaged Navier-Stokes and continuity equations and the calibration of empirical coefficients and constants in standard k-ε turbulence model was recommended to improve accuracy. The modified k-ε model adopted in a 180 degree curved trapezoidal channel by Mosalman et al. [13] demonstrated good agreement with the experimental data. Also, k-ε turbulence model has been used in other different numerical approaches, such as the grid-less models [14,15]. The investigations on the water-air flow properties were conducted by Peng et al. [16–20]. The secondary flow and vortex structures can also be investigated via the direct numerical simulations (DNS) [21–23].

The slope coefficient and radius-to-width ratio are the essential parameter for the investigation on hydraulic characteristics of the channel [8]. The research on small slope coefficient, which is common in city, is not perfect and the objective of the present work is to bridge the gap in previous literature. Simultaneously, the effect of flow properties on the radius-to-width ratios was performed. The present study verified the feasibility of RNG k-ε turbulence model by the existing experiment [24]. Next, the analysis of the effect of varying radius-to-width ratios and slope coefficients on the WTS-BA, longitudinal velocity, secondary flow, shear stress and TKE was presented.

#### **2. Numerical Model**

A large number of U-shaped channel of trapezoidal cross-sections were modeled using Gambit software [25]. Fluent [26] was used to perform the numerical simulation and the renormalization group (RNG) k-ε turbulence model originated by Yakhot and Orszag [27] was employed. The U-shaped channel with symmetrical arrangement of trapezoidal cross-section is consisted of a half circle bend of 180◦ with different centerline radiuses and two straight cross-sections up and down the bend of 4m in length.

The water depth of the inlet and outlet of the numerical model were 0.203 m and 0.19 m, respectively. 3D view of the numerical model and vertical view were shown in Figure 1a,b, respectively. The slope coefficient of the sidewall in the curved channel remains uniform and the shape of the cross-section or centerline radius were varied for the study. The cross-section with different slope coefficients were presented in Figure 1c. Simultaneously, the area of the cross-section which was 0.08 m2 and the flow rate, 0.038 m3/s, remained constant. And constant equivalent width (*B*0) is 0.4 m, which was defined in Equation (1).

$$B\_0 = \mathbb{Q}/Hv \tag{1}$$

where *H* is water depth, *v* is the velocity and *Q* is the flow rate. In this study, *H* and *v* of the inlet of the U-shaped channel were adopted to calculate the equivalent width.

The slope coefficient (*m*) is the ratio of the projection of the wall on the horizontal plane to the vertical plane, as shown in Figure 1d. In order to distinguish the slope coefficient (*m*) from the unit meter (m), the slope coefficient (*m*) was shown in bold. The radius-to-width ratio (*λ*) can be defined as follows:

$$
\lambda = R / B\_0 \tag{2}
$$

where *R* is centerline radius of the bed in bend, which is defined in Equation (3); *B*<sup>0</sup> is equivalent width.

$$R = 0.5(R\_1 + R\_2) \tag{3}$$

where *R*<sup>1</sup> and *R*<sup>2</sup> are the inner radius and outer radius of the bed in bend, respectively.

The relationships between bed width, width of the surface and equivalent width are as follows:

$$B\_0 = B + \mathfrak{m}H \tag{4}$$

$$B\_{up} = B + 2\mathfrak{m}H \tag{5}$$

where *B* is the bed width, *Bup* is the width of the surface and *m* is slope coefficient.

**Figure 1.** *Cont.*

**Figure 1.** U-shaped channel and investigated cross-sections. (**a**) 3D view of the U-shaped channel; (**b**) vertical view of the channel; (**c**) geometry of the cross-section; (**d**) Definition of the slope coefficient.

The aspect ratio, the ratio of the equivalent width (*B*0) to the inlet water depth (*H*) of U-shaped channel, remained 1.97 in this study. The numerical simulations were performed for four radius-to-width ratios (*λ*) of 1.0, 1.5, 2.0 and 3.0 and seven slope coefficients (*m*) of 0, 0.25, 0.5, 0.75, 1.0, 1.10 and 1.24. The test runs program were summarized in Table 1. The *m*<sup>1</sup> and *m*<sup>2</sup> of case9 were 0.5 and 0, respectively and the m1 and m2 of case10 were 0.5 and 0.5, respectively. For convenience, the slope coefficients (*m*) of case9 and case10 were calculated by the definition (shown in Figure 1d) and the values were 1.10 and 1.24, respectively.

**Table 1.** Summary of the test runs for the simulation, where *m* is slope coefficient; *λ* is radius-to-width ratio; *θ* is the central angle of the bend; *B* and *B*<sup>0</sup> are the width of the bed and equivalent width; *R*, *R*<sup>1</sup> and *R*<sup>2</sup> are centerline radius, inner radius and outer radius of the bend, respectively; *Q* is the flow rate; *Fr* is Froude number; *Re* is Reynolds number; *B*0/*H* is the aspect ratio; and cells are the volume elements.


#### *2.1. Numerical Methods*

The volume of fluid (VOF) method, treating complicated free boundary configurations more flexible and efficient via volume fraction in stationary Euler mesh, was applied to track the air-water interface in the present work [28]. The differential equation in each control volume of the discrete region was integrated by the VOF method firstly and next the integral equation was linearized further to obtain the algebraic equation set of the corresponding parameters. Finally, all unknown quantities were solved. The solver employed pressure-based calculation, absolute velocity formulation and transient time and the SIMPLE algorithm was adopted for the pressure-velocity coupling. The computational domain was discretized using structured grid (shown in Figure 2) and the boundary layer cells were meshed in the vicinity of the wall and bed to ensure that they fell within the y<sup>+</sup> values required in the turbulence model. The number of the single trapezoidal U-shaped channel were 20 × 30 × 300, whereas for the compound channel, were 25 × 35 × 300. The number of volume cells for ten tests investigated in the paper were listed in Table 1.

**Figure 2.** Meshing pattern of the U-shaped channel. (**a**) case3; (**b**) case10.

#### *2.2. Governing Equations*

The continuity equation, momentum equation, TKE budget [23] and the turbulent dissipation rate (*ε*) [29], adopted in the model, are shown below:

Continuity equation:

$$\frac{\partial \rho}{\partial t} + \frac{\partial \rho \tilde{u}\_i}{\partial x\_i} = 0 \tag{6}$$

Momentum equation:

$$\frac{\partial(\rho \tilde{u}\_i')}{\partial t} + \frac{\partial(\rho \tilde{u}\_i' \tilde{u}\_j')}{\partial \mathbf{x}\_j} = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \frac{\partial \tilde{u}\_i'}{\partial \mathbf{x}\_j} + \frac{\partial \tilde{u}\_j'}{\partial \mathbf{x}\_i} \right) \right] + \frac{\partial}{\partial \mathbf{x}\_j} \left( -\rho \overline{u\_i' u\_j'} \right) \tag{7}$$

TKE budget:


$$\mathbf{C} = \frac{\partial \overline{\rho} \vec{u}\_{\dot{j}} \overline{\mathbf{K}}}{\partial x\_{\dot{j}}} \tag{9}$$

$$p = -\overline{\rho} \overline{u\_i'} \overline{u\_j'} \frac{\overline{\partial u\_i}}{\partial x\_j} \tag{10}$$

$$T = -\frac{\partial}{\partial x\_j} \left[ \overline{\rho} \frac{\overline{u\_i' u\_i' u\_j'}}{2} \right] + \overline{p' u\_i'} \sigma\_{ij} \tag{11}$$

$$
\Pi = -\overline{p'\frac{\partial u\_i'}{\partial x\_j}}\tag{12}
$$

$$\mathbf{M} = -\overline{\rho} \overline{u\_i'} \left( \frac{\partial \overline{\sigma\_{ij}}}{\partial x\_j} - \frac{\partial \overline{p}}{\partial x\_i} \right) \tag{13}$$

$$\mathcal{D} = \frac{\partial}{\partial x\_{\dot{j}}} \left( \overline{\sigma\_{i\dot{j}}' \mu\_i'} \right) \tag{14}$$

*k* equation:

$$\frac{\partial\left(\rho\overline{k}\right)}{\partial t} + \frac{\partial\left(\rho\overline{k}\overline{u}\_{\overline{i}}\right)}{\partial x\_{\overline{i}}} = \frac{\partial}{\partial x\_{\overline{j}}} \left[ \left(\mu + \frac{\mu\_{\overline{t}}}{\sigma\_{\overline{k}}}\right) \frac{\partial\overline{k}}{\partial x\_{\overline{j}}} \right] + G\_{\overline{k}} - \rho\varepsilon \tag{15}$$

*ε* equation:

$$\frac{\partial(\rho\overline{\varepsilon})}{\partial t} + \frac{\partial(\rho\overline{\varepsilon u\_i})}{\partial \mathbf{x}\_i} = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_t} \right) \frac{\partial \varepsilon}{\partial \mathbf{x}\_j} \right] + \frac{\mathbf{C}\_{1\varepsilon}^\* \overline{\varepsilon}}{\overline{\mathcal{K}}} \mathbf{G}\_k - \mathbf{C}\_{2\varepsilon} \rho \frac{\overline{\varepsilon}^2}{\overline{\mathcal{K}}} \tag{16}$$

$$
\mu\_t = \rho \mathbb{C}\_{\mu} \frac{\mathbb{Z}^2}{\overline{\varepsilon}} \tag{17}
$$

$$\mathbf{G}\_{k} = \mu\_{t} \left( \frac{\partial \overline{u\_{i}}}{\partial \mathbf{x}\_{j}} + \frac{\partial \overline{u\_{j}}}{\partial \mathbf{x}\_{i}} \right) \frac{\partial \overline{u\_{i}}}{\partial \mathbf{x}\_{j}} \tag{18}$$

$$\mathbf{C}\_{1\varepsilon}^{\*} = \mathbf{C}\_{1\varepsilon} - \frac{\eta \left(1 - \eta / \eta\_0\right)}{1 + \beta \eta^3} \tag{19}$$

$$
\eta = \sqrt{\frac{G\_k}{\rho \mathcal{C}\_{\mu\varepsilon}}} \tag{20}
$$

where *K* is the TKE, *μ* is the dynamic viscosity, *μ<sup>t</sup>* is the turbulence viscosity, *ρ* is the corresponding density, *ui* is the velocity component in the *i*th direction, *t* is the time and *Cμ*, *Cμ*, *C*1*ε*, *C*2*ε*, *β* and *η*<sup>0</sup> are empirical constants. In the RNG k-ε model, the empirical constants are given as *C<sup>μ</sup>* = 0.0845, *C<sup>μ</sup>* = *C<sup>ε</sup>* = 0.7179, *C*1*<sup>ε</sup>* = 1.42, *C*2*<sup>ε</sup>* = 1.68, *β* = 0.012, and *η*<sup>0</sup> = 4.38 [29].

The compressible mass flux contribution term is zero for the incompressible fluid. Reynolds stress independent components can be solved via using Equations (9)–(14) above and ε = *σ ik ∂u i ∂xk* , however, supererogatory components, that is, the trible correlation matrix for fluctuating velocity (*u i u i u j* ) and correlation matrix for fluctuating velocity and pressure (*pu i* ), are added in the turbulence kinetic energy budget, which made the budget is not closed. So K and ε equations are recommended to close the budget.

#### *2.3. Boundary Conditions and Assumptions*

In this study, without considering the dissolution of each other, the open channel flow was employed by VOF sub-models with air-water two phases. Pressure-inlet was adopted for the water surface and the pressure of the air-water interface was standard atmospheric pressure. Both sidewalls and bed of the channel, whose roughness height were 0.25 mm, were stationary wall and the fluent standard wall function approach was adopted. The velocity-inlet and pressure-outlet were adopted for the inlet and outlet of the U-shaped channel, respectively and the mean velocity of the inlet cross-section was 0.462 m/s.

#### *2.4. Model Verification*

It is necessary to validate the numerical model before conducting further research. Therefore, the mesh convergence study were performed via using grid convergence index (GCI) [30] with grid numbers of approximately 0.12 million, 0.20 million and 0.30 million. The effect of the grid sizes on the uncertainty of the computational velocity distribution in vertical direction was shown in Figure 3, which showed that the discretization uncertainties were little in most positions. The maximum of the GCI for the velocity profiles in concave bank, center line and convex bank were 6.3%, 2.4% and 8.2% and corresponding to 0.05 m/s, 0.06 m/s and 0.07 m/s, respectively. The grid number was set as 0.20 million to optimize the computational time.

**Figure 3.** Discretization error bars computed using the GCI index for different grid densities, (**a**) the concave bank, (**b**)center line and (**c**)convex bank of the bend apex.

In addition, the results calculated by RNG k-ε turbulence model were compared with the experimental data obtained from Ma [24]. The geometry of the experimental curved channel was exactly the same as case5 in the present article, except the lengths of two straight cross-sections of the experiment were 6 m. The sidewalls and bed composed of cement had a low roughness and the velocity was obtained by ADV in experiment.

In the present work, the length of two straight cross-sections all were 4 m or 6 m in the U-shaped channel, which were compared with the experiment, respectively. The RNG k-ε model was adopted for simulation with the same experimental conditions and the data of simulation and experiment were dimensionless and compared. The longitudinal velocity in different position of cross-section of C0, C90 and C180, the structure of secondary flow of the cross-section in the bend apex and shear stress in the bend were analyzed to verify the numerical model. In Figure 4, the vertical coordinate is the relative depth, where *Z* is the instantaneous height and *H* is the water depth in corresponding position; the horizontal coordinate is the ratio of the longitudinal velocity (*u*) to the corresponding average velocity (*uave*).

**Figure 4.** Comparison of dimensionless longitudinal velocity of various cross-sections between the numerical simulation and experiment, (**a**–**c**) the concave bank, center line and convex bank of the inlet of the bend, respectively; (**d**–**f**) the concave bank, center line and convex bank of the bend apex, respectively; and (**g**–**i**) the concave bank, center line and convex bank of the outlet of the bend, respectively.

It can be seen from Figure 4 that, in most locations, the longitudinal velocity distribution of the simulation exhibited a good agreement with the experimental in the vertical direction. The maximums of absolute error (the absolute value of the difference between experiment and numerical simulation) and relative error (the ratio of the absolute error to experiment) in different location were presented in Table 2. A well consistency was shown for the structure of secondary flow in Figure 5. Although the simulated shear stress in the bed was slightly larger than the experimental data, shown in Figure 6b, the commutated results were fairly consistent with Ma's investigation [24].

**Table 2.** Summary of the maximums of relative error and absolute error in different location. The locations of (a–i) were introduced in Figure 4.

**Figure 5.** Comparison of the secondary flow between the numerical simulation and experiment, (**a**) experiment [24], (**b**) numerical simulation L = 6 m, (**c**) numerical simulation L = 4 m.

Consequently, the RNG k-ε model could be used for the numerical investigation of the U-shaped channel and ensure the accuracy of the data. To optimize the computational time, the length of 4m were selected for the straight sections. Analyses were performed from the aspects of the WTS-BA, longitudinal velocity, secondary flow, shear stress and TKE in next section.

**Figure 6.** Comparison of the shear stress between the numerical simulation and experiment, (**a**) experiment [24], (**b**) numerical simulation L = 6 m, (**c**) numerical simulation L = 4 m.

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

#### *3.1. Analysis of the Water Surface Transverse Slope in Bend Apex*

The centrifugal motion of the water in bend is a result of centrifugal force, gravity and the friction with surrounding water. The transverse slope is presented on the water free surface, due to the descending water level in convex bank and rising water level in the concave bank [23]. The WTS varies with radius-to-width ratio, slope coefficient and cross-section shape [31]. The equation for calculating the WTS of cross-section of the bend apex (C90) is [31]:

$$J\_r = \left(1 + \frac{\mathcal{g}}{\kappa^2 \mathcal{C}^2}\right) \frac{\mu^2}{gr} \tag{21}$$

where *C* is the Chezy coefficient, *κ* is the Carmen constant, *g* is the acceleration of gravity, *r* is the radius and *u* is the corresponding velocity. For this investigation, the empirical values <sup>1</sup> *<sup>κ</sup>*<sup>2</sup> calculated by Zhang [32], Liu [31] and Rozovskii [32] were 5.75, 6.25 and 4, respectively. The cross-section in the bend apex (C90) was selected for further investigations, *r* and *u* were the corresponding centerline radius and the average longitudinal velocity of the cross-section of the bend apex, respectively.

The free surface level across the across-section in the bend apex with varying radius-to-width ratios and slope coefficients were shown in Figure 7a,b, respectively. The vertical coordinate is the water level and the horizontal coordinate is the relative distance from a point on the cross-section to the convex bank on the surface. The WTS-BA with varying radius-to-width ratios and slope coefficients were shown in Figure 8a,b, respectively. The vertical coordinate is a dimensionless WTS and the horizontal ordinate is radius-to-width ratio or slope coefficient.

**Figure 7.** The free surface level across the across-section in the bend apex (C90), (**a**) the radius-to-width ratios (*λ*) of 1.0, 1.5, 2.0 and 3.0, (**b**) the slope coefficients (*m*) of 0, 0.25, 0.5, 0.75, 1.0, 1.10 and 1.24.

**Figure 8.** The WTS-BA (C90), (a) the radius-to-width ratios (*λ*) of 1.0, 1.5, 2.0 and 3.0; (**b**) the slope coefficients (*m*) of 0, 0.25, 0.5, 0.75, 1.0, 1.10 and 1.24.

The WTS decreased with increasing radius-to-width ratios, which can be explained by the decrease in the water level in the concave bank and the increased water level in the convex bank, shown in Figures 7a and 8a. It can be considered that the numerical simulation exhibited a good agreement with the theoretical Equation (21) for the maximum error is 8.2%, in spite of the fact that the data obtained via simulation was smaller than calculated by Equation (21) with the varying Carmen constants.

The WTS increased with the increasing slope coefficients, for the water level in the concave bank slightly decreased first and then increased, while decreased in a relative larger extent in the concave bank, shown in Figures 7b and 8b. As can be seen from Figure 8b, the data obtained via numerical simulation was smaller than calculated by Equation (21) except the compound cross-section, the maximum errors of which were 7.2% and 10.1%, respectively.

The reason was that Equation (21) was derived on the basis of a single trapezoidal cross-section, which could not be applicable to the compound cross-section. In the whole, the numerical simulation was consistent with the Equation (21) for the varying slope coefficients. As can be seen from Figures 7b and 8a, the effect of the radius-to-width ratio was greater than the slope coefficient on the WTS in the bend apex.

#### *3.2. Analysis of the Longitudinal Velocity*

Chang [33] believes that the vertical distribution of the longitudinal velocity in the bend conforms to the Equation (22):

$$\frac{u}{u\_{\text{ave}}} = \left(1 + \frac{1}{n}\right) \left(\frac{Z}{H}\right)^{\frac{1}{n}}\tag{22}$$

where, *n* is obtained by equation *n* = *κ* -8/ *f* , *κ* is the Carmen constant, *f* is the friction constant; *Z* is the instantaneous height and *H* is the corresponding water depth and *u* and *uave* are the longitudinal velocity and the corresponding average in the cross-section. The vertical distributions of the longitudinal velocity calculated by Equation (22) were shown in Figures 9 and 10.

**Figure 9.** Vertical distribution of the dimensionless longitudinal velocity of the bend in different positions with varying radius-to-width ratios, (**a**–c) the concave bank, centerline and convex bank of the inlet of the bend, respectively; (**d**–**f**) the concave bank, centerline and convex bank of the bend apex, respectively; (**g**–**i**) the concave bank, centerline and convex bank of the outlet of the bend, respectively.

**Figure 10.** Vertical distribution of the dimensionless longitudinal velocity of the curved channel in different positions with varying slope coefficients, (**a**–**c**) the concave bank, center line and convex bank of the inlet of the bend, respectively; (**d**–**f**) the concave bank, center line and convex bank of the bend apex, respectively; (**g**–**i**) the concave bank, center line and convex bank of the outlet of the bend, respectively.

Figures 9 and 10 showed the vertical distribution of the dimensionless longitudinal velocity of the bend in different positions with the varying radius-to-width ratios and slope coefficients, respectively. In Figures 9 and 10, the vertical coordinate is the relative depth, where *Z* is instantaneous height, *H* is the corresponding water depth and the horizontal coordinate is the ratio of the longitudinal velocity (*u*) to the corresponding average velocity (*uave*). The velocity locating close to the sidewall at 5%*B*<sup>0</sup> was extracted to represent the velocity at the sidewall, as the velocity at the sidewall is zero.

As shown in Figure 10, the vertical distribution of longitudinal velocity calculated by Equation (22) was relatively consistent with the distribution of case5 (whose cross-section was rectangular). The longitudinal velocity in the bend was affected by the shape of cross-section, radius-to-width ratio and other factors [34–37]. The differences of the velocity distribution in the both banks were bigger than the center line, which was shown in Figures 9 and 10. The longitudinal velocity extracted via the numerical simulation conformed to the Equation (13) in the center line, yet the obvious errors were exhibited in the both banks. In this study, Equation (23) was proposed, which was consistent with the vertical distribution of dimensionless longitudinal velocity for the varying operating conditions in any position:

$$\frac{u}{u\_{\text{ave}}} = A\left(\frac{Z}{H}\right) + \frac{B}{\left(\frac{Z}{H}\right)} + \text{C} \tag{23}$$

where, *A* is equal to −0.09 or −0.03 for the sidewall or the center line, respectively, *B* is equal to −0.05 and the range of *C* is from 1 to 2.

Corresponding optimal fitting curves (the red continuous line) calculated by Equation (16) in different positions were shown in Figures 9 and 10. It can be seen that the vertical distribution of the longitudinal velocity in the bend with varying operating conditions were consistent with the optimal fitting curves in the center line but deviations are presented in the both banks. The deviations for the *λ* = 1.0 were conspicuous for the obvious gradient of the velocity in both banks and the vertical distribution for the *λ* > 1.0 were consistent with the optimal fitting curve. The deviations were lesser than others when the slope coefficient was 0.5. As shown in Figures 9 and 10, the influence of the slope coefficient on the distribution of longitudinal velocity was more conspicuous than the radius-to-width ratio.

The vertical distribution of the longitudinal velocity in both banks at the inlet of the bend were different for the varying operating conditions and demarcation were each *Z*/*H* = 0.5. As shown in Figure 9, the longitudinal velocity in the concave bank gradually decreased or increased for the *Z*/*H* < 0.5 or *Z*/*H* > 0.5 with increasing radius-to-width ratios but did not show the obvious trend for the varying radius-to-width ratios in the convex bank. As shown in Figure 10, the longitudinal velocity gradually increased or decreased for the *Z*/*H* < 0.5 or *Z*/*H* > 0.5 with increasing slope coefficients in the both banks. The difference of the velocity distribution between two types of compound cross-sections was slight.

The vertical distributions in the concave bank of the bend apex failed to show regular law for the varying operating conditions and were similar to the distributions in the convex bank of the inlet of the bend in the convex bank, whose demarcations were each *Z*/*H* = 0.6. The distributions of the longitudinal velocity of the outlet in the bend were similar to the distribution of the cross-section in the bend apex and the demarcations were *Z*/*H* = 0.4 and *Z*/*H* = 0.5.

#### *3.3. Analysis of the Secondary Flow*

The secondary flow is occurred for the momentum exchange [7]. The force of weeny water column are as follow [32]:

$$
\begin{array}{ccccc}
\Delta F &=& \mathbf{g} &+& \rho ghI\_{\mathcal{T}} &-& a\mathbf{u}\frac{\mathbf{u}^{2}}{\mathcal{S}^{2}} &+& \mathbf{F}\mathbf{o} \\
& & \text{gravity} & & \text{pressure} & & \text{centriffugal} & & \text{resistance} \\
& & & \text{force} & & \text{force} & & \\
\end{array}
\tag{24}
$$

where, *g* is the acceleration of gravity, *r* is the radius, *h* and *u* is the corresponding height and velocity, *Jr* is the transverse slope, *α*<sup>0</sup> is coefficient and *F*<sup>0</sup> is the summation of other resistance; in this investigation, gravity force is neglected in horizontal direction.

The structure of secondary flow of the cross-section in the bend apex for varying operating conditions was shown in Figure 11. In the present study, the clockwise vortex was defined as positive and counterclockwise vortex is reverse. As shown in Figure 11a–d, a reverse vortex is captured in case1 and the vortex structure do not change significantly from case2 to case4. The reason was that the WTS decreased sharply though the velocity increase in the concave bank with increasing radius-to-width

ratios (shown in Figure 8a). The ΔF was positive for the convex bank and negative for the concave bank only in case1 in horizontal direction.

**Figure 11.** Structure of secondary flow of the cross-section in the bend apex for varying operating conditions, (**a**–**j**) case1–case10, respectively.

As shown in Figure 11c,e–h, a reverse vortex are gradually obvious in the cross-section of the bend apex with increasing slope coefficients. The obvious double vortexes were captured in Figure 11g,h, the positive vortex near the convex bank was close to the bed but the reverse vortex was generated near the water surface in the concave bank. The reason was that velocity decreased in the concave bank (shown in Figure 11), yet the transverse slope increased (shown in Figure 8b) with increasing slope coefficients and the reverse vortex occur for the ΔF was negative at a certain position. Obvious

multi-vortex was captured in Figure 10i,j in the present study, due to the irregular velocity distribution caused by the compound section.

In this study, the equation of the intensity of secondary flow [38] is as follows:

$$
\omega = \frac{\mu\_w^2 + \mu\_v^2}{\mu\_l^2 + \mu\_w^2 + \mu\_v^2} \tag{25}
$$

where, *ω* is the intensity of secondary flow, *uw* is the streamwise velocity, *uw* is the spanwise velocity and *uV* is the vertical velocity.

Figure 12a,b showed the graph of dimensionless intensity of secondary flow for varying radius-to-width ratios and slope coefficients. The vertical or horizontal coordinate were the dimensionless intensity of secondary flow or the angle of the bend, respectively. The distribution of the intensity of secondary flow for *m* = 0 is consistent with the investigation of Ma [24], and Vaghefi et al. [39] also presented that the maximum intensity of secondary flow appears at approximately 60◦~80◦ via various analysis methods, which demonstrated good agreement with the above conclusion.

**Figure 12.** Graph of the dimensionless intensity of secondary flow for varying radius-to-width ratios and slope coefficients, (**a**) radius-to-width ratio, (**b**) slope coefficient.

The intensity of secondary flow decreased firstly, then increased and finally decreased with increasing slope coefficients. The reason why sudden increase was occurred was that the cross-section changes from trapezoid to compound cross-section. It can be considered that the intensity of secondary flow decreased with increasing slope coefficients but the intensity was the weakest when the slope coefficient was 1.0. As the radius-to-width ratio increased, the intensity of secondary flow decreased and the corresponding angle of the maximum intensity of secondary flow transforms from 105◦ to 60◦.

#### *3.4. Analysis of the Shear Stress*

The higher the shear stress on the bed and sidewall generated by the water in the bend, the more obvious erosion or siltation in channel and shear stress was investigated by many scholars [8,9,38,39]. The shear stress cannot be measured directly in actual research and various calculation methods were adopted to calculate the shear stress via other physical parameters, for example, energy slope and hydraulic radius [40]. The average velocity and correlation coefficient were employed via linear calculation to determine the shear stress by Williams et al. but this calculation method was not accurate in local location [41] Fluctuating velocity [42] and TKE [43] can be used to calculate the shear stress. The TKE method less affected by local streamline is more suitable for complex areas in the bend [44]. In this study, TKE method was adopted to calculate the shear stress and the calculation equation is as follows:

$$\tau = \mathbb{C}\_1 \left[ 0.5 \rho \left( \overline{\boldsymbol{\nu}}^2 + \overline{\boldsymbol{\nu}}^2 + \overline{\boldsymbol{\omega}}^2 \right) \right] \tag{26}$$

where *C*<sup>1</sup> is 0.19 [43], *ρ* is corresponding density and *u* , *v and ω* are the streamwise, spanwise and vertical fluctuating velocity, respectively.

The distributions of the shear stress in the bend were shown in Figure 13. The maximum shear stress and the corresponding angle for varying operating conditions were shown in Figures 14 and 15, respectively.

**Figure 13.** *Cont.*

**Figure 13.** Shear stress distribution of the riverbed and sidewall for varying operating conditions, (**a**)–(**j**) case1–case10, respectively.

**Figure 14.** Maximum shear stress for varying operating conditions, (**a**) radius-to-width ratio, (**b**) slope coefficient.

**Figure 15.** Corresponding angle of the maximum shear stress for varying operating conditions, (**a**) radius-to-width ratio, (**b**) slope coefficient.

According to Figure 13a–d, the transformation law of the shear stress in the convex bank, bed, concave bank decreased firstly and then increased with increasing radius-to-width ratios and the extent of maximum value variation were 27%, 15% and 12%, respectively (shown in Figure 14a). The position of the maximum shear stress moves to the upstream and the corresponding angle moved from 140◦ to 37◦ or 180◦ to 170◦ in the convex bank or bed (shown in Figure 15a) but did not move in the concave bank. Why shear stress was significantly higher than others for *λ* = 1 was that the momentum exchange and turbulence transfer violently for the dramatic change of the water in the bend.

The shear stress increased in the convex bank and concave bank with the increasing slope coefficients, which was consistent with the conclusion inferred by Ansar et al. [9], who revealed that the shear stress of sidewall increased with increasing slope coefficients in straight trapezoidal channel. Hence, the law of shear stress in both banks increased with increasing slope coefficients can be adopted in both straight and cured trapezoidal channels. The shear stress in the bed failed to present the regular variations with increasing slope coefficients and was higher than the trapezoidal cross-section in the convex bank of the compound cross-section in the same position. The position of maximum shear stress did not show the obvious law in both banks and bed with increasing slope coefficients.

#### *3.5. Analysis of the TKE*

Blanckaert K. [45] and Ma [24] investigated the distribution of TKE in the cross-section of the bend. In this study, the equation of TKE calculation is as follows:

$$K = \frac{1}{2} \left( \mu'^2 + \upsilon'^2 + \omega'^2 \right) \tag{27}$$

where *K* presents the TKE, *u* , *v* , *ω* are the streamwise, spanwise and vertical fluctuating velocity, respectively.

The graphs of the dimensionless TKE with the radius-to-width ratios and slope coefficients at various locations of the bend were shown in Figure 16. The relational graphs of the change extent of the dimensionless TKE with the radius-to-width ratios and the slope coefficients were shown in Figures 17 and 18, respectively. The TKE located at 5%*B*<sup>0</sup> was extracted to represent the sidewall. The turbulent kinetic energy was nondimensionalized via the mean turbulent kinetic energy (TKEave).

**Figure 16.** *Cont.*

**Figure 16.** Graphs of the dimensionless TKE with the radius-to-width ratios and slope coefficients in different positions of the bend, (**a**,**b**) the inlet of the bend; (**c**,**d**) the apex of the bend; (**e**,**f**) the outlet of the bend.

**Figure 17.** Relation graphs of the change extent of the dimensionless TKE with radius-to-width ratios, (**a**) the difference between inlet and apex of the bend, (**b**) the difference between apex and outlet of the bend.

**Figure 18.** Relation graphs of the change extent of the dimensionless TKE with slope coefficients, (**a**) the difference between inlet and apex of the bend, (**b**) the difference between apex and outlet of the bend.

According to Figure 16, the TKE on both sidewalls was larger than the center line in the bend. The TKE in the convex bank and central line gradually decreased and increased slightly in the concave bank with increasing radius-to-width ratios. The reason was that the turbulence and pulsation were weaker in the convex bank and central line and more violent in the concave bank with the increasing radius-to-width ratios. The TKE of all single trapezoidal cross-section increased with increasing slope coefficients.

As the flow of water in the bend, the TKE obviously increased in the convex bank and center line and obviously reduced for the concave bank in the cross-section of the bend apex, as shown in Figure 16a–d. The amplitudes decreased with increasing radius-to-width ratios and increased with increasing slope coefficients, which were shown in Figures 17a and 18a.

In the outlet of the bend, the TKE significantly increased in the concave bank and center line and decreased in the convex bank, as shown in Figure 16c–f. The amplitudes decreased with increasing radius-to-width ratios and decreased in the convex bank and center line with increasing slope coefficients. An increase was observed in the concave bank for the single trapezoidal cross-section but an opposite trend was presented for the compound cross-section with the increasing slope coefficients.

#### **4. Conclusions**

In this study, the numerical model was verified via the existing experimental data and then numerical simulation was performed to investigate the hydraulic characteristics of the U-shaped channel with trapezoidal cross-section under various conditions, that is, four types of radius-to-width ratio and seven types of slope coefficient with the aspect ratio being constant. The conclusions are as follow:


**Author Contributions:** Conceptualization and Methodology, J.Z.; Numerical Simulation, R.H.; Investigation, R.H.; Writing-Original Draft Preparation, R.H.; Writing-Review & Editing, J.Z.; Supervision, J.Z.; Project administration, J.Z.; Funding acquisition, J.Z.

**Funding:** This research was funded by the National Key Research and Development Program of China (No. 2016YFC0401707), National Science Fund for Distinguished Young Scholars (No. 51625901).

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

#### **References**


© 2018 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 (http://creativecommons.org/licenses/by/4.0/).
