A weakly compressible smoothed particle hydrodynamics (WCSPH) method was developed to model open-channel 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 mixed-length-based subparticle scale (SPS) turbulence model is adopted to simulate uniform flow in open channel, and this model is compared with the Smagorinsky-based 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 Smagorinsky-based 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 semi-implicit (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 fluid-structure interaction [35], underwater explosion [6, 7], local scour [8, 9], other transport phenomena [10], and related applications [1113]. Despite advances in SPH methodology, the SPH method has not yet been widely adopted for the simulation of open-channel flows in natural rivers and artificial channels. Thus, the focus of this study was to address the use of SPH for simulation of open-channel flows.

Boundary conditions can directly affect the calculation accuracy and efficiency for the accurate simulation of open-channel 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 [1517] 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 open-channel flow simulation [1822]. In some works, the bottom boundary is roughly treated as slip or no-slip wall [1820]. López et al. [21] proposed that uneven distribution of Lennard-Jones 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 triangle-type, trapezoid-type, and sinusoid-type 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 open-channel 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, time-averaged pressure, Reynolds stress, and turbulent energy. The slip/no-slip treatment on the bottom boundary directly affects the overall accuracy of the simulation results for the open-channel 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 open-channel flow. One approach is to express the effect of the rough bed on the fluid by adding a source term in the fluid-governing equation [24, 25]. For example, Khayyer and Gotoh [24] investigated the effect of bed friction on a dam-break flow by applying a frictional drag force between a near-bed 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 drag-based formulation to account for the large-scale roughness bed surface in depth-limited turbulent open-channel 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 small-sized 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 open-channel 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 open-channel flows. The turbulence phenomena are quite complex, preventing full resolution through full-proof 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 mesh-free 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: large-scale vortex and small-scale vortex. Large-scale vortex can be directly solved using the Navier–Stokes (N-S) equation, and small-scale 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 average-based 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 near-wall particle distribution is chaotic. Recently, De Padova et al. [35] investigated undular hydraulic jumps over small-scale rough bed for bed roughness of 0.02H, where H denotes the water depth. The SPH mixing-length turbulence model exhibits satisfactory performance for simulating hydraulic jump as long as strong turbulent rollers are not present and is applicable for open-channel flow simulation [25, 36]. Overall, there have been insufficient studies of the open-channel flow over rough bed using SPH methods. Therefore, the aim of this work was to further investigate rough-bed boundary treatment and the turbulent model and extend this model to simulate open-channel 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 small-scale rough bed. Next, a numerical open channel is established with a revised turbulence model for open-channel flow in Section 3. In Section 4, turbulent open channel flows over a small-scale 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 ra and rb, respectively; f(ra) and f(rb) 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 ΔVb denotes the volume of particle b. ΔVb = mb/ρb, where mb 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πh2 for the 2D case and q denotes the dimensionless distance between particles.

2.2. Governing Equations

The Reynolds-averaged 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. ua,pa, , and refer to the velocity, pressure, shear stress, and gravity acceleration of particle a, respectively. uab and rab refer to the relative velocity and location between particles a and b and are expressed as uab = uaub and rab = ra – rb. The ∇aWab refers to the gradient of smoothing kernel W(ra − rb, h) for particles a and b.

To close the governing equations (4) and (5), the subparticle-scale (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, Sij denotes an element of the SPS strain tensor, , k denotes the turbulence kinetic energy, CI = 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  = (CsΔl)2|S|, where Cs 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| = (2SijSij)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 lm 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 c0 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 high-frequency 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 delta-SPH 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 inflow domain, internal fluid domain, and outflow domain (Figure 1). The inflow and outflow domains are considered part of the buffer layer to apply appropriate boundary conditions to the internal fluid. The width of inflow and outflow 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 inflow, internal, and outflow domains are called inflow particles, internal fluid particles, and outflow 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 Gomez-Gesteira 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 open-channel flows over a small-scale 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

Open-channel 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, ks 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 open-channel flow over the small-scale 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 near-wall mesh. In the present study, the velocity of the solid-wall 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 near-wall fluid particle is corrected by adjusting the virtual velocity of the solid-walled 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 small-scale rough bed can be adopted. The scale of roughness in the SPH method is a relative concept, where small-scale rough is if , and large-scale rough is if , where is equal to 0.25ks [44]. The wall function method is only applicable to small-scale roughness because the horizontal current velocity in the range does not meet the log-profile 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 cross-sectional velocity distribution is closest to the theoretical flow velocity distribution, the ghost particle velocity corresponds to the equivalent roughness.

2.4. Solution Algorithm

A two-stage predictor-corrector algorithm is used for the numerical integration scheme in this study [45], with the fluid properties of particle a, i.e., velocity (ua), position (ra), and density (ua) 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 fa/ma denotes the fluid force per unit mass exerted on particle a.

3. Simulation of Open-Channel 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 ks = 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 (dp) 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 octa-core Intel® Core™ i7-6700 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.

3.2. Optimization of the Calculation Domain

The fluid boundary layer and open-channel 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 open-channel 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.

3.3. Optimization of Time-Average Operation

The development of open-channel 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.

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 time-average operation is adopted to update the physical quantities of the boundary particles. To determine the time interval of the time-average 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 (dp = 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, ks = 0.00005 m), case II (H = 0.08 m, s = 0.001, ks = 0.01 m), case III (H = 0.08 m, s = 0.002, ks = 0.002 m), and case IV (H = 0.1 m, s = 0.001, ks = 0.01 m), where H denotes water depth, s denotes bed slope, and ks 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.

3.4. Smagorinsky-Based and Mixed-Length-Based SPS Model

The original SPS model [31], known as the Smagorinsky-based SPS turbulence model, has been successfully applied to some coastal applications such as fluid-structure interaction [4, 47, 48]. However, the performance of this model to simulate open-channel flow is not satisfactory. Figure 8 shows the comparison of current velocity and turbulent shear stress with changing water depth for the Smagorinsky-based SPS turbulence model with different Smagorinsky constants (Cs). None of the simulated results were consistent with the theoretical result. This is likely because the tested Cs values are suitable for nonconstant flow problems and should be revised for a case of constant uniform flow.

For cases with Cs value being of similar magnitude, the velocity profiles for different Cs values are almost the same (Figure 8(a)). The value of Cs affects the current velocity distribution of the open-channel flow if the Cs value has a different order of magnitude. In general, the larger the Cs 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 Cs value also has a significant influence on the fluid shear stress. A larger Cs value corresponds to larger fluid shear stress (Figure 8(c)), but even the largest shear stress obtained in the Smagorinsky-based model is less than the theoretical value (Figures 8(d)8(f)).

The CsΔl value shows a certain degree of influence on the flow field distribution of the open-channel flow. Figure 9 shows the effect of particle spacing on the current velocity and shear stress for the Smagorinsky-based SPS turbulence model. A more uniform velocity profile can be obtained by using a relatively larger CsΔl value (Figure 9(a)). The shear stress increases with increasing CsΔ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 Smagorinsky-based turbulence model is attributed to the loss of the turbulent stress by spatial filtering during the long-term 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.

Mayrhofer et al. [49] simulated open-channel flow using a Smagorinsky-based eddy viscosity model and found that the resolution of the particle-based method needed to be at least 16 times higher than that of the grid-based method in order to correctly redistribute the energy between the Reynolds stress components in the SPH-LES frame. This indicates that the open-channel flow cannot be well resolved if the resolution size is not sufficiently dense. For the Smagorinsky-based 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 Smagorinsky-based model in the near-wall region. The original turbulence eddy viscosity νt = (CsΔl)2S was modified to νt = [min (CsΔl, κ·dwall)]2S, where Cs denotes the Smagorinsky constant, Δl denotes the particle spacing, denotes the von Karman constant, dwall 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 mixed-length-based SPS model (equation (7)) is adopted to simulate the open-channel 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.

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 open-channel 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 open-channel 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 open-channel 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.25ks.

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, ks = 0.01 m, dp = 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 near-bottom 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.

Figure 12 shows the effect of particle spacing on the wall function method. The particle spacings shown in Figure 12 are dp1, dp2, and dp3, respectively, where dp1 < dp2 < dp3. 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 near-surface zone (y/h > 0.8). There is only a small relative average error in the mainstream and near-surface 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 near-surface 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.

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 data-mining method, the evolutionary polynomial regression with multiobjective genetic algorithm (MOGA-EPR) 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 open-channel flow. In this study, six potential dimensionless candidate variables were selected in terms of relative water depth (hdp = H/dp), relative bed surface roughness (ksdp = ks/dp), channel bed slope (s), roughness scale (hks = H/ks), dimensionless water depth (), dimensionless particle resolution (), and dimensionless bed surface roughness (). The nonnegative least-squares 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 EPR-MOGA in terms of hydraulic smooth flow, hydraulic rough flow, and hydraulic transition flow and were written aswhere , hdp, ksdp, hks, , 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.

5. Conclusions

In this study, an in-depth study on the simulation of the open-channel 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 Smagorinsky-based SPS turbulence model performs poorly near the rough bed surface in the open-channel flow, showing deviation of the vertical velocity and turbulent shear stress distribution from the theoretical values, but the improved mixed-length-based 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 small-scale 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 near-surface 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 open-channel 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.


The improved WCSPH code used in this paper was developed from the open source SPHysics code developed by Prof. M. Gomez-Gesteira 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 Hydro-Science and Engineering, Tsinghua University, Beijing, China (grant no. 2017-ky-04); 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 open-sourced for a potential researcher. (Supplementary Materials)