*Article* **Numerical Investigation on the Hydraulic Properties of the Skimming Flow over Pooled Stepped Spillway**

#### **Shicheng Li and Jianmin Zhang \***

State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu 610065, China; tyutscu@163.com

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

Received: 21 September 2018; Accepted: 15 October 2018; Published: 19 October 2018

**Abstract:** Pooled stepped spillway is known for high aeration efficiency and energy dissipation, but the understanding for the effects of pool weir configuration on the flow properties and energy loss is relatively limited, so RNG *k* − *ε* εturbulence model with VOF method was employed to simulate the hydraulic characteristics of the stepped spillways with four types of pool weirs. The calculated results suggested the flow in the stepped spillway with staggered configuration of' two-sided pooled and central pooled steps (TP-CP) was highly three dimensional and created more flow instabilities and vortex structures, leading to 1.5 times higher energy dissipation rate than the fully pooled configuration (FP-FP). In FP-FP configuration, the stepped spillway with fully pooled and two-sided pooled steps (FP-TP) and the spillway with fully pooled and central pooled steps (FP-CP), the pressure on the horizontal step surfaces presented U-shaped variation, and TP-CP showed the greatest pressure fluctuation. For FP-TP and FP-CP, the vortex development in the transverse direction presented the opposite phenomenon, and the maximum vortex intensity in TP-CP occurred at Z/W = 0.25, while FP-FP illustrated no significant change in the transverse direction. The overlaying flow velocity distribution in the spanwise direction demonstrated no obvious difference among FP-FP, FP-TP, and FP-CP, while the velocity in TP-CP increased from the axial plane to the sidewalls, but the maximum velocity for all cases were approximately the same.

**Keywords:** numerical simulation; pooled stepped spillway; pool weir; flow property

#### **1. Introduction**

Stepped spillway, an important water-transferring hydraulic infrastructure, has long been used in water related projects owing to its high efficiency of energy dissipation and air entrainment, making it possible to release massive and high speed flood flows while reducing and preventing dangerous scour and cavitation. Additionally, economic benefits will be enjoyed due to the reduction in the construction of the downstream energy dissipation basin [1,2]. Due to the advantages of stepped spillway, a growing number of this particular hydraulic structure is being applied to engineering projects.

For the past decades, stepped spillway has been studied globally by the means of laboratory experiments and numerical simulations [3–7], providing significant insights into the complicated air-water flow properties, momentum exchange process, and energy loss mechanism. For example, model tests performed by Felder and Chanson [8] to investigate the influence of uniform and non-uniform step height on the hydraulic performance illustrated that there was no significant change between the two configurations regarding the residual head. Channel slope also plays a vital role in the air-water flow properties, such as velocity profile and turbulence intensity reported by Felder and Chanson [9]. The time-averaged pressure on the step faces was experimentally measured by Zhang et al. [10], drawing the conclusions that the pressure distributions on the horizontal step surfaces exhibited S-shaped variations with the maximum and minimum values occurring at the downstream

end and upstream start, respectively. Physical models were established by Zhang and Chanson [11] to examine the effect of step edge and cavity geometry on the air entrainment and energy loss performance in stepped chutes. Apart from experimental studies, thanks to the advances in the computational performance of computers, numerical methods have also been successfully used to study the detailed flow properties in complex contexts [12–15]. Baylar et al. [16] numerically examined the energy loss and aeration efficiency of the stepped chutes and found that both of them in nappe flow regime were greater than that of the skimming flow regime. Numerical simulations by researchers [17–19] demonstrated its high performance in investigating the detailed flow field and pressure distribution over stepped spillway.

In recent years, an innovative design has been made to stepped spillway to improve the energy dissipation performance. Unlike the conventional stepped spillway whose horizontal step faces are flat, this novel kind of stepped spillway, pooled stepped spillway, is equipped with pool weirs at the step edges. A few of the researches have been carried out to investigate the air-water structures in this unique type of stepped spillway. For instance, the differences concerning the air entrainment processes of flat, pooled, and a combination of flat and pooled stepped spillway have been clarified by André [20] and Kökpinar [21]. Thorwarth [22] reported self-induced instabilities over pooled stepped spillways with channel slopes of 8.9◦ and 14.6◦. The detailed air-water flow properties over several kinds of stepped pooled spillway were also physically studied by Felder [7] and Felder and Chanson [9].

