Open access peer-reviewed chapter

A Numerical Simulation of the Shallow Water Flow on a Complex Topography

Written By

Alexander Khoperskov and Sergey Khrapov

Submitted: March 13th, 2017 Reviewed: September 15th, 2017 Published: December 20th, 2017

DOI: 10.5772/intechopen.71026

Chapter metrics overview

1,122 Chapter Downloads

View Full Metrics


In current chapter, we have thoroughly described a numerical integration scheme of nonstationary 2D equations of shallow water. The scheme combines the smoothed particle hydrodynamics (SPH) and the total variation diminishing (TVD) methods, which are sequentially used at various steps of the combined SPH-TVD algorithm. The method is conservative and well balanced. It provides stable through calculations in presence of nonstationary “water-dry bottom” boundaries on complex irregular bottom topography including the transition of such a boundary between wet and dry bottom through the computational boundary. Multifarious tests demonstrate the effectiveness of the combined SPH-TVD scheme application for a solution of diverse problems of the engineering hydrology.


  • computational fluid dynamics
  • numerical methods
  • hydrology
  • shallow water
  • digital elevation model
  • numerical experiment
  • GPU

1. Introduction

The usage of the shallow water models (SWMs) allows to solve a wide range of engineering tasks related to the dynamics of surface waters (a seasonal flooding [1], a drain and rain flows [2]), the emergence and expansion of the marine nonlinear waves [3] (problems of tsunami impact on a shore, nonlinear waves formation due to earthquakes, and meteorological waves generation by an open ocean resonance [4]), and the inundation in a coastal area by storm surge [5]. The SWM modifications are effective for studying various geophysical problems such as the dynamics of the pyroclastic flows [6] and the riverbed processes including the sediment dynamics and the diffusion of pollutant particles in reservoirs. The multilayer models utilization significantly expands the opportunities of the shallow water approach [7]. The tasks associated with flooding research of river valleys or interfluves [1] are stood out among the hydrological problems. In the framework of the SWM, the important results have also been obtained for the atmospheric phenomena, meteorological forecasts, and global climate models [8, 9].

It should be noted that the SWM is also actively applied and developed for theoretical research of various cosmic gas flows. The shallow water approximation allows describing a whole series of astrophysical objects: the protoplanetary and circumstellar disks [10, 11], the accretion disks around the compact relativistic objects [12], the cyclonic movements in the giant planets’ atmospheres [13], and the spiral galaxies gas disk components [14]. The gravitational fields play the topography role in such astrophysical problems.

A lot of numerical methods have been proposed for the shallow water dynamics modeling employed for diverse tasks and conditions. Despite the fact that conservative methods of finite volume poorly describe stationary states, they allow correctly calculating shock waves and contact discontinuities. The latter problem could be overcome by the so-called well-balanced (WB) circuits [15, 16, 17].

Our main aim is through calculation for a flow with various Froude number within 0 ≤ Fr < 100 (Fr = u/cg, where cg=gH is an analog of the sound speed, H is the depth, u is the flow velocity, g is the specific gravitational force) in order to simulate subcritical (Fr <  < 1), transcritical (Fr~1), and supercritical (Fr >  > 1) flows. Highly heterogeneous terrain topography including vertical discontinuities and small-scale inhomogeneities at the computational domain boundary makes the calculations noticeably complicated. The latter leads to special quality requirements in numerical algorithms. Hence, the modern numerical schemes should simulate fluid movements along the dry bottom and correctly describe the interfaces between wet and dry bottom [1, 18, 19].

Among the numerous numerical methods solving shallow water equations (SWEs), the following methods should be mentioned: the discontinuous Galerkin method based on triangulation [8], the weighted surface-depth gradient method for the MUSCL scheme [18], and the modified finite difference method [20]. The so-called constrained interpolation profile/multimoment finite-volume method utilizing the shallow water approximation is developed to simulate geophysical currents on a rotating planet in spherical coordinate system ([9] and see the references there). The particle-mesh method demonstrates good opportunities for calculation of rotating shallow water [21]. As a rule, numerical schemes of the second-order accuracy give satisfactory results and allow to correctly solve a wide range of tasks for diverse applications [17]. Special attention should be focused on the numerical way of a source term setting, since in the case of discontinuous topography, it may induce an error at the shock wave front [22].

We consider the original CSPH-TVD (combined smoothed particle hydrodynamics—total variation diminishing) algorithm of numerical integration of the Saint-Venant equations. It accounts for the new modifications improving the computational properties of the scheme. A detailed description of the numerical scheme is the main aim of the chapter.


