Next Article in Journal
Study on the Pricing of Water Rights Transaction between Irrigation Water Users Based on Cooperative Game in China
Previous Article in Journal
Mechanistic Model and Optimization of the Diclofenac Degradation Kinetic for Ozonation Processes Intensification
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Application of Smooth Particle Hydrodynamics to Particular Flow Cases Solved by Saint-Venant Equations

1
Al-Hikma College of Science and Technology, Khartoum, Sudan
2
Banat Water Administration, Romanian Waters, 300173 Timisoara, Romania
3
Faculty of Civil Engineering, Politehnica University of Timisoara, 300223 Timisoara, Romania
4
Department of Hydroinformatics and Socio-Technical Innovation, IHE Delft, 2611 DA Delft, The Netherlands
*
Author to whom correspondence should be addressed.
Water 2021, 13(12), 1671; https://doi.org/10.3390/w13121671
Submission received: 30 April 2021 / Revised: 7 June 2021 / Accepted: 12 June 2021 / Published: 16 June 2021
(This article belongs to the Section Hydraulics and Hydrodynamics)

Abstract

:
Smoothed particle hydrodynamics (SPH) is a Lagrangian mesh free particle method which has been developed and widely applied to different areas in engineering. Recently, the SPH method has also been used to solve the shallow water equations, resulting in (SPH-SWEs) formulations. With the significant developments made, SPH-SWEs provide an accurate computational tool for solving problems of wave propagation, flood inundation, and wet-dry interfaces. Capabilities of the SPH method to solve Saint-Venant equations have been tested using a SPH-SWE code to simulate different hydraulic test cases. Results were compared to other established and commercial hydraulic modelling packages that use Eulerian approaches. The test cases cover non-uniform steady state profiles, wave propagation, and flood inundation cases. The SPH-SWEs simulations provided results that compared well with other established and commercial hydraulic modeling packages. Nevertheless, SPH-SWEs simulations experienced some drawbacks such as loss of inflow water volume of up to 2%, for 2D flood propagation. Simulations were carried out using an open source solver, named SWE-SPHysics.

1. Introduction

Smoothed particle hydrodynamics (SPH) is a mesh-free Lagrangian particle method originally designed for continuum scale applications. The method was first introduced in the late seventies independently by [1,2] for astrophysical problems in the three-dimensional open space. Initially, the method’s main feature uses statistical techniques to describe the physical variables quantities from a known distribution of fluid elements. In the 1990s, the method was extended to the simulation of free surface flow by [3]. The method has been improved throughout the years, and now it is attracting more researchers and has been successfully used in many engineering areas [4,5]. With the development of the SPH method and its application to a wide range of problems in hydraulics [6,7], more attractive features have been demonstrated while some drawbacks have also been identified. Among the advantages listed in [8], SPH can handle very large deformations of a continuum body given that the connectivity between particles describing the continuum are updated every computational time step. Furthermore, the resolution can be changed with particle position and time.
The SPH method obtains approximate numerical solution of the fluid flow equations using the notion of representing the fluid by unconnected and unrestrained particles carrying fluid properties such as mass, density, pressure, velocity, position, etc. The particle properties can change over time due to the interaction with the other surrounding particles using a smoothing function. In SPH method, the function and its spatial derivatives are approximated by the interaction between particles using the so-called kernel and particle approximation. Kernel approximation is obtained by choosing a kernel smoothing function to define the interaction with the surrounding particles and to determine the extent of the influence area of a particle. The kernel smoothing function introduced and published in literature [4] is radial symmetric (gaussian, quadratic, spline, etc.); however, the choice of the smoothing function and the smoothing length has a significant importance in accuracy and speed. Most often, the smoothing function is chosen as a compact function because of the good performance [3,4,5]. The particle approximation is the integration over the nearest neighbor particles, such that the properties variables on a particle are approximated by a summation of the values over the particles within the concern particle support domain. The particle position and the magnitude of the individual fluid properties are updated every computational time step. The particle approximation for functions and their derivatives reduces the set of PDEs to ODEs discretized only with respect to time and solved using explicit numerical schemes [9].
The numerical simulations of different natural phenomena such as flood waves due to dam breaks [10], river flood waves, and tidal flows are important for flood risk analysis [4]. These types of simulations have been conducted using solution of shallow water equations (SWEs) on classical grid-based discretization of the computational domain. Recently, SPH has been used to approximate the SWEs with the so-called SPH-SWEs numerical scheme. The SPH-SWEs scheme shows promising results using different formulations [11,12]. The Lagrangian formulation of the problem where no mesh generation is needed enables the SPH method to describe the wet-dry interface without any special treatment, making SPH particularly suitable for SWEs.
In the SPH-SWEs solution, special consideration is given to aspects such as closed and open boundary conditions, stabilization terms, source terms, and convergence where poor resolution is present. Different methods have been proposed to simulate the closed boundary condition, with each method having its own advantages and drawbacks [13]. An open boundary algorithm is proposed by [14] based on a simplified version of the characteristic method, to simulate both subcritical and supercritical inflow and outflow. In order to avoid numerical oscillation in the presence of shock waves, viscosity is added as a stabilization term to the momentum equation formulation.
A new stabilization term approach based on the idea of Lax–Friedrichs flux is introduced by [11], with the main advantage of having no parameters. This term is to be calibrated in comparison with the artificial viscosity term, which is based on the necessary level of viscosity. As an alternative way to handle shock waves, Riemann solvers are used. A two-shock Riemann solver for the SPH-SWEs scheme is defined by [15]. To reduce the overall numerical viscosity error using both the Lax–Friedrichs and the two-shock Riemann solver, the Monotone Upwind-Centred Scheme for Conservation Laws (MUSCL) is used as shown in [16]. The MUSCL scheme is a non-upwind procedure to reconstruct the viscosity term of the Lax–Friedrichs flux, as well as velocities and water depths in the two-shock Riemann solver. In the momentum equation, the bed gradient source term is addressed in two ways; [17] describes the bed gradient using an analytical function whereas [15] introduces an approach to address irregular bathymetries. In this last technique, the bed is discretized into a new set of interpolation points called bottom particles over which an SPH-based interpolation is performed to calculate the bed gradient and tensor. The method uses two different formulations of the momentum equation to deal with the imbalance associated with the bed discontinuity. The first one obtains a fully conservative formulation using the variational formulation of the SWEs while the second one is based on a non-conservative form of the free-surface gradient. The issue of poor resolution when the flow enters an area with very shallow depths is addressed by using a conservative particle splitting procedure, which improves the method convergence [18]. In this procedure, the resolution is increased by splitting the particle in the region of poor resolution to seven daughter particles given that the mass and the momentum are conserved.
The SPH-SWEs has been validated in literature using different test cases like dam break [9], Tsunami [19], Thacker basin, and flow in a parabolic channel. Except for dam break cases, most of the validation cases were restricted to small scale simulations. Research presented herein examined the applicability of the SPH-SWEs in hydraulic engineering using different practical cases varying from simple non-uniform flow water profile and wave propagation to flood inundation cases with complex geometries in practical scales. The simulations were done using the recent SPH-SWEs improvements. This research includes comparison with semi-analytical solutions and results obtained by using several Eulerian modeling approaches. The simulations were conducted using the SWE-SPHysics solver [20,21,22,23], which incorporates most of the up-to-date improvements. SPHysics is a platform of open source SPH codes to simulate free-surface flows developed by a group of researchers from several universities around the world: the Johns Hopkins University (U.S.A.), the University of Vigo (Spain), the University of Manchester (U.K.), and the University of Rome La Sapienza (Italy). The SWE-SPHysics solver is based on the SPH-SWEs formulations with the option of particle splitting, in which the closed boundary conditions are represented using modified virtual boundary particle method (MVBP), and the open boundary conditions are simulated [14,15]. SWE-SPHysics is an open source and freely downloadable code from SPHysics website, http://www.sphysics.org (accessed on 1 February 2021).
This paper is structured in five parts. After this introduction, the method theoretical formulations are presented followed by the SPH-SWEs models descriptions. Subsequently, the models’ results are analyzed and discussed, and finally, conclusions are drawn.

