# Convergence Improved Lax-Friedrichs Scheme Based Numerical Schemes and Their Applications in Solving the One-Layer and Two-Layer Shallow-Water Equations

Authors: Xinhua  Lu; Bingjiang Dong; Bing Mao; Xiaofeng Zhang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The first-order Lax-Friedrichs (LF) scheme is commonly used in conjunction with other schemes to achieve monotone and stable properties with lower numerical diffusion. Nevertheless, the LF scheme and the schemes devised based on it, for example, the first-order centered (FORCE) and the slope-limited centered (SLIC) schemes, cannot achieve a time-step-independence solution due to the excessive numerical diffusion at a small time step. In this work, two time-step-convergence improved schemes, the C-FORCE and C-SLIC schemes, are proposed to resolve this problem. The performance of the proposed schemes is validated in solving the one-layer and two-layer shallow-water equations, verifying their capabilities in attaining time-step-independence solutions and showing robustness of them in resolving discontinuities with high-resolution.


---

## Body

## 1. Introduction

The Lax-Friedrichs (LF) scheme, also called the Lax method [1], is a classical explicit three-point scheme in solving partial differential equations in, for example, aerodynamics, hydrodynamics, and magnetohydrodynamics [2–4]. This scheme has many advantages, for example, being stable and monotone up to Courant-Friedrichs-Lewy (CFL) number approaching unity in solving the linear convection problem [4] and being computationally efficient and simple to implement (e.g., compared with a Godunov-type scheme). Due to the excellent monotone and stable properties, the LF scheme, which is of the first-order accuracy in space, has been commonly adopted for a combination use with some higher-order schemes to predict monotone and stable solutions with lower numerical diffusion [5]. For instance, Toro and Billett [6] derived the first-order centered (FORCE) scheme, by combining the LF scheme with the two-step, second-order Lax-Wendroff scheme [7]. When solving the linear convection problem, compared with the LF scheme, the numerical viscosity of the FORCE scheme is reduced by half [4]. The FORCE scheme and its second-order (both in time and in space) extension, the slope-limited centered scheme (SLIC) [8], have been used and verified to perform well in solving the one-layer shallow-water equations (1LSWEs), the Euler equations of compressible gas dynamics, and the magnetohydrodynamic equations [8, 9].Nevertheless, the LF scheme and the schemes devised based on it (e.g., the FORCE and SLIC schemes) have a common shortcoming, which has raised less attention so far. We demonstrate this shortcoming by using the LF scheme to solve the following linear convection problem in the time-space domain[0,T]×[0,L]:(1a)PDE:∂ϕ∂t+∂Fϕ∂x=0,Fϕ=λϕ,(1b)Initial  condition:ϕx,0=ϕ0x,(1c)Boundary  conditions:ϕ0,t=ϕ1t,ϕL,t=ϕ2t.Here, λ is the wave propagation speed and, without loss of generality, is assumed to be positive. An explicit finite-volume discretization for (1a) gives(2)ϕin+1=ϕin+ΔtΔxF^i-1/2-F^i+1/2,where the superscripts “n” and “n+1” denote the values at the old and new time steps, respectively; Δt and Δx denote the time step and space interval, respectively; F^ denotes the numerical flux and, for the LF scheme, it reads [1](3)F^i-1/2=12Fϕi-1n+Fϕin-Δx2Δtϕin-ϕi-1n.A truncation error analysis shows that the LF scheme in solving (1a), (1b), and (1c), with the discretizations given in (2) and (3), leads to a numerical viscosity of [4](4)αLF=1-Cr22CrλΔx,where Cr is the CFL number and, for the linear convection problem,(5)Cr=λΔtΔx.In Figure 1, αLF is plotted against Cr. For clear view purpose, the relation in the region Cr⩽0.01 is plotted on a semilog scale (see Figure 1(a)). For convenience, in Figure 1, λ=Δx=1 is assumed so that Cr=Δt (see (5)). As seen in Figure 1, αLF increases as Δt decreases. When Δt→0, αLF→∞, suggesting that the numerical diffusion is absolutely dominant over the physical convection; this is, of course, unrealistic. In numerical practice, just as usually done via implementing a grid-independence study to choose an appropriate spatial resolution (can be done for the LF scheme as when Δx→0, αLF→0), a Δt-independence study is required to choose a reasonable time step. However, due to the rapidly increased numerical diffusion when Δt decreases, the solution computed by using the LF scheme does not converge and thus a Δt-independence study cannot be implemented. Failing to achieve a Δt-independence solution introduces uncertainty in a numerical solution, and the time step can only be determined empirically. Note that, in some situations, for instance, when rapid wet-dry transitions exist (e.g., for simulations in the nearshore regions), or when some other constraints on the time step (e.g., those due to surface tension and fluid viscosity effects [10]) dominate the CFL condition, one may have to use a rather small time step compared with that determined by the CFL condition. In these situations, the LF-based schemes are inappropriate due to the excessive numerical diffusion at a small value of Δt.Figure 1
Numerical viscosity versusCr for the LF and LLF schemes in solving the linear convection problem: (a) Cr⩽0.01 and (b) 0.01⩽Cr⩽1. Note that the result for the LLF scheme is plotted on the right y-axis. (a) is plotted on a semilog scale.
(a)
(b)We observe that the problem, related to the LF scheme being unable to achieve aΔt-independence solution, can be resolved by the local LF (LLF hereafter) scheme [11]. The LLF scheme calculates the numerical flux as(6)F^i-1/2=12Fϕi-1n+Fϕin-λ~i-1/22ϕin-ϕi-1n,where λ~ denotes the local bound of the wave propagation speed and is determined by(7)λ~i-1/2=max⁡λi-1,λi,and, in solving (1a), (1b), and (1c),  λ~i-1/2≡λ. In the Appendix, we show that the LLF scheme can be viewed as a first-order accurate characteristics scheme when solving the linear convection problem.Inserting (6) and (7) into (2), one attains the following when solving the linear convection problem:(8)ϕin+1=Crϕi-1n+1-Crϕin,which introduces a numerical viscosity of(9)αLLF=1-Cr2λΔx.From (9) (also see Figure 1), it is seen that αLLF increases with decreasing Cr. However, αLLF has a maximum and limiting value of 0.5λΔx when Cr→0. Therefore, the uncontrolled excessive numerical viscosity in the LF scheme is resolved and a Δt-independence study can be performed. It is worth remarking that, besides the advantage in achieving a Δt-independence solution, the LLF scheme also has the merits (e.g., being stable and monotone) as the LF scheme. To date, the LLF scheme has been widely used to solve the aerodynamic and magnetohydrodynamic problems [2, 3, 11] but has not raised much attention in computational hydrodynamics (e.g., in solving the shallow-water equations).The motivation of this work is simple. As the LF-based schemes (e.g., the FORCE and SLIC schemes) share the same drawback as the LF scheme to be unable to achieve aΔt-independence solution (note that the numerical viscosity of the FORCE scheme is just half of that defined in (4)), while the LLF scheme overcomes this drawback, it is meaningful to improve the LF-based schemes by a substitution of the LF scheme with the LLF scheme. In this work, the FORCE and SLIC schemes are improved following this idea. For convenience, the modified FORCE (SLIC) scheme will be named as C-FORCE (C-SLIC) scheme hereafter, where the prefix “C” denotes convergence ensured. For illustration purpose and as representatives, the performance of the C-FORCE and C-SLIC schemes will be assessed in solving the 1LSWEs and the two-layer shallow-water equations (2LSWEs), in 1D space and without considering the interfacial shear stress (for the 2LSWEs) and bed friction.The rest of the paper is organized as follows. In Section2, a comparison of the performances between the C-FORCE (C-SLIC) scheme and the original FORCE (SLIC) scheme is made in solving the 1LSWEs, and a similar comparison in solving the 2LSWEs is presented in Section 3. Finally, conclusions and discussions are given in Section 4.


