A meshless Lagrangian (particle) method based on the weakly compressible moving particle semi-implicit formulation (WC-MPS) is developed and analysed for simulation of flow over spillways. To improve the accuracy of the model for pressure and velocity calculation, some modifications are proposed and evaluated for the inflow and wall boundary conditions implementation methods. The final model is applied for simulation of flow over the 45° and 60° ogee spillways (with different inflow rates) and also shallow flow over a spillway-like curved bed channel. To evaluate the model, the numerical results of free surface profile and velocity and pressure field are compared with the available experimental measurements. Comparisons show the results’ accuracy of the developed model and proposed improvements. The results of this study will not only provide a reliable numerical tool for modelling of flow over spillways, but also provide an insight for better understating flow pattern over these hydraulic structures.

## INTRODUCTION

Spillways are important structures in dams and are used for controlled release of flow from the upstream reservoir to the downstream channel. Proper design of spillways is a key factor for safe passage of floods from the reservoir pool. The type selection of the spillway depends on design discharge, location of dam and its structure. The ogee-crest spillways, due to their effective and safe passage of flow and proper performance in a wide range of discharges, are the most commonly applied spillways. These spillways transfer the upstream excess water to downstream through a smooth slope following an ogee curve crest (shaped to conform to the lower nappe of a water sheet flowing over an aerated sharp crest). The shape of these spillways has been determined by inspiration from trajectory of flow over the sharp-crest rectangular weirs. One of the effective factors in design and proper performance of a spillway structure is formation of proper flow pattern in different parts of the structure. The United States Army Corps of Engineers Waterways Experiment Station (USACE 1995) studied the behaviour of water flow over spillways and developed a series of design charts. Physical models have been used extensively (e.g*.*, Savage & Johnson 2001; Chatila & Tabbara 2004). However, the physical models are known to be costly and associated with errors due to the scale (ratio of the prototype to model) effects.

Alternatively, numerical modelling can be a reliable and cost-effective tool for understanding of the complexity of hydraulic phenomena, such as flow over spillways. Over the years, the numerical modelling of free-surface flows, like flow over spillways, has mainly been on the basis of the Eulerian mesh-based formulations such as finite element (FE) and finite volume (FV) (e.g*.*, Chatila & Tabbara 2004; Kim & Park 2005; Bhajantri *et al.* 2007; Chanel & Doering 2008; Kirkgoz *et al.* 2009; Yazdi *et al.* 2010). Mesh-based methods have difficulties dealing with applications with large interfacial deformations (due to mesh adaptability and connectivity). Even using an advanced interface tracking scheme such as volume-of-fluid (VOF) methods (Hirt & Nichols 1981), the mesh-based methods have a problem in maintaining sharp interfaces and are associated with numerical diffusions (Koshizuka *et al.* 1998; Shakibaeinia & Jin 2011a).

A newer generation of numerical methods, namely, meshless particle (Lagrangian) methods, such as smoothed-particle hydrodynamics (SPH) (Gingold & Monaghan 1977) and moving particle semi-implicit (MPS) (Koshizuka & Oka 1996) methods, has provided a promising opportunity for modelling of highly deformative flows. Through the past few years, MPS has proved its capabilities for a number of hydraulic engineering problems, such as dam break (Koshizuka *et al.* 1995, 1998; Shakibaeinia & Jin 2010, 2011b), hydraulic jump formation (Shakibaeinia & Jin 2010; Nazari *et al.* 2012), flow over sills and trenches (Shakibaeinia & Jin 2011a; Jabbari Sahebari *et al.* 2011) and multiphase flows (Shakibaeinia & Jin 2012; Liu *et al.* 2005).

Ataei-Ashtiani & Farhadi (2006) investigated the effects of different kernel functions on stability and accuracy of the MPS. The artificial pressure fluctuation in the MPS was addressed by Khayyer & Gotoh (2009), Shakibaeinia & Jin (2010) and Kondo & Koshizuka (2011) and some modification in pressure-gradient calculation was introduced. Shakibaeinia & Jin (2010) proposed the weakly compressible MPS method (WC-MPS) for modelling of incompressible fluids and showed that this method not only improved the unphysical fluctuations, but also slightly increased the model efficiency compared to standard MPS method. They also proposed a scheme for implementation of the inflow/outflow boundaries in MPS. Shakibaeinia & Jin (2011a) and Jabbari Sahebari *et al.* (2011) evaluated the WC-MPS model for various open-channel problems.