2. Theoretical Background of the SPH-SWEs Formulations

The SWEs are written in Lagrangian form as defined by Equations (1) and (2):
D d D t = d . v
D v D t = g d + g ( b + S f )
where v is the horizontal depth-averaged velocity vector, d is the water depth, b is the bottom elevation, g is the acceleration due to gravity, and S f is the bed friction source term.
These equations are formulated and implemented as a system of equations of density and momentum.

2.1. Density Formulation

In SPH-SWEs, the density ρ is defined as the mass of fluid per unit area (two-dimensional density), as detailed in Equation (3).
ρ = ρ w d
where ρ w is the constant water density. The density evaluation in SPH-SWEs for a particle i is given by Equation (4) as:
ρ i = j m j W i ( x j , h i )
where m j is the particle density, and W i is the kernel, expressed as a function of the vector particles distance x j and the smoothing length h i . In order to maintain the solution accuracy, SPH-SWEs uses a varying smoothing length h i related to the density, as shown in Equation (5).
h i = h 0 ( ρ 0 ρ i ) 1 / D m
In Equation (5), D m is the number of the space dimensions, and h 0 and ρ 0 are the initial smoothing length and density, respectively. This implicit problem is solved by using a Newton–Raphson iteration [17]. The stopping criteria are the non-dimensional density error threshold and/or the number of iterations.

2.2. Momentum Formulation

The momentum formulation presented herein follows the detailed derivation in [14,15,16,17], where the particle acceleration a i can be calculated by Equation (6) as
a i = g + v i . k i . v i + t i b i 1 + b i b i b i t i
where b i is the bed gradient, k i is the tensor curvature, given by k i = ( b i ) , and t i = T i / m i , the internal force. The expression of T i is given by the Equation (7)
T i = j m i m j g 2 ρ w [ 1 α j W j ( x i , h j ) 1 α i W i ( x j , h i ) ]
The bed friction term S i , j for each fluid particle is determined using Equation (8)
S f , i = v i   g n i 2 | v i | R i
where n i denotes the Manning coefficient, and R i is the hydraulic radius. In the case of very shallow water depth, where R i is too small, the friction term becomes unphysical, and hence, a minimum value of R i is introduced.

2.3. Time Stepping

In order to march the particle in time, particle positions and velocities are integrated in time using an explicit Leap-frog scheme [20], where the time step must satisfy a Courant–Friedrichs–Levy (CFL) condition given by Equation (9)
Δ t C CFL   . min ( h i c i + | v i | )
where c is the wave propagation speed, c = g d .