2. Mathematical models

In this section, the mathematical models for terrestrial hydrology, which accounts for the maximum number of physical and meteorological factors and can be described by the shallow water equations, are considered. The SWE or the Saint-Venant equations are fairly simple model describing the free surface dynamics of incompressible fluid (see the details in the review of [23]). The model assumes that vertical equilibrium in a medium exists at every moment of time. The conservation laws of mass and momentum in the integral form for a thin layer of moving substance with additional sources and forces are:


where Q is the sources function, ∇ = ex∂/∂x + ey∂/∂y is the nabla operator, S(t) is the cross-sectional area of the “liquid particle”, r = xex + yey is the radius vector, v = uex + vey is the velocity vector, η(rt) = H(rt) + b(r) is the free surface elevation, b(r) is the bottom profile, f(rt) = fb + ffr + fCor + fw + fs is the sum of the external forces, accounting for the bottom friction fb, the internal friction (viscosity) ffr, the Coriolis force fCor, the effect of atmospheric wind fw, and the force fs determined by the liquid momentum due to the action of the sources Q. The surface density of sources is σ(rt) = σ(+) + σ(-) = dQ/dS, where σ(+) and σ(−) are the liquid sources (rain, melting snow, flows through the hydro-constructions, groundwater, etc.) and sinks (infiltration, evaporation, etc.), respectively.

The digital elevation model (DEM) provides the quality of numerical simulations of real hydrological objects to a significant extent. The DEM is determined by the height matrix bij = b(xiyj) on numerical grid {xiyj} (i = 1, …, Nx, j = 1, …, Ny). The DEM elaboration utilizes diverse geoinformation methods for the processing of spatial data obtained from various sources. A matrix of heights has been built in several stages. At the beginning stage, the remote sensing data are accounted by the function b(xiyj). The river sailing directions and the actual water depth measurements allow to construct the DEM for large river beds (for example, for the Volga River and the Akhtuba River). To improve our DEM, the data on various small topography objects such as small waterways, roads, small dams, etc., should be included into the consideration. The numerical simulation results of the shallow water dynamics reveal flooding areas which may be compared with real observational data (both from the remote methods and our own GPS measurements). Such approach qualitatively improves the DEM as a result of the iterative topography refinement.


3. Numerical method СSPH-TVD

In current paragraph, the CSPH-TVD numerical method is thoroughly examined. The method combines classical TVD and SPH methods at various stages of numerical integration of hyperbolic partial differential equations and uses the particular benefits of both. The advantages of graphical processing units (GPUs) with CUDA acceleration for hydrological simulations are also discussed.

The computational domain has been covered by a fixed uniform grid with a spatial step h, while mobile “liquid particles” (hereinafter particles) are placed in the centers of cells. The time layers tn have a nonuniform step τn = tn − tn − 1, where n denotes the index of the time layer. The vector space index i = (ij) characterizes the radius vectors of the fixed centers of the cells ri0=xi0yj0. The total number of cells and particles is Np = Nx ⋅ Ny. At the initial time moment, the particles are placed in the centers of cells. The main characteristics of the particles are the volume Vin=SitHdS, momentum Pin=SitHvdS, coordinates of their centers rin, and linear sizes in=Sin.

The CSPH-TVD method main stages are the following:

  1. The Lagrangian’s stage. It includes an application of the modified CSPH approach [1, 24]. At this stage, the changes of the integral characteristics and particles positions ‍(VinVin+1, PinPin+1, and rinrin+1) due to the hydrodynamic and external forces acting on them are calculated. Figure 1 demonstrates the deformations of particles occurred during their motion, while the positions of particles are shifted relatively to the cells centers (it is assumed that in+1h). The difference between the CSPH-TVD method and the traditional SPH approach [25] consists in the fact that the particles at each time interval [tntn + 1] of the numerical time integration at the moment of time tn are found in the initial state (rinri0 and in=h).

  2. The Euler’s stage. At the Euler’s stage, a fixed grid is used to calculate the mass and momentum fluxes through the cells boundaries at time moment tn + 1/2 = tn + τn/2. The corresponding changes of the particles integral characteristics are proportional to the difference between the inflow and outflow fluxes. If the region shaded at Figure 1 is inside of the cell, then the substance flows into the cell, otherwise it flows out. The modified TVD approach [1, 24] and approximate solution of the Riemann’s problem [26] are applied to calculate the flows. At the end of the Euler’s stage, the particles return to their initial state.

Figure 1.