This study aims to investigate the application of the WC-MPS method for simulation of flow characteristics (i.e., velocity, pressure and free surface) in spillways (with the focus on ogee-crest spillways). Furthermore, to improve the accuracy of the method, this study will propose some modifications in solid and inflow boundary implementation methods and evaluate their effect on the simulation results. The capabilities of the model for accurate prediction velocity and pressure fields and free surface profile are comprehensively evaluated by comparing the results with the experimental measurement in a broad range of flow and geometrical conditions.

## GOVERNING EQUATIONS

**is velocity vector,**

*u**t*is time,

*ρ*is fluid density,

*p*is pressure,

**is gravitational force,**

*g***is the position vector and**

*r**μ*is fluid viscosity. Note that in Lagrangian frame there is no convective acceleration term in the mass and momentum conservation and the motion of particle is simply calculated by D

**/D**

*r**t*=

**. An equation of state is used to calculate the pressure field using the density field.**

*u*### MPS method

*i*with the position vector

*r*

_{i}, interacts with its neighbouring particle

*j*located within its radius of influence

*r*

_{e}using a kernel (weight) function,

*W*(

*r*,

_{ij}*r*); where

_{e}*r*= |

_{ij}

*r**-*

_{j}

*r**| is distance between two particles*

_{i}*i*and

*j*.

*n*, is defined as (Shakibaeinia & Jin 2011a):This parameter is a measure for the presence of particles around a certain particle. Assuming equal mass

*m*for all particles,

*n*is related to the fluid density by:where

*n*

_{0}is the average initial particle number density, operator is kernel approximation. The kernel function used in this study is a third-order polynomial function, proposed by Shakibaeinia & Jin (2010) as:The basic MPS approximation for the spatial operations: interpolation, gradient, divergence, and vector Laplacian of a scalar

*ϕ*and vector

**are, respectively, defined by (Koshizuka & Oka 1996):where**

*ψ**d*is the number of space dimensions,

*e**=*

_{ij}

*r**is the direction vector from particle*

_{ij}/r_{ij}*i*to particle

*j*and λ is corrective parameter defined by:These MPS spatial integration formula can be used to approximate the spatial derivatives in governing equations of flow. The MPS approximation of the momentum equation then can be expressed by: is used in place of

*p*to stabilize the solution by guaranteeing the positivity of the pressure gradient (Koshizuka

_{i}*et al.*1998). Note that in a uniform-viscosity flow field (e.g., single phase Newtonian flows)

*μ*and

_{i}*μ*are equal.

_{j}*γ*

*=*7 and

*c*

_{0}is a numerical sound speed. Since using real sound speed of the fluid (i.e., water) requires very small time intervals, a value smaller than the actual sound speed is selected based on the maximum acceptable compressibility (density variation) for the application case. For instance, to keep variations of fluid density (Δ

*ρ*/

*ρ*) less than 1% of the reference density, a sound speed larger than ten times of the maximum fluid velocity, |

*u*|

*, must be selected.Since an explicit time division is used, the stability condition (CFL conditions) should be satisfied. The CFL conditions are as follows:where Δ*

_{max}*l*is initial particle distance (particle size) and 0 <

*C*≤ 1 is the Courant number. The sub-particle-scale (SPS) large eddy simulation (LES) turbulence closure (Gotoh & Sakai 2006) is used to model the effect of unresolved small-scale (smaller than particle size) turbulence fluctuations. Further details of SPS model implication for WC-MPS are given in Shakibaeinia & Jin (2011a).

### Time integration

*u**, and a correction velocity

*u*′. The velocity of the prediction step is calculated using the viscous terms, body forces and part of the pressure gradient force as:where

*α*is a relaxation factor and

*m*denotes the previous time step. This velocity is then used to calculate the predicted particle positions

*** and particle number**

*r**n** density and the pressure in new time step (m + 1) as:The velocity correction is calculated from the rest of the pressure gradient term in the momentum equation. The particle velocity and position in the new time step (

*m*+ 1) is then calculated by:

### Modified boundary conditions

This section presents the methods (and proposed modifications) of implementing the free-surface, solid and open boundaries.

#### Solid boundaries

*et al.*1995). These particles are arranged within a distance covered by interaction radius

