Computational Methods and Applications to Simulate WaterRelated Natural Hazards
View this Special IssueResearch Article  Open Access
A Meshless WCSPH Boundary Treatment for OpenChannel Flow over SmallScale Rough Bed
Abstract
A weakly compressible smoothed particle hydrodynamics (WCSPH) method was developed to model openchannel flow over rough bed. An improved boundary treatment is proposed to quantitatively characterize bed roughness based on the ghost boundary particles (GBPs). In this model, the velocities of GBPs are explicitly calculated by using evolutionary polynomial regression with a multiobjective genetic algorithm. The simulation results show that the proposed boundary treatment can well reflect the influence of wall roughness on the vertical flow structure. A fully developed open channel is established, and its flume length, processing time, and turbulent model are discussed. The mixedlengthbased subparticle scale (SPS) turbulence model is adopted to simulate uniform flow in open channel, and this model is compared with the Smagorinskybased one. For the modified WCSPH model, the results show that the calculated vertical velocity and turbulent shear stress distribution are in good agreement with experimental data and fit better than the calculations obtained from the traditional Smagorinskybased model.
1. Introduction
Computational fluid dynamics (CFD) is widely used in hydraulic engineering simulations with good computational efficiency and accuracy. As an important branch of CFD, meshless methods such as smoothed particle hydrodynamics (SPH) method, moving particle semiimplicit (MPS) method, and discrete element method (DEM) have received wide attention due to their ability to consider large deformations of free interface and multiple interfaces without mesh distortion. The SPH method is a widely used Lagrangian meshless method that was originally developed by Gingold and Monaghan [1] and Lucy [2] to solve astrophysics problems and later used in fluid dynamics applications such as fluidstructure interaction [3–5], underwater explosion [6, 7], local scour [8, 9], other transport phenomena [10], and related applications [11–13]. Despite advances in SPH methodology, the SPH method has not yet been widely adopted for the simulation of openchannel flows in natural rivers and artificial channels. Thus, the focus of this study was to address the use of SPH for simulation of openchannel flows.
Boundary conditions can directly affect the calculation accuracy and efficiency for the accurate simulation of openchannel flows. However, it is quite challenging to achieve strict boundary conditions in the SPH framework due to the interpolation procedures. Previous studies of numerical simulations of free surface flows mainly focused on boundary treatments of inflow and outflow. A variety of nonreflecting boundary schemes such as open boundary condition [6, 14] and sponge layer [15–17] have been applied in the field of free surface flow to minimize undesirable boundary effects. However, the influence of bed surface roughness on the flow field has not been explicitly stated or well considered in the most openchannel flow simulation [18–22]. In some works, the bottom boundary is roughly treated as slip or noslip wall [18–20]. López et al. [21] proposed that uneven distribution of LennardJones repulsive force exerting on particles near the solid wall can be considered as numerical resistance caused by roughness, which may suffer from spurious behavior [23]. Chern and Syamsuri [22] used corrugated solid walls to represent rough bed surface resistance, in which different bed boundaries such as triangletype, trapezoidtype, and sinusoidtype are characterized by discrete particles. Although these simple approaches are reasonable approximations, the roughness of the bed surface can only be qualitatively expressed in these works, and thus, the effect of the bed roughness on the fluid structure cannot be accurately described.
An additional complication is that openchannel beds in nature have different forms and different scales of roughness, such as fine sand bed, pebble river bed, or vegetation river bed. The bed roughness can exert an important influence on the flow field characteristics, such as average flow rate, timeaveraged pressure, Reynolds stress, and turbulent energy. The slip/noslip treatment on the bottom boundary directly affects the overall accuracy of the simulation results for the openchannel flow. Therefore, the bottom boundary should be carefully described. To the best of our knowledge, there are two major types of approaches to accurately account for the effects of bed roughness on the openchannel flow. One approach is to express the effect of the rough bed on the fluid by adding a source term in the fluidgoverning equation [24, 25]. For example, Khayyer and Gotoh [24] investigated the effect of bed friction on a dambreak flow by applying a frictional drag force between a nearbed fluid particle and its neighboring wall particle. However, the proposed drag force is a purely artificial drag force that lacks physical meaning such as bed roughness. Kazemi et al. [25] also proposed a dragbased formulation to account for the largescale roughness bed surface in depthlimited turbulent openchannel flow. In that work, the thickness of the roughness zone ranged from 0.009 m to 0.011 m, which refers to gravel according to river sediment classification [26]. It is not clear whether this drag force model could be applied to a smallsized rough bed surface, such as silt and sand with particle size that is at least an order of magnitude lower than that of gravel. Another approach to account for the effects of bed roughness on the openchannel flow is to describe the shearing effect of the rough boundary on the ambient flow by indirectly modifying the velocity gradient of the bed particles and ambient fluid particles [27, 28]. Violeau and Issa [27] used the wall function model to reproduce the logarithmic velocity distribution near a solid wall, but the particle disorder can be clearly found near the bed surface (see Figure 9 in their paper). The velocity of ghost particles in the bottom boundary was revised by Fu and Jin [28] to force the vertical velocity distribution of water flow to logarithmic form. The ghost particle velocity is thought to be only related to water depth and roughness height but independent of the slope of the bed surface. Similar to the second approach [27, 28], Barreiro et al. [29] modified the artificial viscosity value for interactions between fluid and bottom to mimic the Manning roughness coefficient of vegetation [26]. However, there is still no universal roughness boundary of free surface flow in the current SPH framework for different flow patterns such as hydraulic smooth flow, hydraulic rough flow, or hydraulic transition flow. Therefore, a focus of this work is the proposal of a more comprehensive rough bed boundary treatment method for the SPH framework.
In addition to boundary conditions, the turbulence model is equally important in simulating openchannel flows. The turbulence phenomena are quite complex, preventing full resolution through fullproof theory. Existing turbulence models are primarily based on semiempirical assumptions. For example, the Reynolds stress caused by momentum exchange can be related to the average flow field by the eddy viscosity assumption [30]. One of the earliest and commonly used turbulent models in meshfree method is the subparticle scale (SPS) turbulence model, proposed by Gotoh et al. [31]. The main idea of the SPS turbulence model is to apply a certain filtering function to decompose the instantaneous motion of the turbulent flow into two parts: largescale vortex and smallscale vortex. Largescale vortex can be directly solved using the Navier–Stokes (NS) equation, and smallscale vortex is set by SPS formulation. Subsequently, the SPS turbulence model was adopted to applications such as dam break [32], hydraulic jump [22], and sediment transport [33], although most cases are mainly based on a smooth bed surface. In addition to the spatial averagebased SPS model, other turbulence models are based on the Reynolds average of the flow field physical quantities in the time domain, such as the kε model [23, 34] and the explicit algebraic Reynolds stress model (EARSM) [27]. These two models can successfully predict complex free surface with strong distortion and rotation. However, these models [27, 34] do not consider the effect of bed roughness on fluid turbulence, and the nearwall particle distribution is chaotic. Recently, De Padova et al. [35] investigated undular hydraulic jumps over smallscale rough bed for bed roughness of 0.02H, where H denotes the water depth. The SPH mixinglength turbulence model exhibits satisfactory performance for simulating hydraulic jump as long as strong turbulent rollers are not present and is applicable for openchannel flow simulation [25, 36]. Overall, there have been insufficient studies of the openchannel flow over rough bed using SPH methods. Therefore, the aim of this work was to further investigate roughbed boundary treatment and the turbulent model and extend this model to simulate openchannel flows in natural rivers and artificial channels with wide hydraulic conditions.
The remainder of this paper is organized as follows: After the introduction, the methodology of the SPH model is presented in detail and an improved bottom boundary treatment is proposed for smallscale rough bed. Next, a numerical open channel is established with a revised turbulence model for openchannel flow in Section 3. In Section 4, turbulent open channel flows over a smallscale rough bed surface are simulated, and the error is analyzed. In Section 5, the conclusions are summarized.
2. Materials and Methods
2.1. SPH Approximation
SPH is a purely Lagrangian method, using a series of particles to discretize the continuum. For the fluid dynamics, the value of field functions such as mass, density, velocity, and pressure at each discrete particle can be approximated based on the physical properties of surrounding particles:where the subscripts a and b denote the particles at locations of r_{a} and r_{b}, respectively; f(r_{a}) and f(r_{b}) denote the field function of the target particle a and its neighboring particle b, respectively; W denotes the kernel function; h denotes the influencing domain of the kernel function or smoothing length; and ΔV_{b} denotes the volume of particle b. ΔV_{b} = m_{b}/ρ_{b}, where m_{b} and ρ_{b} denote the mass and density of particle b, respectively.
A commonly used smoothing kernel is a quintic function, as proposed by Wendland [37], and can be expressed aswhere α_{D} is taken as 7/4πh^{2} for the 2D case and q denotes the dimensionless distance between particles.
2.2. Governing Equations
The Reynoldsaveraged Navier–Stokes (RANS) equations can be expressed aswhere ρ denotes the fluid density, u denotes the velocity vector of fluid, p denotes the fluid pressure, ν_{0} denotes the kinetic viscosity, denotes the turbulence stress tensor, and denotes the gravitational acceleration.
In the SPH notation, the governing equations can be rewritten as [38,39]where the subscripts a and b denote the index of the target particle and its neighbor in the support domain, respectively. u_{a,}p_{a}, , and refer to the velocity, pressure, shear stress, and gravity acceleration of particle a, respectively. u_{ab} and r_{ab} refer to the relative velocity and location between particles a and b and are expressed as u_{ab} = u_{a}–u_{b} and r_{ab} = r_{a} – r_{b}. The ∇_{a}W_{ab} refers to the gradient of smoothing kernel W(r_{a} − r_{b}, h) for particles a and b.
To close the governing equations (4) and (5), the subparticlescale (SPS) model is adopted to model the turbulence stress tensor with its elements aswhere subscripts i and j denote the index of spatial coordinates, respectively. ν_{t} denotes the eddy viscosity coefficient, S_{ij} denotes an element of the SPS strain tensor, , k denotes the turbulence kinetic energy, C_{I} = 0.0066, Δl denotes particle spacing, and δ_{ij} denotes the Kronecker delta function.
In the original SPS turbulence model, the eddy viscosity coefficient is taken as = (C_{s}Δl)^{2}S, where C_{s} denotes the Smagorinsky constants with a range of 0.1∼0.2 and S denotes the local strain rate and can be calculated as S = (2S_{ij}S_{ij})^{1/2}. In this study, the eddy viscosity coefficient is modified with reference to mixed length theory [36] and can be expressed as follows:where l_{m} denotes the mixed length, denotes the local strain rate, κ denotes the von Karman constant, and ζ denotes the vertical relative position and is equal to y/H. П denotes the Coles parameter, and when it approaches 0, equation (8) degenerates to
Following Dalrymple and Rogers [38], the fluid in the SPH formalism is considered weakly compressible with fluid pressure determined via the Tait state equation as follows:where c_{0} denotes the speed of sound at the reference density ρ_{0} and γ denotes the isentropic expansion factor and is taken as 7 for the case of water. The sound speed should be at least ten times larger than the maximum fluid velocity to restrict density variation to less than 0.01. Fluid pressure is explicitly expressed by the state equation without solving the Poisson equation iteratively, so the problem of poor convergence and low computational efficiency of the linear system can be avoided.
To limit spurious highfrequency noise in the density field, the dissipative term proposed by Molteni and Colagrossi [40] should be included in the right side of equation (4) and can be written aswhere δ denotes the deltaSPH coefficient and is taken as 0.1. Note that this is a numerical artifact that is used only in SPH notation.
2.3. Boundary Treatment
2.3.1. Inlet and Outlet Conditions
The computing domain is divided into three parts of inﬂow domain, internal fluid domain, and outﬂow domain (Figure 1). The inflow and outﬂow domains are considered part of the buffer layer to apply appropriate boundary conditions to the internal fluid. The width of inflow and outﬂow domains should be at least as wide as the smoothing radius to guarantee full function of kernel interpolation for fluid particles in the internal fluid region. In this study, the buffer length of the inlet/outlet is taken as two times the smoothing length (or three times the particle diameter), as used in Wang et al. [6] and Tafuni et al. [14]. The particles in the inﬂow, internal, and outﬂow domains are called inflow particles, internal fluid particles, and outﬂow particles, respectively. The inflow and outflow flow particles are treated as fluid particles, as are the internal fluid particles.
The properties of the internal fluid particles are determined by the governing equations, and the properties of the inflow and outflow particles are controlled by the specified boundary conditions. The velocities of particles in the inflow and outflow domains are controlled by a specified flow velocity, which is analogous to the physical pumping process. The density and pressure conditions for the inflow and outflow particles are both determined by the hydrostatic pressure of the free surface flow. Because the velocity and pressure of particles in the buffer layer are forcibly allocated, particle instability and diverging velocity field near the interface between the internal fluid and the buffer layer could arise. Therefore, the flume length should be long enough to avoid undesirable boundary effects and to reduce any overconstrained issues in inflow and outflow regions.
Similar to the periodic boundary condition proposed by GomezGesteira et al. [41], the outflow particles that pass through the end of the outflow region will immediately flow into the computational domain from the front of the inflow region (Figure 1). Unlike most open boundary treatments like Verbrugghe et al. [42], there are no particle additions and deletions within the calculation domain, so the total number of particles remains the same and mass conservation of fluid is implicitly satisfied. The particle positions are adjusted according to the periodic boundary, and the particle properties such as velocity, density, and pressure in the buffer layer of the inlet and outlet are artificially specified. Thus, it is not necessary to extend the support domain of outflow particles into the inlet domain.
2.3.2. Free Surface Condition
No special attention is paid to the particles at the free surface because the continuous density method (equation (4)) is adopted to solve particle density. The kinematic condition is intrinsically satisfied at the free surface by the Lagrangian solver. In the simulation of openchannel flows over a smallscale rough bed, when the flow in the water flume tends to stabilize, the water elevation in the streamline direction is almost constant.
2.3.3. Bottom Boundary Condition
Openchannel water flow is always affected by fluid resistance and energy dissipation, especially for rough bed. The roughness of the bed surface significantly influences the formation and development of the water flow structure, i.e., the vertical flow velocity distribution. To quantitatively describe the bed surface roughness, equivalent roughness is generally introduced, also known as Nikuradse’s equivalent roughness [43]. The flow state can be divided into hydraulic smooth flow, hydraulic rough flow, and hydraulic transition flow according to the bed surface roughness. For a hydraulic smooth flow (, where denotes friction Reynolds number, is the friction current velocity, k_{s} is the equivalent roughness, and ν_{0} is the kinematic viscosity), the bed surface roughness height is smaller than the adhesive bottom layer thickness, the coarse elements are submerged in the viscous bottom layer, and the bulk flow does not perceive bed roughness. For this condition, the flow velocity distribution in the logarithmic zone and the outer zone is not affected by the bed surface roughness but is affected by the viscosity of the fluid. With regard to hydraulic rough flow (), the bed surface roughness height is much greater than the viscous bottom layer thickness. The flow velocity distribution in the logarithmic zone and the outer zone is affected by the roughness of the bed surface but is not affected by the viscosity of the fluid. For hydraulic transition flow, the roughness of the bed surface is equivalent to the thickness of the viscous bottom layer, and the flow velocity distribution can be affected by both surface roughness and fluid viscosity.
The treatment of bed surface roughness in the SPH method has always been a complex issue. If the governing equations of equations (1) and (2) are solved using a direct numerical simulation (DNS) method, the solid boundary and the nearby turbulence field can be directly resolved in the computational mesh. However, DNS methods are computationally intensive and thus are not usable for most engineering problems. Therefore, a generalized wall function model is proposed to simulate turbulent openchannel flow over the smallscale rough bed surface. The wall function in the particle method is similar to that in the grid method, but this function does not directly modify the velocity of the first layer of the nearwall mesh. In the present study, the velocity of the solidwall boundary particle is modified to create a velocity gradient with the ambient fluid particles. Based on the equivalent roughness of the bed surface, the bottom shear stress of the nearwall fluid particle is corrected by adjusting the virtual velocity of the solidwalled particle, thereby modifying the flow velocity profile according to the following expression:where denotes the dimensionless current velocity (), denotes the dimensionless vertical distance (), and κ denotes the von Karman constant and is equal to 0.41. denotes the dimensionless distance between the top of the rough element and the reference bed surface, where the logarithmic curve zero is located. denotes the dimensionless vertical distance offset of velocity log profile (Figure 2). B denotes the logarithmic integration constant and can be written aswhere denotes the friction Reynolds number.
For the hydraulic smooth flow and the hydraulic transition flow, the wall function method can be adopted. For the hydraulic rough flow, the rough scale should be evaluated, and the smallscale rough bed can be adopted. The scale of roughness in the SPH method is a relative concept, where smallscale rough is if , and largescale rough is if , where is equal to 0.25k_{s} [44]. The wall function method is only applicable to smallscale roughness because the horizontal current velocity in the range does not meet the logprofile distribution. The shear stress can be correctly transferred to the fluid in the log rate region only if the distance between the inner water particles and the side wall is greater than the thickness of the bottom boundary layer (Figure 2).
The ghost particle velocity of the bottom boundary must be determined by trial calculation. When the numerical result of the crosssectional velocity distribution is closest to the theoretical flow velocity distribution, the ghost particle velocity corresponds to the equivalent roughness.
2.4. Solution Algorithm
A twostage predictorcorrector algorithm is used for the numerical integration scheme in this study [45], with the fluid properties of particle a, i.e., velocity (u_{a}), position (r_{a}), and density (u_{a}) initially predicted at the half time step n + 1/2, and can be written asand corrected at the time step n + 1 with the following expressions:
As suggested by Monaghan and Kos [46], the time step is then adjusted according to the Courant–Friedrichs–Lewy (CFL) condition, and the forcing terms and the viscous diffusion term can be written aswhere CFL denotes the Courant–Friedrichs–Lewy number and f_{a}/m_{a} denotes the fluid force per unit mass exerted on particle a.
3. Simulation of OpenChannel Flows
3.1. Model Setup and Verification
To validate the proposed model, a case with water depth H = 0.05 m, bed slope s = 0.002, and bed roughness k_{s} = 0.01 m was tested. The bed roughness of the open channel was calculated in accordance with Section 2.3, with the velocity of ghost particles in bottom boundary being determined by trial calculation and taken as −0.025 m/s. The height of the model calculation domain is equal to the water depth, and the length is 100 times the water depth (Figure 3). The initial particle spacing values (d_{p}) were set as 0.002 m, 0.003 m, or 0.004 m. For the case with particle spacing of 0.002 m, the whole computational domain is modeled by using 3753 solid particles and 31275 water particles. A computer with an octacore Intel® Core™ i76700 processor (clock speed of 3.40 GHz and 16.0 GB RAM) was used to carry out the simulation. The simulated time of 8 s costs nearly 8.8 h of CPU time.
In the test scenario, the numerical sensors were evenly arranged in the middle section to monitor the sectional water level process (Figure 3). The wave gauges WG1–4 were positioned in the streamline direction at x = 50H, x = 60H, x = 70H, and x = 80H, respectively. The water surface is considered stable as long as the error between the sensor value and the target water level is less than 2%. As described in Section 2.3, inflow and outflow particles in the buffer layer were assigned hydrostatic pressure and specific flow velocity. The specific flow velocity was obtained by trial calculation based on the following principle. When the water levels determined by the four sensors are consistent with the target water level (Figure 3), the flow velocity for the particles in the buffer layer is considered reasonable. The water surface can be stabilized when the specific current velocity in the buffer layer is equal to 0.312 m/s. The maximum error between the sensor value of water elevation and the target value is 0.08% (Table 1).

