Listen

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
Listen

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.

Listen
The governing flow equations including mass conservation and motion in Lagrangian frame for a weakly compressible flow are given by (Shakibaeinia & Jin 2010):
1
where u is velocity vector, t is time, ρ is fluid density, p is pressure, g is gravitational force, r is the position vector and μ 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 Dr/Dt = u. An equation of state is used to calculate the pressure field using the density field.

MPS method

Listen
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 ri, interacts with its neighbouring particle j located within its radius of influence re using a kernel (weight) function, W(rij, re); where rij = |rj-ri| is distance between two particles i and j.
Figure 1
Figure 1. MPS kernel interpolation conceptual model.
View largeDownload slide

MPS kernel interpolation conceptual model.

Figure 1
Figure 1. MPS kernel interpolation conceptual model.
View largeDownload slide

MPS kernel interpolation conceptual model.

Close modal
A dimensionless parameter, namely particle number density, n, is defined as (Shakibaeinia & Jin 2011a):
2
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:
3
where n0 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:
4
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):
5
6
7
where d is the number of space dimensions, eij = rij/rij is the direction vector from particle i to particle j and λ is corrective parameter defined by:
8
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:
9
10
is used in place of pi to stabilize the solution by guaranteeing the positivity of the pressure gradient (Koshizuka et al. 1998). Note that in a uniform-viscosity flow field (e.g., single phase Newtonian flows) μi and μj are equal.
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:
11
where γ= 7 and c0 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.
12
Since an explicit time division is used, the stability condition (CFL conditions) should be satisfied. The CFL conditions are as follows:
13
where Δ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

Listen
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:
14
where α 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:
15
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:
16

Modified boundary conditions

Listen

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

Solid boundaries

Listen
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 re (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.
Figure 2
Figure 2. Particle arrangement near the solid boundary (uθ, uη are tangential and normal velocity components, respectively).
View largeDownload slide

Particle arrangement near the solid boundary (uθ, uη are tangential and normal velocity components, respectively).

Figure 2
Figure 2. Particle arrangement near the solid boundary (uθ, uη are tangential and normal velocity components, respectively).
View largeDownload slide

Particle arrangement near the solid boundary (uθ, uη are tangential and normal velocity components, respectively).

Close modal
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:
17
where 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

Listen
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)).
Figure 3
Figure 3. Open boundary implementation method: (a) for outflow boundary, (b) for inflow (original method) and (c) for inflow (modified method of this study).
View largeDownload slide

Open boundary implementation method: (a) for outflow boundary, (b) for inflow (original method) and (c) for inflow (modified method of this study).

Figure 3
Figure 3. Open boundary implementation method: (a) for outflow boundary, (b) for inflow (original method) and (c) for inflow (modified method of this study).
View largeDownload slide

Open boundary implementation method: (a) for outflow boundary, (b) for inflow (original method) and (c) for inflow (modified method of this study).

Close modal
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:
18
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

Listen
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):
19
A 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

Listen

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

  1. Input initial particle positions ri and assignment of field variables (e.g., ui and pi).

  2. Time integration:

    • Prediction of the velocities u* from (14) and calculation of r* and 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 (rim+1, uim+1, pim+1) to the output; preparation for next time step (t+ Δt).

  3. Repeating steps 2 for the next time step.

