Corrected SPH methods for solving shallow-water equations*

Shan-qun CHEN （陈善群）， Bin LIAO （廖斌）， Tao HUANG （黄涛）

College of Architecture and Civil Engineering， Anhui Polytechnic University， Wuhu 241000， China，

E-mail: chenshanqun@126.com



Corrected SPH methods for solving shallow-water equations*

Shan-qun CHEN （陈善群）， Bin LIAO （廖斌）， Tao HUANG （黄涛）

College of Architecture and Civil Engineering， Anhui Polytechnic University， Wuhu 241000， China，

E-mail: chenshanqun@126.com

The artificial viscosity in the traditional smoothed particle hydrodynamics （SPH） methodology concerns some empirical coefficients， which limits the capability of the SPH methodology. To overcome this disadvantage and further improve the accuracy of shock capturing， this paper introduces two other ways for numerical viscosity， which are the Lax-Friedrichs flux and the twoshock Riemann solver with MUSCL reconstruction to provide stability. Six SPH methods with different kinds of numerical viscosity are tested against the analytical solution for a 1-D dam break with a wet bed. The comparison shows that the Lax-Friedrichs flux with MUSCL reconstruction can capture the shock wave more accurate than other five methods. The Lax-Friedrichs flux and the artificial viscosity with MUSCL reconstruction are finally both applied to a 2-D dam-break test case in a L-shaped channel and the numerical results are compared with experimental data. It is concluded that this corrected SPH method can be used to solve shallow-water equations well.

smoothed particle hydrodynamics （SPH） methodology， artificial viscosity， Lax-Friedrichs flux， two-shock Riemann solver， MUSCL reconstruction， shallow water equations

Introduction

The shallow-water equations （SWEs） are widely used for hydrodynamic simulations in coastal regions，bays， estuaries and lakes， to predict tsunamis， dam breaks， storm surges， floods and other natural disasters［1］. Due to the obvious nonlinearity of the SWEs， the analytical solutions can only be obtained in rare special circumstances， so numerical simulation methods are required in actual projects.

The grid-based classical Euler methods are now widely applied to solve the SWEs， such as the finite difference method and the finite volume method. However， due to restrictions of grid， grid-based methods suffer many limitations in simulating multi-phase effects， most importantly， the debris flows in flood modeling. On the other hand， the particle method requires no grid， therefore， the grid distortion and reconstruction problems can be avoided， with a natural advantage in dealing with large deformation for free interface. This feature makes particle methods promising in solving the SWEs.

The smoothed particle hydrodynamics（SPH） is a purely Lagrangian meshless method originally introduced to simulate astrophysical problems by Lucy［2］in 1997. The method was then applied to solve the Navier-Stokes equations， and it now becomes increasingly popular as a technique to study a range of applications， including wave breaking， impact-fracture problems， and bio-medical problems. The SWEs are based on the incompressible Navier-Stokes equations with the assumption of the hydrostatic pressure and the Boussinesq approximation， which provides a theoretical basis for the use of the SPH method in their solutions. Ata and Soulaïmani［3］obtained some results in the wet bed test case by improving the formulation of the stabilization term. Rodriguez-Paz and Bonet［4］presented a SPH formulation for shallow water， based on the variational formulation， which can conserve boththe total mass and momentum. De Leffe et al.［5］solved the nonlinear SWEs by an SPH method and presented coastal flow simulations.

In the process of solving SWEs by the SPH methods， virtual numerical oscillations are produced in the vicinity of the shock wave. Traditionally， an artificial viscosity was added to the SPH momentum equation to suppress the non-physical oscillation. Nevertheless， the artificial viscosity contains some empirical coefficients［6］， which is different in different test cases and so the conventional SPH methods suffer some limitations in solving the SWEs.

To improve the accuracy and the generality of the conventional SPH method in modeling the SWEs， this paper introduces two schemes of numerical viscosities in the SPH method， and also uses the MUSCL reconstruction to reduce the level of numerical viscosity. Then， the corrected SPH methods are employed to solve classic shallow-water test cases and the results are compared with exact solutions. Finally， the ability of the shock capturing of the corrected SPH methods is verified by a more complicated numerical experiments.

1. SPH for shallow water

1.1 Lagrangian formulation of SWEs

The SWEs are the depth-integrated equations of mass and momentum conservations and are written in the Lagrangian form as