While the existing researches have gained understanding into this particular type of stepped spillway, little has been known about the impact of pool weir configurations on the flow properties. Therefore, in order to expand the limited insights into the effect of pool weir's geometrical properties on the flow characteristics, numerical investigations of pooled stepped spillways with different pool weir configurations were carried out.

#### **2. Numerical Simulation**

In the present study, ANSYS 16.0 (ANSYS®, Canonsburg, PA, USA) was utilized to investigate the flow characteristics over stepped spillway with different pool weir configurations. The numerical model consists of 4 parts: the upstream water tank, the broad-crested weir with an upstream rounded corner (radius = 0.08 m), the stepped spillway and the downstream channel. The length, *L*crest, and width, *W*, of the broad-crested weir are 1.01 m and 0.52 m, respectively. The spillway has 10 steps with vertical step height, *h*, and horizontal step length, *l*, of 0.1 m and 0.2 m, respectively. The channel slope is 26.6◦ and constant in each configuration. Moreover, the height, *d*, and the thickness, *lw*, of pool weir are, respectively, 0.09 m and 0.015 m. The schematics of the spillways are shown in Figure 1.

**Figure 1.** *Cont*.

**Figure 1.** Sketch of pooled stepped spillway configurations. (**a**) Fully pooled steps (FP-FP); (**b**) fully pooled and two-sided pooled steps (FP-TP); (**c**) fully pooled and central pooled steps (FP-CP); and (**d**) two-sided pooled and central pooled steps (TP-CP).

#### *2.1. Volume of Fluid Method*

Developed by Hirt and Nichols [23], the volume of fluid (VOF) model is believed to be a highly effective method to determine the position of the interface of two or more immiscible flows [19,24–26]. For the purpose of tracking the air-water interface, the continuity equation for the volume fraction of water is adopted in the following format:

$$
\frac{\partial \alpha\_w}{\partial t} + u\_i \frac{\partial \alpha\_w}{\partial x\_i} = 0 \tag{1}
$$

where *α<sup>w</sup>* is the volume fraction of water, *ui* and *xi* are the velocity components and coordinates (*i* = 1, 2, 3), respectively. Hence, by solving the continuity equation, the free surface of water can be determined. In addition, it is clear that the sum of their volume fractions, *α<sup>w</sup>* and *αa*, respectively, is unity in every computational cell, which can be described as:

$$
\kappa\_w + \kappa\_a = 1 \tag{2}
$$

Obviously, the maximum and minimum value for the volume fraction of water, *αa*, is 1 and 0, indicating the given cell is full of water or air, respectively, between which there must be a mixture of them containing the interface.

As for the density, *ρ*, and molecular viscosity, *μ*, they can be given as

$$
\rho = \alpha\_w \rho\_w + (1 - \alpha\_w)\rho\_a \tag{3}
$$

$$
\mu = \mathfrak{a}\_w \mu\_w + (1 - \mathfrak{a}\_w) \mu\_a \tag{4}
$$

where *ρ<sup>w</sup>* and *ρ<sup>a</sup>* represents the density of water and air, respectively, while *μ<sup>w</sup>* and *μ<sup>a</sup>* are the molecular viscosity of water and air, respectively. Therefore the values of *ρ* and *μ* can be known by iterating the Equations (3) and (4).

#### *2.2. Turbulence Model*

Proposed by Yakhot and Orszag [27], the Renormalization Group (RNG) *k* − *ε* closure model has been employed to simulate flow over stepped spillways in many studies [28–30], yielding reliable results compared with laboratory experiments. The continuity equation, momentum equation, and transport equations of turbulent kinetic energy, and the turbulent dissipation rate, used in the model, are shown below.

Continuity equation:

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

Momentum equation:

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

*k* equation:

$$\frac{\partial(\rho k)}{\partial t} + \frac{\partial(\rho\_i k)}{\partial x\_i} = \frac{\partial}{\partial x\_i} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_k} \right) \frac{\partial k}{\partial x\_i} \right] + G\_k - \rho \varepsilon \tag{7}$$

*ε* equation:

$$\frac{\partial(\rho\varepsilon)}{\partial t} + \frac{\partial(\rho u\_i \varepsilon)}{\partial x\_i} = \frac{\partial}{\partial x\_i} \left[ \left( \mu + \frac{\mu\_l}{\sigma\_\varepsilon} \right) \frac{\partial \varepsilon}{\partial x\_i} \right] + C\_{1t} \rho \frac{\varepsilon}{k} G\_k - C\_{2t} \rho \frac{\varepsilon^2}{k} \tag{8}$$

where *ρ* and *μ* respectively represent the average density of the volume fraction and the molecular viscosity. *p* is the modified pressure, and *μ<sup>t</sup>* is the turbulent viscosity, which can be calculated from the following equations consisting of *k* and *ε*. *μ<sup>t</sup>* = *C<sup>μ</sup> k*2 *<sup>ε</sup>* , *<sup>C</sup>*1*<sup>ε</sup>* <sup>=</sup> 1.42 <sup>−</sup> *<sup>η</sup>* <sup>1</sup><sup>−</sup> *<sup>η</sup> η*0 <sup>1</sup>+*βη*<sup>3</sup> , *<sup>η</sup>* <sup>=</sup> *Sk <sup>ε</sup>* , *S* = 2*SijSij*. The constants in the above expressions are given in the Table 1.

**Table 1.** Constants of governing equations.


#### *2.3. Boundary Conditions*

Combined with VOF method, RNG *k* − *ε* turbulence model was applied to numerically study the flow characteristics over stepped spillway with different pool weirs and the entire simulation domain was meshed with structured grids, as shown in Figure 2. The boundary conditions for the mathematical models were set as follows:


**Figure 2.** Meshing pattern in the computational domain.

#### *2.4. Grid Testing and Model Verification*

Before applying the numerical model to further analysis, it is necessary to validate its reliability. Therefore, the accuracy of the calculation results were tested by utilizing grid convergence index (GCI) [31] with grid numbers of approximately 0.31 million, 0.55 million and 0.76 million. Figure 3 presents the impact of the grid sizes on the uncertainty of the computational velocity distribution at the axial plane. The discretization uncertainties at most locations were quite small. The maximum GCI values for the velocity profiles over step 6 and step 7 were 7.7% and 5.6%, corresponding to ±0.19 m/s and ±0.15 m/s, respectively. Considering the simulation efficiency and accuracy, the grid number set as 0.55 million is reasonable.

In addition, experimental data obtained from Felder et al. [32] were used to verify the numerical results by comparing the velocity distribution over the pool at steps 8, 9, and 10, as shown in Figure 4. The schematic of the physical model used for validation was the same as the FP-FP configuration in the present article, except the pool weir height. In the research of Felder et al. [32], the pool weir height was 0.031 m. It can be seen from Figure 4 that, at most locations, the calculated results of velocity distribution were fairly consistent with that of laboratory tests. The maximums of the relative error for steps 8, 9 and 10 were 14.8%, 7.1% and 6.5%, respectively, indicating that the numerical simulations produced reliable and acceptable results.

**Figure 3.** Fine-grid solution with discretization error bars computed using the GCI index.

**Figure 4.** Comparison of numerical and experimental data. Notes: Vc: critical flow velocity; *y* represents the distance normal to the pseudo-bottom formed by the pool edges with *y* = 0 at the pool edges.

#### **3. Results and Analysis**

#### *3.1. Flow Pattern*