2.4. Boundary Conditions

The SPH method was initially designed for astrophysics and galaxy simulation, an open space where there are no boundaries; however, for fluid flow simulation, there are either open or closed boundaries. In SPH, the main problem in representation of a boundary is associated with the truncation of the kernel function. Closed boundary is supposed to block and prevent the particles to penetrate through the wall without changing the fluid physical properties. Closed boundaries are presented in detail in [16]. The virtual boundary particle method (VBP) for closed boundary is proposed by [21], which considers the notion of mirrored ghost particles. In the VBP, virtual particles are placed along the closed boundary at a distance from the boundary, which is less than its influence domain. A line of fictious interior particles is generated by using local point symmetry. The new fictitious particles carry the same physical properties of the fluid particle. Modified virtual boundary particle (MVBP) method is introduced by [14] in order to minimize the errors associated with kernel truncation in the original VBP. MVBP provides a stencil that is closer to the interior fluid particles with a complete kernel support. For all test cases presented herein, hydrograph boundary conditions were transformed into velocity using Manning formula for discharge computation.

3. Materials and Methods

In order to test the suitability of the SPH-SWEs approach to solve Saint-Venant equations describing hydraulic problems, a number of theoretical and practical cases were considered: the non-uniform water profiles, wave propagation, and flooding. The first set of test cases are the non-uniform flow water profile (M1 and M2) for which the SPH-SWEs model results were compared with the semi-analytical results using the direct step method. The second set of tests simulate the wave propagation. Results of second tests were analyzed and discussed based on the general understanding of the wave behavior. Finally, two cases of flood inundation over an initially dry bed were conducted, and the results were compared with the results obtained using several Eulerian modeling approaches.
A common setup is adopted for all the SPH-SWEs considered models. The Lax–Friedrichs flux with MUSCL reconstruction is used as a stabilization term. The initial smoothing length is taken to be 1.2 times the initial particles spacing; however, the simulations were marched using variable time step with Courant number equal to 0.4. The inflow and outflow open boundary conditions were chosen as discharge hydrographs. In SWE-SPHysics, particles can only have properties in terms of velocities and water depths; therefore, Manning formula was used to transform the discharge hydrographs into velocities and water depths.

3.1. Non-Uniform Steady State Profiles

Two simple shallow water cases, water surface profiles M1 and M2, for mild slope channels are represented and solved using SPH-SWEs formulation. The results obtained were compared with the ones obtained by applying the semi-analytical direct step method. In order to evaluate the effect of wall boundary in a 2D simulation domain, both 1D and 2D simulations were conducted. Details of the considered models in terms of channel physical properties (length, width, slope, roughness), boundary conditions, and numerical parameters (number of particles and their spacing) are summarized in Table 1. The models are run until the simulation reaches the steady state. The maximum non-dimensional errors were calculated using Equations (10)–(12).
E r r o r   d = m a x ( d d s a d s a )
E r r o r   v x = m a x ( v x v x s a g d s a )
E r r o r   v y = m a x ( v y g d s a )
where the superscript sa is for the result obtained with the semi-analytically solution using the direct step method, and v x and v y are velocities in x and y directions, respectively.

3.2. Wave Propagation

Wave propagation behavior was studied by simulating its behavior in a straight channel. The main objective is to assess the method’s ability to represent the attenuation and translation of the wave.

3.2.1. Wave Attenuation

Wave attenuation takes place when the channel has the capacity to reduce the peak of the wave, as it is translated towards downstream. This is checked in the downstream part of a river, in regions with flat topography. Such a wave is referred to as diffusive wave. To simulate wave attenuation, the considered model consists of a channel with a relatively small slope of 0.001. The channel has a rectangular cross-section of 60 m width with a uniform Manning coefficient of 0.03. An initial discharge of 80 m3/s is applied, and the upstream boundary condition represents an inflow flood hydrograph with a flood peak of 500 m3/s over 2000 s as shown in Figure 1. A constant outflow discharge of 80 m3/s was set as the downstream boundary condition.

3.2.2. Wave Translation

The considered wave translation was the one of kinematic wave when the flood hydrograph travels downstream with its shape unchanged. This phenomenon occurs when the gravitational force is the dominant one as compared to the frictional force; therefore, the wave propagates downstream with almost no attenuation. This situation is seen during flash floods in a steep channel. The SPH-SWEs test model is set-up as a rectangular cross-sectional channel with 60 m width and a steep slope of 0.01. The friction parameter uses a uniform Manning coefficient 0.02. The initial discharge is 150 m3/s, and the inflow hydrograph is applied upstream with 800 m3/s peak discharge over a time base of 1500 s (Figure 1). A constant outflow discharge of 150 m3/s was set as the downstream boundary condition.
Detailed model set-up for the wave propagation tests is given in Table 2.

3.3. Flooding

