Analysis of a new Kolgan-type scheme motivated by the shallow water equations

Applied Numerical Mathematics 62 (2012) 489–506

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Analysis of a new Kolgan-type scheme motivated by the shallow water equations M. Elena Vázquez-Cendón a,∗ , Luis Cea b a b

Departamento de Matemática Aplicada, Universidade de Santiago de Compostela, Spain Grupo de Ingeniería del Agua y del Medio Ambiente ETSICCP, Universidade da Coruña, Spain

a r t i c l e

i n f o

Article history: Available online 15 June 2011 Keywords: Shallow water equations Finite volume method Upwind sources C -property High-order schemes Kolgan scheme

a b s t r a c t This paper presents the analysis of two different finite volume schemes for hyperbolic conservation laws: the Kolgan high-resolution scheme, and a new Kolgan-type scheme in which the high-order extrapolation of the conservative variables is used just in the upwind contribution of the numerical flux and source terms. Both schemes are compared in terms of the local truncation error, the stability conditions and the C -property. The schemes are applied to different hyperbolic conservation equations, including the one-dimensional scalar transport equation, the Burgers equation and the 2D shallow water equations, in order to compute the observed order of accuracy and to verify the C -property. When applied to the 2D shallow water equations, the new approach avoids spurious oscillations in the solution without the need of using high-order corrections in the definition of the bed slope source term. © 2011 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The shallow water equations are widely used to model free surface flow in rivers and coastal regions. The numerical solution of this system of hyperbolic conservation equations has motivated the development of a large number of highorder finite volume schemes for hyperbolic conservation laws [1,3,14,15]. The discretisation of non-conservative source terms in the shallow water equations, and specially the bed slope source term, has been also subject of study in previous publications [10,6,19]. This paper presents the analysis of two different finite volume discretisation schemes for hyperbolic conservation laws: the Kolgan high-resolution scheme,1 which uses a second-order spatial discretisation and a first-order time discretisation [18]; and a new Kolgan-type scheme in which the high-order extrapolation of the conservative variables is used just in the upwind contribution of the numerical flux and source terms. Following the appraisal of the ideas of V. Kolgan outside the USSR (as it has been done by van Leer [18] and Toro [17]), the new scheme will be called Kolgantype scheme. This new scheme was proposed in [4] in order to discretise the 2D shallow water equations with porosity, showing a well-balanced behaviour, in the sense that it reproduces the exact hydrostatic solution without the need of using high-order corrections in the source terms discretisation. The scheme proposed in [4] can be applied to the classic shallow water equations without porosity by just setting the porosity value equal to one in the whole spatial domain. This way of proceeding allows us to achieve a high-order scheme for the classic shallow water equations without the need of highorder corrections in the bed slope discretisation, as it is needed in other numerical discretisations [5,8], which simplifies the numerical implementation of the schemes.

*

Corresponding author. E-mail addresses: [email protected] (M.E. Vázquez-Cendón), [email protected] (L. Cea). 1 Recently B. van Leer in [18] presented a historical oversight to honour Kolgans’ contribution to CFD, previous to the English translation of the 1972 original paper [18]. 0168-9274/$30.00 © 2011 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2011.06.002

490

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

The idea of the new discretisation is to use high-order extrapolations just in the upwind contribution of the numerical flux and source terms. This is justified because a centred discretisation of the convective flux, which is numerically unstable, is of second-order accuracy in space. It is the numerical diffusion introduced by the upwind contribution which stabilises the numerical method and reduces the accuracy of the scheme to the first order. The high-order extension aims to reduce the numerical diffusion and therefore, it only needs to be applied to the upwind contribution in the numerical flux and source terms. As it is shown numerically in [4], in the case of the 2D shallow water equations this approach avoids spurious oscillations in the free surface elevation without the need of using high-order corrections in the definition of the bed slope source term. In this paper several test cases are used in order to analyse the accuracy and numerical stability of the proposed scheme. These test cases include the scalar transport equation and the Burgers equation with source terms, as classical reference problems, in order to prove that this numerical scheme preserves the second order. A specific analysis of the local truncation error, stability conditions and C -property of different discretisations of the flux and source term are also considered in order to state the properties of the numerical method considered. The outline of the paper is as follows. In Section 2 the transport equation without source terms is considered as a reference equation in order to analyse the numerical flux of the Kolgan method and the new Kolgan-type scheme. The transport equation with source terms is used to analyse the numerical discretisation of the source terms in Section 3. A non-linear problem defined by Burgers equation with source term, is solved using both numerical schemes in Section 4. Finally, in Section 5 the schemes are used to solve the 2D shallow water equations including viscous stresses, as a motivation test case for the future developments of the scheme. 2. Analysis of a new Kolgan-type method for a scalar model without source terms As it has been mentioned, the new second-order in space numerical scheme which is analysed in this paper was presented in [4] for the 2D shallow water equations with porosity. The scheme uses the high-order extrapolation of the conservative variables just in the upwind contribution of the numerical flux, but not in the centred contribution. The order of accuracy of the resulting scheme was not proven in [4] using numerical analysis, due the difficulty which implies the demonstration for the 2D shallow water equations with porosity. In this work, we will study the order of accuracy and stability of the proposed scheme following the methodology which is normally used in the numerical analysis of hyperbolic equations. 2.1. Data reconstruction and numerical flux for a scalar problem without sources The well-known scalar transport equation is given by:

∂ f (w) ∂w (x, t ) + (x, t ) = 0, ∂t ∂x

f (w) = λw,

(1)

where w is the conservative variable, and f ( w ) is the physical flux of w. A conservative numerical scheme for Eq. (1), can be written in terms of a numerical flux φ as:

w nj +1 = w nj −

  t   n n  φ w j , w j +1 − φ w nj−1 , w nj . x

(2)

We will start comparing a family of first-order upwind schemes, the so-called Q -schemes, with the new proposed scheme. The numerical flux of the Q -schemes can be written as:

φ(u , v ) =

f (u ) + f ( v ) 2

 1 −  Q (u , v )( v − u ), 2





Q (u , v ) = A w Q (u , v ) ,

(3)

where A for a scalar conservation law is A ( w ) = f  ( w ), and w Q is the mean state of the corresponding Q -scheme, which depends on the variables u , v (the arithmetic mean for the van Leer Q -scheme and the Roe-mean for the Roe scheme). In the linear scalar case Q (u , v ) = λ. Assuming, without loss of generality, that λ  0, the numerical flux can be written as:

φ(u , v ) = λ

u+v 2

1

− λ( v − u ) = λu . 2

(4)

Using Eq. (2), the first-order upwind scheme for the scalar transport equation is obtained as:

w nj +1 = w nj −

 λt  n w j − w nj−1 . x

(5)

The Kolgan scheme (which will be referred to as 2OS) can be written in conservative form as:

w nj +1 = w nj −

   t   n φ w L , j +1/2 , w nR , j +1/2 − φ w nL, j −1/2 , w nR , j −1/2 x

(6)

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

491

Fig. 1. Data reconstruction of the states obtained by linear interpolation.

where w L , j −1/2 , w R , j −1/2 , w L , j +1/2 , w R , j +1/2 are the data reconstruction (see Toro [17]) of the conservative variables at the two boundaries of a finite volume, x j −1/2 and x j +1/2 , which are obtained using linear extrapolation as:

w L , j −1/2 = w j −1 + w L , j +1/2 = w j +

1 2

1 2

∗  Lj − 1,

 Lj ∗ ,

w R , j −1/2 = w j + w R , j +1/2 = w j +1 +

1 2

1 2

 Rj ∗ ,

(7)

 Rj +∗ 1

(8)

where  Lj ∗ , and  Rj ∗ are the left and right limited slopes at the node x j

  Lj ∗

  Rj ∗

=

w j − w j −1 , w j +1 2 w j − w j −1 min[0, max( , w j +1 2

max[0, min(

=

− w j )] if w j +1 − w j > 0,

w j +1 − w j , w j −1 2 w j +1 − w j min[0, max(− , w j −1 2

max[0, min(−

(9)

− w j )] if w j +1 − w j < 0, − w j )] if w j −1 − w j > 0,

(10)

− w j )] if w j −1 − w j < 0.

(See Fig. 1.) In order to analyse the order of convergence of the scheme, we will consider the case in which no slope or flux limiters are needed. In this case, the slopes used to reconstruct the variables are given by:

 Lj ∗ =

w j − w j −1 2

,

 Rj ∗ = −

w j +1 − w j 2

(11)

.

Using Eq. (11) into Eq. (6), the second-order upwind scheme for the linear scalar transport equation, without limiters, reduces to:

w nj +1 = w nj −

λt  2 x



w nj−2 − 4w nj−1 + 3w nj ,

λ > 0,

(12)

this expression is the same that Kolgan referred in Eq. (23) of [9] as second branch. The new Kolgan-type upwind scheme (which will be referred to as 2NS) was introduced in [4] as a scheme of secondorder accuracy in space and first-order accuracy in time, and it is defined by means of the numerical flux φ(u , v , u L , v R ). It computes the centred discretisation of the convective flux from the values of the variables at the nodes, u and v, since this is a second-order accuracy interpolation. The scheme uses the data reconstruction of the variables at the cell faces, u L and v R , only for the numerical viscosity part of the numerical flux. In this way, the numerical flux is computed as:

φ(u , v , u L , v R ) =

f (u ) + f ( v ) 2

 1 −  Q (u L , v R )( v R − u L ).

(13)

2

After some simple mathematical manipulation the 2NS scheme for the linear transport equation, in the case in which no slope or flux limiters are needed, can be written as:

w nj +1

=

w nj

  λt 1 n 1 n n n n − − 3w j −1 + 3w j − w j +1 + w j +2 , w 2 x 2 j −2 2

λ > 0.

(14)

492

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Table 1 Sufficient conditions for linear stability (SCLS),

α = ωx.

Scheme

 g (ω, x, t )2  1

2OS

μ((1 − cos(α )) + 4 sin(α ) )  2

2NS

μ[(1 − cos(α )) + 1 − cos (α )] − 2(1 − cos(α ))  0

2

4

SCLS

μ

2

2

2

μ

1 2 1 2

2.2. Stability analysis In order to compare the previous schemes we will start by a linear stability analysis for λ > 0. The amplification factor for a wave number ω of the Kolgan scheme (2OS) given by Eq. (12), g (ω, x, t )2OS , and for the new Kolgan-type scheme (2NS) given by Eq. (14), g (ω, x, t )2NS , are given by:



2



2

g (ω, x, t )2OS = 1 − μ 1 − cos(ωx)

g (ω, x, t )2NS = 1 − μ 1 − cos(ωx)

  − i2μ sin(ωx) 1 − cos(ωx) , − i μ sin(ωx)

(15)

where μ is the Courant number. The stability condition over the norm of the amplification factor  g (ω, x, t )2  1 for any wave number ω is shown in Table 1, together with the sufficient conditions of linear stability (SCLS) over the Courant number for the two second-order methods. The stability region for both second-order schemes is reduced in comparison with the first-order upwind scheme for which the stability condition is μ  1. 2.3. Local truncation error In order to analyse their accuracy, the local truncation error τ n of the two schemes has been obtained. The local truncation error τ n is computed by replacing in Eqs. (12) and (14) the approximated solution w nj by the exact solution w (x j , t n ), and doing its Taylor-series expansion. For the 2OS and 2NS we obtain respectively:

 

τ nj τ nj

 2OS

    1      1 = t w tt x j , t n + O t 2 − x2 λ w xxx x j , t n + O x3 ,

(16)

2NS

    = t w tt x j , t n + O t 2 +

(17)



2 1

3 1

2

6

     x2 λ w xxx x j , t n + O x3 .

As it can be observed, (τ nj )2OS and (τ nj )2NS are dominated by an O (x2 ) term which confirms that the methods are second-order accuracy in space, for the Kolgan scheme it is also mentioned in [9]. The fact that the dominant error term depends on w xxx (x j , t n ) reflects the dispersive nature of both methods. The ratio between the coefficients of the dominant spatial error terms is −

1 λ w xxx (x j ,t n ) 3 1 λ w xxx (x j ,t n ) 6

= −2. This will be confirmed numerically in the following section, in Tables 3–4 which

report the global error of both schemes. 2.4. Numerical convergence error for a test case The local truncation error and the formal order of accuracy analysed in the previous section were determined by substitution of the Taylor-series expansion into the discretized equations. Alternatively, the observed order of accuracy of a numerical scheme can be determined by examination of the numerical solution obtained for a given test case. The observed order of accuracy will not necessarily be equal to the formal order of accuracy, especially when the solutions are outside of the asymptotic grid convergence range, as mentioned in [16]. In order to compute the observed order of accuracy in a test case for the one-dimensional scalar transport equation, we will consider as initial condition w (x, 0) = sin(2π x) on the range [a, b] = [0, 1], a value of λ = 1, and periodic boundary conditions. The exact solution for this test case is w (x, t ) = sin(2π (x − λt )). In the implementation of the numerical methods it was not used any flux limiter, in order to follow the algebra used to obtain the local truncation error in the previous section. The observed order of accuracy of each method is estimated by looking at the behaviour of various norms of the spatial error, εx . We will consider the L 1 and L ∞ norms, and the spatial errors at a fixed time t f = N t in a uniform mesh with −a , where M is the number of nodes. The spatial error εx is defined for each norm as: spacing x = bM

 

1 f = x ε x t