After ensuring the accuracy of the calculation, the effects of pool weir geometrical parameters on the hydraulic behaviors were numerically examined. Figure 5 illustrates the flow pattern at the axial plane in the calculated cases. It is noticeable that the free water surface of FP-FP configuration was parallel to the pseudo-bottom, and showed the least instability among the simulated configurations. FP-TP and FP-CP showed the similar free surface profile, but slight flow turbulences caused by the partially pooled steps can be observed. An obvious feature of the flow in TP-CP was the highly three-dimensional flow motion in the spanwise direction, so TP-CP highlighted strong fluctuations in free water surface.

**Figure 5.** Flow pattern at the axial plane in different configurations. (**a**) FP-FP; (**b**) FP-TP; (**c**) FP-CP; and (**d**) TP-CP.

Figure 6 highlights the streamlines on the representative steps in each configuration. As can be seen, like the traditional flat stepped spillway, the flow surface over the pool weir was approximately parallel to the pseudo-bottom formed by the pool edges, although TP-CP presented some instabilities caused by the staggered configurations of central pooled and two-sided pooled steps. For FP-FP configuration, vortexes with the similar size and position were formed in every pool and there was no significant change in the transverse direction, which was different from other configurations, that is, the flow pattern in FP-TP, FP-CP, and TP-CP was three-dimensional. In FP-TP, the scale of the vortex on the two-sided pooled steps increased from the axial plane to the sidewalls, while the scale of the vortex on the fully pooled steps decreased from the axial plane to the sidewalls. FP-CP presented the opposite phenomenon compared with FP-TP. In TP-CP configuration, the maximum vortex intensity occurred at Z/W = 0.25 and no eddy can be found on the central pooled steps at the sidewalls and on the two-sided pooled steps at the axial plane.

**Figure 6.** Instantaneous streamlines in the transeverse direction. (**a**) FP-FP; (**b**) FP-TP; (**c**) FP-CP; and (**d**) TP-CP.

For the identification of turbulence structures, *Q* criterion [33] is introduced to examine the effects of pool weir geometry on the vortex flow. The *Q* criterion is defined as follows:

$$Q = \frac{1}{2} \left( \left\| \Omega\_{\hat{i}\hat{j}} \right\|^2 - \left\| S\_{\hat{i}\hat{j}} \right\|^2 \right) \tag{9}$$

where

$$
\Omega\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} - \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{10}
$$

represents the rate of rotation tensor corresponding to the pure rotational motion, and

$$S\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \tag{11}$$

represents the rate of strain tensor corresponding to the irrotational motion.

The parameter *Q* is an effective approach to reveal the main motion pattern of fluid element. It can been seen from Equation (9) that the rotation motion will dominate the region where the value of *Q* is positive, while the deformation motion will dominate when *Q* is negative. However, it should be noted that the vortex structures cannot be accurately represented in the region with *Q* = 0, due to experimental error and viscous diffusion.

Figure 7 demonstrates the iso-surface of *Q* = 100. In FP-FP configuration, the number of lumps of the *Q* iso-surface distributing on each step was approximately the same and there was no apparent change in the transverse direction. The distribution of lumps in FP-TP and FP-CP exhibited the similar pattern. The major difference was that small lumps can be found near the sidewalls on the two-sided pooled steps in FP-TP, while small lumps were located near the axial plane in FP-CP, meaning turbulent structures existed at the two locations. It is noticeable that distribution of lumps in TP-CP configuration was much denser than that of other cases, indicating greater vortex and more turbulent structures.

**Figure 7.** Distributions of lumps of the iso-surface of *Q* = 100. (**a**) FP-FP; (**b**) FP-TP; (**c**) FP-CP; and (**d**) TP-CP.

#### *3.2. Velocity Distribution*

Figure 8 illustrates the velocity vectors at the axial plane. It implies that, in all configurations, the skimming flow transferring from the upstream to the downstream over the steps was partly trapped in the recirculating zones under the pseudo-bottom. However, it is evident that stable vortices can only be found on some steps at the axial plane in FP-TP, FP-CP, and TP-CP, which is consistent with the analysis from flow pattern. Furthermore, the velocity magnitude of the overlaying flow was significantly greater than that of the vortex.

**Figure 8.** Velocity vectors at the axial plane. (**a**) FP-FP; (**b**) FP-TP; (**c**) FP-CP; and (**d**) TP-CP.