In 2009, UK Environment Agency (EA) carried out a benchmarking study of several 2D hydraulic modeling packages [18]. Selected studies cover flooding cases with detailed discussion and models’ evaluations and provides an insight of the latest 2D hydraulic modeling tools capabilities. An updated report was published in 2013 [24], and the data are available to use for comparison. Same benchmarking tests were used by [23,24,25,26] to investigate the capability of XBeach software. These cases offer an opportunity to weigh SPH method against the other grid-based methods implemented in the hydraulic modeling packages used in the EA report. Two test cases of the EA report were chosen for the simulation of floods using SPH-SWA: momentum conservation over a hump and filling of floodplain depressions. The comparisons were carried out using InfoWorks RS 2D, Flood Modeller 2D component, and XBeach software. The InfoWorks RS 2D algorithm is based on solving the shallow water equations using finite volume scheme while Flood Modeller 2D and XBeach use different finite difference schemes.

3.3.1. Momentum Conservation over a Hump

The main objective of this test case is to examine SPH-SWEs capability to conserve momentum over a hump. It also shows the method’s ability to simulate advancing of a wave front over an initially dry bed and to handle disconnected water bodies simulation. Momentum conservation is an essential property in the case of sewer or pluvial flooding in urban floodplain areas.
As described in [24], the test topography consists of two depressions separated by a hump, with a longitudinal profile as shown in Figure 2. The domain is initially fully dry where a varying inflow hydrograph is applied at the left boundary. The inflow hydrograph has a peak flow of 65.5 m3/s and a time base of 30 s (Figure 2). The flow travels downhill with a steep slope of 1:200. The total volume of the inflow hydrograph is just sufficient to fill in the left depression, and some of the volume is expected to overtop the hump as a result of the flow inertia. The total simulation time is 15 min to allow the water to settle. The channel friction is represented by a uniform Manning coefficient value of 0.01.
Five cases were simulated with different values of minimum friction depth, inflow velocity, water depth, and particle spacing as presented in Table 3. Particle splitting procedure was applied in Case 4. The domain was initially empty, and therefore, starting to split the particle as it enters the domain causes some of the daughter particles to be located outside. Thus, the splitting process was set to take place at enough distance from the boundary; however, this resulted in having the particle travel without split until it reaches the splitting region.

3.3.2. Filling of Floodplain Depressions

This test investigates the method capability to determine the final inundation extent and the water depth of low momentum flood over a complex geometry. It also assesses the ability to handle wetting and drying of a floodplain and to simulate disconnected water bodies. The test domain is a square area of 2000 × 2000 m2 floodplain with 16 flattened egg-shape depressions of 0.5 m deep. The general slope of 1:1500 is applied from the north to south direction and the one of 1:3000 from the west to east direction, with a 2 m drop in elevation between the top left corner to the bottom right corner. The Digital elevation model (DEM) of this area is shown in Figure 3a. The Manning coefficient for bed friction is 0.03.
The inflow boundary is applied along a 100 m line running south from the north-western corner of the test domain. The hydrograph applied on this boundary has a peak flow of 20 m3/s and a time base of 85 min, as shown in Figure 3b. Except for the inflow boundary, all other boundaries are closed boundaries. The total simulation time is 48 h, set to reach an inundation state over the whole domain. The water level in the middle center point of each depression is observed together with the final inundation extent.
Four cases were carried out to investigate the effect of having different inflow momentum and smaller resolution, as detailed in Table 4.
Total inflow volume of the inflow hydrograph in Figure 3 is 97,200 m3 for the time base of 85 min.

4. Results and Discussion

4.1. Water Profile for Non-Uniform Flow

The results obtained using 1D and 2D SPH-SWEs were compared with the semi-analytical results (direct step method) for both velocities and water levels. In general, there are very small differences towards the downstream end; however, these deviations almost diminish when the water level is at normal depth state at the upstream boundary (see Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8).
The results of the 1D simulation were better than the 2D simulation, showing no disturbances and small errors as presented in Table 5. The small differences between semi-analytical solution and SPH-SWE in the 1D simulation results are due to the coarse resolution used in the model representation (100 m).
In the case of the 2D simulation, there were approximately 2% differences in results both along the longitudinal profile and in the cross-sections. This is due to the fact that the wall boundary in MVBP method is not well represented. The disturbance along the cross-sections created a small velocity component in the transverse direction. As shown in Figure 5 and Figure 8, the 2D simulation results improved significantly when a high resolution (10 m) was used, because the negative effect of the wall boundary was reduced by having smaller particle spacing. The velocity field was found to be more sensitive to the effect of the wall boundary. The cross-sectional view in Figure 6 and Figure 9 shows the noise created by the wall boundary.

4.2. Wave Propagation

For the considered wave propagation cases, the hydrograph peak shows attenuation as the flow travels downstream (Figure 10). The wave attenuation cases are in subcritical conditions where the flow information is also received from downstream. In order to properly assess the propagation, the selected presented results are on locations far from the effect of the downstream boundary.
The typical rating curve loop pattern for the mild slope channel is presented in Figure 11. It shows that in the same station and at the same stage, more flow passes through the channel during the rising stage than during the falling stage of the hydrograph.
In the wave translation case, the considered flow is supercritical, and the information is mainly received from upstream; therefore, the backwater effect of the downstream boundary is limited. The results in Figure 12 show hydrographs, at different positions along the channel, confirming a full translation of the wave downstream. The hydrographs are slightly tilted back due to the fact that the wave celerity is proportional to the water depth, and therefore, the part in the hydrograph with high depth travels faster than the one with low depth. The obtained rating curve in Figure 13 represents the behavior of the kinematic wave. In the kinematic wave, the rating curve loop is very narrow, and therefore, the flow can be approximated using the normal steady flow formulas.