M    w x j , t f − w N , j

j =1



 





∞ f = max w x j , t f − w Nj . ε x t j

(18) (19)

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

493

Table 2 Number of nodes and spacing of the different grids used in the 1D test case without source terms.

Number of nodes M i , i = 1, 4 Spacing xi , i = 1, 4

G1

G2

G3

G4

40 0.025000

80 0.012500

160 0.006250

320 0.003125

Table 3 Spatial error of the Kolgan scheme (2OS) with the maximum norm used in the 1D test case without source terms. ∞ (1), 2OS ε x

t

G1

G2

G3

G4

0.00125 0.000625 0.0003125 0.00015625 0.000078125 3.90625E−05 1.95313E−05 9.76563E−06 4.88281E−06

0.027401 0.025831 0.025617 0.025577 0.025644 0.025695 0.025720 0.025733 0.025739

0.013677 0.008691 0.007001 0.006552 0.006457 0.006445 0.006444 0.006447 0.006451

0.012475 0.006350 0.003444 0.002201 0.001769 0.001649 0.001620 0.001614 0.001614

0.012414 0.006195 0.003109 0.001589 0.000865 0.000554 0.000444 0.000413 0.000405

Table 4 Spatial error of the new Kolgan-type scheme (2OS) with the maximum norm. ∞ (1), 2NS ε x

t

G1

G2

G3

G4

0.00125 0.000625 0.0003125 0.00015625 0.000078125 3.90625E−05 1.95313E−05 9.76563E−06 4.88281E−06

0.015953 0.013258 0.012901 0.012970 0.013081 0.013136 0.013164 0.013178 0.013185

0.012442 0.006638 0.004215 0.003431 0.003250 0.003228 0.003233 0.003240 0.003244

0.012389 0.006190 0.003146 0.001699 0.001084 0.000875 0.000820 0.000808 0.000807

0.012408 0.006184 0.003089 0.001550 0.000791 0.000429 0.000275 0.000221 0.000206

A common approach to analyse grid convergence is to compute the ratios of spatial errors on grids of different size. In this case we will also consider different time steps in order to report the order of accuracy on time. The details of the four grids G i , i = 1, 4, considered are shown in Table 2. For all the numerical computations the final time is t f = 1 s and all t satisfy, for the four grids considered, the sufficient linear stability condition over the Courant number for the two second-order schemes (see Table 1). x The selected grids verify x i = 2 and therefore, if the scheme is of order of accuracy p in space the following relation i +1

applies:

Kk =

εxk (t f ) (xk ) p  = 2p εxk+1 (t f ) (xk+1 ) p

⇒

log(Kk ) = p log(2).

(20)

Therefore, in terms of Kk , for two grids (G k , G k+1 ), the order of accuracy in space p can be estimated as

pk =

log(Kk ) log(2)

,

k = 1, 3.

(21)

Tables 5–6 have been obtained computing the observed order of accuracy pk , k = 1, 3, for the same time steps tl , in rows l = 1, 9, considered in Tables 3–4 (t 1 = 0.00125, tl = tl−1 /2, 1 < l  9). Both methods, 2OS and 2NS, are secondorder accurate in space, as it was estimated in the previous section. Notice that, in order to observe the second-order behaviour in Tables 5–6, it is necessary to reduce the time step, so that the error due to the time discretisation is much smaller than the error due to the spacial discretisation. The order of accuracy on time can also be reported by ratios of errors for different time steps with an equivalent procedure as the one used to obtain the order of accuracy on space. The time steps tl , l = 1, 9, are also defined with the t ratios t l = 2. Therefore, the scheme is of order of accuracy q in time if the following relation applies: l +1

Ll =

εxk (t l N l ) εxk

(t l+1 N l+1 )



(t l )q = 2q , (t l+1 )q

l = 1, 8

where the pairs (t l+1 , N l ) and (t l+1 , N l+1 ) verify:

t f = t l N l = t l+1 N l+1 .

⇒

log(Lk ) = q log(2),

(22)

494

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Table 5 Observed spatial order of accuracy p with L ∞ norm for schemes 2OS and 2NS. Kolgan scheme 2OS

New Kolgan-type scheme 2NS

log(K1 ) log(2)

log(K2 ) log(2)

log(K3 ) log(2)

log(K1 ) log(2)

log(K2 ) log(2)

log(K3 ) log(2)

1.00 1.57 1.87 1.96 1.99 2.00 2.00 2.00 2.00

0.13 0.45 1.02 1.57 1.87 1.97 1.99 2.00 2.00

0.01 0.04 0.15 0.47 0.90 1.57 1.87 1.96 1.99

0.36 1.00 1.61 1.92 2.01 2.02 2.03 2.02 2.02

0.01 0.10 0.42 1.01 1.58 1.88 1.98 2.00 2.01

0.00 0.00 0.03 0.13 0.45 1.03 1.58 1.87 1.97

Table 6 Observed spatial order of accuracy p with L 1 norm for schemes 2OS and 2NS. Kolgan scheme 2OS

New Kolgan-type scheme 2NS

log(K1 ) log(2)

log(K2 ) log(2)

0.394 1.042 1.645 1.947 2.034 2.054 2.056 2.055 2.054

0.014 0.114 0.440 1.034 1.602 1.897 1.993 2.016 2.021

log(K3 ) log(2)

−0.001 0.003 0.030 0.138 0.463 1.035 1.585 1.878 1.976

log(K1 ) log(2)

log(K2 ) log(2)

log(K3 ) log(2)

1.041 1.607 1.900 1.997 2.020 2.025 2.024 2.023 2.023

0.145 0.470 1.042 1.591 1.883 1.980 2.007 2.012 2.014

0.009 0.039 0.154 0.478 0.916 1.582 1.873 1.972 1.998

Table 7 Observed time order of accuracy q with L 1 norm for schemes 2OS and 2NS. Kolgan scheme 2OS log(L1 ) log(2) log(L2 ) log(2) log(L3 ) log(2)

New Kolgan-type scheme 2NS

G1

G2

G3

G4

0.08

0.65

0.97

1.00

0.02

0.31

0.88

0.99

0.00

0.09

0.64

0.97

G1 log(L1 ) log(2) log(L2 ) log(2) log(L3 ) log(2)

0.25

G2

G3

G4

0.90

1.00

1.00

0.05

0.65

0.97

1.00

−0.01

0.29

0.89

0.99

Consequently, the order of accuracy on time q is estimated as:

ql =

log(Ll ) log(2)

,

l = 1, 8.

(23)