*r*(Figure 2). Boundary values are enforced to these ghost particles. The normal velocity of ghost particles is set to zero. For the free-slip boundary conditions, the tangential velocity of a ghost particle is equal to that of corresponding fluid particle, while in no-slip conditions, the velocity of a ghost particle is the opposite of fluid particle tangential velocity (Figure 2). Pressure is calculated for the wall particles and is extrapolated to other ghost particle layers to repel the fluid particles from the boundary and to avoid particle penetration of the solid boundary.

_{e}*u*is the friction velocity (shear velocity),

_{τ}*C*is a constant (≈0.5 for smoothed bed),

^{+}*y*is the distance from the solid boundary and

*κ*is the Von Kármán constant (≈0.41). The wall law velocity is applied to the near wall particle (

*y*< Δ

*l*).

#### Open boundaries

*x*-direction inflow, in each depth

*y,*a new particle is added at each

*k*time step based on the inflow velocity

*u*(

*y*) as:The algorithm proposed by Shakibaeinia & Jin (2010, 2011a) successfully prescribes the inflow velocity and controls the outflow depth. However, it also introduces some unphysical pressure pulses at the inflow. This is due to the sudden adding of the new particles. These particles are reflected by the inflow ghost particles creating the unphysical pressure pulses. Here, in this study, we use an alternative approach which damps the sudden adding of the new particles. In this method, the inflow static ghost particles are replaced by dynamics fluid-type particles in a transition zone. The new particles are added to the outer layer of the transition zone. The particles in the transition zone are prescribed a hydrostatic pressure and a velocity equal to the inflow velocity. They are moved according to this velocity. In this way, the new particles not only are not reflected by the ghost particles (as there is no ghost particle at inflow), but also their sudden adding pressure plus is damped as they are passed through the transition zone (Figure 3).

#### Free surface boundary

### Final algorithm

The algorithm used here for each time step can be summarized as:

Input initial particle positions

*r*and assignment of field variables_{i}**(**e.g*.*,*u*)._{i}and p_{i}Time integration:

Prediction of the velocities

* from (14) and calculation of*u** and*r**n**.Pressure calculation from equation of state.

Calculation of the velocity correction (15).

Calculation of the velocity, particle position (moving particle by their velocity) and particle number density.

Sending the new results (

*r*_{i}^{m+1}**,***u*) to the output; preparation for next time step (_{i}^{m+1}, p_{i}^{m+1}*t*+ Δ*t*).

Repeating steps 2 for the next time step.

## STUDY PROBLEMS

*h*:

*l*) of 1:1 and 1:2, respectively (Figure 4(a)). The third case (case III) includes shallow flow over a symmetric curved bed in an open channel (Figure 4(b)). Table 1 summarizes the geometrical and hydraulic conditions including: the inflow discharge

*q*, upstream and downstream length

*L*,

_{u}*L*, spillway height

_{d}*H*, water level above the spillway crest

_{p}*H*, and design head

_{s}*H*(summation of

_{d}*H*and the velocity head in approach canal for design discharge). The available free surface, pressure, velocity experimental profiles taken from Chow (1988), Kim & Park (2005), Chatila & Tabbara (2004) and Sivakumaran

_{s}*et al.*(1983) are used to validate and evaluate the simulation results. Note that the unit-width discharge values used as inflow in the numerical model are calculated for desired experimental

*H*/

_{s}*H*.

_{d}Geometry . | q (m^{2}/s)
. | H (m)
. _{s} | H
. _{s/}H_{d} | H (m)
. _{p} | H (m)
. _{d} | L (m)
. _{u} | L (m)
. _{d} | h/l
. | Experiments ref. . |
---|---|---|---|---|---|---|---|---|---|

45 ° ogee spillway (Case I) | 0.008 | 0.025 | 0.50 | 0.075 | 0.050 | 0.25 | 0.25 | 1:1 | Chow (1988); Kim & Park (2005) |

0.031 | 0.050 | 1.00 | |||||||

0.039 | 0.067 | 1.33 | |||||||

60 ° ogee spillway (Case II) | 0.015 | 0.038 | 0.75 | 0.38 | 0.051 | 0.60 | 1.0 | 1.73:1 | Chatila & Tabbara (2004) |

0.024 | 0.058 | 1.00 | |||||||

0.047 | 0.077 | 1.50 | |||||||