The OX-projection of the main stages scheme of the CSPH-TVD method. The dash-dotted lines xi(t) correspond to the trajectories of the particles due to their displacement xinxin+1, the dashed black lines show the boundaries of cells. The particles boundaries deformed during the motion are shown by solid lines. The shaded regions correspond to a change in the particles integral characteristics caused by a flow of matter through the cells boundaries.

It should be noted that when the dry bottom regions are considered, the particles with zero depth and momentum are placed in the corresponding cells. For such particles, the Lagrangian’s stage is skipped. At the Euler’s stage, the mass and momentum fluxes flowing into the corresponding cell with zero depth are calculated in case of the water presence in the neighboring cells. The changes in the particles integral characteristics are determined too. Such a method permits to carry out an effective through calculation in the presence of nonstationary “water-dry bottom” boundaries in the computational domain.

3.1. The Lagrangian’s stage

Let us define the following quantities


being the analogous of the function mean values φ = {HHuHv} in cells in the finite volume approximation on a fixed grid. Taking into the account Eq. (2), the integral characteristics of the particles can be written in the form:


After substituting relations (2) and (3) into the system of Eq. (1), we obtain


where Ui=HiHvi, Φi=σigHiiη+Hifi, iη=Hr=ri+br=ri0, σi = σ(ri), fi = f(ri). The law of motion of a particle is


where vi = (Hv)i/Hi is the velocity of the i-th particle.

For the numerical integration of the system of Eq. (4), an approximation for the spatial derivatives is required. In accordance with the SPH-approach [25], any characteristic of the medium φ and its derivatives ∇φ are replaced by their smoothed values in the flow area Ω:


where W is the smoothing kernel function, h is the smoothing length, k is the vector index similar to the spatial index i. Accounting for Eq. (2), the quantity ψ can be determined as


Our modification of the SPH-method consists in the generalization of Eq. (7). In case of the traditional SPH-approach for the SWE [27] expression ψtrk=Vktφ̂kt/Ĥkt is used instead of (7). Substituting the SPH-approximation (6) and (7) in the right-hand side of the expressions for system (4), we obtain


where W¯ik=h2Wrirkh, W¯ik0=h2Wri0rk0h are the functions computed in the cells centers and afterward used to approximate the gradient of the fixed bottom relief b(r). The approximation (8) ensures the conservatism of the CSPH-TVD method at the Lagrangian’s stage and the fulfillment of the WB-condition on the inhomogeneous topography.

There are three conditions superimposed on the smoothing kernels W:

  1. the kernel is finiteness;

  2. the normalization condition is fulfilled ΩWrrhdr=1;

  3. the zero spatial step limit limh0Wrrh=δrrh is the Dirac delta-function.

For the smoothing kernel W, spline functions of different orders or Gaussian distribution are used in the SPH hydrodynamic simulations [25, 28, 29]. The Monaghan’s cubic spline


is the most commonly applied approximation, where q =  ∣ ri − rk ∣ /h is the relative distance between the particles, Aw=23h107πh21πh3 is the normalization constant for the 1D, 2D and 3D cases, respectively. In the CSPH-TVD method, the smoothing kernel W approximates only the gradients of physical quantities in Eq. (8). Therefore, the value of the normalization constant may be corrected to increase the accuracy of the numerical solution. For example, it was done by the authors of [1] when for the two-dimensional case, the normalization constant was chosen to be Aw0.987107πh2. The smoothing kernel (9) can lead to unphysical clustering of particles in the simulation of hydrodynamic flows by the SPH method; therefore, an approximation for the smoothing kernel of the pressure gradient is applied in form [29]:


where Bw=18h516πh21564πh3 is the normalization constant. Figure 2 shows the comparison between two smoothing kernels (9) and (10). As seen from Figure 2b for q < 0.7, the characteristic feature of the smoothing kernel (9) is the impetuous decrease of the derivative value to zero as the particles approach each other. At small q values, this feature leads to a significant weakening of the hydrodynamic repulsion force between the particles and appearance of nonphysical clusters of particles.

Figure 2.

The spatial distributions of smoothing kernels (a) and their derivatives (b). The smoothing kernels (9) and (10) are shown by solid and dashed lines, correspondingly.

For the numerical integration of the set of differential Eqs. (4) and (5), the method of a predictor-corrector type of second-order accuracy (so-called leapfrog method) is employed. The main steps of the leapfrog method for the Lagrange stage of CSPH-TVD method are:

  1. At the predictor step, the water depth Hi and velocity vi are calculated at time moment tn + 1 = tn + τn:


  2. The Particles’ spatial positions ri are determined at time tn + 1:


  3. During the corrector step, the water depth Hi and velocity vi values are recalculated (tn + 1):