Figure 4 presents the simulated velocity and pressure fields for the model including the wall function. The instantaneous flow field is captured at the moment of 6 s. The results show that fully developed uniform flow is obtained at the middle of water flume without numerical noise.
(a)
(b)
3.2. Optimization of the Calculation Domain
The fluid boundary layer and openchannel turbulence can be fully developed under a sufficiently long computational domain in the streamline direction. However, it is necessary to minimize the length of the computing domain due to limited computing resources. The suitable flume length must be determined by trial calculation and is initially taken as up to 100 times the water depth. The length of the front, middle, and rear sections of the numerical water flume is taken as 50, 20, and 30 times the water depth, respectively (Figure 3). The development of the openchannel flow directly depends on the flow distance. Figure 5 presents the simulated results of current velocity versus the variation of flow distance for measuring points located at different water depths. Due to speed constraints in the inflow and outflow regions, the upstream uniform velocity fluid gradually stratifies until the velocity profile along the water depth becomes consistent with the logarithmic one. The downstream fluid is also affected by the outflow region, and the average flow velocity of this fluid will converge to 0.312 m/s. For a particle size of 0.002 m, the current velocity of the upper fluid with y = H increases with increasing flow distance and tends to stabilize to 0.36 m/s at x ≈ 35H. However, the fluids in the lower layer show different trends. For example, at y = 1/5H, the current velocity gradually decreases to 0.25 m/s at x ≈ 25H as flow distance increases due to the friction from the bottom boundary. In Figure 5(b), it can be seen that the profile of the flow velocity along the water depth convergences to a stable value until the flow distance reaches 40 times the water depth. The flow velocity constraint in the outflow region will nonphysically affect the flow field structure. At a distance of about 5H from the end of the water flume, the flow velocity of each monitoring section significantly changes to 0.312 m/s. The results show little effect of particle size on the influence range of the inflow and outflow boundaries (Figure 5(a)). The velocity development range of inflow particles and the influence range of the outlet boundary are about 40H and 10H, respectively, for different particle sizes. This indicates that the fully developed flow field structure including the flow velocity profile can be achieved with a flume length of 50H. Thus, to reduce computing resources, the computing domain is taken as 50 times the water depth in the following study.
(a)
(b)
3.3. Optimization of TimeAverage Operation
The development of openchannel flow requires a certain amount of time, and the simulation time is directly related to the boundary and initial conditions. Therefore, for a specific model, it is crucial to determine the appropriate simulation time. Figure 6 presents the evolution process of the streamline component of current velocity at the position of x = 40H. The flow velocities at different water depths are developed from an initial value of 0.312 m/s towards stable values so that the vertical velocity distribution is in accordance with logarithm law. In general, a large particle spacing corresponds to a shorter flow time or convergence time (Figure 6(a)). For fine particle spacing of 0.002 m, the current velocities at all depths become stable after a flow time of 6 s. The flow time required for the fluid near the free surface to reach a steady velocity is larger than that close to the flume bottom, indicating that the boundary layer is developing upward from the flume bottom.
(a)
(b)
Since fluid turbulence can cause spatiotemporal fluctuation in fluid parameters such as current velocity and water pressure, especially for fluid near the flume bottom, time and space averaging are required to reduce the oscillation of flow parameters. The particle approximation process (equation (1)) in the SPH method is essentially spatial interpolation [6], where the physical quantities of a target particle are obtained from its neighbor particles in the supporting domain. A timeaverage operation is adopted to update the physical quantities of the boundary particles. To determine the time interval of the timeaverage operation, it is necessary to carry out time convergence analysis of the main flow rate at different depths. Figure 6(b) presents the relative standard deviation (RSD) values of the current velocities in the streamline direction for measuring points at different water depths. The RSD values of coarse particles (d_{p} = 0.004 m) are relatively smaller than those of finer ones, and the convergence time is relatively shorter for coarser particles. For cases with different particle spacings, the flow velocities at all depths can reach convergence within 1.75 s∼2 s, with an RSD value in the range of ±0.5%. Except for the current velocity near the bottom wall, current velocities at other water depths can be considered to converge by 1 s.
The average current velocities under four other flow conditions were analyzed in the same way, with case I (H = 0.05 m, s = 0.001, k_{s} = 0.00005 m), case II (H = 0.08 m, s = 0.001, k_{s} = 0.01 m), case III (H = 0.08 m, s = 0.002, k_{s} = 0.002 m), and case IV (H = 0.1 m, s = 0.001, k_{s} = 0.01 m), where H denotes water depth, s denotes bed slope, and k_{s} denotes equivalent roughness. The simulation results show that the average main flow velocities under the four flow conditions with different water depths, bed slopes, and particle spacings can be stabilized after 6 s (Figure 7(a)) and that the time interval used for time averaging is less than 2 s (Figure 7(b)). Therefore, in this analysis, the processing time of 2 s was adopted to average fluid quantities.
(a)
(b)
3.4. SmagorinskyBased and MixedLengthBased SPS Model
The original SPS model [31], known as the Smagorinskybased SPS turbulence model, has been successfully applied to some coastal applications such as fluidstructure interaction [4, 47, 48]. However, the performance of this model to simulate openchannel flow is not satisfactory. Figure 8 shows the comparison of current velocity and turbulent shear stress with changing water depth for the Smagorinskybased SPS turbulence model with different Smagorinsky constants (C_{s}). None of the simulated results were consistent with the theoretical result. This is likely because the tested C_{s} values are suitable for nonconstant flow problems and should be revised for a case of constant uniform flow.
(a)
(b)
(c)
(d)
(e)
(f)
For cases with C_{s} value being of similar magnitude, the velocity profiles for different C_{s} values are almost the same (Figure 8(a)). The value of C_{s} affects the current velocity distribution of the openchannel flow if the C_{s} value has a different order of magnitude. In general, the larger the C_{s} value, the greater the numerical error, with overestimation of current velocity near the water surface and underestimation near the channel bottom (Figure 8(b)). The C_{s} value also has a significant influence on the fluid shear stress. A larger C_{s} value corresponds to larger fluid shear stress (Figure 8(c)), but even the largest shear stress obtained in the Smagorinskybased model is less than the theoretical value (Figures 8(d)–8(f)).
The C_{s}Δl value shows a certain degree of influence on the flow field distribution of the openchannel flow. Figure 9 shows the effect of particle spacing on the current velocity and shear stress for the Smagorinskybased SPS turbulence model. A more uniform velocity profile can be obtained by using a relatively larger C_{s}Δl value (Figure 9(a)). The shear stress increases with increasing C_{s}Δl (Figure 9(b)), while even the largest value is still lower than the theoretical shear stress. The low turbulent shear stress obtained by the Smagorinskybased turbulence model is attributed to the loss of the turbulent stress by spatial filtering during the longterm development of the water flow. In this model, the particle spacing has almost no effect on the flow velocity at free surface due to the loss of the surface shear force.
(a)
(b)
Mayrhofer et al. [49] simulated openchannel flow using a Smagorinskybased eddy viscosity model and found that the resolution of the particlebased method needed to be at least 16 times higher than that of the gridbased method in order to correctly redistribute the energy between the Reynolds stress components in the SPHLES frame. This indicates that the openchannel flow cannot be well resolved if the resolution size is not sufficiently dense. For the Smagorinskybased turbulence model, it is difficult to reconcile the simulated velocity profile with the theoretical value by adjusting particle resolution for this particle scale.
To obtain a reasonable shear force, some researchers have tried to modify the original SPS model based on other perspectives. For example, Gotoh et al. [48] modified the Smagorinskybased model in the nearwall region. The original turbulence eddy viscosity ν_{t} = (C_{s}Δl)^{2}S was modified to ν_{t} = [min (C_{s}Δl, κ·d_{wall})]^{2}S, where C_{s} denotes the Smagorinsky constant, Δl denotes the particle spacing, denotes the von Karman constant, d_{wall} denotes the minimum distance of particles from the wall, and S denotes the local strain rate. The treatment proposed by Gotoh et al. [48] is quite similar to the simplified mixed length formula proposed by Nezu et al. [50]. In this study, a mixedlengthbased SPS model (equation (7)) is adopted to simulate the openchannel flow. As shown in Figure 10, the turbulent shear stress is obviously improved by using the SPS model based on mixed length. The vertical distribution of the mixed length and Reynolds stress obtained by numerical simulation are quite similar to the theoretical values. Since this method allows determination of the correct turbulent shear stress, the vertical velocity profile can be well matched with the theoretical curve.
(a)
(b)
(c)
(d)
4. General Boundary Treatment for Rough Bed Surface
In Section 3, the improved wall function model and turbulence model in the WCSPH framework are successfully employed to represent the openchannel flow over a rough bed with a roughness height of 0.001. To apply the wall function model to various types of rough bed, a general boundary treatment for a rough bed surface is proposed in this section.
4.1. Investigation of Bed Surface with Wide Roughness
To investigate openchannel flow with wide bed roughness, a series of numerical cases were designed with reference to the actual size of river sediment particles [26] (Table 2). The values of equivalent roughness were selected as 0.00005 m, 0.001 m, 0.002 m, and 0.01 m; the values of the water depth are 0.05 m, 0.08 m, and 0.1 m; the values of bed slope are 0.001, 0.002, and 0.004; and the values of particle spacing are 0.002 m, 0.003 m, and 0.004 m. As described in Section 2.3, bed surface roughness is achieved by exerting proper wall shear stress on the bottom boundary. The shear stress can be correctly transferred to the fluid only if the first layer of fluid particles near the lower wall is in the log rate region of openchannel flow. This means that the particle spacing should be larger than the dimensionless distance between the top of the rough element and the reference bed surface, or , where equal to 0.25k_{s}.