Curved bed (Case III) | 0.112 | 0.15 | – | 0.20 | – | 2.10 | 1.5 | – | Sivakumaran et al. (1983) |

0.036 | 0.07 | – |

Geometry . | q (m^{2}/s)
. | H (m)
. _{s} | H
. _{s/}H_{d} | H (m)
. _{p} | H (m)
. _{d} | L (m)
. _{u} | L (m)
. _{d} | h/l
. | Experiments ref. . |
---|---|---|---|---|---|---|---|---|---|

45 ° ogee spillway (Case I) | 0.008 | 0.025 | 0.50 | 0.075 | 0.050 | 0.25 | 0.25 | 1:1 | Chow (1988); Kim & Park (2005) |

0.031 | 0.050 | 1.00 | |||||||

0.039 | 0.067 | 1.33 | |||||||

60 ° ogee spillway (Case II) | 0.015 | 0.038 | 0.75 | 0.38 | 0.051 | 0.60 | 1.0 | 1.73:1 | Chatila & Tabbara (2004) |

0.024 | 0.058 | 1.00 | |||||||

0.047 | 0.077 | 1.50 | |||||||

Curved bed (Case III) | 0.112 | 0.15 | – | 0.20 | – | 2.10 | 1.5 | – | Sivakumaran et al. (1983) |

0.036 | 0.07 | – |

## RESULTS AND DISCUSSION

### Case I

*.*, particle size Δ

*l*

*=*0.002 m) was selected. This particle distance was selected in order to have a reasonable number of particles in depth at the shallowest areas (e.g., spillway crest). Sensitivity analysis of the model with respect to the particle size shows that the larger particle size may not illustrate the flow features appropriately, and a smaller size will increase the computational time. At the beginning, the particles were uniformly distributed and the velocity and pressure values were set to zero and hydrostatic pressure, respectively.

#### Effect of boundary modifications

*L*, is not long enough (only about twice the depth) to damp the fluctuations, this can affect the accuracy of the results significantly. Replacing the fixed ghost particles with transition particles helps to damp the pressure pulse created by adding the new particles within the transition region (Figure 6(b) and 6(c)).

_{u}*u*) at the spillway crest (

_{max}*x*

*=*0) for two hydraulic conditions (

*H*/

_{s}*H*= 1.00 and 1.33). The numerical results are achieved with and without applying the logarithmic wall law for the near wall velocities. The numerical results are achieved by spatial and temporal interpolation of particles instantaneous velocities on the equally distributed depths at

_{d}*x*

*=*0. As the results show, the numerical velocity profiles underestimate the experiments when no wall law is applied for the near wall velocity. The numerical profiles become more compatible with the experimental ones when the near wall velocities are modified applying a wall law. Figure 8 shows the simulated versus experimentally measured velocity values at the spillway crest for different hydraulic conditions with and without applying the logarithmic wall law for the near wall velocities (Figure 8(a) and 8(b), respectively). The figure also shows the root-mean-square error (RMSE) values. As the figure shows, the proposed modifications of the near-bed velocities will decrease the RSME for all the hydraulic conditions (

*H*/

_{s}*H*= 0.50, 1.00 and 1.33) from 0.124, 0.063 and 0.059 to 0.052, 0.036 and 0.015, respectively.

_{d}#### Simulation results

*H*= 0.5. As the simulation starts, the flat water (of the initial condition) begins to collapse under the gravity force and flows down slope. As the time goes on a gradually varied flow with a subcritical constant depth flow at upstream, a critical flow at the crest and supercritical flow after the crest down the slope forms. The simulation continues until the relatively steady state condition is reached (around

_{s}/H_{d}*t*

*=*3.5 sec). Figure 10 shows the final (

*t*

*=*3.5) velocity vectors.

*H*= 0.50, 1.00 and 1.33. It shows a good compatibility with the RSME of 0.005, 0.027 and 0.048, respectively. Figure 12 shows the pressure head (normalized using

_{s}/H_{d}*H*) on the spillway for

_{d}*H*= 0.50, 1.00. The numerical values are achieved by spatial and temporal averaging of particles instantaneous pressures on the equally distributed distances. As the results show, the pressure head drops from the hydrostatic pressure of the reservoir to near zero on the spillway in the high velocity low-depth area. The numerical results are in relative agreement with the experiments with a RSME around 0.06. The discrepancies, mostly observed at the upstream of the crest, can be the result of MPS unphysical pressure fluctuations (especially at the reservoir) due to the weak compressibilities or MPS approximation errors.