4.3. Flooding

The results obtained for the nine considered flooding cases were compared with InfoWorks RS 2D, Flood Modeller 2D component, and XBeach software. These software results of the considered cases are taken as presented and available in the literature in [24,25,26].
The loss in volume and computer running time results for the nine flooding cases are presented in Table 6.

4.3.1. Momentum Conservation over a Hump

Analysis of SPH-SWEs models’ results shows an overestimation of the wave front velocity, where the water arrived at the center of first depression earlier by ~10 s. As a result of faster advancing wave front, water spills out of the depression. Therefore, there is a small underestimation of the water level at the center of the first depression and a small overestimation of water levels at the center of second depression with the water arriving earlier by ~30 s, as compared to the other available models in the literature (see Figure 14 and Figure 15). in the transverse direction of the domain, a small oscillatory behavior was observed which is due to the wall boundary representation. A small amount of particles penetrated through the boundary wall out of the domain.
Case 1 and Case 2 showed a limited effect of the minimum water depth when considering the friction parameter. This is due to the high values of velocity at the open boundary, and therefore, the velocity becomes a dominant factor in the friction term calculation. Despite the loss in the total inflow volume shown in Table 6, the results of the considered cases overestimated the water level at the center of second depression. In Case 3, due to smaller spacing, more particles were used, which shows a loss in the inflow water volume of only 0.5%, which resulted in higher water level at the center of second depression. The particle splitting method was used in Case 4 using the same velocity of Cases 1, 2, and 3. Case 4 presented better results particularly at the center of second depression, but after t ≈ 300 s, the water level decreased as the split particles close to the wall boundary started to go out through the wall boundary and leaved the domain; this is attributed to the poor wall boundary representation. The open boundary parameters, velocity, and water depth may play a major role in shaping the output results. In this test case, the domain of the first depression is not long enough to influence the importance of the particles velocity at the inflow boundary. This is noticed in Case 5, where the maximum velocity at the boundary was reduced to 2.0 m/s using the same resolution as for Case 1 and Case 2. Results show less water overtopping the hump for the second depression. Case 5 provided the best results in comparison to the other 2D hydraulic model’s results.
In the case of a finite inflow hydrograph, the exact volume of inflow water is difficult to achieve as it depends on the open particles spacing, the shape of the hydrograph, and the inflow velocity. Large open particle spacing, or steep hydrographs give high probability not to achieve the exact inflow volume, which can be noticed in Table 6 where only 0.5% loss in volume was observed when using a high particle resolution. Table 6 also shows that SPH-SWEs method is computationally expensive because it requires long running times, when run on a CPU. When particle splitting method is used, simulation time is longer because solution engine contains more loops to calculate the physical properties of the daughter particles.

4.3.2. Filling of Floodplain Depressions

The SPH-SWEs cases show similar results to those obtained by the other software modelling tools (see Figure 16 and Figure 17). In general, SPH-SWEs cases’ results are the same in the region of high momentum close to the inflow location; however, significant underestimation of the water level is noticed in the depressions far from the inflow location. Arrival times on most of the computational domain depressions are observed to be the same for all four considered SPH-SWE cases. Most of the models built with other software packages use a resolution of 20 m; however, SPH-SWEs had difficulties to predict the correct inundation extent using 20 m resolution. When higher resolution of 10 m was used, as in Case 9, obtained results show similar values as the ones in other Eulerian approaches. This shows that, for overflow simulation with low inflow momentum, high particle resolution is needed to achieve better solutions. The use of high resolution of the open boundary particles, as in Case 9, caused less loss in the water volume.

5. Conclusions

The SPH approximation of the shallow water equations, SPH-SWEs, proved to have good results for river flow simulations. There are two different aspects to address in concluding remarks: limitation of the method to solve the posed problem and the limitations due to implementation.
Obtained results for the cases of water profiles in non-uniform flow and wave propagation are good, with a convergence to semi-analytical solution of approximate 1.4 × 10−3. These subcritical and supercritical flows were handled easily with no restrictions regarding bed complexity and steepness, which shows that the method can be used to solve such hydraulic problems. The momentum conservation property was verified, and the results were compared to other grid-based software results with an agreement of up to 0.5%. The overflow cases associated with low-momentum flow required high particle resolution to achieve better results.
Apart from the abovementioned advantages in using the SPH method of solution, some drawbacks are recognized, including the ones due to the implementation of the code. The existing approaches for closed boundary simulations cause some oscillatory behavior and decrease the method accuracy order. The closed boundary development needs to be extended to include an option for permeable close boundary for further coupling with sub-surface flow models. In terms of implementation, the open boundary simulation requires to associate velocities and water depth values to the open boundary particles. This requires the correct calculations of the velocities and water depths from the discharge hydrographs particularly when the open boundary is close to the area of interest. The effect of this is a loss in the inflow water volume especially when using coarse open boundary particle resolution.
The method requires a good definition of the friction term, which is very important, particularly in the region of very shallow water as it might overcome the gravitational force and causes particle to move opposite to the flow direction. This problem was noticed in cases with wave front advancing in an initially dry bed. Therefore, the determination of the minimum depth for friction is important for calculating the momentum for wave front particles. Furthermore, regarding the implementation, the distribution pattern used in the particle splitting procedure cannot be applied for particles close to the boundary as some of the daughter particles will be located outside of the domain. An adaptive particle splitting procedure that takes into consideration the particle locations when defining the daughter particles’ distribution pattern is welcomed.
The implementation of the SPH-SWEs allows for higher time step computation as compared to the SPH formulations for Navier–Stokes equations (SPH-NS), which might lead to less simulation computational time. Coupling of SPH-SWEs with other grid-based method can be done to utilize SPH-SWEs to solve the issue of wet-dry interfaces.
Finally, this research shows the SPH-SWEs method capabilities of handling hydraulic simulations while it opens up the possibility for the method to be used for both fluvial and coastal areas. SWE-SPHysics, as an open source code, gives an insight of how the method does work and provides a platform for researchers and users to further improve the method and add more functionalities to it.