## 2. Performance Comparisons in Solving the 1LSWEs

We consider solving the 1LSWEs in this section. A sketch of the one-layer shallow-water system is illustrated in Figure2(a).Figure 2
Sketch of the shallow-water systems: (a) one-layer and (b) two-layer. Not to scale.
(a)
(b)


### 2.1. Governing Equations and Numerical Methods for the 1LSWEs

The 1LSWEs, in a conservative vector form, can be written as(10)∂U∂t+∂F∂x=S,with vectors defined by(11)U=zqx,F=qxuqx+g2h2,S=0-gh∂zb∂x,where t denotes time; x denotes the streamwise coordinate; z and zb denote the water surface and bed elevations above a horizontal reference level zr; h=z-zb denotes the water depth; and qx=uh denotes the discharge per unit width with u being the depth-averaged velocity in the x-direction.An explicit finite-volume discretization for (10) gives(12)Uin+1-UinΔt=F^i-1/2-F^i+1/2Δx+S¯i.Here, the overbar on a term denotes that this term is evaluated by the temporally evolved values (described later in this section); the source term vector is discretized by using a central difference scheme as [12](13)S¯i=0-gh¯i-1/2R+h¯i+1/2L2zbi+1/2-zbi-1/2Δx;the flux vector F^ is evaluated by the FORCE or SLIC scheme, which reads uniformly(14)F^i-1/2=ωFi-1/2LW+1-ωFi-1/2LF,where(15a)Fi-1/2LW=FUi-1/2LW,(15b)Ui-1/2LW=12U-i-1/2L+U-i-1/2R-12ΔtΔxFU-i-1/2R-FU-i-1/2L,(16)Fi-1/2LF=12FU-i-1/2L+FU-i-1/2R-12λ~i-1/2U-i-1/2R-U-i-1/2L.Here, the superscripts “LW” and “LF” denote that the variables are related to the two-step Lax-Wendroff and Lax-Friedrichs fluxes, respectively; the superscripts “L” and “R” denote the reconstructed values at the left and right sides of the cell-interface, respectively. For both the FORCE and SLIC schemes, ω=0.5 and λ~≡Δx/Δt; for both the C-FORCE and C-SLIC schemes, ω is set to be 0.3 empirically and λ~ is set to be the maximum local characteristic speed as(17)λ~i-1/2=max⁡u+ghi-1/2L,u+ghi-1/2R.For face value reconstructions, first, the monotonic upstream-centered scheme for conservation laws (MUSCL) [13] is applied; for example,(18a)zi-1/2L=zin+12Ψrizin-zi-1n,zi-1/2R=zi+1n-12Ψri+1zi+1n-zin,(18b)hi-1/2L=hin+12Ψrihin-hi-1n,hi-1/2R=hi+1n-12Ψri+1hi+1n-hin,where r denotes the ratio of successive gradients (e.g., ri=zi+1n-zin/zin-zi-1n) and Ψ is the flux limiter. For both the FORCE and C-FORCE schemes, Ψ≡0, and, for both the SLIC and C-SLIC schemes, a minmod slope limiter [14] is employed to ensure numerical stability. Second, the hydrostatic reconstruction approach proposed in [12] is employed to achieve the well-balanced property [15] and nonnegative water depth solutions; in this approach, a single-valued bed elevation at the cell-interface is defined as(19)zbi-1/2=0.5zi-1/2L-hi-1/2L+zi-1/2R-hi-1/2R.Third, the reconstructed face values are temporally evolved as(20)U¯i-1/2R=Ui-1/2R+ΔUi,U¯i+1/2L=Ui+1/2L+ΔUi.For both the FORCE and C-FORCE schemes, ΔU≡0, and, for both the SLIC and C-SLIC schemes, ΔU is calculated according to(21)ΔUi=Δt2ΔxFUi-1/2R-FUi+1/2L+Δt2Sin.Here, the source term vector Sin is evaluated similarly to (13) but calculated with the values before being temporally evolved (i.e., Ui-1/2L,R instead of U¯i-1/2L,R).