To reduce the rounding errors during the numerical implementation of the algorithm in recurrent relations (11)(13), the value of the relative displacement of particle ξit=ritri0 (here ri0 is the initial position) should be used instead of ri(t).

3.2. The Euler’s stage

At this stage, the relationship between the values of U˜in+1 and Uin+1 is determined at time moment tn + 1. According to the CSPH-TVD approach, the changes in the particles integral characteristics at their return to the initial state ri0 are computed. In order to do it, the following expression


should be integrated over time from tn to tn + 1 . Taking into the account Eq. (2), we obtain


where Fi±1/2,jn+1/2 and Gi±1/2,jn+1/2 are the average values (in the interval [tntn + 1]) of mass and momentum fluxes through the cell boundaries in x and y directions, respectively. To determine these fluxes, we use the TVD-approach [30] and approximate methods of the Riemann’s Problem (RP) solving for an arbitrary discontinuity decay (Lax-Friedrichs (LF) method, Harten-Lax-van Leer (HLL) method [31]). The Riemann’s problem is solved separately for each boundary of the Euler’s cells.

The values of the flow parameters to the left UL and right UR of the considered boundary are the initial conditions for RP. These quantities define magnitude of the jump and can be obtained by a piecewise polynomial reconstruction of the function U(tr). We limit the consideration by the piecewise linear reconstruction providing the second order of accuracy by space coordinate for a numerical scheme [32]. The modification of the TVD-approach used in the CSPH-TVD method is concerned with the reconstruction of the value of U(tr) relatively to the position of the particles at time tn + 1/2 at distance ξi(tn + 1/2) from the centers of the cells. Thus, the piecewise linear profile of the function U(tr) inside of the i-th cell at time tn + 1/2 in the projection onto the Ox axis (Figure 3) has the following form


Figure 3.

The projection onto the Ox axis of the piecewise linear reconstruction of the function U(t, r) with the respect to the mass centers of the particles for the CSPH-TVD scheme at time tn + 1/2.

where Θi,jn+1/2 is the vector of slopes (angular coefficients), xin+1/2=xi0+ξin+1/2. Taking into the account Eq. (15), the expressions for ULand URat the boundary xi+1/20 can be written as


The slopes of the piecewise linear distribution (16) should satisfy the TVD-condition [30]. To fulfill this condition, the limiter function


has to be used. The analysis shows that the limiters of the minmod [33], van Leer [32], and van Albada et al. [34] satisfy the TVD-condition and suppress the nonphysical oscillations of the numerical solution in the vicinity of discontinuities.

3.3. The stability condition and parallel implementation on GPUs

For the stability of the numerical scheme CSPH-TVD the following conditions have to be fulfilled during the integration time τn: (1) at the Lagrangian’s stage the shift of particles’ center of mass should not exceed h/2 relatively to the initial position; (2) at the Euler’s stage the perturbations should not propagate over a distance greater than the cell size h. These requirements provide the numerical stability condition for the time step τn of the CSPH-TVD algorithm:


where 0 < K < 1 is the Courant number.

The CSPH-TVD method demonstrates a good ability for parallelization on graphics processors [35]. Figure 4 shows the computer implementation of the CSPH-TVD method on GPUs. It indicates the execution sequence of CUDA kernels at various stages of our algorithm.

Figure 4.

The flow diagram for each of the calculation module: K1 determines the presence of water in the CUDA block; K2 calculates the forces at the time moment tn at the Lagrangian’s stage; K3 calculates the time step tn + 1; K4 calculates the new positions of the particles and their integral characteristics at time tn + 1/2; K5 defines the forces on the time layer tn + 1/2; K6 calculates the positions of the particles and their integrated characteristics for the next time layer tn + 1; K7 calculates the flux of physical quantities through the cells boundaries at time moment tn + 1/2; and K8 determines the final hydrodynamic parameters at time moment tn + 1 (see details in Ref. [35]).


4. The set of test problems verifying the numerical models