_{s}/H_{d}*C*, is an important design parameter of ogee crest spillways, and defined as:The experimental and numerical discharge ratios of spillway Case I are provided in Table 2. The result shows that the numerical model has accurately predicted the experimental discharge coefficient with a small relative error of around 0.4%.

_{d}. | Discharge ratio, C_{d}. | . | |
---|---|---|---|

H_{s}/H_{d}
. | Numerical . | Experimental . | Relative error (%) . |

0.5 | 0.463 | – | – |

1.0 | 0.5047 | 0.501 | 0.435 |

1.33 | 0.5252 | 0.523 | 0.418 |

. | Discharge ratio, C_{d}. | . | |
---|---|---|---|

H_{s}/H_{d}
. | Numerical . | Experimental . | Relative error (%) . |

0.5 | 0.463 | – | – |

1.0 | 0.5047 | 0.501 | 0.435 |

1.33 | 0.5252 | 0.523 | 0.418 |

### Case II

*H*

_{s}/

*H*= 0.75, 1.0 and 1.5. The domain is initially represented by equal-distance particles with a resolution of 200 particles per unit length (i.e., particle distance of 0.005 m). Figure 13 compares the experimental (Chatila & Tabbara 2004) and numerical (WC-MPS) free surfaces for Case II for different hydraulic conditions. As the results show the numerical model has predicted the water surface upstream and downstream of the crest, with a good accuracy. The RMSE values of predicted water level are 0.013, 0.017, 0.024 for

_{d}*H*

_{s}/

*H*= 0.75, 1.0 and 1.5, respectively.

_{d}### Case III

The third case (Case III) of this study is a bench mark test case, a shallow flow over a curved bed (acting like a spillway) in an open channel. This case is relatively different from the others as the spillway is curved both at the upstream and downstream. Laboratory experiments were conducted by Sivakumaran *et al.* (1983) in different flow conditions. There are various analytical solutions for this case based on frictionless shallow flow assumption. Here, two discharge values resulting in *H _{s}* of 0.15 and 0.07 m were used as inflow boundary condition. The outflow boundary is set to be critical depth (i.e., the reflection from outflow is brought to zero by removing outflow ghost particles). The simulations continued until the relatively steady state condition is achieved.

*et al.*1983) and an analytical solution (based on the formulation provided by Delestre

*et al.*2013). The particle resolution of 1/100 (particle size of 0.010 m) overestimates the water surface elevation. This is due to the fact that there are not enough particles to represent the depth especially in shallower areas. For instance, on the crest with a critical depth of around 0.05 m (

*H*= 0.07 m), there will be only five particles in depth for the resolution of 1/100. Increasing the number of the resolution by two times improves the compatibility (with the experimental measurements) considerably, although it can increase the computational time by around eight times.

_{s}*Hs*value, respectively. The highest velocities occur downstream of the curvature. The steady state condition is achieved after around 4 seconds. Figure 16 shows the final pressure field and velocity vectors. The pressure at the upstream reservoir is relatively hydrostatic and it decreases as the flow is accelerated on the crest. Figure 17 compares the experimental and numerical pressure head on the curved bed. The pressure drops from a hydrostatic pressure at the upstream to minimum pressure right after the crest and then it increases slightly further downstream. The figure shows a good compatibility of experimental and numerical results with a RMSE value of around 0.03.

## CONCLUSIONS

A meshless Lagrangian (particle) model based on the WC-MPS formulation was developed, analysed and evaluated for the simulation flow over ogee-crest spillways with different geometrical and flow conditions. Considering its meshless nature, the developed model was capable of dealing with free surface deformations/fragmentations without using any surface tracking scheme. The inflow boundary condition implementation method was modified by considering a transition region at inflow. This new method improved the near-boundary unphysical pressure fluctuations. Moreover, the implemented standard wall law for near boundary velocity improved the predicted velocity profile. The available experimental measurements for the free surface profile, and pressure and velocity fields were used for validation purposes. The comparisons and error analysis showed a good compatibility of the numerical results with the measurements. The results showed the capabilities of the developed model for accurate modelling flow over ogee-crest spillways with a broad range of conditions. The model of this study can also be applied for other spillway types.