### 2.2. Numerical Test for the 1LSWEs

The benchmark dam-break problem is considered here, for which the analytic solution is available [16]. The simulation domain is 10 m long (x∈[0,10m]). The initial flow is at rest and the water surface elevation is initially prescribed as(22)zx,t=0=0.2m,x⩽5m,1.8m,x>5m.Open boundary conditions (i.e., zero gradients of variables at boundaries) are imposed at both ends of the simulation domain. The simulation domain is discretized by equally spaced 1600 mesh cells, determined based on a grid-independence study using the FORCE scheme with Cr=0.5. Here, Cr is defined as(23)Cr=min⁡λ~i-1/2ΔtΔx,where λ~i-1/2 is calculated by (17); the minimum is defined over the whole simulation domain. Note that, for this particular test, Cr is commonly set to be a large value (e.g., Cr≈0.5) in literatures [8] and thus the numerical diffusion in the predictions is insignificant. We choose this test here as it is good to demonstrate the different performances of the LF- and LLF-based numerical schemes in achieving Δt-independence solutions.In Figures3(a) and 3(b), the computed water surface profiles at t=0.5s by using the FORCE and C-FORCE schemes with different values of Cr are displayed. Also shown is the analytic solution of Stoker [16]. As seen, for the FORCE scheme, as Cr decreases, the numerical diffusion increases, causing a smearing of the wave profile; when Cr⩽0.05, the numerical diffusion becomes significantly large and gradually dominates other processes. It is surprising that, compared with the analytic solution, a larger value of Cr contributes to a more accurate solution when the FORCE scheme is employed for flux evaluations. However, when the C-FORCE scheme is used for flux evaluations, as can be seen in Figure 3(b), a Δt-independence solution is successfully achieved, which matches the analytic solution with reasonable accuracy.Figure 3
Computed water surface profiles by using the FORCE scheme (a) and the C-FORCE scheme (b) att=0.5s, for the one-layer dam-break problem test. The analytic solution of Stoker [16] is also shown.
(a)
(b)In Figures4(a) and 4(b), the computed water surface profiles at t=0.5s by using the SLIC and C-SLIC schemes with different values of Cr are illustrated. As seen, the solution does not converge when the SLIC scheme is used for flux evaluations, though the numerical diffusion is significantly reduced compared with that simulated by the FORCE scheme (cf. Figures 3(a) and 4(a)). The C-SLIC scheme successfully predicts a Δt-independence solution which matches the analytic solution almost exactly (see Figure 4(b)).Figure 4
Same as Figure3, but computed by using the SLIC scheme (a) and C-SLIC scheme (b), respectively.
(a)
(b)