In Table 7 the values of ql are presented in rows l = 1, 8 for the four grids, for both norms and for the two numerical schemes studied. In the first row of Table 7, which corresponds to the coarsest mesh on time, the numerical solution exhibits a first-order accuracy on time. This can be observed on the grids G 3 and G 4 , where the error in space is less log(L ) relevant than the error in time, in order to simplify the information we only consider the rows corresponding to log(2l) , l = 1, 3. If the time step is reduced, the error in time is reduced and hidden by the spatial error. Similar results are obtained with the norm L 1 . 3. Analysis of the new Kolgan-type scheme for a linear scalar equation with source terms In this section we will consider the following linear scalar equation with a source term g (x, w ):

∂w ∂w (x, t ) + λ (x, t ) = g (x, w ). ∂t ∂x

(24)

3.1. An analytical solution with a geometry dependent scalar source term Some kinds of numerical models consider a source term function g that depends on a geometric function (as for example the bed gradient or the breath gradient in the 1D shallow water equations), as well as on the conservative variable (as for example the water depth or the water discharge). A source term which commonly appears in shallow water models can be written as:

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

g (x, w ) = − w (x, t )σx (x),

495

(25)

where g (x, w ) leads to exponential decay (if σx (x) > 0) or growth (if σx (x) < 0) along the characteristic lines, as it is mentioned in [11]. Moreover, for any initial conditions w 0 (x, t ) the solution of the Cauchy problem is given by: 1

w (x, t ) = w 0 (x − λt )e λ (σ (x−λt )−σ (x)) .

(26)

In order to analyse the flux and source term discretisations it should be noticed that this model [11] has an equilibrium solution w e , which can be obtained for the boundary problem on the range [a, b] with inflow boundary condition w (a, 0) = η as: 1

w e (x, t ) = ηe − λ (σ (x)−σ (a)) .

(27)

The linear source term g (x, w ) = − w (x, t )β can be obtained in terms of the geometry function as:

σ (x) = β x + γ ,

(28)

and the transient and steady solutions of this model are:

w (x, t ) = w 0 (x − λt )e −β t , e

w (x, t ) = ηe

− βλ (x−a)

(29) (30)

.

This test case was also studied in [2] and [13]. We will show in this section that it is a very relevant way in which the geometric part of the source term is discretized in this linear case. 3.2. The numerical sources for the Kolgan-type scheme The numerical source functions ψL and ψR have been introduced in [2] for non-linear hyperbolic systems as an extension of the Q -schemes in order to define an upwind discretisation of the source terms as:





    ψL x j −1 , x j , w nj−1 , w nj + ψR x j , x j +1 , w nj , w nj+1 ,  

 ψL (x, y , u , v ) = X ( w Q ) I + Λ( w Q )Λ−1 ( w Q ) X −1 ( w Q ) gˆ (x, y , u , v ),  

 ψR (x, y , u , v ) = X ( w Q ) I − Λ( w Q )Λ−1 ( w Q ) X −1 ( w Q ) gˆ (x, y , u , v ), g x j , t n ≈ g nj =

1

2

(31) (32) (33)

where X and Λ are respectively the eigenvector and the diagonal eigenvalues matrix of the Jacobian A, w Q is the mean state of the corresponding Q -scheme, and gˆ (x, y , u , v ) is a centred approximation of g. We recall the system formulation because this is the discretisation considered in Section 5 for the geometric source terms. For the linear scalar problem, the numerical source functions are simplified to:

 1

1 + sign(λ) gˆ (x, y , u , v ), 2  1

ψR (x, y , u , v ) = 1 − sign(λ) gˆ (x, y , u , v ), 2 ψL (x, y , u , v ) =

(34) (35)

where sign(x) is the signum function. The extension of the Kolgan scheme 2OS to the upwind source term discretisation is obtained with the data reconstruction values w nL, j ±1/2 and w nR , j ±1/2 defined in Section 2.1 as:

g nj =

 1

     1

1 + sign(λ) gˆ x j −1 , x j , w nL, j −1/2 , w nR , j −1/2 + 1 − sign(λ) gˆ x j , x j +1 , w nL, j +1/2 , w nR , j +1/2 . 2 2

The way in which the centred discretisation gˆ The numerical source functions defined by position is related in [4] with the definition of the values of the variables at the nodes, while

g nj =

1 2

(36)

is done will be detailed in Section 3.3. Eqs. (32)–(33) have a centred part and an upwind component. This decomthe new Kolgan-type scheme 2NS: the centred component is computed from the upwind component is evaluated from the data reconstruction values:

  1   ψL x j −1 , x j , w nj−1 , w nj , w nL, j −1/2 , w nR , j −1/2 + ψR x j , x j +1 , w nj , w nj+1 , w nL, j +1/2 , w nR , j +1/2 . 2

(37)

Therefore, the numerical sources of the 2NS scheme for a hyperbolic non-linear system can be written as:

ψL (x, y , u , v , u L , v R ) = gˆ (x, y , u , v ) + X ( w QLR )Λsign ( w QLR ) X −1 ( w QLR ) gˆ (x, y , u L , v R ), ψR (x, y , u , v , u L , v R ) = gˆ (x, y , u , v ) − X ( w QLR )Λsign ( w QLR ) X

−1

( w QLR ) gˆ (x, y , u L , v R )

(38) (39)

where Λsign is the diagonal matrix with the signum of the eigenvalues of A, that is, (Λsign )ii = sign(λi ). After some simple mathematical manipulation, the numerical sources for the linear transport equation can be written as:

496

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

ψL (x, y , u , v , u L , v R ) = gˆ (x, y , u , v ) + sign(λ) gˆ (x, y , u L , v R ),

(40)

ψR (x, y , u , v , u L , v R ) = gˆ (x, y , u , v ) − sign(λ) gˆ (x, y , u L , v R ).

(41)

Therefore, the source discretisation used in the 2NS scheme is given by:

g nj =

   1  gˆ x j −1 , x j , w nj−1 , w nj + sign(λ) gˆ x j −1 , x j , w nL, j −1/2 , w nR , j −1/2 2    1  + gˆ x j , x j +1 , w nj , w nj+1 − sign(λ) gˆ x j , x j +1 , w nL, j +1/2 , w nR , j +1/2 . 2

(42)

3.3. Centred source discretisation The centred discretisation gˆ is very relevant in order to prove the second-order accuracy of the numerical schemes. We will start with a remark about the function σ used in Eq. (28). In this case the source term is g (x, w ) ≡ g ( w ) = −β w. Therefore, we can avoid the existence of σ , that is, the spatial dependence of g (from now on CS1 option), or we can solve our problem as a particular case and we define gˆ having into account this dependence (from now on CS2 option). We have also the possibility to introduce a data reconstruction of this geometry function σ . 3.3.1. CS1 option for 1OS, 2OS and 2NS schemes If we consider that g ( w ) = −β w, the source term is linear on w, so a natural election for gˆ is

gˆ (x, y , u , v ) =

1 2

   u+v 1 g (u ) + g ( v ) = g = −β (u + v ). 2

(43)

2