The numerical cases were divided into hydraulic smooth flow, hydraulic rough flow, and hydraulic transition flow based on calculation of the friction Reynolds number. Here, a total of 159 sets of cases were performed. Among them, 51 sets made up the hydraulic smooth flow experimental group, 39 sets made up the hydraulic transition flow experimental group, and 69 sets made up the hydraulic rough flow experimental group (see Supplemental data (available here)). The test case is named according to the water depth, slope, equivalent roughness, and particle spacing. For example, for the case mentioned in Section 3.1 (H = 0.05 m, s = 0.002, k_{s} = 0.01 m, d_{p} = 0.002 m), the test case name is H50s2ks10(2).
The distribution of current velocity in the entire vertical section should be consistent with that of the theoretical curve by specifying the appropriate ghost particle velocity, so the ghost particle velocity can be matched with flow condition and bed surface roughness. Figure 11 presents the simulated results of the dimensionless velocity distribution for the model with wall function. Satisfactory results were obtained regardless of the flow regime. However, there was a significant deviation at the free surface and near the bottom surface, especially for the nearbottom area. The deviation at the free surface is mainly due to the absence of Coles’ wake function [51]. However, the apparent deviation near the bottom is the inevitable result of using a meshless method with the wall function model. The wall function itself cannot resolve particle velocity inside the transition zone and the viscous bottom layer, only improving the flow velocity distribution in and/or outside the logarithmic region. In addition, since the computational domain is composed of particles, the shear stress of the fluid must be transmitted through interactions between particles. It is necessary to exert reasonable force on the internal fluid particles through the ghost particles. However, strictly speaking, the shear force is exerted at the positions of the upper fluid particles yet not directly on the upper fluid particle itself, making the external force subject to the particle size.
(a)
(b)
Figure 12 shows the effect of particle spacing on the wall function method. The particle spacings shown in Figure 12 are d_{p1}, d_{p2}, and d_{p3}, respectively, where d_{p1} < d_{p2} < d_{p3}. The shear stress required by the bottom fluid particles for each case is equal to the theoretical shear stress corresponding to the height of the fluid particle. That is, the smaller the particle size, the larger the required shear stress of the fluid and the smaller the required ghost particle velocity.
We next performed an error analysis of the velocity distribution of all the test cases (Figure 13). The average relative error (MRE) describes the error of the WCSPH model coupled with the wall function method. In the vertical direction, the fluid region can be roughly divided into three layers, the inner zone (y/h < 0.2), the main zone (0.2 ≤ y/h ≤ 0.8), and the nearsurface zone (y/h > 0.8). There is only a small relative average error in the mainstream and nearsurface regions, roughly below 2%, while the relative average error in the inner region varies widely, mainly due to the influence of particle size. As shown in Figure 13, common laws can be obtained from the error distribution of typical cases with different water flow conditions: (i) the error in the mainstream area is generally smaller than that in the inner area; (ii) the error in the mainstream area and the nearsurface area is hardly affected by the test conditions but significantly affected by the test conditions in the inner area; (iii) for the case with the same particle size, the deeper the water depth, the smaller the data error in the inner zone; (iv) the error is positively correlated with bed roughness, especially for hydraulic rough flow; and (v) the bed slope has almost no effect on the error of the inner zone for hydraulic smooth flow, but the error in the inner zone increases with the increase of the channel slope for hydraulic rough flow.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
According to the above findings, the bed slope and the equivalent roughness are positively correlated determinants of the numerical error. Obviously, for the same hydraulic condition, the result will be more accurate for smaller particle spacing, but more computing resources will be needed.
4.2. Explicit Rough Wall Treatment
From the analysis presented above, a reasonable current velocity distribution can be obtained by adjusting the velocity of the ghost particles on the bed surface. However, this method would be cumbersome if it was necessary to determine the ghost particle velocity every time by trial and error. Therefore, it is necessary to establish relationships between ghost particle velocity and case condition such as bed slope, bed surface roughness, water depth, and particle resolution.
As a powerful datamining method, the evolutionary polynomial regression with multiobjective genetic algorithm (MOGAEPR) technique is widely used to construct nonlinear polynomial expression [52, 53]. It is interesting that the variable relationship and structure parameters of polynomial expression do not need to be identified in prior, making this approach especially suitable to solve multiobjective value optimization problems with unclear relationship structure.
Based on the simulation results presented in Section 4.1, a total of 159 sets of training data including hydraulic smooth flow, hydraulic rough flow, and hydraulic transition flow were used to determine the ghost particle velocity, essential to exert reasonable turbulence shear stress in openchannel flow. In this study, six potential dimensionless candidate variables were selected in terms of relative water depth (h_{dp} = H/d_{p}), relative bed surface roughness (k_{sdp} = k_{s}/d_{p}), channel bed slope (s), roughness scale (h_{ks} = H/k_{s}), dimensionless water depth (), dimensionless particle resolution (), and dimensionless bed surface roughness (). The nonnegative leastsquares method was adopted for global regression, with logarithmic function as the interfunction of the Pareto front [54]. The maximum number of terms of the Pareto sets was taken as 4. Three empirical formulae were derived from the EPRMOGA in terms of hydraulic smooth flow, hydraulic rough flow, and hydraulic transition flow and were written aswhere , h_{dp}, k_{sdp}, h_{ks}, , and denote the dimensionless ghost particle velocity (), the relative water depth, relative bed surface roughness, dimensionless water depth, dimensionless particle resolution, and dimensionless bed surface roughness, respectively.
From Table 3, it can be seen that (i) numerical boundary condition is indeed related to particle resolution as described in Section 4.1; (ii) returned equations have satisfactory accuracy and also conform to physical insights. For example, the velocity of the ghost particle on the bed surface is positively correlated with bed slope and dimensionless particle resolution (see Figure 12 for detail); and (iii) some potential information can be inferred from the empirical formula. For example, the ghost particle velocity is almost irrelevant to the relative water depth.
 