## 2.1. Governing Equations and Numerical Methods for the 1LSWEs

The 1LSWEs, in a conservative vector form, can be written as(10)∂U∂t+∂F∂x=S,with vectors defined by(11)U=zqx,F=qxuqx+g2h2,S=0-gh∂zb∂x,where t denotes time; x denotes the streamwise coordinate; z and zb denote the water surface and bed elevations above a horizontal reference level zr; h=z-zb denotes the water depth; and qx=uh denotes the discharge per unit width with u being the depth-averaged velocity in the x-direction.An explicit finite-volume discretization for (10) gives(12)Uin+1-UinΔt=F^i-1/2-F^i+1/2Δx+S¯i.Here, the overbar on a term denotes that this term is evaluated by the temporally evolved values (described later in this section); the source term vector is discretized by using a central difference scheme as [12](13)S¯i=0-gh¯i-1/2R+h¯i+1/2L2zbi+1/2-zbi-1/2Δx;the flux vector F^ is evaluated by the FORCE or SLIC scheme, which reads uniformly(14)F^i-1/2=ωFi-1/2LW+1-ωFi-1/2LF,where(15a)Fi-1/2LW=FUi-1/2LW,(15b)Ui-1/2LW=12U-i-1/2L+U-i-1/2R-12ΔtΔxFU-i-1/2R-FU-i-1/2L,(16)Fi-1/2LF=12FU-i-1/2L+FU-i-1/2R-12λ~i-1/2U-i-1/2R-U-i-1/2L.Here, the superscripts “LW” and “LF” denote that the variables are related to the two-step Lax-Wendroff and Lax-Friedrichs fluxes, respectively; the superscripts “L” and “R” denote the reconstructed values at the left and right sides of the cell-interface, respectively. For both the FORCE and SLIC schemes, ω=0.5 and λ~≡Δx/Δt; for both the C-FORCE and C-SLIC schemes, ω is set to be 0.3 empirically and λ~ is set to be the maximum local characteristic speed as(17)λ~i-1/2=max⁡u+ghi-1/2L,u+ghi-1/2R.For face value reconstructions, first, the monotonic upstream-centered scheme for conservation laws (MUSCL) [13] is applied; for example,(18a)zi-1/2L=zin+12Ψrizin-zi-1n,zi-1/2R=zi+1n-12Ψri+1zi+1n-zin,(18b)hi-1/2L=hin+12Ψrihin-hi-1n,hi-1/2R=hi+1n-12Ψri+1hi+1n-hin,where r denotes the ratio of successive gradients (e.g., ri=zi+1n-zin/zin-zi-1n) and Ψ is the flux limiter. For both the FORCE and C-FORCE schemes, Ψ≡0, and, for both the SLIC and C-SLIC schemes, a minmod slope limiter [14] is employed to ensure numerical stability. Second, the hydrostatic reconstruction approach proposed in [12] is employed to achieve the well-balanced property [15] and nonnegative water depth solutions; in this approach, a single-valued bed elevation at the cell-interface is defined as(19)zbi-1/2=0.5zi-1/2L-hi-1/2L+zi-1/2R-hi-1/2R.Third, the reconstructed face values are temporally evolved as(20)U¯i-1/2R=Ui-1/2R+ΔUi,U¯i+1/2L=Ui+1/2L+ΔUi.For both the FORCE and C-FORCE schemes, ΔU≡0, and, for both the SLIC and C-SLIC schemes, ΔU is calculated according to(21)ΔUi=Δt2ΔxFUi-1/2R-FUi+1/2L+Δt2Sin.Here, the source term vector Sin is evaluated similarly to (13) but calculated with the values before being temporally evolved (i.e., Ui-1/2L,R instead of U¯i-1/2L,R).