In this case we use the same arithmetic mean for gˆ (x, y , u L , v R ). After some simple mathematical manipulation of Eqs. (6), (36) and (43), the Kolgan scheme for the linear transport equation without limiters can be written as:

 βt   λt  n − w nj−2 + 3w nj−1 + 3w nj − w nj+1 , w j −2 − 4w nj−1 + 3w nj − 2 x 4 In this case the local truncation error (τ j )2OS reduces to w nj +1 = w nj −

λ > 0.

(44)

    1 1 (τ j )2OS = t w tt x j , t n + xβ w x x j , t n 2

2

1



+ (x) β w xx x j , t

2

n



4

  1 − (x)2 λ w xxx x j , t n + O(t )2 + O(x)3 . 3

(45)

The dominant term xβ 12 w x (x j , t n ) confirms that this method is first-order accurate in space. Moreover, the factor β in this term shows that the lack of second-order accuracy is due to the presence of the source term. In the case of the 2NS scheme, using the same centred discretisation of the source gˆ , after some mathematical manipulation of Eqs. (14), (37) and (43), the 2NS scheme can be written as:



λt 1 n 1 w j −2 − 3w nj−1 + 3w nj − w nj+1 + w nj+2 2 x 2 2

βt 1 1 − − w nj−2 + 3w nj−1 + 2w nj − w nj+1 + w nj+2 .

w nj +1 = w nj −

4

2

2

(46)

The local truncation error in this scheme (τ j )2NS is:

    1 1 (τ j )2NS = t w tt x j , t n + xβ w x x j , t n 2

1



− (x)2 β w xx x j , t 4

2

 n

  1 + (x)2 λ w xxx x j , t n + O(t )2 + O(x)3 . 6

(47)

Therefore the order of accuracy is the same that for the 2OS, and the lack of second-order accuracy in space is also related with β . The numerical convergence error is computed with the same procedure reported in Section 2.4. In order to analyse the dependence on the source term, we will try to reproduce the same functions and parameters of the test case solved in Section 2.4. Therefore, the grids and initial conditions are the same that have been considered in Table 2. The exact solution now is given by Eq. (29). We will take β = 1. Table 8 shows the observed order of accuracy pk (k = 1, 3), for the same time steps tl (in rows l = 1, 9) considered in Tables 3–4. Both schemes 2OS and 2NS show a first-order accuracy behaviour. Nevertheless the values of pk are larger with 2NS scheme than with 2OS scheme. This behaviour appears with the two norms L 1 and L ∞ .

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

497

Table 8 Observed spatial order of accuracy p with L 1 norm for schemes 2OS and 2NS, β = 1. 2OS, CS1 option β = 1

2NS, CS1 option β = 1

log(K1 ) log(2)

−0.08 0.10 0.08 0.14 0.21 0.21 0.21 0.72 0.72

log(K2 ) log(2)

log(K3 ) log(2)

log(K1 ) log(2)

log(K2 ) log(2)

log(K3 ) log(2)

0.45 0.67 0.66 0.70 0.71 0.71 0.71 0.87 0.87

0.27 0.67 0.68 0.86 0.91 0.88 0.88 0.94 0.94

1.24 1.25 1.22 1.21 1.23 1.23 1.23 1.12 1.11

0.99 1.13 1.12 1.13 1.13 1.13 1.12 1.06 1.06

0.49 0.91 0.92 1.07 1.09 1.07 1.07 1.04 1.03

3.3.2. CS2 option for 1OS, 2OS and 2NS schemes In this case we will define gˆ having into account the dependency on σ . This option is related to the ideas considered in [13] and the model formulation used in [11]. The centred approximation of gˆ is defined as:

gˆ (x, y , u , v )1 = −

σ ( y ) − σ (x) u + v

(48)

, x 2 σ ( y ) R − σ (x)L u + v gˆ (x, y , u , v )2 = − , x 2   x+ y u+v gˆ (x, y , u , v )3 = −σx . 2

(49) (50)

2

• Discretisation gˆ (x, y , u , v )1 is the classical one used in first-order schemes. It is equivalent to CS1 option if σ is a linear function. We consider this for 1OS, 2OS, and for the centred source discretisation of 2NS.

• Discretisation gˆ (x, y , u , v )2 uses the data reconstruction proposed in Section 2.1 for the geometry function σ . From this point of view it seems an optimal election for second-order accuracy in space. We have tried this option for the 2OS scheme and for the upwind part of the 2NS scheme. We will show that this choice is worse for 2OS than gˆ (x, y , u , v )1 . • Finally, discretisation gˆ (x, y , u , v )3 is equivalent to CS1 option, if σ is given by Eq. (28). The analysis of this alternative discretisation for gˆ will be done in this section for the schemes 2OS and 2NS. From Eqs. (6), (36) and (49), the 2OS scheme for the linear transport equation (λ > 0) without limiters, can be written as:

w nj +1 = w nj −



n n  σ R , j−1/2 − σL , j−1/2 w L , j−1/2 + w R , j−1/2 λt  n w j −2 − 4w nj−1 + 3w nj − t 2 x x 2

(51)

where σ L , j −1/2 and σ R , j −1/2 are the linear interpolation values for the function σ obtained as in Eq. (8). It should be noticed that in the cases where σ is a linear function the source discretisation is null, which can be easily shown by replacing in it the approximated solution w nj by the exact solution w (x j , t n ), and doing its Taylor-series expansion n n σ R , j−1/2 − σL , j−1/2 w L , j−1/2 + w R , j−1/2

x

2

  1 ≈ − (x)2 σxxx (x j ) w x j , t n . 2

(52)

This is equivalent to neglect the source term discretisation, and the local truncation error (τ j )2OS reduces to:

    1     (τ j )2OS = w t x j , t n + λ w x x j , t n + t w tt x j , t n + O t 2 2

 1   1   − x2 λ w xxx x j , t n + (x)2 σxxx (x j ) w x j , t n + O(x)3 . 

3

2

(53)

1 Thus, we conclude that this scheme has zero-order accuracy. This is shown in Table 9, which shows ε x for different spatial discretisations (see Table 2) and for three time steps. In the case of the 2NS scheme we will use the discretisations gˆ (x, y , u , v )1 and gˆ (x, y , u , v )2 for the centred and upwind parts respectively. After some mathematical manipulation of Eqs. (14) and (37), using de discretisations gˆ (x, y , u , v )1 and gˆ (x, y , u , v )2 for the centred and upwind parts of the source discretisation, the 2NS scheme can be written as:

w nj +1 = w nj −



λt 1 n 1 w j −2 − 3w nj−1 + 3w nj − w nj+1 + w nj+2 2 x 2 2

498

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Table 9 Spatial error of Kolgan scheme 2OS with L 1 norm. 1 ε x (1), CS2 option σx (x) = β = 1