X_{i} and A_{j} denote the number of independent variables and term number, respectively. SSE and CoD denote the sum of squared errors and determining coefficient, respectively. The symbols +, −, and / denote positive correlation, negative correlation, and nocorrelation between the dependent variable and the independent variables, respectively. 
5. Conclusions
In this study, an indepth study on the simulation of the openchannel flow over roughness bed based on WCSPH was performed. Through the numerical test, the influence of bed surface roughness on the flow structure was clarified. Numerical tests showed that the Smagorinskybased SPS turbulence model performs poorly near the rough bed surface in the openchannel flow, showing deviation of the vertical velocity and turbulent shear stress distribution from the theoretical values, but the improved mixedlengthbased SPS turbulence model achieved satisfactory fluid structure.
Based on the principle of the wall function model in the grid method, an improved meshless wall function was developed for smallscale rough boundary and the sensitivity of this method was determined and error analysis was performed. There was relatively small numerical error in the mainstream and nearsurface regions, but the error varied widely for the inner region. This error is positively correlated with the channel bed slope and the equivalent roughness. In addition, for the convenience of model application, empirical formulae are proposed to calculate the ghost particle velocity under different hydraulic conditions using data mining. The proposed rough bed model is expected to be useful for openchannel turbulence simulation.
Data Availability
The training data of the multiobjective genetic algorithm used to support the findings of this study are included within the supplementary information file.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The improved WCSPH code used in this paper was developed from the open source SPHysics code developed by Prof. M. GomezGesteira at University of Vigo, Prof. R.A. Dalrymple at Johns Hopkins University, and others. This research was supported by the National Key Research and Development Project (grant nos. 2017YFC0403600 and 2016YFE0201900), the State Key Laboratory of HydroScience and Engineering, Tsinghua University, Beijing, China (grant no. 2017ky04); and the China Postdoctoral Science Foundation (grant no. 2018M641372).
Supplementary Materials
Supplemental test results including average currently in inflow and outflow regions, friction velocity, Reynolds number, Froude number, friction Reynolds number, and the velocity of ghost particle in the wall function model can be accessed from the online version of the paper. Totally, 159 sets of training data in terms of hydraulic smooth flow, hydraulic rough flow, and hydraulic transition flow are opensourced for a potential researcher. (Supplementary Materials)
References
 R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamics: theory and application to nonspherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, no. 3, pp. 375–389, 1977. View at: Publisher Site  Google Scholar
 L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,” The Astronomical Journal, vol. 82, pp. 1013–1024, 1977. View at: Publisher Site  Google Scholar
 B. Ren, M. He, Y. Li, and P. Dong, “Application of smoothed particle hydrodynamics for modeling the wavemoored floating breakwater interaction,” Applied Ocean Research, vol. 67, pp. 277–290, 2017. View at: Publisher Site  Google Scholar
 B. Ren, H. Wen, P. Dong, and Y. Wang, “Numerical simulation of wave interaction with porous structures using an improved smoothed particle hydrodynamic method,” Coastal Engineering, vol. 88, no. 2, pp. 88–100, 2014. View at: Publisher Site  Google Scholar
 A.M. Zhang, P.N. Sun, F.R. Ming, and A. Colagrossi, “Smoothed particle hydrodynamics and its applications in fluidstructure interactions,” Journal of Hydrodynamics, vol. 29, no. 2, pp. 187–216, 2017. View at: Publisher Site  Google Scholar
 P. Wang, A.M. Zhang, F. Ming, P. Sun, and H. Cheng, “A novel nonreflecting boundary condition for fluid dynamics solved by smoothed particle hydrodynamics,” Journal of Fluid Mechanics, vol. 860, pp. 81–114, 2018. View at: Publisher Site  Google Scholar
 F. R. Ming, A. M. Zhang, Y. Z. Xue, and S. P. Wang, “Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions,” Ocean Engineering, vol. 117, pp. 359–382, 2016. View at: Publisher Site  Google Scholar
 D. Wang, S. Li, T. Arikawa, and H. Gen, “ISPH simulation of scour behind seawall due to continuous tsunami overflow,” Coastal Engineering Journal, vol. 58, no. 3, Article ID 1650014, p. 1, 2016. View at: Publisher Site  Google Scholar
 D. Wang, S. Shao, S. Li, Y. Shi, T. Arikawa, and H. Zhang, “3D ISPH erosion model for flow passing a vertical cylinder,” Journal of Fluids and Structures, vol. 78, pp. 374–399, 2018. View at: Publisher Site  Google Scholar
 Y. Shi, S. Li, H. Chen, M. He, and S. Shao, “Improved SPH simulation of spilled oil contained by flexible floating boom under wavecurrent coupling condition,” Journal of Fluids and Structures, vol. 76, pp. 272–300, 2018. View at: Publisher Site  Google Scholar
 H. Gotoh and A. Khayyer, “On the stateoftheart of particle methods for coastal and ocean engineering,” Coastal Engineering Journal, vol. 60, no. 1, pp. 79–103, 2018. View at: Publisher Site  Google Scholar
 D. Violeau and B. D. Rogers, “Smoothed particle hydrodynamics (SPH) for freesurface flows: past, present and future,” Journal of Hydraulic Research, vol. 54, no. 1, pp. 1–26, 2016. View at: Publisher Site  Google Scholar
 M. S. Shadloo, G. Oger, and D. Le Touzé, “Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: motivations, current state, and challenges,” Computers & Fluids, vol. 136, pp. 11–34, 2016. View at: Publisher Site  Google Scholar
 A. Tafuni, J. M. Domínguez, R. Vacondio, and A. J. C. Crespo, “A versatile algorithm for the treatment of open boundary conditions in smoothed particle hydrodynamics GPU models,” Computer Methods in Applied Mechanics and Engineering, vol. 342, pp. 604–624, 2018. View at: Publisher Site  Google Scholar
 C. Altomare, J. M. Domínguez, A. J. C. Crespo et al., “Longcrested wave generation and absorption for SPHbased DualSPHysics model,” Coastal Engineering, vol. 127, pp. 37–54, 2017. View at: Publisher Site  Google Scholar
 S. J. Lind, R. Xu, P. K. Stansby, and B. D. Rogers, “Incompressible smoothed particle hydrodynamics for freesurface flows: a generalised diffusionbased algorithm for stability and validations for impulsive flows and propagating waves,” Journal of Computational Physics, vol. 231, no. 4, pp. 1499–1523, 2012. View at: Publisher Site  Google Scholar
 K. Gong, H. Liu, and B.L. Wang, “Water entry of a wedge based on SPH model with an improved boundary treatment,” Journal of Hydrodynamics, vol. 21, no. 6, pp. 750–757, 2009. View at: Publisher Site  Google Scholar
 A. J. Sahebari, Y.C. Jin, and A. Shakibaeinia, “Flow over sills by the MPS meshfree particle method,” Journal of Hydraulic Research, vol. 49, no. 5, pp. 649–656, 2011. View at: Publisher Site  Google Scholar
 A. Shakibaeinia and Y.C. Jin, “MPSbased meshfree particle method for modeling openchannel flows,” Journal of Hydraulic Engineering, vol. 137, no. 11, pp. 1375–1384, 2011. View at: Publisher Site  Google Scholar
 I. Federico, S. Marrone, A. Colagrossi, F. Aristodemo, and M. Antuono, “Simulating 2D openchannel flows through an SPH model,” European Journal of Mechanics—B/Fluids, vol. 34, pp. 35–46, 2012. View at: Publisher Site  Google Scholar
 D. López, R. Marivela, and L. Garrote, “Smoothed particle hydrodynamics model applied to hydraulic structures: a hydraulic jump test case,” Journal of Hydraulic Research, vol. 48, no. S1, pp. 142–158, 2010. View at: Publisher Site  Google Scholar
 M.J. Chern and S. Syamsuri, “Effect of corrugated bed on hydraulic jump characteristic using SPH method,” Journal of Hydraulic Engineering, vol. 139, no. 2, pp. 221–232, 2013. View at: Publisher Site  Google Scholar
 M. Ferrand, D. R. Laurence, B. D. Rogers, D. Violeau, and C. Kassiotis, “Unified semianalytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method,” International Journal for Numerical Methods in Fluids, vol. 71, no. 4, pp. 446–472, 2013. View at: Publisher Site  Google Scholar
 A. Khayyer and H. Gotoh, “On particlebased simulation of a dam break over a wet bed,” Journal of Hydraulic Research, vol. 48, no. 2, pp. 238–249, 2010. View at: Publisher Site  Google Scholar
 E. Kazemi, A. Nichols, S. Tait, and S. Shao, “SPH modelling of depthlimited turbulent open channel flows over rough boundaries,” International Journal for Numerical Methods in Fluids, vol. 83, no. 1, pp. 3–27, 2017. View at: Publisher Site  Google Scholar
 S. Dey, Fluvial Hydrodynamics: Hydrodynamic and Sediment Transport Phenomena, Springer, Berlin, Germany, 2014.
 D. Violeau and R. Issa, “Numerical modelling of complex turbulent freesurface flows with the SPH method: an overview,” International Journal for Numerical Methods in Fluids, vol. 53, no. 2, pp. 277–304, 2007. View at: Publisher Site  Google Scholar
 L. Fu and Y.C. Jin, “A meshfree method boundary condition technique in open channel flow simulation,” Journal of Hydraulic Research, vol. 51, no. 2, pp. 174–185, 2013. View at: Publisher Site  Google Scholar
 A. Barreiro, J. M. Dominguez, A. J. C. Crespo, H. GonzalezJorge, D. Roca, and M. GomezGesteira, “Integration of UAV photogrammetry and SPH modelling of fluids to study runoff on real terrains,” PLoS One, vol. 9, no. 11, pp. 1–11, 2014. View at: Publisher Site  Google Scholar
 J. Boussinesq, Theorie Analytique de la Chaleur, GauthierVillars, Paris, France, 1903.
 H. Gotoh, T. Shibahara, and T. Sakai, “Subparticlescale turbulence model for the MPS method Lagrangian flow model for hydraulic engineering,” Advanced Methods for Computational Fluid Dynamics, vol. 9, pp. 339–347, 2001. View at: Google Scholar
 R. Issa, D. Violeau, and D. Laurence, “A first attempt to adapt 3D large eddy simulation to the smoothed particle hydrodynamics gridless method,” in Proceedings of the International Conference on Computational & Experimental Engineering Sciences, 1st Symposium on Meshless Methods, StaraLesna, Slovakia, June 2005. View at: Google Scholar
 G. Pahar and A. Dhar, “On force consideration in coupled ISPH framework for sediment transport in presence of freesurface flow,” Environmental Fluid Mechanics, vol. 18, no. 3, pp. 555–579, 2018. View at: Publisher Site  Google Scholar
 D. Violeau, “One and twoequations turbulent closures for smoothed particle hydrodynamics,” Hydroinformatics, vol. 2, pp. 87–94, 2004. View at: Publisher Site  Google Scholar
 D. De Padova, M. Mossa, and S. Sibilla, “SPH modelling of hydraulic jump oscillations at an abrupt drop,” Water, vol. 9, no. 10, p. 790, 2017. View at: Publisher Site  Google Scholar
 E. Gabreil, S. J. Tait, S. Shao, and A. Nichols, “SPHysics simulation of laboratory shallow free surface turbulent flows over a rough bed,” Journal of Hydraulic Research, vol. 56, no. 5, pp. 727–747, 2018. View at: Publisher Site  Google Scholar
 H. Wendland, “Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree,” Advances in Computational Mathematics, vol. 4, no. 1, pp. 389–396, 1995. View at: Publisher Site  Google Scholar
 R. A. Dalrymple and B. D. Rogers, “Numerical modeling of water waves with the SPH method,” Coastal Engineering, vol. 53, no. 23, pp. 141–147, 2006. View at: Publisher Site  Google Scholar
 E. Y. Lo and S. Shao, “Simulation of nearshore solitary wave mechanics by an incompressible SPH method,” Applied Ocean Research, vol. 24, no. 5, pp. 275–286, 2002. View at: Publisher Site  Google Scholar
 D. Molteni and A. Colagrossi, “A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH,” Computer Physics Communications, vol. 180, no. 6, pp. 861–872, 2009. View at: Publisher Site  Google Scholar
 M. GomezGesteira, B. D. Rogers, A. J. C. Crespo, R. A. Dalrymple, M. Narayanaswamy, and J. M. Dominguez, “SPHysics—development of a freesurface fluid solver—Part 1: theory and formulations,” Computers & Geosciences, vol. 48, pp. 289–299, 2012. View at: Publisher Site  Google Scholar
 T. Verbrugghe, J. M. Domínguez, C. Altomare et al., “Nonlinear wave generation and absorption using open boundaries within DualSPHysics,” Computer Physics Communications, vol. 240, pp. 46–59, 2019. View at: Publisher Site  Google Scholar
 H. Schlichting and K. Gersten, BoundaryLayer Theory, Springer, New York, NY, USA, 2016.
 L. C. Van Rijn, Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas, Aqua Publications, Netherlands, Amsterdam, 1993.
 J. J. Monaghan, “On the problem of penetration in particle methods,” Journal of Computational Physics, vol. 82, no. 1, pp. 1–15, 1989. View at: Publisher Site  Google Scholar
 J. J. Monaghan and A. Kos, “Solitary waves on a Cretan beach,” Journal of Waterway, Port, Coastal, and Ocean Engineering, vol. 125, no. 3, pp. 145–155, 1999. View at: Publisher Site  Google Scholar
 S. Shao and H. Gotoh, “Simulating coupled motion of progressive wave and floating curtain wall by SPHLES model,” Coastal Engineering Journal, vol. 46, no. 2, pp. 171–202, 2004. View at: Publisher Site  Google Scholar
 H. Gotoh, S. Shao, and T. Memita, “SPHLES model for numerical investigation of wave interaction with partially immersed breakwater,” Coastal Engineering Journal, vol. 46, no. 1, pp. 39–63, 2004. View at: Publisher Site  Google Scholar
 A. Mayrhofer, D. Laurence, B. D. Rogers, and D. Violeau, “DNS and LES of 3D wallbounded turbulence using smoothed particle hydrodynamics,” Computers & Fluids, vol. 115, pp. 86–97, 2015. View at: Publisher Site  Google Scholar
 I. Nezu, H. Nakagawa, and G. H. Jirka, “Turbulence in openchannel flows,” Journal of Hydraulic Engineering, vol. 120, no. 10, pp. 1235–1237, 1994. View at: Publisher Site  Google Scholar
 D. Coles, “The law of the wake in the turbulent boundary layer,” Journal of Fluid Mechanics, vol. 1, no. 2, pp. 191–226, 1956. View at: Publisher Site  Google Scholar
 O. Giustolisi and D. A. Savic, “Advances in datadriven analyses and modelling using EPRMOGA,” Journal of Hydroinformatics, vol. 11, no. 34, pp. 225–236, 2009. View at: Publisher Site  Google Scholar
 O. Giustolisi and D. A. Savic, “A symbolic datadriven technique based on evolutionary polynomial regression,” Journal of Hydroinformatics, vol. 8, no. 3, pp. 207–222, 2006. View at: Publisher Site  Google Scholar
 V. Pareto, Cours d’économie Politique, Librairie Droz, Rouge, Lausanne, 1964.
Copyright
Copyright © 2019 Yang Shi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.