## 2.2. Numerical Test for the 1LSWEs

The benchmark dam-break problem is considered here, for which the analytic solution is available [16]. The simulation domain is 10 m long (x∈[0,10m]). The initial flow is at rest and the water surface elevation is initially prescribed as(22)zx,t=0=0.2m,x⩽5m,1.8m,x>5m.Open boundary conditions (i.e., zero gradients of variables at boundaries) are imposed at both ends of the simulation domain. The simulation domain is discretized by equally spaced 1600 mesh cells, determined based on a grid-independence study using the FORCE scheme with Cr=0.5. Here, Cr is defined as(23)Cr=min⁡λ~i-1/2ΔtΔx,where λ~i-1/2 is calculated by (17); the minimum is defined over the whole simulation domain. Note that, for this particular test, Cr is commonly set to be a large value (e.g., Cr≈0.5) in literatures [8] and thus the numerical diffusion in the predictions is insignificant. We choose this test here as it is good to demonstrate the different performances of the LF- and LLF-based numerical schemes in achieving Δt-independence solutions.In Figures3(a) and 3(b), the computed water surface profiles at t=0.5s by using the FORCE and C-FORCE schemes with different values of Cr are displayed. Also shown is the analytic solution of Stoker [16]. As seen, for the FORCE scheme, as Cr decreases, the numerical diffusion increases, causing a smearing of the wave profile; when Cr⩽0.05, the numerical diffusion becomes significantly large and gradually dominates other processes. It is surprising that, compared with the analytic solution, a larger value of Cr contributes to a more accurate solution when the FORCE scheme is employed for flux evaluations. However, when the C-FORCE scheme is used for flux evaluations, as can be seen in Figure 3(b), a Δt-independence solution is successfully achieved, which matches the analytic solution with reasonable accuracy.Figure 3
Computed water surface profiles by using the FORCE scheme (a) and the C-FORCE scheme (b) att=0.5s, for the one-layer dam-break problem test. The analytic solution of Stoker [16] is also shown.
(a)
(b)In Figures4(a) and 4(b), the computed water surface profiles at t=0.5s by using the SLIC and C-SLIC schemes with different values of Cr are illustrated. As seen, the solution does not converge when the SLIC scheme is used for flux evaluations, though the numerical diffusion is significantly reduced compared with that simulated by the FORCE scheme (cf. Figures 3(a) and 4(a)). The C-SLIC scheme successfully predicts a Δt-independence solution which matches the analytic solution almost exactly (see Figure 4(b)).Figure 4
Same as Figure3, but computed by using the SLIC scheme (a) and C-SLIC scheme (b), respectively.
(a)
(b)