Figure 9 displays the velocity distributions above the center of the horizontal step surfaces. In FP-FP, the velocity showed the same distribution in the transverse direction; the velocity distributions on each step presented the similar pattern, and no significant change between different steps can be observed for *y*/*y*max < 0.8, but in the flow surface region, the velocity above the downstream step surfaces was greater than that of the upstream ones; the minimum and maximum velocity above each step occurred at the same position: *y*/*y*max = 0.3 and *y*/*y*max = 0.9, respectively. In FP-TP configuration, at the axial plane, the velocity distributions above the fully pooled and two-sided pooled steps were obviously different in the bottom region, with the former greater than the latter; However, at Z/W = 0.25, there was no major difference in the velocity distribution above different steps except the free surface region. In FP-CP, at Z/W = 0.5, the velocity at 0.5 < *y*/*y*max < 0.8 for different steps indicated the same distribution; different from Z/W = 0.5, the velocity over different steps at Z/W = 0.25 exhibited similar distributions. In TP-CP, the velocity distributions presented a totally different pattern from other configurations; the velocity demonstrated the different distribution on two-sided pooled steps and central pooled steps, and velocity fluctuated greatly in the spanwise direction, proving its highly three dimensional characteristics.

**Figure 9.** Velocity distributions above the step surfaces in each configuration. (**a**) Z/W = 0.5; and (**b**) Z/W = 0.25.

In order to better compare the velocity field in the transverse direction, velocity contour maps at different cross-sections are shown in Figure 10. It can be noticed that, in all cases, the velocity distribution along the transverse direction was symmetrical and the overlaying flow over the pool weir showed the maximum velocity magnitude, while the minimum velocity value occurred near the bottom, but both the minimum and maximum indicated no apparent change in different configurations; In FP-FP, FP-TP and FP-CP, the velocity of the mainstream illustrated no significant difference along the spanwise direction, but the velocity in TP-CP fluctuated greatly with its value increasing from the axial plane to the sidewalls.

**Figure 10.** *Cont*.

**Figure 10.** Velocity countour maps at the cross-sections. (**a**) X = 1.91 m; and (**b**) X = 2.11 m.

#### *3.3. Pressure Distribution*

Figure 11 highlights the pressure distributions on the horizontal step surfaces of two representative steps (steps 6 and 7). In FP-FP configuration, the minimum and maximum pressure on different horizontal step surfaces showed no significant change, as well as the locations where the minimum and maximum pressure occurred. The pressure distribution patterns on the fully pooled steps in FP-TP and on the central pooled steps in FP-CP were similar with the extreme pressure occurring at the axial plane near the downstream end of the step face and the minimum pressure occurring in the middle the step surface; the pressure distributions on the two-sided pooled steps in FP-TP and on the fully pooled steps in FP-CP demonstrated the similar patterns with the extreme pressure occurring at the sidewalls near the downstream end of the step face. As for TP-CP configuration, the extreme pressures on the two-sided pooled and central pooled steps occurred near the sidewalls and axial plane, respectively.

**Figure 11.** *Cont*.

**Figure 11.** Pressure distributions on the horizontal step surfaces. (**a**) Step 6; (**b**) Step 7.

Figure 12 illustrates the pressure distributions on the vertical step surfaces of two representative steps (steps 6 and 7). As can been seen, in FP-FP, FP-TP and FP-CP, the pressure distributions on the vertical step surfaces highlighted the same pattern: the value of pressure decreased with the increasing in water depth, but the pressure for TP-CP presented some variations in the vertical direction. The maximum pressures on odd-number step surfaces in different configurations were in the following order: FP-CP > TP-CP > FP-TP > FP-FP, while for the even-number steps, it was TP-CP > FP-FP > FP-TP > FP-CP.

**Figure 12.** *Cont*.

**Figure 12.** Pressure distributions on the vertical step surfaces. (**a**) Step 6; and (**b**) Step 7.