Equations （1） are in the form the same as the Euler equations if we redefine the densityas the amount of fluid per unit of area in a 2-D domain； with this new definition of， it can be related to the depth of wateras

1.2 Density evaluation

The SPH approximation for the density of theparticleis

In the shallow-water approximation， the fluid will follow the terrain and its projected 2-D density will expand or contract according to the height of the water column as shown by Eq.（2）. A variable smoothing length is therefore needed in order to maintain the accuracy of the solution. In general，must vary according to［4］

According to Eq.（5）， the smoothing lengthof theparticle is related to the densityas

Differentiating Eq.（4） and using the chain rule for the kernel leads to

The derivative of the kernel function with respect tois obtained as

Substituting Eqs.（8） and （9） andinto Eq.（7） leads to

Converting the directional derivative of Eq.（11）into the derivative of the density， we have

The aforementioned Eq.（6） is implicit because the densityis a function ofas in Eq.（4）， a Newton-Raphson iteration is adopted to solve Eqs.（4）and （6）.

The root of Eq.（7） can be found by using the Newton-Raphson iterative formula

The derivative of the residual is calculated by differentiating Eq.（13） and using the chain rule for the kernel function

Substituting the derived results of Eq.（9） into the above equation and remembering that， we have

Substituting Eq.（16） into Eq.（14） gives the final iterative formula for

where

The Newton-Raphson iterations can be conducted independently for each particle and will be stopped when. Then entering into the iteration process of particlesand until all particles are covered. Since now we focus on the problem of the poor ability in shock capturing in solving the SWEs and the accuracy requirements are not very important， we let the coefficient

1.3 Momentum equation

The Lagrangian equation of motion for a particle iis

where the Lagrangian functionalis defined in terms of the kinetic energyand the potential energyasis a function of particle positions but not velocities. The kinetic energy for a system of particles can be calculated as the sum of the energies of all particles

Fig.1 Flow with a free surface under the effect of gravity

According to Newton’s second law， Eq.（19） is equivalent to

Substituting the kinetic energyinto the inertial forcegives

The total internal energy stored in the group of particles is

Substituting the compression ratioand the pressureinto Eq.（25） for the equivalent transformation， the directional derivative ofis （see Ref.［7］）

Substituting the derivative of（Eq.（11）） in the above equation and rearranging the summations gives

The comparison of Eq.（26） with Eq.（28） gives the internal force

Substituting in Eq.（29） the pressureobtained by means of the hydrostatic law:， the final formulation foris

By substituting Eqs.（23）， （24） and （30） into （22）and taking into account also the effect of the friction source term， the particle accelerationcan be finally obtained as

1.4 Time integration scheme

To integrate in time the particle positions and velocities， we use the leap-frog time integration scheme［8］defined as:

As an explicit method， the time step must satisfy a Courant-Friedrichs-Lewy （CFL） condition［9］. In the SPH， this condition is imposed with the smoothing length as a reference length

2. Stabilizing treatments

2.1 Numerical viscous improvement

In the Von Neumann stability analysis system［10］，the SPH method can be interpreted as a central finite difference scheme and some viscosity is needed to avoid numerical oscillations in the presence of shock waves. Therefore， Eq.（30） should be modified as follows

In the original SPH formulation introduced by Monaghan［6］，is an artificial viscosity activated when two particles are approaching. The main drawback of this formulation is that it needs to be tuned according to the necessary numerical viscosity， which is different in different test cases. In order to overcome this drawback， the paper introduces two modified schemes.

（1） Lax-Friedrichs flux

According to Ata and Soulaïmani［3］， the centred fluxin the Lax-Friedrichs flux is replaced as

Fig.2 Initial condition of 1D dam break flow with wet bed

After some algebraic operations， the following expression of the stabilizingis obtained

Fig.3 Water depth for 1-D dam break with wet bed at time 0.05 s

（2） Two-shock Riemann solvers

The Riemann solvers are widely used in finite volume schemes for hyperbolic equations［11，12］and there were some attempts to introduce them in the SPH formalism［13］. Comparing these approaches with the artificial viscosity method， the advantage of the Riemann solvers is that no extra numerical dissipation is introduced.

In this work we introduce the two-shock Riemann solver［14］into our Shallow Water models. The main idea is to consider each interaction between theandparticles as a Riemann problem and therefore to replace the pressuresin Eq.（29） with the resultant pressure