## 3. Performance Comparisons in Solving the 2LSWEs

In this section, the proposed scheme in solving the 2LSWEs is considered. A sketch of the two-layer shallow-water system is illustrated in Figure2(b).


### 3.1. Governing Equations and Numerical Methods for the 2LSWEs

The 2LSWEs are generally written as (10) with vectors defined by [17](24)U=h2qx2h1qx1,F=qx2u2qx2+g2h22qx1u1qx1+g2h12,S=0-gh2∂z1∂x0-gh1∂zb∂x-χgh1∂h2∂x,where hk(k=1,2) denote the layer depths, with k=1 and k=2, respectively, referring to the lower- and upper-layers (see Figure 2(b)); zb is the bed elevation above a reference level zr; zk(k=1,2) denote the layer surface elevations with respect to zr and have the relations z2=z1+h2 and z1=zb+h1; qxk(k=1,2) denote the layer discharges per unit width in the streamwise direction; uk=qxk/hk(k=1,2) denote the layer-averaged velocities; and ρk(k=1,2) are the layer densities and χ=ρ2/ρ1 is the density ratio.-
χ
g
h
1
∂
h
2
/
∂
x in (24) is a nonconservative product term, which has no definition around shocks and discontinuities and is tough to process in numerical discretizations. Following Spinewine et al. [18], and Lu et al. [17], the 2LSWEs are reformulated to make the nonconservative product term vanish. This is done by replacing the equations governing the lower-layer flow motions by a combination of the equations for the two-layer flow system as a whole; namely, the 2LSWEs are recast to (10) with vectors defined by(25)U=h2qx2hwqxw,F=qx2u2qx2+g2h22qxwFw,S=0-gh2∂z1∂x0-gH+1-χχh1∂zb∂x,where(26)H=h1+h2,hw=h2+1χh1,qxw=qx2+1χqx1,Fw=u2qx2+1χu1qx1+g2H2+1-χχg2h12.As the 1LSWEs and 2LSWEs are written in the same conservative form (see (10)), the 2LSWEs are solved similarly as that for the 1LSWEs presented in Section 2.1, with a difference in calculating λ~ in (17). In solving the 2LSWEs, for both the FORCE and SLIC schemes, λ~≡Δx/Δt, and, for both the C-FORCE and C-SLIC schemes, λ~ is set to be the approximately maximum local characteristic speed as [17, 18](27)λ~i-1/2=max⁡λ~1,i-1/2,λ~2,i-1/2,where(28a)λ~1,i-1/2=max⁡u1+gh1+χh2i-1/2L,u1+gh1+χh2i-1/2R,(28b)λ~2,i-1/2=max⁡u2+gh2i-1/2L,u2+gh2i-1/2R.


### 3.2. Numerical Test for the 2LSWEs