t

G1

G2

G3

G4

1.25E−03 6.25E−04 3.13E−04

0.416849589 0.408746527 0.404731751

0.418115988 0.410089856 0.40611303

0.418280476 0.410274543 0.406307922

0.418308825 0.410311374 0.406347245

Table 10 Observed order of accuracy p with L 1 norm for scheme 2NS, with β = 1. 2NS, CS2 option

σx (x) = β = 1

log(K1 ) log(2)

log(K2 ) log(2)

log(K3 ) log(2)

−0.02 0.06 0.37 0.99 1.59 1.90 1.98

0.27 0.94 1.59 1.91 1.99 2.00 2.00

−0.01 −0.01 0.00 0.09 0.41 1.02 1.60

Table 11 1 Equilibrium solution. Spatial error ε x (1) of 1OS, 2OS, 2NS schemes with L 1 norm. 1 ε x (1), σ (x) = sin(x), t = 10

1OS 2OS 2NS

G1

G2

G3

G4

0.026507 3.689322 0.147283

0.006731 2.071982 0.017223

0.001689 1.105094 0.001939

0.000373 0.539237 0.000325



n n n n σ R , j−1/2 − σL , j−1/2 w L , j−1/2 + w R , j−1/2 t σ j − σ j −1 w j −1 + w j + 2 x 2 x 2

n n n n t σ j +1 − σ j w j + w j +1 σ R , j +1/2 − σ L , j +1/2 w L , j +1/2 + w R , j +1/2 − . − 2 x 2 x 2

−

(54)

The local truncation error for this scheme (τ j )2NS is:

      1 1 (τ j )2NS = t w tt x j , t n + O t 2 + (x)2 λ w xxx x j , t n 2 6

      (x)2 7 σxxx (x j ) w x j , t n + σxx (x j ) w x x j , t n + σx (x j ) w xx x j , t n + 2

+

12

(x)3

4













σxxxx (x j ) w x j , t n + σxxx (x j ) w x x j , t n + O x4 .

(55)

Therefore, the 2NS is second-order accurate in space. This is consistent with the observed order of accuracy pk , k = 1, 3, for the scheme 2NS, which is shown in Table 10. Notice in Eq. (55) that the discretisation of the upwind part of the numerical sources is third-order accurate. If σ is a linear function this term is null, and the scheme is equivalent (to the second-order accuracy) to a centred discretisation of the source term given by:

w nj +1

=

w nj



 n  n n βt w j −1 + 2w j + w j +1 λt 1 n 1 n n n n . − − 3w j −1 + 3w j − w j +1 + w j +2 − w 2 x 2 j −2 2 2 2

(56)

3.4. Solving an equilibrium solution, analysis of the C -property A stationary solution is considered here in order to compare the different numerical methods. In order to illustrate the behaviour of both methods with a non-linear σ function we will consider the function σ (x) = sin(x). The equilibrium solution is given by Eq. (27). This problem is a good test to validate the C -property [2] and to prove if the numerical schemes are well-balanced. (See Figs. 2 and 3.) Tables 11 and 12 show that the schemes 1OS and 2NS verify the approximated C -property, because the equilibrium solution is obtained with at least second-order accuracy.

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 2. Numerical solution at t f = 1 for

Fig. 3. Absolute error at t f = 1 for

499

σ (x) = sin(x), x = x4 = 0.029411764, t = 0.002941176 and λ = 1.

σ (x) = sin(x), x = x4 = 0.029411764, t = 0.002941176, β = 1 and λ = 1.

Table 12 Observed order of accuracy p, t = 10 with L 1 norm. pk , k = 1, 3,

σ (x) = sin(x)

1OS 2OS 2NS

1.98 0.83 3.10

1.99 0.91 3.15

2.18 1.04 2.57

3.5. Solving a discontinuous solution The numerical results presented in previous sections are related to the analysis of the order of accuracy of the numerical schemes with continuous solutions. However, hyperbolic problems are often related to discontinuous solutions. A wellknown problem with source terms, which was introduced in [12], is defined by Eq. (24) with a non-linear source term defined by:

500

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 4. Numerical solution at t f = 0.3 for x = 0.00625, t = 0.0028125,



g (x, w ) = −η w (x, t ) w (x, t ) − 1



 w (x, t ) −

1 2

η = 1000, λ = 1 and μ = 0.45.

 .

(57)

If the previous source term is linearised, the parameter η is equivalent to the parameter β in Eq. (28). If the source term is zero for the initial conditions, the solution is the same as the solution of the transport equation w (x, t ) = w (x − λt , 0). In [12] the following initial conditions are considered:

w (x, 0) =



1 if x  0.3, 0 if x > 0.3.

(58)

The solution will be obtained with the three schemes studied: 1OS, 2OS, 2NS. The centred discretisation of the source term considered for all of them is given by:

gˆ (x, y , u , v ) =

g (u ) + g ( v ) 2

(59)

.

To solve this problem flux limiters are considered. In Fig. 4 it is clearly shown the dispersive effect of the solution, which was predicted in the analysis of the local truncation error (Section 2.3) for the 2OS and 2NS schemes. The different signum and factor related with the third space derivative for the two schemes are also shown in Fig. 4. Fig. 5 shows for the three schemes the average speed or numerical velocity defined in [12] as: ∗

λ =

x  tn

w ni

−

i



w 0i

.

(60)

i

4. Analysis of the new Kolgan-type scheme for a non-linear scalar equation with source terms The main aim of this section is to solve the difficulties related to the discretisation of a non-linear flux combined with source terms. We will consider the Burgers equation with the source term introduced in Section 3.1, which is given by:

∂w ∂w (x, t ) + w (x, t ) = − w (x, t )σx (x), ∂t ∂x

(61)

2

where f ( w ) = w2 is the non-linear flux, and g (x, w ) = − w (x, t )σx (x) is the same source term which appears in Eq. (25). This problem has a family of steady solutions which are given by:

w e (x) = K − σ (x),

K = const.

(62)

This steady solution allows us to verify the C -property of the numerical methods analysed in this paper. Since w Q is equal to the arithmetic mean, which is equal to the Roe-mean for the Burgers equation, the numerical flux for the 2OS and 2NS schemes in this problem are given by:

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 5. Numerical velocity at t f = 0.3 for x = x3 = 0.001, t = t 1 = 0.00125,

501

η = 1000, λ = 1 and μ = 0.15.

Table 13 1 Equilibrium solution. 1OS, 2OS, 2NS schemes, ε x (1) with L 1 norm. 1 ε x (2)

1OS 2OS 2NS

G1

G2

G3

G4

2.22044e−17 0.278209 5.24566e−14

9.71445e−18 0.1637840 4.03844e−08

2.72351e−17 0.0985740 7.20693e−09