where， according to the two-shock Riemann solver，are the left and right water depthsrepresent the normal velocities in the left and right statesandis an estimate of the water depth that can be obtained from some other direct Riemann Solvers.

2.2 MUSCL-type reconstruction

To reduce the level of the numerical viscosity and improve the accuracy of the numerical calculation，a monotone upwind-centred scheme for the conservation law （MUSCL） non-upwind procedure［15］is used to reconstruct a generic physical quantityin the left and right states of the Riemann problem

The MUSCL reconstruction used on two amendments in above two sections is applied to reconstruct the velocitiesand the water depthsin theterm of the Lax-Friedrichs flux （see Eq.（36））， whereis replaced with， and in the two-shock Riemann solver （see Eq.（39））.

3. Test cases

The dam break flow can cause disasters in the downstream， which would propagate in rivers in the form of standing wave， and the wave crest would generate a sudden rise of water level along its path［16-18］. The SWEs are widely used to simulate the water level in dam break flows， which are difficult to capture exactly. Hence we choose two dam break cases to test the corrected SPH methods.

3.1 1-D dam break flow with wet bed

To verify the effect of the shock capturing capability by using the stabilization schemes in solving shallow water equations， a 1-D dam break case is considered in this section. The initial conditions are shown in Fig.2， where the water depthin the upstream part of the domainandin the downstream part of the domain. There are 150 particles scattered nonuniformly inside the domain according to the water depth. In the upstream part of the reservoir， we letin the six different simulations， for the particles placed downstream，is twice the initial particle spacing for the upstream part. In the test， no source terms are considered and the initial velocity is 0.

Fig.4 Water depth for 1-D dam break with wet bed at time 0.05 s

Figure 3 shows the comparison between the analytical solution and the SPH results obtained by using different kinds of numerical viscosity. The results in Figs.3（a）， 3（c） and 3（e） show that the three kinds of numerical viscosity， the artificial viscosity， the Lax-Friedrichs flux and the two-shock Riemann solver can all capture the shock wave with a certain degree of accuracy， but the additional numerical viscosity can cause unnecessary oscillations with a significant deviation of the water line in corners. In order to improve the shock capturing ability， the reconstruction technique is required to introduce into the three kinds of numerical viscosity mentioned above to prevent from producing the rarefaction wave. The comparison isshown in Figs.3（a）-3（f）. The results indicate that the viscosity terms with the MUSCL reconstruction in comparison to those without reconstruction can reproduce the sharper shock without introducing unnecessary oscillatios in the rarefaction wave. Finally， the comparison is shown in Figs.3（b）， 3（d）and 3（f）. The results show that the artificial viscosity and the twoshock Riemann solver with the MUSCL reconstruction both overpredict the water depth in the initial part of the rarefaction wave， but with the Lax-Friedrichs flux with reconstruction， more accurate results are obtained.

Table 1non-dimensional norm of water depth error calculated for 1-D dam break with wet bed atconsidering 6 different stabilization terms

Table 1non-dimensional norm of water depth error calculated for 1-D dam break with wet bed atconsidering 6 different stabilization terms

x∆AV 　LF 　TS 　AV+MUSCL 　LF+MUSCL 　TS+MUSCL 0.0100 m 　1.52×10-2　1.57×10-2　1.56×10-2　1.03×10-2　1.00×10-2　1.02×10-20.0050 m 　9.29×10-3　9.81×10-3　9.81×10-3　6.01×10-3　5.88×10-3　6.01×10-30.0025 m 　5.82×10-3　6.41×10-3　6.41×10-3　4.12×10-3　3.87×10-3　4.12×10-3

Figure 4 shows the comparison of water levels between the exact solution and the three viscosity schemes with MUSCL reconstruction， as well as， the local amplification of the shock wave front. From 4（b）and 4（c）， it is seen that the Lax-Friedrichs flux method is more accuate than other two in the positon of a sudden drop of the water level.

In order to illustrate the shock capturing capability of the six terms， a convergence analysis is also performed by using three different initial particle spacings0.01 m， 0.005 m and 0.0025 m， respectively， theerror norm of the nondimensional water depth is defined as

3.2 2-D dam-break flow in a L-shaped channel

In order to validate the shock wave capturing ability of the Lax-Friedrichs flux with the MUSCL reconstruction， the case of a 2-D dam-break flow in a channel with abend［19］is taken for simulation.

Fig.5 Geometry of the reservoir and L-shaped channel: plane and profile views and positions of the gauges of the experimental setup （m）

