3. Numerical method СSPH-TVD

also been obtained for the atmospheric phenomena, meteorological forecasts, and global climate

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

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

Our main aim is through calculation for a flow with various Froude number within 0 ≤ Fr < 100

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

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

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

detailed description of the numerical scheme is the main aim of the chapter.

gH p is an analog of the sound speed, H is the depth, u is the flow

play the topography role in such astrophysical problems.

models [8, 9].

238 Numerical Simulations in Engineering and Science

(WB) circuits [15–17].

(Fr <sup>=</sup> <sup>u</sup>/cg, where cg <sup>¼</sup> ffiffiffiffiffiffi

wet and dry bottom [1, 18, 19].

2. Mathematical models

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 τ<sup>n</sup> = tn � tn � 1, where n denotes the index of the time layer. The vector space index i = (i, j) characterizes the radius vectors of the fixed centers of the cells r0 <sup>i</sup> <sup>¼</sup> <sup>x</sup><sup>0</sup> <sup>i</sup> ; y<sup>0</sup> j � �. The total number of cells and particles is Np <sup>=</sup> 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 V<sup>n</sup> <sup>i</sup> <sup>¼</sup> ÐÐ <sup>S</sup>ið Þ<sup>t</sup> HdS, momentumP<sup>n</sup> <sup>i</sup> <sup>¼</sup> ÐÐ <sup>S</sup>ið Þ<sup>t</sup> <sup>H</sup>vdS, coordinates of their centers <sup>r</sup><sup>n</sup> <sup>i</sup> , and linear sizes ℓ n <sup>i</sup> <sup>¼</sup> ffiffiffiffiffi Sn i p .

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

Sið Þt

dUi

� �, <sup>∇</sup>i<sup>η</sup> <sup>¼</sup> <sup>∇</sup>H<sup>j</sup>

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

dri

being the analogous of the function mean values φ = {H, Hu, Hv} in cells in the finite volume approximation on a fixed grid. Taking into the account Eq. (2), the integral characteristics of the

, PiðÞ¼ t ð Þ Hv <sup>i</sup>

φð Þ t; x; y dxdy (2)

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

http://dx.doi.org/10.5772/intechopen.71026

241

dt <sup>¼</sup> <sup>Φ</sup><sup>i</sup> (4)

<sup>r</sup>¼ri <sup>þ</sup> <sup>∇</sup>bj<sup>r</sup>¼r<sup>0</sup>

dt <sup>¼</sup> vi (5)

ð Þ <sup>N</sup> Xx;Ny

k¼ð Þ 1;1

ð Þ<sup>t</sup> <sup>h</sup><sup>2</sup> (3)

i

ψð Þ t;rk ∇Wð Þ jr � rkj; h (6)

, σ<sup>i</sup> = σ(ri), fi = f

"water-dry bottom" boundaries in the computational domain.

φi ðÞ¼ t 1 h2 ðð

<sup>V</sup>iðÞ¼ <sup>t</sup> <sup>H</sup>ið Þ<sup>t</sup> <sup>h</sup><sup>2</sup>

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

�gH<sup>i</sup> ∇iη þ Hifi

<sup>ψ</sup>ð Þ <sup>t</sup>;rk <sup>W</sup>ð Þ <sup>j</sup><sup>r</sup> � rkj; <sup>h</sup> , <sup>∇</sup>φbð Þ¼ <sup>t</sup>;<sup>r</sup>

3.1. The Lagrangian's stage

Let us define the following quantities

particles can be written in the form:

ð Þ Hv <sup>i</sup>

<sup>φ</sup>bð Þ¼ <sup>t</sup>;<sup>r</sup>

(ri). The law of motion of a particle is

� �, <sup>Φ</sup><sup>i</sup> <sup>¼</sup> <sup>σ</sup><sup>i</sup>

where vi = (Hv)i/H<sup>i</sup> is the velocity of the i-th particle.

ð Þ <sup>N</sup> Xx;Ny

k¼ð Þ 1;1

where Ui <sup>¼</sup> <sup>H</sup><sup>i</sup>

The CSPH-TVD method main stages are the following:


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 x<sup>n</sup> <sup>i</sup> ! <sup>x</sup><sup>n</sup>þ<sup>1</sup> <sup>i</sup> , 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

r0 <sup>i</sup> <sup>¼</sup> <sup>x</sup><sup>0</sup> <sup>i</sup> ; y<sup>0</sup> j � �

volume V<sup>n</sup>

(V<sup>n</sup>

sizes ℓ n <sup>i</sup> <sup>¼</sup> ffiffiffiffiffi Sn i p .

<sup>i</sup> <sup>¼</sup> ÐÐ

<sup>i</sup> ! <sup>V</sup><sup>n</sup>þ<sup>1</sup>

the initial state (r<sup>n</sup>

<sup>i</sup> , P<sup>n</sup>

240 Numerical Simulations in Engineering and Science

<sup>S</sup>ið Þ<sup>t</sup> HdS, momentumP<sup>n</sup>

The CSPH-TVD method main stages are the following:

<sup>i</sup> , and r<sup>n</sup>

<sup>i</sup> and ℓ n <sup>i</sup> ¼ h).

<sup>i</sup> ! <sup>P</sup><sup>n</sup>þ<sup>1</sup>

<sup>i</sup> � <sup>r</sup><sup>0</sup>

to the trajectories of the particles due to their displacement x<sup>n</sup>

cells centers (it is assumed that ℓ

. The total number of cells and particles is Np = Nx � Ny. At the initial time moment,

<sup>S</sup>ið Þ<sup>t</sup> <sup>H</sup>vdS, coordinates of their centers <sup>r</sup><sup>n</sup>

<sup>i</sup> ) due to the hydrodynamic and external forces

<sup>i</sup> ! <sup>x</sup><sup>n</sup>þ<sup>1</sup> <sup>i</sup> , the dashed black lines show the boundaries of cells.

<sup>i</sup> 6¼ h). The difference between the CSPH-TVD method

<sup>i</sup> , and linear

the particles are placed in the centers of cells. The main characteristics of the particles are the

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

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

and the traditional SPH approach [25] consists in the fact that the particles at each time interval [tn, tn + 1] of the numerical time integration at the moment of time tn are found in

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

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.

<sup>i</sup> <sup>¼</sup> ÐÐ

<sup>i</sup> ! r nþ1

nþ1

Let us define the following quantities

$$\varphi\_{\mathbf{i}}(t) = \frac{1}{h^2} \iint\_{S\_{\mathbf{i}}(t)} \varphi(t, \mathbf{x}, y) d\mathbf{x} dy \tag{2}$$

being the analogous of the function mean values φ = {H, Hu, Hv} 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:

$$V\_{\mathbf{i}}(t) = H\_{\mathbf{i}}(t)h^2, \mathbf{P}\_{\mathbf{i}}(t) = (H\mathbf{v})\_{\mathbf{i}}(t)h^2 \tag{3}$$

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

$$\frac{d\mathbf{U}\_i}{dt} = \mathbf{O}\_i \tag{4}$$

where Ui <sup>¼</sup> <sup>H</sup><sup>i</sup> ð Þ Hv <sup>i</sup> � �, <sup>Φ</sup><sup>i</sup> <sup>¼</sup> <sup>σ</sup><sup>i</sup> �gH<sup>i</sup> ∇iη þ Hifi � �, <sup>∇</sup>i<sup>η</sup> <sup>¼</sup> <sup>∇</sup>H<sup>j</sup> <sup>r</sup>¼ri <sup>þ</sup> <sup>∇</sup>bj<sup>r</sup>¼r<sup>0</sup> i , σ<sup>i</sup> = σ(ri), fi = f (ri). The law of motion of a particle is

$$\frac{d\mathbf{r\_i}}{dt} = \mathbf{v\_i} \tag{5}$$

where vi = (Hv)i/H<sup>i</sup> 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 Ω:

$$\hat{\boldsymbol{\varphi}}(t,\mathbf{r}) = \sum\_{\mathbf{k}=(1,1)}^{\left(N\_{\mathbf{r}}N\_{\mathbf{y}}\right)} \boldsymbol{\psi}(t,\mathbf{r\_{k}}) \, \mathcal{W}(|\mathbf{r}-\mathbf{r\_{k}}|,h) \, \nabla \hat{\boldsymbol{\varphi}}(t,\mathbf{r}) = \sum\_{\mathbf{k}=(1,1)}^{\left(N\_{\mathbf{r}}N\_{\mathbf{y}}\right)} \boldsymbol{\psi}(t,\mathbf{r\_{k}}) \nabla \mathcal{W}(|\mathbf{r}-\mathbf{r\_{k}}|,h) \tag{6}$$

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

$$\psi(t, \mathbf{r\_k}) = \iint\_{\mathbf{S\_k}(t)} \varphi(t, \mathbf{r}) d\mathbf{S} = \varphi\_\mathbf{k}(t) h^2 \tag{7}$$

W qð Þ¼ Bw

<sup>8</sup><sup>h</sup> ; <sup>5</sup> <sup>16</sup><sup>π</sup> <sup>h</sup><sup>2</sup> ; <sup>15</sup> 64π h<sup>3</sup>

nonphysical clusters of particles.

tn + 1 = tn + τn:

(tn + 1):

where Bw <sup>¼</sup> <sup>1</sup>

ð Þ 2 � q 3 , 0 ≤ q ≤ 2 q > 2

, W<sup>0</sup>

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

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:

I. At the predictor step, the water depth H<sup>i</sup> and velocity vi are calculated at time moment

τn

III. During the corrector step, the water depth H<sup>i</sup> and velocity vi values are recalculated

Figure 2. The spatial distributions of smoothing kernels (a) and their derivatives (b). The smoothing kernels (9) and (10)

<sup>i</sup> ð Þ tnþ<sup>1</sup> 2 þ

<sup>2</sup> við Þþ tn <sup>v</sup><sup>∗</sup>

τn

ð Þ¼� q

� � is the normalization constant. Figure 2 shows the comparison

Bw h

8 < :

3 2ð Þ � q 2 , 0 ≤ q ≤ 2

http://dx.doi.org/10.5772/intechopen.71026

q > 2 , (10)

243

0,

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

<sup>i</sup> ð Þ¼ tnþ<sup>1</sup> Uið Þþ tn τnΦið Þ rið Þ tn ; Uið Þ tn (11)

<sup>i</sup>ð Þ tnþ<sup>1</sup>

<sup>2</sup> <sup>Φ</sup><sup>i</sup> rið Þ tnþ<sup>1</sup> ; <sup>U</sup><sup>∗</sup>

� � (12)

<sup>i</sup>ð Þ tnþ<sup>1</sup>

� � (13)

0,

U∗

<sup>U</sup><sup>~</sup> <sup>i</sup>ð Þ¼ tnþ<sup>1</sup>

are shown by solid and dashed lines, correspondingly.

II. The Particles' spatial positions ri are determined at time tn + 1:

rið Þ¼ tnþ<sup>1</sup> rið Þþ tn

Uið Þþ tn <sup>U</sup><sup>∗</sup>

8 < :

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 <sup>ψ</sup>ð Þ¼ <sup>t</sup>;rk <sup>V</sup>kð Þ<sup>t</sup> <sup>φ</sup><sup>b</sup> <sup>k</sup>ð Þ<sup>t</sup> <sup>=</sup>H<sup>b</sup> <sup>k</sup>ð Þ<sup>t</sup> 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

$$\mathbf{OP}\_{\mathbf{i}} \approx \left( -gH\_{\mathbf{i}} \sum\_{\mathbf{k}} \left( H\_{\mathbf{k}} \nabla\_{\mathbf{i}} \overline{W}\_{\mathbf{i}\mathbf{k}} + b\_{\mathbf{k}} \nabla\_{\mathbf{i}} \overline{W}\_{\mathbf{i}\mathbf{k}}^{0} \right) + H\_{\mathbf{i}} \mathbf{f}\_{\mathbf{i}} \right) \tag{8}$$

where <sup>W</sup>ik <sup>¼</sup> <sup>h</sup><sup>2</sup> <sup>W</sup>ð Þ <sup>j</sup>ri � rkj; <sup>h</sup> , <sup>W</sup><sup>0</sup> ik <sup>¼</sup> <sup>h</sup><sup>2</sup> <sup>W</sup> <sup>j</sup>r<sup>0</sup> <sup>i</sup> � <sup>r</sup><sup>0</sup> <sup>k</sup>j; <sup>h</sup> � � 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:


3. the zero spatial step limit lim<sup>h</sup>!<sup>0</sup> W jr � r<sup>0</sup> ð Þ¼ j; h δ jr � r<sup>0</sup> ð Þ j; h 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

$$\mathcal{W}(q) = A\_w \begin{cases} 1 - 1.5q^2 + 0.75q^3, & 0 \le q \le 1 \\ 0.25 \left( 2 - q \right)^3, & 1 \le q \le 2, \mathcal{W}'(q) = -\frac{A\_w}{h} \\ 0, & q > 2 \end{cases} \begin{cases} 3q - 2.25q^2, & 0 \le q \le 1 \\ 0.75 \left( 2 - q \right)^2, & 1 \le q \le 2 \\ 0, & q > 2 \end{cases} \tag{9}$$

is the most commonly applied approximation, where q = ∣ ri � rk ∣ /h is the relative distance between the particles, Aw <sup>¼</sup> <sup>2</sup> <sup>3</sup><sup>h</sup> ; <sup>10</sup> <sup>7</sup><sup>π</sup> <sup>h</sup><sup>2</sup> ; <sup>1</sup> π h<sup>3</sup> � � 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 Aw ≈ 0:987 <sup>10</sup> 7π h<sup>2</sup> � �. 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]:

A Numerical Simulation of the Shallow Water Flow on a Complex Topography http://dx.doi.org/10.5772/intechopen.71026 243

$$\mathcal{W}(q) = B\_w \begin{cases} \left(2 - q\right)^3, & 0 \le q \le 2 \\ 0, & q > 2 \end{cases} \\ \mathcal{W}'(q) = -\frac{B\_w}{h} \begin{cases} 3(2 - q)^2, & 0 \le q \le 2 \\ 0, & q > 2 \end{cases} \tag{10}$$

where Bw <sup>¼</sup> <sup>1</sup> <sup>8</sup><sup>h</sup> ; <sup>5</sup> <sup>16</sup><sup>π</sup> <sup>h</sup><sup>2</sup> ; <sup>15</sup> 64π h<sup>3</sup> � � 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.

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:

I. At the predictor step, the water depth H<sup>i</sup> and velocity vi are calculated at time moment tn + 1 = tn + τn:

$$\mathbf{U\_i^\*(t\_{n+1}) = \mathbf{U\_i(t\_n)} + \tau\_n \mathbf{O\_i(r\_i(t\_n), \mathbf{U\_i(t\_n)})} \tag{11}$$

II. The Particles' spatial positions ri are determined at time tn + 1:

ψð Þ¼ t;rk

expressions for system (4), we obtain

242 Numerical Simulations in Engineering and Science

where <sup>W</sup>ik <sup>¼</sup> <sup>h</sup><sup>2</sup> <sup>W</sup>ð Þ <sup>j</sup>ri � rkj; <sup>h</sup> , <sup>W</sup><sup>0</sup>

1. the kernel is finiteness;

W qð Þ¼ Aw

3. the zero spatial step limit lim<sup>h</sup>!<sup>0</sup>

2. the normalization condition is fulfilled

<sup>1</sup> � <sup>1</sup>:5q<sup>2</sup> <sup>þ</sup> <sup>0</sup>:75q<sup>3</sup>,

3 ,

<sup>3</sup><sup>h</sup> ; <sup>10</sup> <sup>7</sup><sup>π</sup> <sup>h</sup><sup>2</sup> ; <sup>1</sup> π h<sup>3</sup> � �

0:25 2ð Þ � q

7π h<sup>2</sup> � �

0,

8 >>>><

>>>>:

between the particles, Aw <sup>¼</sup> <sup>2</sup>

chosen to be Aw ≈ 0:987 <sup>10</sup>

Φ<sup>i</sup> ≈

0 @

�gH<sup>i</sup>

X k

ik <sup>¼</sup> <sup>h</sup><sup>2</sup> <sup>W</sup> <sup>j</sup>r<sup>0</sup>

stage and the fulfillment of the WB-condition on the inhomogeneous topography.

There are three conditions superimposed on the smoothing kernels W:

ðð

<sup>φ</sup>ð Þ <sup>t</sup>;<sup>r</sup> dS <sup>¼</sup> <sup>φ</sup>kð Þ<sup>t</sup> <sup>h</sup><sup>2</sup> (7)

Skð Þt

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 <sup>ψ</sup>ð Þ¼ <sup>t</sup>;rk <sup>V</sup>kð Þ<sup>t</sup> <sup>φ</sup><sup>b</sup> <sup>k</sup>ð Þ<sup>t</sup> <sup>=</sup>H<sup>b</sup> <sup>k</sup>ð Þ<sup>t</sup> is used instead of (7). Substituting the SPH-approximation (6) and (7) in the right-hand side of the

σi

<sup>i</sup> � <sup>r</sup><sup>0</sup>

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

ðð

Ω

used in the SPH hydrodynamic simulations [25, 28, 29]. The Monaghan's cubic spline

0 ≤ q ≤ 1 1 ≤ q ≤ 2 q > 2

For the smoothing kernel W, spline functions of different orders or Gaussian distribution are

, W<sup>0</sup>

is the most commonly applied approximation, where q = ∣ ri � rk ∣ /h is the relative distance

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

particles in the simulation of hydrodynamic flows by the SPH method; therefore, an approxi-

mation for the smoothing kernel of the pressure gradient is applied in form [29]:

ð Þ¼� q

<sup>H</sup>k∇iWik <sup>þ</sup> <sup>b</sup>k∇iW<sup>0</sup>

� �

ik

W jr � r<sup>0</sup> ð Þ j; h dr<sup>0</sup> ¼ 1;

W jr � r<sup>0</sup> ð Þ¼ j; h δ jr � r<sup>0</sup> ð Þ j; h is the Dirac delta-function.

Aw h

8 >>>><

>>>>:

. The smoothing kernel (9) can lead to unphysical clustering of

<sup>3</sup><sup>q</sup> � <sup>2</sup>:25q<sup>2</sup>,

2 , 0 ≤ q ≤ 1

1 ≤ q ≤ 2

(9)

q > 2

0:75 2ð Þ � q

0,

is the normalization constant for the 1D, 2D and 3D

þ Hifi

<sup>k</sup>j; <sup>h</sup> � � are the functions computed in the cells

1

A (8)

$$\mathbf{r\_i(t\_{n+1})} = \mathbf{r\_i(t\_n)} + \frac{\tau\_n}{2} \left( \mathbf{v\_i(t\_n)} + \mathbf{v\_i^\*(t\_{n+1})} \right) \tag{12}$$

III. During the corrector step, the water depth H<sup>i</sup> and velocity vi values are recalculated (tn + 1):

$$
\hat{\mathbf{U}}\_{\mathbf{i}}(t\_{n+1}) = \frac{\mathbf{U}\_{\mathbf{i}}(t\_n) + \mathbf{U}\_{\mathbf{i}}^\*(t\_{n+1})}{2} + \frac{\tau\_n}{2} \mathbf{O}\_{\mathbf{i}}\left(\mathbf{r}\_{\mathbf{i}}(t\_{n+1}), \mathbf{U}\_{\mathbf{i}}^\*(t\_{n+1})\right) \tag{13}
$$

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.

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 <sup>ξ</sup>iðÞ¼ <sup>t</sup> rið Þ� <sup>t</sup> <sup>r</sup><sup>0</sup> i (here r<sup>0</sup> <sup>i</sup> 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~ <sup>n</sup>þ<sup>1</sup> <sup>i</sup> and U<sup>n</sup>þ<sup>1</sup> <sup>i</sup> 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 r<sup>0</sup> <sup>i</sup> are computed. In order to do it, the following expression

$$\frac{d}{dt} \iint\_{\mathbf{S}(t)} \mathbf{U} dx dy = \iint\_{\mathbf{S}(t)} \left( \frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} \right) dx dy, \text{ where } \mathbf{F} = \left( \begin{array}{c} Hv \\ Hu^2 \\ Hvu \end{array} \right), \text{ } \mathbf{G} = \left( \begin{array}{c} Hv \\ Hw \\ Hv^2 \end{array} \right)$$

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

$$\mathbf{U\_{i}^{n+1}} = \mathbf{\tilde{U}\_{i}^{n+1}} - \frac{\tau\_{n}}{h} \left( \mathbf{F}\_{i+1/2,j}^{n+1/2} - \mathbf{F}\_{i-1/2,j}^{n+1/2} + \mathbf{G}\_{i,j+1/2}^{n+1/2} - \mathbf{G}\_{i,j-1/2}^{n+1/2} \right) \tag{14}$$

UL

<sup>i</sup>þ1=<sup>2</sup> <sup>¼</sup> <sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup>

i,j þ

To fulfill this condition, the limiter function

Θ<sup>n</sup>þ1=<sup>2</sup> i,j ¼ Λ

numerical solution in the vicinity of discontinuities.

where 0 < K < 1 is the Courant number.

h <sup>2</sup> � <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup> i � �Θ<sup>n</sup>þ1=<sup>2</sup>

the mass centers of the particles for the CSPH-TVD scheme at time tn + 1/2.

i,j , <sup>U</sup><sup>R</sup>

Figure 3. The projection onto the Ox axis of the piecewise linear reconstruction of the function U(t,r) with the respect to

The slopes of the piecewise linear distribution (16) should satisfy the TVD-condition [30].

<sup>i</sup>þ<sup>1</sup> � <sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup> i

> <sup>i</sup>þ<sup>1</sup> � <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup> i

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

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 τ<sup>n</sup> of the CSPH-TVD algorithm:

> ; <sup>h</sup> max<sup>i</sup> <sup>j</sup>v<sup>n</sup>

!

<sup>i</sup> j þ ffiffiffiffiffiffiffiffiffi gH<sup>n</sup> i

� � p

<sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup>

0 B@

3.3. The stability condition and parallel implementation on GPUs

<sup>τ</sup><sup>n</sup> <sup>¼</sup> <sup>K</sup>min <sup>h</sup>

2 max<sup>i</sup> ∣v<sup>n</sup> i ∣

<sup>h</sup> <sup>þ</sup> <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>i</sup>þ1=<sup>2</sup> <sup>¼</sup> <sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup>

, <sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup>

<sup>h</sup> <sup>þ</sup> <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>i</sup>þ1,j � <sup>h</sup>

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

http://dx.doi.org/10.5772/intechopen.71026

245

<sup>i</sup> � <sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup> i�1

<sup>i</sup> � <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup> i�1

<sup>2</sup> <sup>þ</sup> <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup> iþ1 � �Θ<sup>n</sup>þ1=<sup>2</sup>

> 1 CA

<sup>i</sup>þ1,j (16)

(17)

where F<sup>n</sup>þ1=<sup>2</sup> <sup>i</sup>�1=2,j and <sup>G</sup><sup>n</sup>þ1=<sup>2</sup> <sup>i</sup>�1=2,j are the average values (in the interval [tn, tn + 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 U<sup>L</sup> and right U<sup>R</sup> 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(t, r). 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(t, r) 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(t, r) inside of the i-th cell at time tn + 1/2 in the projection onto the Ox axis (Figure 3) has the following form

$$\mathbf{U}\{t\_{n+1/2},\mathbf{x}\} = \tilde{\mathbf{U}}\_{i,j}^{n+1/2} + \left(\mathbf{x} - \mathbf{x}\_i^{n+1/2}\right) \Theta\_{i,j}^{n+1/2}, \quad \mathbf{x} \in \left[\mathbf{x}\_{i-1/2}^0, \mathbf{x}\_{i+1/2}^0\right] \tag{15}$$

where Θ<sup>n</sup>þ1=<sup>2</sup> i,j is the vector of slopes (angular coefficients), x nþ1=2 <sup>i</sup> <sup>¼</sup> <sup>x</sup><sup>0</sup> <sup>i</sup> <sup>þ</sup> <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup> <sup>i</sup> . Taking into the account Eq. (15), the expressions for U<sup>L</sup> and URat the boundary x<sup>0</sup> <sup>i</sup>þ1=<sup>2</sup> can be written as

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.

$$\mathbf{U}\_{i+1/2}^{L} = \tilde{\mathbf{U}}\_{i,j}^{n+1/2} + \left(\frac{\hbar}{2} - \xi\_{i}^{n+1/2}\right) \Theta\_{i,j}^{n+1/2}, \\ \mathbf{U}\_{i+1/2}^{R} = \tilde{\mathbf{U}}\_{i+1,j}^{n+1/2} - \left(\frac{\hbar}{2} + \xi\_{i+1}^{n+1/2}\right) \Theta\_{i+1,j}^{n+1/2} \tag{16}$$

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

$$
\boldsymbol{\Theta}\_{i,j}^{n+1/2} = \Lambda \left( \frac{\boldsymbol{\tilde{\mathbf{U}}}\_{i+1}^{n+1/2} - \boldsymbol{\tilde{\mathbf{U}}}\_{i}^{n+1/2}}{h + \boldsymbol{\xi}\_{i+1}^{n+1/2} - \boldsymbol{\xi}\_{i}^{n+1/2}}, \frac{\boldsymbol{\tilde{\mathbf{U}}}\_{i}^{n+1/2} - \boldsymbol{\tilde{\mathbf{U}}}\_{i-1}^{n+1/2}}{h + \boldsymbol{\xi}\_{i}^{n+1/2} - \boldsymbol{\xi}\_{i-1}^{n+1/2}} \right)
$$

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 τ<sup>n</sup> of the CSPH-TVD algorithm:

$$\tau\_{\boldsymbol{\pi}} = \text{Kmin}\left(\frac{h}{2\max\_{\mathbf{i}}|\mathbf{v}\_{\mathbf{i}}^{\boldsymbol{n}}|}, \frac{h}{\max\_{\mathbf{i}}\left(|\mathbf{v}\_{\mathbf{i}}^{\boldsymbol{n}}| + \sqrt{gH\_{\mathbf{i}}^{\boldsymbol{n}}}\right)}\right) \tag{17}$$

where 0 < K < 1 is the Courant number.

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 <sup>ξ</sup>iðÞ¼ <sup>t</sup> rið Þ� <sup>t</sup> <sup>r</sup><sup>0</sup>

moment tn + 1. According to the CSPH-TVD approach, the changes in the particles integral

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

<sup>i</sup>þ1=2,j � <sup>F</sup><sup>n</sup>þ1=<sup>2</sup>

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

The values of the flow parameters to the left U<sup>L</sup> and right U<sup>R</sup> 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(t, r). 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(t, r) 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(t, r) inside of the i-th cell at time tn + 1/2 in the

> nþ1=2 i � �

Θ<sup>n</sup>þ1=<sup>2</sup>

and URat the boundary x<sup>0</sup>

i,j , x<sup>∈</sup> <sup>x</sup><sup>0</sup>

nþ1=2 <sup>i</sup> <sup>¼</sup> <sup>x</sup><sup>0</sup>

<sup>i</sup>�1=<sup>2</sup>; <sup>x</sup><sup>0</sup>

iþ1=2 h i

<sup>i</sup>þ1=<sup>2</sup> can be written as

<sup>i</sup> . Taking into the

<sup>i</sup> <sup>þ</sup> <sup>ξ</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>h</sup> <sup>F</sup><sup>n</sup>þ1=<sup>2</sup>

dxdy, where F ¼

<sup>i</sup>�1=2,j <sup>þ</sup> <sup>G</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>i</sup>�1=2,j are the average values (in the interval [tn, tn + 1]) of mass and

� �

<sup>i</sup> and U<sup>n</sup>þ<sup>1</sup>

Hv Hu<sup>2</sup>

0

BB@

Hvu

i,jþ1=<sup>2</sup> � <sup>G</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>i</sup> are computed. In order to do it, the following

1

CCA, <sup>G</sup> <sup>¼</sup>

i,j�1=2

<sup>i</sup> is the initial position) should be used instead of ri(t).

At this stage, the relationship between the values of U~ <sup>n</sup>þ<sup>1</sup>

∂U ∂t þ ∂F ∂x þ ∂G ∂y

<sup>i</sup> � <sup>τ</sup><sup>n</sup>

projection onto the Ox axis (Figure 3) has the following form

i,j þ x � x

i,j is the vector of slopes (angular coefficients), x

<sup>U</sup> tnþ1=<sup>2</sup>; <sup>x</sup> � � <sup>¼</sup> <sup>U</sup><sup>~</sup> <sup>n</sup>þ1=<sup>2</sup>

account Eq. (15), the expressions for U<sup>L</sup>

� �

characteristics at their return to the initial state r<sup>0</sup>

ðð

Sið Þt

<sup>i</sup> <sup>¼</sup> <sup>U</sup><sup>~</sup> <sup>n</sup>þ<sup>1</sup>

U<sup>n</sup>þ<sup>1</sup>

(here r<sup>0</sup>

expression

d dt ðð

where F<sup>n</sup>þ1=<sup>2</sup>

where Θ<sup>n</sup>þ1=<sup>2</sup>

Sið Þt

Udxdy ¼

244 Numerical Simulations in Engineering and Science

<sup>i</sup>�1=2,j and <sup>G</sup><sup>n</sup>þ1=<sup>2</sup>

boundary of the Euler's cells.

3.2. The Euler's stage

i

<sup>i</sup> is determined at time

Hv Huv Hv<sup>2</sup>

1

CCA

(14)

(15)

0

BB@

correct description of the moving boundary between the water and dry bottom. It also

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

http://dx.doi.org/10.5772/intechopen.71026

247

allows comparison of the numerical result with exact solution of the RP [37].

rarefaction wave [36].

been studied.

(Fr > 1) [24, 39].

basin b(x, y) = bL(x<sup>2</sup> + y<sup>2</sup>

water volume) [18].

3D numerical solutions.

dry bed.

c. 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]. d. The two shock waves collision is described in the terms of the RP and contains a

e. 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–18, 38] or triangular obstacle with a break (see line 1) [24] have

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

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

h. 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 twodimensional analog of such problem is the oscillations in the frictionless paraboloid

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

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

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

simple SWM due to the formation of negative wave propagation [42].

a. Dam-break flow over the initially dry bed with a bottom obstacle is a tough test for a

, for which the free surface remains flat [18, 41]. The

absence of a numerical storm produced by the scheme [40].

)/L<sup>2</sup>

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

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]).

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.