3.24393e−17 0.053392 1.57018e−15

  1  uL + v R  ( v R − u L ), −   2 2 2   f (u ) + f ( v ) 1  uL + v R  ( v R − u L ). φ(u , v , u L , v R )2NS = −   φ(u L , v R )2OS =

f (u L ) + f ( v R )

2

2

2

(63) (64)

The numerical source term discretisation is given by:





g nj 2OS

=

1 2

+ 

1 + sign 1 2

 wn

L , j −1/2

1 − sign

+ w nR , j −1/2    gˆ x j −1 , x j , w nL, j −1/2 , w nR , j −1/2 2

 wn

L , j +1/2

+ w nR , j +1/2    gˆ x j , x j +1 , w nL, j +1/2 , w nR , j +1/2 , 2

 n

 1   1  g j 2NS = gˆ x j −1 , x j , w nj−1 , w nj + gˆ x j , x j +1 , w nj , w nj+1 2 2 1

+ sign 2 1

− sign 2

 wn

L , j −1/2

 wn

+ w nR , j −1/2    gˆ x j −1 , x j , w nL, j −1/2 , w nR , j −1/2 2

+ w nR , j −1/2    L , j +1/2 gˆ x j , x j +1 , w nL, j +1/2 , w nR , j +1/2 . 2

(65)

Finally, we will use the definition of gˆ given by Eq. (48) in the centred and upwind parts of the source term for schemes 1OS and 2OS, and also in the centred part of the source term for scheme 2NS. Eq. (49) will be used in the upwind discretisation of the source term for the scheme 2NS. The equilibrium solution is solved exactly by the schemes 1OS and 2NS, their accuracy being independent of the mesh size. Therefore, both schemes verify the exact C -property. This is not the case for the scheme 2OS. Figs. 6, 7 and Table 13 show the absolute error for grid G 1 . It is clear that the scheme 2OS does not verify the C -property and that it is not a well-balanced scheme.

502

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 6. Absolute error. The schemes 1OS and 2NS verify the exact C -property. x = 0.05.

Fig. 7. Approximated solution (σ (x j ) + w n (x j )), j = 1, M, with 1OS, 2ON and 2NS schemes. Test exact C -property. x = 0.05.

4.1. Local truncation error and analysis of the C -property for the new Kolgan-type scheme 2NS The local truncation error of the scheme 2NS is given by:

(τ j )2NS =

t 2

+ +









w tt x j , t n + O t 2 +

(x)2

4



6

 

f xxx w x j , t n







  + σxxx (x j ) w x j , t n

σxx (x j ) w x x j , t n + σx (x j ) w xx x j , t n +

(x)3

4



(x)2



w xxx x j , t

 n

(x)3

4

     + σxxx (x j ) w x x j , t n + O x4 .





 

w xxxx x j , t n + σxxxx (x j ) w x j , t n

 (66)

The third-order dominant terms are linked with the discretisation of the upwind part of both the numerical flux and the numerical source terms. Therefore, a centred discretisation of the source term verifies the approximated C -property. Nevertheless, it is necessary to use an upwind discretisation of the flux in order to obtain numerical stability. Moreover,

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

503

for the equilibrium equation, in which the term w (x) + σ (x) and all its derivatives of the same order vanish, this thirdorder dominant terms are zero. The second-order dominant terms are also zero, considering that for Burgers equation f xxx ( w ) = 3w x w xx + w w xxx and for the equilibrium solution σx (x j ) = − w x (x j ), σxx (x j ) = − w xx (x j ). Therefore, the following relations apply:

(x)2

6 (x)2

6



f xxx ( w j ) + σxxx (x j ) w (x j ) =

(x)2 2



w x w xx ,

σxx (x j ) w x (x j ) + σx (x j ) w xx (x j ) = −

(67)

2

(x) 2

w x w xx

(68)

and then, from the local truncation error of 2NS we prove at least fourth order. It is straightforward to prove the exact C -property if we detail the expression of the 2NS for the equilibrium solution. It is only necessary to consider the following relations:

w i + σi = const,

σ R , j±1/2 − σL , j±1/2 = w R , j±1/2 − w L , j±1/2 ,

j = 1, M .

(69)

The equilibrium solution is verified at each node if the same linear interpolation for w and σ is considered. This is the same conclusion that was achieved in [4] in order to verify the C -property for the 2D shallow water equations with porosity. 4.2. Embid’s problem The equation considered in this section was chosen in order to analyse the behaviour of the new scheme 2NS, in terms of the well-balanced property. The Embid’s problem [7] also considers other difficulties related with the hyperbolic equations with source terms: it has two steady solutions which satisfy the entropy condition. This problem was presented as a simple scalar approximation to the one-dimensional equations that model the flow of a gas through a duct of variable cross-section. The source term is given by

g (x, w ) = −(6x − 3) w (x, t ) = −σx (x) w (x, t ),

σ (x) = 3x2 − 3x,

(70)

where two different ways of writing the source term are considered. We will consider the same notation introduced in Section 3.3: without σx (CS1 option) and with σ (CS2 option). The numerical schemes are tested to approach the stable entropy satisfying steady solution, with a standing shock at xs = 0.18 given by



w (x) =

1 + 3x2 − 3x, −0.1 + 3x2 − 3x,

x < xs , x > xs .

(71)

The boundary conditions are w (0, t ) = 1 and w (1, t ) = −0.1. As in [13], we consider initial conditions with a jump at xs = 0.18. It is clear in Figs. 8–9 that the continuous part of the solution is approximated better with the option CS2 than with the option CS1, while the shock is computed in a similar way with both alternatives. The results are obtained with the same parameters proposed in [13]. The scheme 2NS has a similar behaviour as the well-balanced scheme (WB2) defined in [13] where, in order to curb the oscillations near the shock their authors propose a hybrid scheme. 5. 2D flow in a square cavity In this last section we will consider a motivation test case for future developments of the scheme in which we solve the 2D shallow water equations with viscous stresses in a square cavity. We will consider the case of a flat bottom and the case with a slanted bottom. 5.1. Flat bottom In this test case, the 2D shallow water equations including viscous stresses are used to compute the velocity field in a square cavity of 1 × 1 m (Fig. 10) with a horizontal bed elevation. The equations solved are given by:

∂ h ∂ qx ∂ q y + + = 0, ∂t ∂x ∂y  2        ∂ qxq y ∂ ∂ qx gh2 ∂ zb ∂Ux ∂Ux ∂ qx ∂ + = − gh + , + + + νh νh ∂t ∂x h 2 ∂y h ∂x ∂x ∂x ∂y ∂y         ∂q y ∂U y ∂U y ∂ qxq y ∂ ∂ q2y ∂ gh2 ∂ zb + + + νh νh + = − gh + ∂t ∂x h ∂y h 2 ∂y ∂x ∂x ∂y ∂y