The internal dam-break problem is considered here, which is widely adopted to test the robustness and accuracy of numerical schemes in solving the 2LSWEs [19–21]. The configuration of this test is the same as that in the one-layer dam-break test presented in Section 2.2, but with an upper-layer superimposed on the lower-layer. The density ratio between the layers is set to χ=0.7. The initial condition for the lower-layer flow is the same as that defined in (22), and the initial surface level of the upper-layer is z2(x,t=0)=2m. The spatial resolution is set to be the same as that in the one-layer dam-break problem test.As the analytic solution for this test is unavailable, the numerical result computed by using the C-SLIC scheme withCr=0.5 and nx=6400 is used as a reference solution. Here, Cr is defined by (23) with λ~i-1/2 calculated by (27); nx denotes the number of the mesh cells. In Figure 5, the interface profile of the reference solution at t=1s is shown, along with the numerical results computed by Bouchut and de Luna [19] using a time-splitting method to decouple the 2LSWEs into two 1LSWEs with an entropy inequality for controlling the numerical instability (the relevant scheme is abbreviated as BD08 scheme hereafter), by Bouchut and Zeitlin [21] using a source-centered hydrostatic reconstruction scheme (BZ10 scheme), and by Dudzinski and Lukáčová-Medvid’ová [20] using a bicharacteristics theory based finite-volume evolution Galerkin scheme (DL13 scheme). A zoom view in the region x∈[4.9m,5.7m] is given on the bottom-right of Figure 5. As seen, spurious oscillations can be easily predicted in simulations for this test. The solution computed by the C-SLIC scheme matches the predictions computed by the BD08 and DL13 schemes well, but without oscillations (see the zoom view in Figure 5).Figure 5
Computed interface profile by using the C-SLIC scheme withCr=0.5 and nx=6400 at t=1s, along with those predicted by the BC08, BZ10, and DL13 schemes, for the two-layer internal dam-break problem test. A zoom view is shown on the bottom-right.In Figures6(a) and 6(b), the computed layer surface profiles at t=1.0s by using the FORCE and C-FORCE schemes with different values of Cr are shown. The relevant results for the SLIC and C-SLIC schemes are presented in Figure 7. Note that, in Figures 6 and 7, the reference solution is obtained by using the C-SLIC scheme with Cr=0.5 and nx=6400. As we compare Figure 3 with Figure 6 and Figure 4 with Figure 7, it is found that the effects of Cr on the results of the two-layer internal dam-break problem and the one-layer dam-break problem are quite similar: the numerical predictions do not converge when the FORCE or SLIC scheme is used, while a Δt-independence solution is achieved by using the C-FORCE or C-SLIC scheme.Figure 6
Computed layer surface profiles by using the FORCE scheme (a) and the C-FORCE scheme (b) att=1s, for the two-layer internal dam-break problem test. The reference solution is obtained by using the C-SLIC scheme with Cr=0.5 and nx=6400.
(a)
(b)Figure 7
Same as Figure6, but computed by using the SLIC scheme (a) and C-SLIC scheme (b), respectively.
(a)
(b)


## 3.1. Governing Equations and Numerical Methods for the 2LSWEs

The 2LSWEs are generally written as (10) with vectors defined by [17](24)U=h2qx2h1qx1,F=qx2u2qx2+g2h22qx1u1qx1+g2h12,S=0-gh2∂z1∂x0-gh1∂zb∂x-χgh1∂h2∂x,where hk(k=1,2) denote the layer depths, with k=1 and k=2, respectively, referring to the lower- and upper-layers (see Figure 2(b)); zb is the bed elevation above a reference level zr; zk(k=1,2) denote the layer surface elevations with respect to zr and have the relations z2=z1+h2 and z1=zb+h1; qxk(k=1,2) denote the layer discharges per unit width in the streamwise direction; uk=qxk/hk(k=1,2) denote the layer-averaged velocities; and ρk(k=1,2) are the layer densities and χ=ρ2/ρ1 is the density ratio.-
χ
g
h
1
∂
h
2
/
∂
x in (24) is a nonconservative product term, which has no definition around shocks and discontinuities and is tough to process in numerical discretizations. Following Spinewine et al. [18], and Lu et al. [17], the 2LSWEs are reformulated to make the nonconservative product term vanish. This is done by replacing the equations governing the lower-layer flow motions by a combination of the equations for the two-layer flow system as a whole; namely, the 2LSWEs are recast to (10) with vectors defined by(25)U=h2qx2hwqxw,F=qx2u2qx2+g2h22qxwFw,S=0-gh2∂z1∂x0-gH+1-χχh1∂zb∂x,where(26)H=h1+h2,hw=h2+1χh1,qxw=qx2+1χqx1,Fw=u2qx2+1χu1qx1+g2H2+1-χχg2h12.As the 1LSWEs and 2LSWEs are written in the same conservative form (see (10)), the 2LSWEs are solved similarly as that for the 1LSWEs presented in Section 2.1, with a difference in calculating λ~ in (17). In solving the 2LSWEs, for both the FORCE and SLIC schemes, λ~≡Δx/Δt, and, for both the C-FORCE and C-SLIC schemes, λ~ is set to be the approximately maximum local characteristic speed as [17, 18](27)λ~i-1/2=max⁡λ~1,i-1/2,λ~2,i-1/2,where(28a)λ~1,i-1/2=max⁡u1+gh1+χh2i-1/2L,u1+gh1+χh2i-1/2R,(28b)λ~2,i-1/2=max⁡u2+gh2i-1/2L,u2+gh2i-1/2R.