In the paragraph, the basic set of various tests for one- and two-dimensional problems suitable for revelation of positive and negative sides of the numerical shallow water models is reviewed. According to Toro [26], there are four main test types appropriate to evaluate the numerical solutions (Figure 5).

  1. The comparison with the exact solutions for several 1D problems.

    1. Dam-break waves propagating on the wet bed [16, 36] with the initial distribution of the water depth H(x ≤ 0) = H01 and H(x > 0) = H02, (H01 > 0, H02 > 0). This solution is an analog of the pressure jump decay in an ordinary hydrodynamics. It also contains a hydraulic shock (an analog of a shock wave) and a rarefaction wave for which the exact solutions of the RP for a self-similar type f(xt) = f(ξ) ≡ f(x/t) exist [37]. The bora (hydraulic jump) front for the CSPH-TVD scheme is smeared into three cells and does not contain dispersion oscillations typical for sufficiently good numerical schemes.

    2. Dam-break waves propagating on the dry bed [16]: H(x ≤ 0) = H01 and H(x > 0) = 0. This is an important test imposing special requirements on a numerical scheme for a correct description of the moving boundary between the water and dry bottom. It also allows comparison of the numerical result with exact solution of the RP [37].

    3. The dam-break problem in a sloping channel is more complicated test comparing the dynamics of numerical wetting and drying fronts with the analytic solutions [18].

    4. The two shock waves collision is described in the terms of the RP and contains a rarefaction wave [36].

    5. It is possible to check the coincidence of the numerical and exact solutions for the flow over the bed profiles of different shapes (Figure 5). The transcritical mode with and without a shock for the flow over bed profile in the form of a parabolic bump (see line 3, Figure 5) [16, 17, 18, 38] or triangular obstacle with a break (see line 1) [24] have been studied.

    6. The more complicated bottom profile (see line 2 in Figure 5) helps identifying the best advantages of numerical schemes [22]. We have exact stationary solutions on a stepped bottom for subcritical (Froude number Fr < 1) and supercritical regimes (Fr > 1) [24, 39].

    7. To verify the properties of well-balanced numerical model, the oscillating lake with a heterogeneous bottom [15] or the basin with a stationary fluid with an island in the middle of the numerical domain are suitable to modeling [19]. The study of a propagation of small perturbation in stationary water with a nonuniform bottom allows to determine the quality of the well-balanced scheme [16, 19] and guarantees the absence of a numerical storm produced by the scheme [40].

    8. The same problem of assessing the quality of the WB-scheme is solved by studying the oscillations of the liquid level on the 1D parabolic bottom [24], when the profile of the oscillating free surface of the liquid should remain flat at any time. The two-dimensional analog of such problem is the oscillations in the frictionless paraboloid basin b(xy) = bL(x2 + y2)/L2, for which the free surface remains flat [18, 41]. The analytical solution for the Thacker test with curved water surface in the paraboloidic bottom basin is described in [41]. The initial circular paraboloidic water volume oscillates periodically changing the shape of the free surface in a complex manner (see Figure 5, lines 5a and 5b corresponds to the bottom profile and the initial shape of the water volume) [18].

  2. The second approach is based on the tests with reliable numerical solutions for 1D equations. A simple and effective two-dimensional test is circular dam-break [16, 36], since for the exact solution can be obtained from the one-dimensional radial equation in the polar coordinate system [18]. It is necessary to distinguish the circular dam break on the wet and dry bed.

  3. The third direction is represented by tests using other numerical schemes.

  4. Eventually, the numerical SWM solutions can be compared with the observables data or 3D numerical solutions.

    1. Dam-break flow over the initially dry bed with a bottom obstacle is a tough test for a simple SWM due to the formation of negative wave propagation [42].

    2. The experimental results of the dam-break on a dry bed channel with varying width are used to examine numerical models [16]. The laboratory measurements of flow parameters for the dam-break in the nonprismatic converging-diverging channel are described in Ref. [38].

    3. The so-called right-angled open channel junction test gives a complex two-dimensional flow structure which is highly inhomogeneous in the cross section of the main channel [43]. Such a kind of comparison with essentially 2D currents is an important stage in testing of the numerical model [16].

    4. The dam-break flow over a triangular obstacle contains a large number of characteristic features [38], including the downstream over the dry channel bed, the smooth rarefaction wave, the upstream reflected shock, and the downstream shock after passing the triangular obstacle. The formation of a reflected shock wave from the boundary wall passing back through the hump causes complex nonlinear oscillations in such a basin.

    5. An analytical solution for steady flow case with spatially varying width, bed, and friction slope of special kind is known [38, 44]. There is a singular point in the flow in the transition from the inclined profile to the shallow bottom. A fine structure of the solution in the vicinity of this point is a good test for numerical models.

    6. The use of various software, for example, FLOW-3D, FLOW-3D-SWM, MIS2D, is a fast and effective way for the verification of new numerical models of shallow water [45].

Figure 5.

The bottom profiles with known exact solutions of SWE.