Author Contributions

Conceptualization, methodology, writing—review and editing, analysis, I.P.; methodology, model development, original draft preparation, preliminary analysis, S.A.M.F.-E.; software testing, analysis, review and editing, C.M.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lucy, L.B. A numerical approach to the testing of the fission hypothesis. Astron. J. 1977, 1013–1024. [Google Scholar] [CrossRef]
  2. Gingold, R.A.; Monaghan, J.J. Smoothed Particle Hydrodynamics: Theory and application to non-spherical stars. Mon. Not. R. Astron. Soc. 1977, 181, 375–389. [Google Scholar] [CrossRef]
  3. Monaghan, J.J. An introduction to SPH. Comput. Phys. Commun. 1988, 48, 89–96. [Google Scholar] [CrossRef]
  4. Liu, M.B.; Liu, G.R. Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments. Arch. Comput. Methods Eng. 2010, 17, 25–76. [Google Scholar] [CrossRef] [Green Version]
  5. Monaghan, J.J. Smooth Particle Hydrodynamics. Rep. Prog. Phys. 2005, 68, 1703–1754. [Google Scholar] [CrossRef]
  6. Mirauda, D.; Albano, R.; Sole, A.; Adamowski, J. Smoothed Particle Hydrodynamics Modeling with Advanced Boundary Conditions for Two-Dimensional Dam-Break Floods. Water 2020, 12, 1142. [Google Scholar] [CrossRef] [Green Version]
  7. Novak, G.; Tafuni, A.; Domínguez, J.M.; Cetina, M.; Žagar, D. A Numerical Study of Fluid Flow in a Vertical Slot Fishway with the Smoothed Particle Hydrodynamics Method. Water 2019, 11, 1928. [Google Scholar] [CrossRef] [Green Version]
  8. Kazemi, E.; Tait, S.; Shao, S.; Nichols, A. Potential Application of Mesh-Free SPH Method in Turbulent River Flows. In Hydrodynamic and Mass Transport at Freshwater Aquatic Interfaces; Rowiński, P., Marion, A., Eds.; GeoPlanet: Earth and Planetary Sciences; Springer: Cham, Switzerland, 2016. [Google Scholar] [CrossRef] [Green Version]
  9. Tran-Duc, T.; Meylan, M.H.; Thamwattana, N.; Lamichhane, B.P. Wave Interaction and Overwash with a Flexible Plate by Smoothed Particle Hydrodynamics. Water 2020, 12, 3354. [Google Scholar] [CrossRef]
  10. Gu, S.; Zheng, X.; Ren, L.; Xie, H.; Huang, Y.; Wei, J.; Shao, S. SWE-SPHysics Simulation of Dam Break Flows at South-Gate Gorges Reservoir. Water 2017, 9, 387. [Google Scholar] [CrossRef] [Green Version]
  11. Ata, R.; Soulaïmani, A. A stabilized SPH method for inviscid shallow water flows. Int. J. Numer. Methods Fluids 2005, 47, 139–159. [Google Scholar] [CrossRef]
  12. Meister, M.; Rauch, W. Modelling aerated flows with smoothed particle hydrodynamics. J. Hydroinform. 2015, 17, 493–504. [Google Scholar] [CrossRef] [Green Version]
  13. Gan, B.S.; Nguyen, D.K.; Han, A.; Alisjahbana, S.W. Proposal for fast calculation of particle interactions in SPH simulations. J. Comput. Fluids 2016, 104, 20–29. [Google Scholar] [CrossRef]
  14. Vacondio, R.; Rogers, B.; Stansby, P.; Mignosa, P. SPH Modeling of Shallow Flow with Open Boundaries for Practical Flood Simulation. J. Hydraul. Eng. 2012, 138, 530–541. [Google Scholar] [CrossRef]
  15. Vacondio, R.; Rogers, B.D.; Stansby, P.K.; Mignosa, P. A correction for balancing discontinuous bed slopes in two-dimensional smoothed particle hydrodynamics shallow water modeling. Int. J. Numer. Methods Fluids 2013, 71, 850–872. [Google Scholar] [CrossRef]
  16. Vacondio, R.; Rogers, B.D.; Stansby, P.K. Accurate particle splitting for smoothed particle hydrodynamics in shallow water with shock capturing. Int. J. Numer. Methods Fluids 2012, 69, 1337–1410. [Google Scholar] [CrossRef]
  17. Rodriguez-Paz, M.; Bonet, J. A corrected smooth particle hydrodynamics formulation of the shallow-water equations. Comput. Struct. 2005, 83, 1396–1410. [Google Scholar] [CrossRef]
  18. Ferrand, M.; Violeau, D.; Mayhoffer, A.; Mahmood, O. Correct boundary conditions for turbulent SPH. In Advances in Hydroinformatics; Gourbesville, P., Cunge, J., Caignaert, G., Eds.; Springer: Singapore, 2015; pp. 245–258. [Google Scholar]
  19. Sarfaraz, M.; Pak, A. SPH Numerical Simulation of Tsunami Wave Forces Impinged on Bridge Superstructures. Coast. Eng. 2017, 121, 145–157. [Google Scholar] [CrossRef]
  20. SWE-SPHysics SWE-SPHysics Code v1.0. Available online: http://www.sphysics.org (accessed on 20 March 2021).
  21. Monaghan, J.J.; Kajtar, J.B. SPH particle boundary forces for arbitrary boundaries. Comput. Phys. Commun. 2009, 180, 1811–1820. [Google Scholar] [CrossRef]
  22. Ferrari, A.; Dumbser, M.; Toro, E.F.; Armanini, A. A new 3D parallel SPH scheme for free surface flows. Comput. Fluids 2009, 38, 1203–1217. [Google Scholar] [CrossRef]
  23. Napoli, E.; De Marchis, M.; Vitanza, E. Panormus-SPH. 2015, A new Smoothed Particle Hydrodynamics solver for incompressible flows. Comput. Fluids 2015, 106, 185–195. [Google Scholar] [CrossRef]
  24. Néelz, S.; Pender, G. Benchmarking the Latest Generation of 2D Hydraulic Modelling Packages. 2013. Available online: https://publications.environment-agency.gov.uk/ms/BXoKPi (accessed on 1 March 2021).
  25. Hartanto, I.M.; Beevers, L.; Popescu, I.; Wright, N.G. Application of a coastal modelling code in fluvial environments. Environ. Model. Softw. 2011, 26, 1685–1695. [Google Scholar] [CrossRef] [Green Version]
  26. Beevers, L.; Popescu, I.; Pan, Q.; Pender, D. Applicability of a coastal morphodynamic model for fluvial environments. Environ. Model. Softw. 2016, 80, 83–99. [Google Scholar] [CrossRef]