Listen
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 Lu, Ld, spillway height Hp, water level above the spillway crest Hs, and design head Hd (summation of Hs 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 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 Hs/Hd.
Table 1

Geometrical and hydraulic characteristics of the study problems

Geometryq (m2/s)Hs (m)Hs/HdHp (m)Hd (m)Lu (m)Ld (m)h/lExperiments 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 – 
Geometryq (m2/s)Hs (m)Hs/HdHp (m)Hd (m)Lu (m)Ld (m)h/lExperiments 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 – 
View Large
Figure 4
Figure 4. Geometries of the study problems: (a) type-I standard ogee-crested spillway (Case I and Case II); (b) symmetric curved bed in an open channel (Case II).
View largeDownload slide

Geometries of the study problems: (a) type-I standard ogee-crested spillway (Case I and Case II); (b) symmetric curved bed in an open channel (Case II).

Figure 4
Figure 4. Geometries of the study problems: (a) type-I standard ogee-crested spillway (Case I and Case II); (b) symmetric curved bed in an open channel (Case II).
View largeDownload slide

Geometries of the study problems: (a) type-I standard ogee-crested spillway (Case I and Case II); (b) symmetric curved bed in an open channel (Case II).

Close modal
Listen

Case I

Listen
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.
Figure 5
Figure 5. Initial position of particles for ogee spillway of Case I (fluid-type, ghost-type and wall-type particles).
View largeDownload slide

Initial position of particles for ogee spillway of Case I (fluid-type, ghost-type and wall-type particles).

Figure 5
Figure 5. Initial position of particles for ogee spillway of Case I (fluid-type, ghost-type and wall-type particles).
View largeDownload slide

Initial position of particles for ogee spillway of Case I (fluid-type, ghost-type and wall-type particles).

Close modal

Effect of boundary modifications

Listen
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, Lu, 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 6
Figure 6. Snapshot of simulated pressure field (at t = 0.95sec) for models: (a) without transition region (Shakibaeinia & Jin 2010) and (b) modified inflow boundary method with transition region. (c) Time history of near-inflow pressure at point (x, y) = (−0.2, −0.05).
View largeDownload slide

Snapshot of simulated pressure field (at t= 0.95sec) for models: (a) without transition region (Shakibaeinia & Jin 2010) and (b) modified inflow boundary method with transition region. (c) Time history of near-inflow pressure at point (x, y) = (−0.2, −0.05).

Figure 6
Figure 6. Snapshot of simulated pressure field (at t = 0.95sec) for models: (a) without transition region (Shakibaeinia & Jin 2010) and (b) modified inflow boundary method with transition region. (c) Time history of near-inflow pressure at point (x, y) = (−0.2, −0.05).
View largeDownload slide

Snapshot of simulated pressure field (at t= 0.95sec) for models: (a) without transition region (Shakibaeinia & Jin 2010) and (b) modified inflow boundary method with transition region. (c) Time history of near-inflow pressure at point (x, y) = (−0.2, −0.05).

Close modal
Figure 7 shows the experimental (Kim & Park 2005) and numerical stream-wise velocity (normalized using its maximum value umax) at the spillway crest (x= 0) for two hydraulic conditions (Hs/Hd = 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 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 (Hs/Hd = 0.50, 1.00 and 1.33) from 0.124, 0.063 and 0.059 to 0.052, 0.036 and 0.015, respectively.
Figure 7
Figure 7. Experimental (Kim & Park 2005) and WC-MPS (with and without proposed modifications) normalized stream-wise velocity at the spillway crest (x = 0) for (a) Hs/Hd = 1.33 and (b) Hs/Hd = 1.00.
View largeDownload slide

Experimental (Kim & Park 2005) and WC-MPS (with and without proposed modifications) normalized stream-wise velocity at the spillway crest (x= 0) for (a) Hs/Hd = 1.33 and (b) Hs/Hd = 1.00.

Figure 7
Figure 7. Experimental (Kim & Park 2005) and WC-MPS (with and without proposed modifications) normalized stream-wise velocity at the spillway crest (x = 0) for (a) Hs/Hd = 1.33 and (b) Hs/Hd = 1.00.
View largeDownload slide

Experimental (Kim & Park 2005) and WC-MPS (with and without proposed modifications) normalized stream-wise velocity at the spillway crest (x= 0) for (a) Hs/Hd = 1.33 and (b) Hs/Hd = 1.00.

Close modal
Figure 8
Figure 8. Simulated (WC-MPS) versus experimentally measured velocity values at x = 0, (a) with and (b) without the modified near-boundary velocity.
View largeDownload slide

Simulated (WC-MPS) versus experimentally measured velocity values at x= 0, (a) with and (b) without the modified near-boundary velocity.

Figure 8
Figure 8. Simulated (WC-MPS) versus experimentally measured velocity values at x = 0, (a) with and (b) without the modified near-boundary velocity.
View largeDownload slide

Simulated (WC-MPS) versus experimentally measured velocity values at x= 0, (a) with and (b) without the modified near-boundary velocity.

Close modal

Simulation results

Listen
Figure 9 shows the time evolution of the simulated flow over Case I spillway for Hs/Hd = 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.
Figure 9
Figure 9. Time history of simulated flow simulated flow over 45° spillway (Case I) for Hs/Hd = 0.5.
View largeDownload slide

Time history of simulated flow simulated flow over 45° spillway (Case I) for Hs/Hd = 0.5.

Figure 9
Figure 9. Time history of simulated flow simulated flow over 45° spillway (Case I) for Hs/Hd = 0.5.
View largeDownload slide

Time history of simulated flow simulated flow over 45° spillway (Case I) for Hs/Hd = 0.5.

Close modal
Figure 10
Figure 10. Simulated velocity vectors after 4.5 seconds (Hs/Hd = 0.5).
View largeDownload slide

Simulated velocity vectors after 4.5 seconds (Hs/Hd = 0.5).

Figure 10
Figure 10. Simulated velocity vectors after 4.5 seconds (Hs/Hd = 0.5).
View largeDownload slide

Simulated velocity vectors after 4.5 seconds (Hs/Hd = 0.5).

Close modal
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 Hs/Hd = 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 Hd) on the spillway for Hs/Hd = 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.
Figure 11
Figure 11. Comparison of experimental (Chow 1988) and simulated (WC-MPS) free-surface profiles for Hs/Hd = 0.50, 1.00 and 1.33: (a) surface profiles; (b) simulated versus measured elevations.
View largeDownload slide

Comparison of experimental (Chow 1988) and simulated (WC-MPS) free-surface profiles for Hs/Hd = 0.50, 1.00 and 1.33: (a) surface profiles; (b) simulated versus measured elevations.

Figure 11
Figure 11. Comparison of experimental (Chow 1988) and simulated (WC-MPS) free-surface profiles for Hs/Hd = 0.50, 1.00 and 1.33: (a) surface profiles; (b) simulated versus measured elevations.
View largeDownload slide

Comparison of experimental (Chow 1988) and simulated (WC-MPS) free-surface profiles for Hs/Hd = 0.50, 1.00 and 1.33: (a) surface profiles; (b) simulated versus measured elevations.

Close modal
Figure 12
Figure 12. Comparison of experimental and simulated (WC-MPS) pressure distribution over the weir crest for Hs/Hd = 0.5 and 1.0: (a) pressure distribution; (b) simulated versus measured values.
View largeDownload slide

Comparison of experimental and simulated (WC-MPS) pressure distribution over the weir crest for Hs/Hd = 0.5 and 1.0: (a) pressure distribution; (b) simulated versus measured values.

Figure 12
Figure 12. Comparison of experimental and simulated (WC-MPS) pressure distribution over the weir crest for Hs/Hd = 0.5 and 1.0: (a) pressure distribution; (b) simulated versus measured values.
View largeDownload slide

Comparison of experimental and simulated (WC-MPS) pressure distribution over the weir crest for Hs/Hd = 0.5 and 1.0: (a) pressure distribution; (b) simulated versus measured values.

Close modal
Discharge coefficient, Cd, is an important design parameter of ogee crest spillways, and defined as:
20
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%.
Table 2

Case I numerical (WC-MPS) and experimental (Chow 1988) discharge ratio

 Discharge ratio, Cd
 
Hs/HdNumericalExperimentalRelative error (%)
0.5 0.463 – – 
1.0 0.5047 0.501 0.435 
1.33 0.5252 0.523 0.418 
 Discharge ratio, Cd
 
Hs/HdNumericalExperimentalRelative error (%)
0.5 0.463 – – 
1.0 0.5047 0.501 0.435 
1.33 0.5252 0.523 0.418 
View Large

Case II

Listen
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 Hs/Hd = 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 Hs/Hd = 0.75, 1.0 and 1.5, respectively.
Figure 13
Figure 13. Comparison of experimental (Chatila & Tabbara 2004) and simulated (WC-MPS) free-surface profiles: (a) Hs/Hd = 0.75, (b) Hs/Hd = 1.0 and (c) Hs/Hd = 1.5.
View largeDownload slide

Comparison of experimental (Chatila & Tabbara 2004) and simulated (WC-MPS) free-surface profiles: (a) Hs/Hd= 0.75, (b) Hs/Hd= 1.0 and (c) Hs/Hd= 1.5.

Figure 13
Figure 13. Comparison of experimental (Chatila & Tabbara 2004) and simulated (WC-MPS) free-surface profiles: (a) Hs/Hd = 0.75, (b) Hs/Hd = 1.0 and (c) Hs/Hd = 1.5.
View largeDownload slide

Comparison of experimental (Chatila & Tabbara 2004) and simulated (WC-MPS) free-surface profiles: (a) Hs/Hd= 0.75, (b) Hs/Hd= 1.0 and (c) Hs/Hd= 1.5.

Close modal

Case III

Listen

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 Hs 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 (Hs = 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 14
Figure 14. Numerical and experimental water surface profile for 1/200 and 1/100 resolutions for (a) Hs = 0.15 and (b) Hs = 0.07 m.
View largeDownload slide

Numerical and experimental water surface profile for 1/200 and 1/100 resolutions for (a) Hs = 0.15 and (b) Hs = 0.07 m.

Figure 14
Figure 14. Numerical and experimental water surface profile for 1/200 and 1/100 resolutions for (a) Hs = 0.15 and (b) Hs = 0.07 m.
View largeDownload slide

Numerical and experimental water surface profile for 1/200 and 1/100 resolutions for (a) Hs = 0.15 and (b) Hs = 0.07 m.

Close modal
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.
Figure 15
Figure 15. Time history of simulated flow over curved bed for Hs = 0.07 m.
View largeDownload slide

Time history of simulated flow over curved bed for Hs = 0.07 m.

Figure 15
Figure 15. Time history of simulated flow over curved bed for Hs = 0.07 m.
View largeDownload slide

Time history of simulated flow over curved bed for Hs = 0.07 m.

Close modal
Figure 16
Figure 16. Snapshot of simulated (a) pressure field and (b) velocity vectors at t = 4sec for Hs = 0.07m.
View largeDownload slide

Snapshot of simulated (a) pressure field and (b) velocity vectors at t = 4sec for Hs = 0.07m.

Figure 16
Figure 16. Snapshot of simulated (a) pressure field and (b) velocity vectors at t = 4sec for Hs = 0.07m.
View largeDownload slide

Snapshot of simulated (a) pressure field and (b) velocity vectors at t = 4sec for Hs = 0.07m.

Close modal
Figure 17
Figure 17. Numerical and experimental pressure head distribution (on curved bed) for (a) Hs = 0.15 and (b) Hs = 0.07 m.
View largeDownload slide

Numerical and experimental pressure head distribution (on curved bed) for (a) Hs = 0.15 and (b) Hs = 0.07 m.

Figure 17
Figure 17. Numerical and experimental pressure head distribution (on curved bed) for (a) Hs = 0.15 and (b) Hs = 0.07 m.
View largeDownload slide

Numerical and experimental pressure head distribution (on curved bed) for (a) Hs = 0.15 and (b) Hs = 0.07 m.

Close modal
Listen

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.

Ataei-Ashtiani
B.
Farhadi
L.
2006
A stable moving–particle semi-implicit method for free surface flows
.
Fluid Dynamics Res.
38
,
241
–
256
.
Google Scholar
Crossref
Search ADS
 
Bhajantri
M. R.
Eldho
T. I.
Deolalikar
P. B.
2007
Modeling hydrodynamic flow over spillway using weakly compressible flow equations
.
J. Hydraul. Res.
45
(
6
),
844
–
852
.
Google Scholar
Crossref
Search ADS
 
Chanel
P. G.
Doering
J. C.
2008
Assessment of spillway modeling using computational fluid dynamics
.
Can. J. Civil Eng.
35
(
12
),
1481
–
1485
.
Google Scholar
Crossref
Search ADS
 
Chatila
J.
Tabbara
M.
2004
Computational modeling of flow over an ogee spillway
.
Comput. Struct.
82
,
1805
–
1812
.
Google Scholar
Crossref
Search ADS
 
Chow
V. T.
1988
Open-Channel Hydraulics
.
McGraw-Hill
,
New York
, pp.
365
–
380
.
Google Scholar
 
Delestre
O.
Lucas
C.
Ksinant
P. A.
Darboux
F.
Laguerre
C.
Vo
T. N.
Cordier
S.
2013
SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies
.
Int. J. Numer. Meth. Fluids
72
(
3
),
269
–
300
.
Google Scholar
Crossref
Search ADS
 
Gingold
R. A.
Monaghan
J. J.
1977
Smoothed particle hydrodynamics: theory and application to non-spherical stars
.
Monthly notices of the Royal Astronomical Society
181
(
3
),
375
–
389
.
Google Scholar
Crossref
Search ADS
 
Gotoh
H.
Sakai
T.
2006
Key issues in the particle method for computation of wave breaking
.
Coastal Eng. J.
53
(
2–3
),
171
–
179
.
Google Scholar
Crossref
Search ADS
 
Hirt
C. W.
Nichols
B. D.
1981
Volume of fluid (VOF) method for the dynamics of free boundaries
.
J. Comput. Phys.
39
(
1
),
201
–
225
.
Google Scholar
Crossref
Search ADS
 
Jabbari Sahebari
A.
Jin
Y. C.
Shakibaeinia
A.
2011
Flow over sills by the MPS meshfree particle method
.
J. Hydraul. Res.
49
(
5
),
649
–
656
.
Google Scholar
Crossref
Search ADS
 
Kim
D. G.
Park
J. H.
2005
Analysis of flow structure over ogee-spillway in consideration of scale and roughness effects by using CFD model
.
KSCE Journal of Civil Engineering
9
(
2
),
161
–
169
.
Google Scholar
Crossref
Search ADS
 
Kirkgoz
M. S.
Akoz
M. S.
Oner
A. A.
2009
Numerical modeling of flow over a chute spillway
.
J. Hydraul. Res.
47
(
6
),
790
–
797
.
Google Scholar
Crossref
Search ADS
 
Kondo
M.
Koshizuka
S.
2011
Improvement of stability in moving particle semi-implicit method
.
Int. J. Numer. Meth. Fluids
65
(
6
),
638
–
654
.
Google Scholar
Crossref
Search ADS
 
Koshizuka
S.
Oka
Y.
1996
Moving particle semi-implicit method for fragmentation of incompressible fluid
.
Nuclear Sci. Eng.
123
(
3
),
421
–
434
.
Google Scholar
Crossref
Search ADS
 
Koshizuka
S.
Tamako
H.
Oka
Y.
1995
A particle method for incompressible viscous flow with fluid fragmentation
.
Comput. Fluid Dyn. J.
4
(
1
),
29
–
46
.
Google Scholar
 
Koshizuka
S.
Nobe
A.
Oka
Y.
1998
Numerical analysis of breaking waves using the moving particle semi-implicit method
.
Int. J. Numer. Meth. Fluids
26
(
7
),
751
–
769
.
Google Scholar
Crossref
Search ADS
 
Savage
B. M.
Johnson
M. C.
2001
Flow over ogee spillway: physical and numerical model case study
.
J. Hydraul. Eng.
127
(
8
),
640
–
649
.
Google Scholar
Crossref
Search ADS
 
Shakibaeinia
A.
Jin
Y. C.
2010
A weakly compressible MPS method for simulation of open-boundary free-surface flow
.
Int. J. Numer. Meth. Fluids
63
(
10
),
1208
–
1232
.
Google Scholar
 
Shakibaeinia
A.
Jin
Y. C.
2011a
MPS-based meshfree particle method for modeling open-channel flows
.
J. Hydraul. Eng.
137
(
11
),
1375
–
1384
.
Google Scholar
Crossref
Search ADS
 
Shakibaeinia
A.
Jin
Y. C.
2011b
A meshfree particle model for simulation of mobile-bed dam break
.
Adv. Water Resour.
34
(
6
),
794
–
807
.
Google Scholar
Crossref
Search ADS
 
Shakibaeinia
A.
Jin
Y. C.
2012
Lagrangian multiphase modeling of sand discharge into still water
.
Adv. Water Resour.
48
,
55
–
67
.
Google Scholar
Crossref
Search ADS
 
Sivakumaran
N. S.
Tingsanchali
T.
Hosking
R. J.
1983
Steady shallow flow over curved beds
.
J. Fluid Mech.
128
,
469
–
487
.
Google Scholar
Crossref
Search ADS
 
United States Army Corps of Engineers (USACE)
1995
Hydraulic design of spillway. Technical Engineering and Design Guides, No. 12, ASCE
.
Yazdi
J.
Sarkardeh
H.
Azamathulla
H. M.
Ghani
A. A.
2010
3D Simulation of flow around a single spur dike with free-surface flow
.
Intl. J. River Basin Manage.
8
(
1
),
55
–
62
.
Google Scholar
Crossref
Search ADS
 
© IWA Publishing 2016