Figure 6 shows the results of the numerical simulations of shear flow instability (Kelvin-Helmholtz instability) in shallow water using the CSPH-TVD method. The characteristic stages of the instability development of the tangential velocity discontinuity in shallow water (see Figure 6a) are obvious: (I) the formation of the eigenmode; (II) the linear stage of instability development, when perturbations increase according to the law ∝ exp(t/T)(T ≈ 9.6 for Fr = 1); (III) the nonlinear mode of perturbations evolution, when the growth of the perturbation amplitude ceases and the characteristic vortex structures are formed (so-called cat eyes, Figure 6b). Thus, in the CSPH-TVD method, there is a lack of important disadvantage of classical SPH algorithm, which does not allow correctly modeling the tangential discontinuity.

Figure 6.

The problem of the tangential velocity discontinuity is represented by: (a) the time dependence of the minimum (shown by solid line) and maximum (shown by dashed line) of the amplitude of the water depth disturbances Hw(t) = H(t) − H0  (H0 ≡ H(0)); (b) the spatial distribution of the vorticity of the velocity vector field ω = ∇ × v/2 at t = 110. The dashed line shown the position of the tangential discontinuity.


5. Engineering applications

We have developed software to solve some engineering applications utilizing the parallel implementation of the CSPH-TVD method on GPUs.

  1. The hydrological regimes of the spring flooding in the Volga-Akhtuba interfluve have been studied thoroughly [1, 35]. The optimal hydrographs of the Volga Hydroelectric Power Station have been constructed. And finally a new approach optimizing the hydrotechnical projects has been developed [46, 47, 48].

  2. Our numerical experiments have reproduced the dynamics of the catastrophic flooding wave in the city of Krymsk area (Russia, Caucasus Mountains) in July 2012, leading to massive losses of life. A number of features of the hydrological regime during the flash flood of 2012 [49], associated with the landscape and the distribution of rainfall has been revealed.

  3. Figure 7 shows the more detailed results of the numerical simulations of the tsunami formation in the Pacific Ocean during an earthquake off the coast of Japan in 2011. A bilateral hydraulic shock (tsunami) is formed due to the displacement of tectonic plates and then the waves propagate to both sides of the discontinuity. The tsunami reaches the shores of Japan in 20–40 min after the earthquake, while the wave height Hw reaches 10–20 m at the coast line.

Figure 7.

Numerical modeling of the tsunami occurred in 2011 off the coast of Japan: (a) the digital terrain model, where the size of the simulation area is 1500 × 1500 cells; the yellow line shows the area of tsunami formation; (b) the 3D wave structure 19 min after the earthquake; frames (c) and (d) are the tsunamis at different moments of time.


6. Conclusion

We have described in detail the numerical method for solving the equations of hydrodynamics in the approximation of the SWM. Our method is a hybrid scheme successfully combining positive properties of Euler’s and Lagrange’s approaches. The CSPH-TVD method allows to model nonstationary flows on a complex heterogeneous bottom relief, containing kinks and sharp jumps of the bottom profile. The numerical scheme is conservative, well balanced and provides a stable through calculation in the presence of non-stationary “water-dry bottom” boundaries on the irregular bottom relief, including the transition through the computing boundary.

On the SPH stage, various smoothing cores can be applied, as well as various TVD-delimiters and methods for the RP solution depending on the features of the problem being solved. The scheme has the second order of convergence on smooth solutions and the first order on discontinuities that corresponds to the accuracy of the Godunov-type schemes. In the case of a non-uniform topography the CSPH-TVD method requires less computational resources than the Godunov’s type schemes, when the WB-condition is necessary to fulfill for the numerical scheme. Comparing with various SPH-method modifications, the numerical CSPH-TVD scheme has higher accuracy and computational speed for the same number of particles; it is less dissipative and better balanced.