The flow domain consists of a square reservoir and the L-shaped channel as shown in Fig.5. The upstream reservoir has dimensions of 2.39 m×2.44 m，the channel is rectangular， 0.495 m wide， the upstream reach is 3.92 m long and the downstream reach， behind thebend， is 2.92 m long. The bottom of the channel is flat and is 0.33 m higher than that of the reservoir. Initially， the water depth is 0.53 m high in the reservoir， which is separated by a gate from the channel and then the gate is suddenly opened to produce a dam-break situation. The water levels are recorded during the experiment in the reservoir and along the channel using 6 gauges， as shown in Fig.5. In the simulation， 2 450 particles are initially placed in the reservoir over an uniform Cartesian grid with size. The channel bed is initially dry， and its Manningʼs friction coefficient is

Fig.6 Water level profiles at typical times of 2-D L-shaped dam break with dry bed （m）

Figure 6 shows some typical water-level profiles at different times in the process of dam break. At3.2 s， the front reaches the bend and a bore forms in the corner， at， the bore travels back to the reservoir， atit disappears， atbecause of the effect of the channel end wall， a bore forms again and travels upstream， at， the second bore disppears and at， the water flow becomes stable.

The variation of the water depth with time by the SPH-SWEs presented in this paper is compared with the experimental data at different gauge positions as shown in Fig.7. The Lax-Friedrich flux with MUSCL reconstruction is applied for its good performance mentioned in the last section. Numerical results obtained from the SPH with the artificial viscosity are presented here for a comparison to those obtained by using the numerical viscosity. Gauge 1 is placed inside the reservoir near the channel， the good agreement of the experimental data with the numerical results means that the discharge entering the channel is correct. Gauges 2， 3， 4 are placed along the channel upstream of the bend， therefore， they can capture the abrupt water level elevations because of the reflected wave travelling to the reservoir. The numerical model can reproduce the water level at Gauges 2， 3 and 4， especially when at 16 s at Gauge 2， at 14 s at Gauge 3 and at 9 s at Gauge 4， where the great change of water level occurs. However， differences are witnessed at Gauge 5. The overall disagreement at Gauge 5 is due to the local head loss caused by thebend which is not taken into account in the numerical computation. At Gauge 6， a good agreement with the experimental data is evidenced. In general， all numerical results are in good agreement with the experimental data， except that the AV+MUSCL method slightly overpredicts the water level at Gauge 2. However， the computational time of the simulation with the LF+MUSCL method is less than that of the AV+MUSCL method by about 1 000 s.

4. Conclusion

In the traditional SPH method， the artificial viscosity is added to the SPH momentum equation foreliminate the non-physical oscillation generated in the vicinity of the shock wave. Nevertheless， the artificial viscosity formulation needs to be tuned according to the necessary numerical viscosity， which is different in different test cases. To improve the accuracy and the generality of the conventional SPH method， this paper introduces two numerical viscosities， which are the two-shock Riemann solver and the Lax-Friedrichs flux. In order to reduce the numerical oscillation in the computational processes， the MUSCL reconstruction is used to reconstruct the velocity and the water depths in the artificial viscosity， the two-shock Riemann solver and the Lax-Friedrichs flux. These improved SPH methods are tested against the analytical solution for the 1-D dam break with wet bed. The results show that the shock waves are simulated accurately by schemes of the numerical viscosity with reconstruction procedures for stability and the best results are obtained by using the Lax-Friedrichs flux with MUSCL reconstruction. Finally， the Lax-Friedrichs flux and the artificial viscosity with MUSCL reconstruction are both applied to the case of a 2-D dam-break flow in a channel with a L-shaped bend. Both methods make good predictions as compared with experimental measurements， but the number of iterations necessary to converge with the LF+MUSCL is less than that for the AV+MUSCL， thus the first method is more efficient. In conclusion， the corrected SPH method can solve shallow water problems with improved accuracy and generality.

Fig.7 Water levels recorded by gauges from G1 to G6

Acknowledgement

This work was supported by the opening fund of key Laboratory of Mechanics， Anhui Polytechnic University （Grant No. 201602）.

References

［1］ ZHAO Zhang-yi. Numerical simulation and application of a Runge-Kutta discontinuous Galerkin scheme for one-dimension shallow water equations［D］. Master Thesis， Tianjin， China: Tianjin University， 2010（in Chinese）.