(72)

where zb is the bed elevation, (U x , U y ) are the two horizontal components of the depth-averaged velocity, |U| is the modulus of the depth-averaged horizontal velocity, h is the water depth, (q x , q y ) are the two components of the unit discharge,

504

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 8. Embid’s problem, 1OS, 2ON and 2NS with CS1 option.

Fig. 9. Embid’s problem, 1OS, 2ON and 2NS, with CS2 option.

Table 14 Size of the meshes used in the square cavity test case.

Element size (m) Number elements

G1

G2

G3

G4

0.05 400

0.025 1600

0.0125 6400

0.00625 25 600

g is the gravity acceleration, and ν is laminar viscosity. Since the bottom is flat, the bed slope source term is zero in this case. The flow enters and exits the domain respectively through a couple of slots of 0.20 m located in the upper-left and upper-right sides of the square. The water depth in the square is 1 m and the inlet discharge is 0.2 m3 /s. The laminar viscosity is ν = 0.02 m2 /s, and a no-slip boundary condition is used at the wall boundaries. For this condition, a large recirculation eddy is generated in the lower-left region of the square, as it is shown in Fig. 10. In order to analyse the order of convergence in space, four numerical meshes with different cell sizes are used, which will be referred to as meshes G 1 to G 4 . The number of elements on each mesh is shown in Table 14. (See Figs. 11 and 12.)

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 10. Velocity field in a square cavity.

Fig. 11. Numerical mesh with 40 × 40 cells.

505

Fig. 12. V x velocity profile. x = 0.5 m.

Fig. 13. Mesh convergence. Square cavity with flat bottom.

The norm L 1 is computed for meshes G 1 , G 2 and G 3 by comparing them with mesh G 4 . The results of the mesh convergence analysis are shown in Fig. 13 for the Kolgan scheme (2OS), the new Kolgan-type scheme (2NS), the first-order scheme (1OS), and the hybrid scheme proposed in [5], which uses a second-order extrapolation for the unit discharges and a first-order extrapolation for the water depth (12HS). The observed order of accuracy is approximately 1.7 for the hybrid scheme and for both second-order schemes, and 0.9 for the first-order scheme. 5.2. Slanted bottom This case is the same as the previous one, but the bed has a slope of 45◦ in the y-direction. Therefore, in this case the bed slope source term is not zero. The results of the mesh convergence analysis are very similar to the ones obtained in the previous case, and are shown in Fig. 14. 6. Conclusions In this paper the properties of a high-order numerical scheme (2NS) for hyperbolic conservation laws have been studied and compared with classic first-order (1OS) and high-order (2OS) discretisation schemes. The analysis has been done in terms of the local truncation error and the stability conditions, and the order of accuracy of the schemes has been computed for different test cases. The idea of the scheme 2NS is to use high-order extrapolation of the conservative variables just in the upwind contribution of the numerical flux and source terms. It has been proved the second-order accuracy of the scheme 2NS for scalar problems without source terms. It has also been pointed out the relevance of the source term discretisation

506

M.E. Vázquez-Cendón, L. Cea / Applied Numerical Mathematics 62 (2012) 489–506

Fig. 14. Mesh convergence. Square cavity with slanted bottom.

in order to preserve the second-order accuracy of the scheme and the C -property [2], which is related with the different equilibrium solutions considered in this paper. The efficiency of the proposed scheme is proved after solving different model problems related with transport and Burgers equations with and without source terms. Acknowledgements The authors are indebted to R. Donat, E.F. Toro and T. Chacón for many valuable discussions. They also want to thank the referees for their suggestions and useful remarks. References [1] E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein, B. Perthame, A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows, SIAM J. Sci. Comput. 25 (2004) 2050–2065. [2] A. Bermúdez, M.E. Vázquez-Cendón, Upwind methods for hyperbolic conservation laws with source terms, Comput. & Fluids 23 (8) (1994) 1049–1071. [3] M.J. Castro, J.M. Gallardo, C. Parés, High-order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems, Math. Comp. 75 (2006) 1103–1134. [4] L. Cea, M.E. Vázquez-Cendón, Unstructured finite volume discretisation of two-dimensional depth averaged shallow water equations with porosity, Internat. J. Numer. Methods Fluids 63 (2010) 903–930. [5] L. Cea, J. Puertas, M.E. Vázquez-Cendón, Depth averaged modelling of turbulent shallow water flow with wet–dry fronts, Arch. Comput. Methods Eng. 14 (3) (2007) 303–341. [6] M. Dumbser, M. Castro, C. Parés, E.F. Toro, ADER schemes on unstructured meshes for non-conservative hyperbolic systems: applications to geophysical flows, Comput. & Fluids 38 (2009) 1731–1748. [7] P. Embid, J. Goodman, A. Majda, Multiple steady states for 1-d transonic flow, SIAM J. Sci. Statist. Comput. 5 (1984) 21–41. [8] M.E. Hubbard, P. García-Navarro, Flux difference splitting and the balancing of source terms and flux gradients, J. Comput. Phys. 165 (2000) 89–125. [9] V.P. Kolgan, Application of the principle of minimizing the derivative to the construction of finite-difference schemes for computing discontinuous solutions of gas dynamics, J. Comput. Phys. 230 (2011) 2384–2390. [10] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wavepropagation algorithm, J. Comput. Phys. 146 (1998) 346–365. [11] R. LeVeque, A well-balanced path-integral f-wave method for hyperbolic problems with source terms, J. Sci. Comput., published online: 21 August, 2010, doi:10.1007/s10915-010-9411-0. [12] R. LeVeque, H.C. Yee, A study of numerical methods for hyperbolic conservation laws with stiff source terms, J. Comput. Phys. 86 (1990) 187–210. [13] A. Martínez-Gavara, R. Donat, A Hybrid second-order scheme for shallow water flows, J. Sci. Comput., published online: 1 December, 2010, doi: 10.1007/s10915-010-9440-8. [14] S. Noelle, N. Pankratz, G. Puppo, J.R. Natvig, Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows, J. Comput. Phys. 213 (2006) 474–499. [15] S. Noelle, Y.L. Xing, C.W. Shu, High-order well-balanced finite volume WENO schemes for shallow water equation with moving water, J. Comput. Phys. 226 (2007) 29–58. [16] C.J. Roy, Grid convergence error analysis for mixed-order numerical schemes, AIAA J. 41 (4) (2003) 595–603. [17] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer, Berlin, 2009. [18] B. van Leer, A historical oversight: Vladimir P. Kolgan and his high-resolution scheme, J. Comput. Phys. 230 (2011) 2378–2383. [19] Y.L. Xing, C.W. Shu, High-order well-balanced finite difference WENO schemes for a class of hyperbolic systems with source terms, J. Sci. Comput. 27 (1– 3) (2006) 477–494.