The first author has been supported by the Ministry of Education and Science of the Russian Federation (government task No. 2.852.2017/4.6). The second author is thankful to the RFBR (grants 16-07-01037, 15-45-02655, 16-45-340152 and 16-48-340147).


  1. 1. Khrapov S, Pisarev A, Kobelev I, Zhumaliev A, Agafonnikova E, Losev A, Khoperskov A. The numerical simulation of shallow water: Estimation of the roughness coefficient on the flood stage. Advances in Mechanical Engineering. 2013 Id. 787016;5:1-11
  2. 2. Singh J, Altinakar MS, Ding Y. Numerical modeling of rainfall-generated overland flow using nonlinear shallow-water equations. Journal of Hydrologic Engineering. 2015;20 Id. 04014089
  3. 3. Kowalik Z. Introduction to Numerical Modeling of Tsunami Waves. Fairbank: University of Alaska; 2012 196 p
  4. 4. Sepic J, Vilibic I, Fine I. Northern Adriatic meteorological tsunamis: Assessment of their potential through ocean modeling experiments. Journal of Geophysical Research: Oceans. 2015;120:2993-3010
  5. 5. Jeong W. A study on simulation of flood inundation in a coastal urban area using a two-dimensional well-balanced finite volume model. Natural Hazards. 2015;77:337-354
  6. 6. Mangeney A, Bouchut F, Thomas N, Vilotte JP, Bristeau MO. Numerical modeling of self-channeling granular flows and of their levee-channel deposits. Journal of Geophysical Research. 2007 id. F02017;112:1-21
  7. 7. Chertock A, Kurganov A, Qu Z, Wu T. Three-layer approximation of two-layer shallow water equations. Mathematical Modelling and Analysis. 2013;18:675-693
  8. 8. Lauter M, Giraldo FX, Handorf D, Dethloff K. A discontiuous galerkin method for the shallow water equations in spherical triangular coordinates. Journal of Computational Physics. 2008;227:10226-10242
  9. 9. Li X, Chen D, Peng X, Takahashi K. A multimoment finite-volume shallow-water model on the yin–yang overset spherical grid. Monthly Weather Review. 2008;136:3066-3086
  10. 10. Khoperskov AV, Khrapov SS, Nedugova EA. Dissipative-acoustic instability in accretion disks at a nonlinear stage. Astronomy Letters. 2003;29(4):246-257
  11. 11. Regaly Z, Vorobyov E. The circumstellar disk response to the motion of the host star. Astronomy and Astrophysics. 2017;601:A24
  12. 12. Donmez O. The effects of different perturbations on the dynamics of the accretion disk around the black hole. International Journal of Modern Physics A. 2007;22:1875-1898
  13. 13. Nezlin MV, Snezhkin EN. Rossby Vortices, Spiral Structures, Solitons. Springer-Verlag; 1993 223 p
  14. 14. Zasov AV, Saburova AS, Khoperskov AV, Khoperskov SA. Dark matter in galaxies. Physics-Uspekhi. 2017;60:3-39
  15. 15. Audusse E, Bouchut F, Bristeau M-O, Perthame LB. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM Journal on Scientific Computing. 2004;25:2050-2065
  16. 16. Amiri SM, Talebbeydokhti N, Baghlani A. A two-dimensional well-balanced numerical model for shallow water equations. Scientia Iranica. 2013;20:97-107
  17. 17. Audusse E, Bouchut F, Bristeau M-O, Sainte-Marie J. Kinetic entropy inequality and hydrostatic reconstruction scheme for the Saint-Venant system. Mathematics of Computation. 2016;85:2815-2837
  18. 18. Aureli F, Maranzoni A, Mignosa P, Ziveri C. A weighted surface-depth gradient method for the numerical integration of the 2D shallow water equations with topography. Advances in Water Resources. 2008;31:962-974
  19. 19. Vater S, Beisiegel N, Behrens J. A limiter-based well-balanced discontinuous galerkin method for shallow-water flows with wetting and drying: One-dimensional case. Advances in Water Resources. 2015;85:1-13
  20. 20. Zarmehi F, Tavakoli A. A simple scheme to solve saint-venant equations by finite element method. International Journal of Computational Methods. 2016 id. 1650001;13:1-22
  21. 21. Frank J, Gottwald G, Reich S. A hamiltonian particle-mesh method for the rotating shallow-water equations. Lecture Notes in Computational Science and Engineering. 2003;26:131-142
  22. 22. Navas-Montilla A, Murillo J. Overcoming numerical shockwave anomalies using energy balanced numerical schemes. Application to the Shallow Water Equations with discontinuous topography. Journal of Computational Physics. 2017;340:575-616
  23. 23. Zeytounian RK. Nonlinear long waves on water and solutions. Physics-Uspekhi. 1995;38:1333-1381
  24. 24. Khrapov SS, Khoperskov AV, Kuz'min NM, Pisarev AV, Kobelev IA. A numerical scheme for simulating the dynamics of surface water on the basis of the combined SPH-TVD approach. Numerical Methods and Programming. 2011;12:282-297
  25. 25. Monaghan J. Smoothed Particle Hydrodynamics. J. Annual Review of Astronomy and Astrophysics. 1992;30:543-574
  26. 26. Toro EF. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. 3rd ed. Berlin Heidelberg: Springer-Verlag; 2009 721 p
  27. 27. Monaghan J. Smoothed particle hydrodynamics. Journal of Report on Progress in Physics. 2005;68:1703-1759
  28. 28. Desbrun M, Cani M-P. Smoothed particles: A new paradigm for animating highly deformable bodies. In: Proceedings of the Eurographics Workshop on Computer Animation and Simulation. 1996. p. 61-76
  29. 29. Muller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications. Proceedings of 2003 ACM SIGGRAPH Symposium on Computer Animation. 2003:154-159
  30. 30. Harten A. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics. 1983;49:357-393
  31. 31. Harten A, Lax P, van B. On upstream differencing and Godunov type methods for hyperbolic conservation laws. SIAM Review. 1983;25:35-61
  32. 32. van B. Towards the ultimate conservation difference scheme V, a second order sequel to Godunov’s method. Journal of Computational Physics. 1979;32:110-136
  33. 33. Roe PL. Some contributions to the modelling of discontinuous flows. Proceedings of the SIAM/AMS Seminar. 1983:163-193
  34. 34. van Albada GD, van B, Roberts WW. A comparative study of computational methods in cosmic gas dynamics. Astronomy and Astrophysics. 1982;108(1):76-84
  35. 35. Dyakonova T, Khoperskov A, Khrapov S. Numerical model of shallow water: The use of NVIDIA CUDA graphics processors. Communications in Computer and Information Science. 2016;687:132-145
  36. 36. LeVeque RJ. Finite Volume Methods for Hyperbolic Problems. UK: Cambridge University Press; 2004 558 p
  37. 37. Kulikovskii AG, Pogorelov NV, Semenov AY. Mathematical Aspects of Numerical Solution of Hyperbolic Systems. Chapman & Hall/CRC; 2000. p. 560
  38. 38. Alias NA, Liang Q, Kesserwani G. A Godunov-type scheme for modeling 1D channel flow with varying width and topography. Computers & Fluids. 2011;46:88-93
  39. 39. Alcrudo F, Benkhaldoun F. Exact solutions to the Riemann problem of the shallow water equations with a bottom step. Computers & Fluids. 2001;30(6):643-671
  40. 40. Noelle S, Pankratz N, Puppo G, Natvig JR. Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows. Journal of Computational Physics. 2005;213:474-499
  41. 41. Thacker WC. Some exact solutions to the nonlinear shallow-water wave equations. Journal of Fluid Mechanics. 1981;107:499-508
  42. 42. Ozmen-Cagatay H, Kocaman S. Dam-break flow in the presence of obstacle: Experiment and CFD simulation. Engineering Applications of Computational Fluid Mechanics. 2011;5:541-552
  43. 43. Weber LJ, Schumate ED, Mawer N. Experiments on flow at a 90° open-channel junction. ASCE Journal of Hydraulic Engineering. 2001;127(5):340-350
  44. 44. Birman A, Falcovitz J. Application of the GRP scheme to open channel flow equations. Journal of Computational Physics. 2007;222:131-154
  45. 45. Wang, T. and Chu, V., Manning friction in steep open-channel flow. Seventh International Conference on Computational Fluid Dynamics, ICCFD7–3303 (2012), 14 p
  46. 46. Vasilchenko A, Voronin A, Svetlov A, Antonyan N. Assessment of the impact of riverbeds depth in the northern part of the Volga-Akhtuba floodplain on the dynamics of its flooding. International Journal of Pure and Applied Mathematics;110(1):183-192
  47. 47. Voronin AA, Vasilchеnko AA, Pisarеv AV, Khrapov SS, Radchеnko YE. Designing mechanisms of the hydrological regime management of the volga-akhtuba floodplain based on geoinformation and hydrodynamic modeling. Science Journal of VolSU. Mathematics. Physics. 2016;1:24-37
  48. 48. Voronin A, Vasilchenko A, Pisareva M, Pisarev A, Khoperskov A, Khrapov S, Podschipkova J. Designing a system for ecological–economical management of the Volga–Akhtuba floodplain on basis of hydrodynamic and geoinformational simulation. Upravlenie bol'simi sistemami. 2015;55:79-102
  49. 49. Agafonnikova EO, Khoperskov AV, Khrapov SS. The problem of forecasting and control the hydrological regime in the mountainous area in the period flash floods on the basis of hydrodynamical numerical experiments. Kibernetika i programmirovanie. 2016;3:35-53

Written By

Alexander Khoperskov and Sergey Khrapov

Submitted: March 13th, 2017 Reviewed: September 15th, 2017 Published: December 20th, 2017