［2］ LUCY L. A numerical approach to the testing of fusion process［J］. Joural of Astronomical， 1977， 8（12）: 1013-1024.

［3］ ATA R.， SOULAÏMANI A. A stabilized SPH method for inviscid shallow water flows［J］. International Journal for Numerical Methods in Fluids，2004， 47（2）: 139-159.

［4］ RODRIGUEZ-PAZ M.， BONET J. A corrected smooth particle hydrodynamics formulation of the shallow-water equations［J］. Computers and Structures， 2005， 83（17-18）: 1396-1410.

［5］ De LEFFE M.， Le TOUZÉ D. and ALESSANDRINI B. Coastal flow simulations using an SPH formulation modeling the nonLinear shallow water equations［C］. Proceedings of the 3th ERCOFTAC SPHERIC workshop on SPH applications. Lausanne， Switzerland， 2008， 48-54.

［6］ MONAGHAN J. J. Smoothed particle hydrodynamics and its diverse applications［J］. Annual Review Fluid Mechanics， 2012， 44: 323-346.

［7］ BONET J.， KULASEGARAM S. Correction and stabilisation of smooth particle hydrodynamics with application in metal forming［J］. International Journal for Numerical Methods in Engineering， 2000， 47（6）: 1189-1214.

［8］ CHEN Shan-qun， LIAO Bin. Numerical simulation of free surface flows based on SPS-SPH Method［J］. Journal of Ship Mechanics， 2013， 17（9）: 969-981.

［9］ De MOURA C. A.， KUBRUSLY C. S. The Courant-Friedrichs-Lewy （CFL） condition［M］. New York， USA: Springer Science+Business Media， 2013.

［10］ YUSTE S. B.， ACEDO L. An explicit finite difference method and a new von Neumann-type stability analysis for fractional diffusion equations［J］. SIAM Journal on Numerical Analysis， 2005， 42（5）: 1862-1874.

［11］ INUTSUKA S.-I. Reformulation of smoothed particle hydrodynamics with Riemann solver［J］. Journal of Computational Physics， 2002， 179（1）: 238-267.

［12］ CHA S. H.， WHITWORTH A. P. Implementations and tests of Godunov-type particle hydrodynamics［J］. Monthly Notice of the Royal Astronomical Society， 2003， 340（1）: 73-90.

［13］ TRICCO T. S.， PRICE D. J.Constrained hyperbolic divergence cleaning for Smoothed Particle Magnetohydrodynamics［J］. Journal of Computational Physics， 2012，231（21）: 7214-7236.

［14］ TORO E. F. Direct Riemann solvers for the time-dependent Euler equations［J］. Shock Waves， 1995， 5（1-2）: 75-80.

［15］ EDWARDS M. G. The dominant wave-capturing flux: A finite-volume scheme without decomposition for systems of hyperbolic conservation laws［J］. Journal of Computational Physics， 2006， 218（1）: 275-294.

［16］ WU Qiao-rui， TAN Ming-yi and XING Jing-tang. An improved moving particle semi-implicit method for dam break simulation［J］. Journal of Ship Mechanics， 2014，18（9）: 1044-1054.

［17］ YUAN Yue， RONG Gui-wen and DAI Hui-chao et al. Simulation of dam-break flow over partially deformed bed based on 2D FEVM-SWEs model［J］. Chinese Journal of Hydrodynamics， 2015， 30（5）: 549-555（in Chinese）.

［18］ ZHANG Ming-liang， XU Yuan-yuan and QIAO Yang et al. Numerical simulation of flow and bed morphology in the case of dam break floods with vegetation effect［J］. Journal of Hydrodynamics， 2016， 28（1）: 23-32.

［19］ SOARES-FRAZAO S.， SILLEN S. and ZECH Y. Dambreak flow through sharp bends: physical model and 2D Boltzmann model validation［C］. Proceedings of the CADAM Meeting. Wallingford， UK， 1998.

May 30， 2014， Revised October 14， 2014）

* Project supported by the National Natural Science Foundation of China （Grant No. 51175001）， the Natural Science Foundation of Anhui Province （Grant No. 1508085QE100） and the Higher Education of Anhui Provincial Scientific Research Project Funds （Grant No. TSKJ2015B03）

Biography: Shan-qun CHEN （1981-）， Female， Ph. D.，

Associate Professor

Bin LIAO，

E-mail: liaobinfluid@126.com