Figure 13 presents the pressure distributions on the horizontal step surfaces in the transverse direction. In FP-FP, FP-TP and TP-CP, the pressure on the upstream step surfaces was greater than that of the downstream ones, but it was quite chaotic and showed no obvious trend in TP-CP. The pressure distributions on horizontal step surfaces of FP-FP, FP-TP, and TP-CP showed U-shaped variation with the maximum pressures occurring at the downstream end of the step surface. In addition, the greatest fluctuation in pressure in the spanwise and streamwise direction can be found in TP-CP, indicating its flow complexities and highly three-dimensional characteristics.

**Figure 13.** *Cont*.

**Figure 13.** Pressure distributions on horizontal step surfaces in the transverse direction.

Figure 14 presents the pressure distributions on the vertical step surfaces in the transverse direction. In FP-FP, FP-TP, and FP-CP, the pressure distributions on the vertical step surfaces displayed no significant change in the transverse direction, but great pressure fluctuation can be seen in TP-CP. In FP-TP and FP-CP, the pressure on the vertical step surfaces of the downstream steps was moderately greater than that of the upstream ones, while FP-FP showed no evident difference.

**Figure 14.** Pressure distributions on vertical step surfaces in the transverse direction.

#### *3.4. Residual Head and Energy Dissipation*

For the conservation of the downstream hydraulic structures, it is critical to quantify the residual head at the downstream of the spillway and the relative energy dissipation rate over the stepped spillway. The residual head, *Hres*, at the last pool weir can be estimated through the following equation:

$$H\_{\rm res} = d\_{\rm t} \times \cos \theta + \frac{\mathcal{U}\_{\rm v}^2}{2 \times g} + d \tag{12}$$

where *de* is the equivalent clear water depth, *θ* represents the channel slope, *Uw* and *d* indicate the flow velocity and the pool weir height, respectively.

Together with the previous studies in Table 2, the dimensionless residual head, *Hres*/*dc*, is shown in Figure 15. The flow discharge of the studies in Table 2 was the same and the flow regime was skimming flow. The residual energy in FP-FP configuration between the present work and the study of Morovati et al. [29] showed little difference, indicating the accuracy of the simulation. For configuration 1~3, the residual head showed no significant change for a given pool height, but the residual data increased with the decrease in pool height. In addition, the residual energy for the spillway with the combination of two-sided pooled steps and central pooled steps (TP-CP) decreased significantly compared with the spillway with only two-sided pooled steps or central pooled steps.



**Figure 15.** Residual energy at the downstream end of the spillway.

The energy dissipation rate, Δ*H*/*H*max, is an important indicator to evaluate the relative energy loss from the start to the end of the stepped spillway. The energy dissipation rate is given as:

$$\frac{\Delta H}{H\_{\text{max}}} = \frac{H\_{\text{max}} - H\_{\text{res}}}{H\_{\text{max}}} = \frac{H\_{\text{dam}} + 1.5d\_{\text{c}} - H\_{\text{res}}}{H\_{\text{dam}} + 1.5d\_{\text{c}}} \tag{13}$$

where *Hdam* represents the dam height. The rate of energy dissipation is presented in Figure 16. It can be seen that TP-CP showed the best energy dissipation performance compared with other configurations, and the energy dissipation rate for TP-CP was almost 1.5 times larger than that of FP-FP, but no obvious change can be observed in FP-FP, FP-TP and FP-CP. That is because there were more turbulent structures in TP-CP caused by the concave-convex pool weirs, leading to greater momentum change between the overlaying flow and the vortex flow, thus creating more energy loss.

**Figure 16.** Energy dissipation rate.

#### **4. Conclusions**

In the present study, four types of pooled stepped spillway with different pool weirs were numerically investigated by using RNG *k* − *ε* turbulence model with VOF method. The numerical results agreed well with the experimental data from Felder et al. [32] and the following conclusions can be drawn:


**Author Contributions:** Conceptualization, S.L.; Methodology, S.L.; Software, S.L.; Formal Analysis, S.L.; Investigation, S.L.; Writing-Original Draft Preparation, S.L.; 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 Natural Science Foundation of China (No. 51579165), and 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/).
