## Abstract

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.

- free surface flow
- meshless particle methods
- ogee-crest spillway
- pressure distribution
- velocity field
- WC-MPS method

## 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

The governing flow equations including mass conservation and motion in Lagrangian frame for a weakly compressible flow are given by (Shakibaeinia & Jin 2010):
1where ** u** is velocity vector,

*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

In MPS the flow domain is represented by a set of particles that possess the flow field variables (such as mass, velocity and pressure) and are free to move in Lagrangian coordinates. The flow governing equations are discretized and solved over these particles. MPS integration method is based on the local kernel interpolation (weighted averaging) of quantities and derivatives over the particle domain. As Figure 1 shows, a target particle *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*.

A dimensionless parameter, namely particle number density, *n*, is defined as (Shakibaeinia & Jin 2011a):
2This 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:
3where *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:
4The 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):
5
6
7where

*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: 8These 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: 9 10 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}The conservation mass (Equation (1)) is automatically satisfied, as mass is represented by the presence of particles. WC-MPS (Shakibaeinia & Jin 2010, 2011a) considers the system as a slightly compressible system, and uses an equation of state to estimate the pressure value for each particle. This study uses the Tait's equation of state (commonly used for high pressure water flow) given by:
11where *γ* *=* 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*|* _{max}*, must be selected.
12Since an explicit time division is used, the stability condition (CFL conditions) should be satisfied. The CFL conditions are as follows:
13where Δ

*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

The time integration method is based on a prediction-correction scheme. Each time step is divided into two pseudo steps of prediction and correction. The velocity in a new time step is the summation of the velocity, calculated in prediction step *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:
14where *α* is a relaxation factor and *m* denotes the previous time step. This velocity is then used to calculate the predicted particle positions ** r*** and particle number

*n** density and the pressure in new time step (m + 1) as: 15The 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: 16

### Modified boundary conditions

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

#### Solid boundaries

In the vicinity of solid boundaries, the particle density decreases. This creates a near-boundary particle deficiency which can cause a problem in the calculation of particle number density and MPS estimation of derivatives. Thus, some layers of so-called ghost particles are placed outside the boundary to prevent this unwanted density reduction (Koshizuka *et al.* 1995). These particles are arranged within a distance covered by interaction radius *r _{e}* (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.

In this study, to better predict the near boundary velocity (in case the particle size is not small enough to resolve the near boundary sub-layer), the standard logarithmic law of wall is applied to the fluid particles in the vicinity of the solid boundaries. As is shown below in the section Case I, this will predict a better near-boundary velocity profile. The logarithmic law of the wall in general form is given by:
17where *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

The method of this study for implementing open-boundaries (inflow and outflow boundaries) is based on a modified version of Shakibaeinia & Jin (2010, 2011a). In their method, some layers of fixed-position ghost particles are defined at the inflow and outflow boundaries to compensate for density deficiency (Figure 3(a) and 3(b)).

At outflow for a known pressure (depth) outflow boundary condition, the boundary pressure is prescribed to the ghost particles. Approaching the first layer of ghost particles, the fluid particles are removed from the domain (Figure 3(a)). The pressure value is transferred to the fluid particle by the repulsive force that has been created by the ghost particles. To control the outflow depth, the height of the ghost particle column is fixed to be the same as the boundary depth. The repulsive force by the ghost particles keeps the fluid depth near the outflow boundary the same as the ghost particle height. At the inflow new particles are added to the distance between the first layer of ghost particles and fluid particles (Figure 3(b)). The boundary condition variable is prescribed to these particles, while the other variables are extrapolated from the fluid particles. For a known velocity inflow boundary condition, particles are added to the inflow based on a rate depending on the velocity distribution at the inflow boundary. For instance, for *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:
18The 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

In the MPS method, particle density is used to track the free water surface. Since there is no particle out of the free water surface, particle density is reduced drastically at the free surface. A particle is recognized as a free surface particle if its particle number density is less than a fraction of the initial average particle number density (Shakibaeinia & Jin 2011a): 19A zero-pressure condition is enforced to the free surface particles. In the MPS method, there is no need to apply further conditions for the free surface.

### 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

Three spillway geometries, each with different hydraulic conditions, are considered for validation and application of the developed model. The first and second cases (cases I, II) are two type-I standard ogee-crested spillways (according to USACE 1995) with straight section slopes (*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}## RESULTS AND DISCUSSION

### Case I

Figure 5 shows the initial position of the particles for Case I (45° ogee spillway). A particle resolution of 500 particles per unit length (i.e*.*, 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

Effects of the proposed modifications in methods of implementing the inflow and solid boundary conditions are investigated here. Figure 6 compares the near inflow pressure field achieved using fixed ghost particles at inflow, as in Shakibaeinia & Jin (2010, 2011a), and using the modified method of this study which uses a transition zone. It also compares the pressure time history for a point near inflow. As Figure 6(a) and 6(c) show, using the fixed ghost particles at inflow creates an unphysical pressure pulse (originated from the inflow) which results in the free surface fluctuations even further from the inflow boundary. As the upstream section of the spillway, *L _{u}*, 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)).

Figure 7 shows the experimental (Kim & Park 2005) and numerical stream-wise velocity (normalized using its maximum value *u _{max}*) at the spillway crest (

*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

Figure 9 shows the time evolution of the simulated flow over Case I spillway for *H _{s}/H_{d}* = 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

*t*

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

*t*

*=*3.5) velocity vectors.

The flow is relatively stagnant at the reservoir and gets accelerated as it is contracted on the crest and driven by the gravity force down the slope. The highest velocities occur on and downstream of the slope, respectively. The figure also shows a relatively hydrostatic pressure at the reservoir, where the flow is slow and stable. The pressure drop flow passes over the crest. Figure 11 compares the experimental and simulated water level for *H _{s}/H_{d}* = 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

*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}Discharge coefficient, *C _{d}*, is an important design parameter of ogee crest spillways, and defined as:
20The 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%.

### Case II

The second case in this study is a 60° ogee-crest spillway (as in Table 1), for which water data of physical modelling are available from Chatila & Tabbara (2004) for three hydraulic conditions *H*_{s}/*H _{d}* = 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

*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.

To evaluate the effect of particle resolution, two different particle resolutions of 1/200 (particle size of 0.005 m) and 1/100 (particle size of 0.010 m) were used. Figure 14 shows the numerical water surface profile with the different particle resolutions. It also compares the results with the experimental measurements (Sivakumaran *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 _{s}* = 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.

Figure 15 illustrates the evolution of particles configuration for flow over the curved bed (particles are coloured by their velocity magnitudes). As the time goes on the depth on and upstream of the crest increase until it reaches to critical depth and *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.

- First received 7 May 2015.
- Accepted in revised form 1 August 2015.

- © IWA Publishing 2016

Sign-up for alerts