## 3.2. Numerical Test for the 2LSWEs

The internal dam-break problem is considered here, which is widely adopted to test the robustness and accuracy of numerical schemes in solving the 2LSWEs [19–21]. The configuration of this test is the same as that in the one-layer dam-break test presented in Section 2.2, but with an upper-layer superimposed on the lower-layer. The density ratio between the layers is set to χ=0.7. The initial condition for the lower-layer flow is the same as that defined in (22), and the initial surface level of the upper-layer is z2(x,t=0)=2m. The spatial resolution is set to be the same as that in the one-layer dam-break problem test.As the analytic solution for this test is unavailable, the numerical result computed by using the C-SLIC scheme withCr=0.5 and nx=6400 is used as a reference solution. Here, Cr is defined by (23) with λ~i-1/2 calculated by (27); nx denotes the number of the mesh cells. In Figure 5, the interface profile of the reference solution at t=1s is shown, along with the numerical results computed by Bouchut and de Luna [19] using a time-splitting method to decouple the 2LSWEs into two 1LSWEs with an entropy inequality for controlling the numerical instability (the relevant scheme is abbreviated as BD08 scheme hereafter), by Bouchut and Zeitlin [21] using a source-centered hydrostatic reconstruction scheme (BZ10 scheme), and by Dudzinski and Lukáčová-Medvid’ová [20] using a bicharacteristics theory based finite-volume evolution Galerkin scheme (DL13 scheme). A zoom view in the region x∈[4.9m,5.7m] is given on the bottom-right of Figure 5. As seen, spurious oscillations can be easily predicted in simulations for this test. The solution computed by the C-SLIC scheme matches the predictions computed by the BD08 and DL13 schemes well, but without oscillations (see the zoom view in Figure 5).Figure 5
Computed interface profile by using the C-SLIC scheme withCr=0.5 and nx=6400 at t=1s, along with those predicted by the BC08, BZ10, and DL13 schemes, for the two-layer internal dam-break problem test. A zoom view is shown on the bottom-right.In Figures6(a) and 6(b), the computed layer surface profiles at t=1.0s by using the FORCE and C-FORCE schemes with different values of Cr are shown. The relevant results for the SLIC and C-SLIC schemes are presented in Figure 7. Note that, in Figures 6 and 7, the reference solution is obtained by using the C-SLIC scheme with Cr=0.5 and nx=6400. As we compare Figure 3 with Figure 6 and Figure 4 with Figure 7, it is found that the effects of Cr on the results of the two-layer internal dam-break problem and the one-layer dam-break problem are quite similar: the numerical predictions do not converge when the FORCE or SLIC scheme is used, while a Δt-independence solution is achieved by using the C-FORCE or C-SLIC scheme.Figure 6
Computed layer surface profiles by using the FORCE scheme (a) and the C-FORCE scheme (b) att=1s, for the two-layer internal dam-break problem test. The reference solution is obtained by using the C-SLIC scheme with Cr=0.5 and nx=6400.
(a)
(b)Figure 7
Same as Figure6, but computed by using the SLIC scheme (a) and C-SLIC scheme (b), respectively.
(a)
(b)


## 4. Conclusions and Discussions

In this work, we have analyzed the problem of the LF-based schemes (e.g., the FORCE and SLIC schemes) regarding achieving time-step-independence solutions; it is found that a time-step-independence study cannot be performed for these schemes as the numerical solutions do not converge with successive finer time steps, owing to the excessive numerical diffusion at a small time step.Two convergence improved schemes, the C-FORCE and C-SLIC schemes, have been proposed to overcome the relevant problem related to the LF-based schemes. As representatives, the C-FORCE and C-SLIC schemes have been applied to solve the one-layer and two-layer shallow-water equations. Results show that the proposed schemes successfully achieveΔt-independence solutions. Also, the proposed schemes are verified to be capable of resolving sharp fronts well without oscillations.The proposed schemes are presented in 1D space in this work. The extension of the proposed schemes to 2D space is straightforward. Also, the use of the C-FORCE and C-SLIC schemes in solving other hyperbolic equation systems is of great interest and is left for future studies.



---

*Source: 379281-2015-11-23.xml*