Figure 1. Inflow hydrograph at the upstream boundary: for wave attenuation and wave translation cases.
Figure 1. Inflow hydrograph at the upstream boundary: for wave attenuation and wave translation cases.
Water 13 01671 g001
Figure 2. Momentum conservation over a hump case: longitudinal profile and inflow hydrograph at the upstream boundary. Total inflow volume of the inflow hydrograph in Figure 2 is 1310 m3.
Figure 2. Momentum conservation over a hump case: longitudinal profile and inflow hydrograph at the upstream boundary. Total inflow volume of the inflow hydrograph in Figure 2 is 1310 m3.
Water 13 01671 g002
Figure 3. Filling of floodplain depressions case: (a) DEM map showing the location of the upstream boundary condition (red line upper left corner) and ground elevation contour lines every 0.05 m; (b) inflow hydrograph on the left upper corner of the domain.
Figure 3. Filling of floodplain depressions case: (a) DEM map showing the location of the upstream boundary condition (red line upper left corner) and ground elevation contour lines every 0.05 m; (b) inflow hydrograph on the left upper corner of the domain.
Water 13 01671 g003
Figure 4. Results of 1D simulation for M1 curve longitudinal profile: (a) velocity and (b) water level.
Figure 4. Results of 1D simulation for M1 curve longitudinal profile: (a) velocity and (b) water level.
Water 13 01671 g004
Figure 5. Results of 2D simulation for M1 curve longitudinal profile: (a) velocity and (b) water level.
Figure 5. Results of 2D simulation for M1 curve longitudinal profile: (a) velocity and (b) water level.
Water 13 01671 g005
Figure 6. Results of 2D simulation for M1 curve cross-sectional view: (a) velocity and (b) water level.
Figure 6. Results of 2D simulation for M1 curve cross-sectional view: (a) velocity and (b) water level.
Water 13 01671 g006
Figure 7. Results of 1D simulation for M2 curve longitudinal profile view: (a) velocity and (b) water level.
Figure 7. Results of 1D simulation for M2 curve longitudinal profile view: (a) velocity and (b) water level.
Water 13 01671 g007
Figure 8. Results of 2D simulation for M2 curve longitudinal profile: (a) velocity and (b) water level.
Figure 8. Results of 2D simulation for M2 curve longitudinal profile: (a) velocity and (b) water level.
Water 13 01671 g008
Figure 9. Results of 2D simulation for M2 curve cross-sectional view: (a) velocity and (b) water level.
Figure 9. Results of 2D simulation for M2 curve cross-sectional view: (a) velocity and (b) water level.
Water 13 01671 g009
Figure 10. Wave attenuation hydrographs at different cross-sections x, along modelled channel.
Figure 10. Wave attenuation hydrographs at different cross-sections x, along modelled channel.
Water 13 01671 g010
Figure 11. Wave attenuation: channel rating curve.
Figure 11. Wave attenuation: channel rating curve.
Water 13 01671 g011
Figure 12. Wave translation hydrographs at different locations along the modelled channel.
Figure 12. Wave translation hydrographs at different locations along the modelled channel.
Water 13 01671 g012
Figure 13. Wave translation: channel rating curve.
Figure 13. Wave translation: channel rating curve.
Water 13 01671 g013
Figure 14. Momentum conservation over a bump: Point 1 water level results.
Figure 14. Momentum conservation over a bump: Point 1 water level results.
Water 13 01671 g014
Figure 15. Momentum conservation over a bump: Point 2 water level results.
Figure 15. Momentum conservation over a bump: Point 2 water level results.
Water 13 01671 g015
Figure 16. Filling of floodplain depressions: water level results for Depression 4 (closest to inflow boundary).
Figure 16. Filling of floodplain depressions: water level results for Depression 4 (closest to inflow boundary).
Water 13 01671 g016
Figure 17. Filling of floodplain depressions: water level results for Depression 6.
Figure 17. Filling of floodplain depressions: water level results for Depression 6.
Water 13 01671 g017
Table 1. Gradually varied flow models description.
Table 1. Gradually varied flow models description.
Model Set-UpM1 Flow CurveM2 Flow Curve
1D
Schematization
2D
Schematization
1D
Schematization
2D
Schematization
Channel length (m)10,00010,00016,65010,000
Channel width (m)14001400
Bed slope (-)0.0010.0010.0010.001
Manning coefficient (m−1/3s)0.010.020.020.02
Specific discharge (m3/s/m)2.004.992.504.99
Normal depth (m)1.262.002.802.00
Critical depth (m)0.741.360.861.36
Upstream velocity (m/s)1.592.490.892.49
Downstream velocity (m/s)0.401.002.503.56
Upstream water level (m)1.252.002.802.00
Downstream water level (m)5.005.001.001.40
Particles spacing (m)100(a) 10100(a) 10
(b) 100(b) 100
Initial number of particles100400100400
40,000 40,000
Table 2. Wave propagation model set-up description.
Table 2. Wave propagation model set-up description.
Model Set-UpDiffusive WaveKinematic Wave
(Attenuation)(Translation)
Channel length (m)15,00015,000
Channel width (m)6060
Bed slope (-)0.0010.01
Manning coefficient (m−1/3s)0.030.02
Specific discharge (m3/s/m)2.002.50
Normal depth (m)1.030.67
Critical depth (m)0.490.86
Upstream BC 1 velocity (m/s)2.307.22
Downstream BC water level (m)3.621.85
Particles spacing (m)1010
Minimum depth for friction (m)0.050.05
Initial number of particles90009000
1 BC—Boundary conditions, for peak inflow conditions.
Table 3. Momentum conservation test cases.
Table 3. Momentum conservation test cases.
Model Set-UpCase 1Case 2Case 3Case 4Case 5
BC 1 maximum velocity (m/s)3.2753.2753.2753.2752.0
BC maximum water depth (m)0.20.20.20.20.3275
Particles spacing (m)2.02.00.552
Split particles NoNoNoYesNo
Min. depth for friction (m)0.050.0050.050.050.05
1 BC—Boundary conditions.
Table 4. Filling of floodplain depressions test cases.
Table 4. Filling of floodplain depressions test cases.
Model Set-UpCase 6Case 7Case 8Case 9
BC 1 maximum velocity (m/s)0.802.670.500.80
B.C. maximum water depth (m)0.250.750.400.25
Particles spacing (m)20202010
Split particles NoNoNoNo
Min. depth for friction (m)0.050.050.050.05
1 BC—Boundary conditions.
Table 5. Maximum reported error.
Table 5. Maximum reported error.
Non-Dimensional Error 11D SPH Model
(Δx = 100 m)
2D SPH Model
(100 m × 100 m)
2D SPH Model
(10 m × 10 m)
M1M2M1M2M1M2
Error d (water depth)1.42 × 10−25.80 × 10−27.91 × 10−26.1 × 10−24.86 ×10−22.43 × 10−2
Error vx4.36 × 10−32.45 × 10−24.21 × 10−25.21 × 10−29.99 × 10−21.99 × 10−2
Error vy--1.19 × 10−24.15 × 10−35.01 × 10−23.16 × 10−2
1 See formulas in Equations (10)–(12).
Table 6. Summary of results for the nine flooding cases.
Table 6. Summary of results for the nine flooding cases.
CaseLoss in Inflow Water Volume (%)Total No. of Particles at the End of SimulationComputer Running Time (s)
Case12.61986262
Case 22.61983244
Case 30.531,73567,956
Case 46.021122246
Case 56.51181228
Case 61.2100020,570
Case 713.7286430,015
Case 80.331423,287
Case 90.16435054,728
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fadl-Elmola, S.A.M.; Ciocan, C.M.; Popescu, I. Application of Smooth Particle Hydrodynamics to Particular Flow Cases Solved by Saint-Venant Equations. Water 2021, 13, 1671. https://doi.org/10.3390/w13121671

AMA Style

Fadl-Elmola SAM, Ciocan CM, Popescu I. Application of Smooth Particle Hydrodynamics to Particular Flow Cases Solved by Saint-Venant Equations. Water. 2021; 13(12):1671. https://doi.org/10.3390/w13121671

Chicago/Turabian Style

Fadl-Elmola, Salman A. M., Cristian Moisescu Ciocan, and Ioana Popescu. 2021. "Application of Smooth Particle Hydrodynamics to Particular Flow Cases Solved by Saint-Venant Equations" Water 13, no. 12: 1671. https://doi.org/10.3390/w13121671

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop