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.
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 :Here, is the wave propagation speed and, without loss of generality, is assumed to be positive. An explicit finite-volume discretization for (1a) giveswhere the superscripts “” and “” denote the values at the old and new time steps, respectively; and denote the time step and space interval, respectively; denotes the numerical flux and, for the LF scheme, it reads [1]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]where is the CFL number and, for the linear convection problem,In Figure 1, is plotted against . For clear view purpose, the relation in the region is plotted on a semilog scale (see Figure 1(a)). For convenience, in Figure 1, is assumed so that (see (5)). As seen in Figure 1, increases as decreases. When , , 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 , ), a -independence study is required to choose a reasonable time step. However, due to the rapidly increased numerical diffusion when decreases, the solution computed by using the LF scheme does not converge and thus a -independence study cannot be implemented. Failing to achieve a -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 .
(a)
(b)
We observe that the problem, related to the LF scheme being unable to achieve a -independence solution, can be resolved by the local LF (LLF hereafter) scheme [11]. The LLF scheme calculates the numerical flux aswhere denotes the local bound of the wave propagation speed and is determined byand, in solving (1a), (1b), and (1c), . 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:which introduces a numerical viscosity ofFrom (9) (also see Figure 1), it is seen that increases with decreasing . However, has a maximum and limiting value of when . Therefore, the uncontrolled excessive numerical viscosity in the LF scheme is resolved and a -independence study can be performed. It is worth remarking that, besides the advantage in achieving a -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 -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 Section 2, 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 Figure 2(a).
(a)
(b)
2.1. Governing Equations and Numerical Methods for the 1LSWEs
The 1LSWEs, in a conservative vector form, can be written aswith vectors defined bywhere denotes time; denotes the streamwise coordinate; and denote the water surface and bed elevations above a horizontal reference level ; denotes the water depth; and denotes the discharge per unit width with being the depth-averaged velocity in the -direction.
An explicit finite-volume discretization for (10) givesHere, 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]the flux vector is evaluated by the FORCE or SLIC scheme, which reads uniformlywhereHere, 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, and ; for both the C-FORCE and C-SLIC schemes, is set to be empirically and is set to be the maximum local characteristic speed asFor face value reconstructions, first, the monotonic upstream-centered scheme for conservation laws (MUSCL) [13] is applied; for example,where denotes the ratio of successive gradients (e.g., ) and is the flux limiter. For both the FORCE and C-FORCE schemes, , 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 asThird, the reconstructed face values are temporally evolved asFor both the FORCE and C-FORCE schemes, , and, for both the SLIC and C-SLIC schemes, is calculated according toHere, the source term vector is evaluated similarly to (13) but calculated with the values before being temporally evolved (i.e., instead of ).
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 (). The initial flow is at rest and the water surface elevation is initially prescribed asOpen 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 . Here, is defined aswhere is calculated by (17); the minimum is defined over the whole simulation domain. Note that, for this particular test, is commonly set to be a large value (e.g., ) 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 -independence solutions.
In Figures 3(a) and 3(b), the computed water surface profiles at by using the FORCE and C-FORCE schemes with different values of are displayed. Also shown is the analytic solution of Stoker [16]. As seen, for the FORCE scheme, as decreases, the numerical diffusion increases, causing a smearing of the wave profile; when , the numerical diffusion becomes significantly large and gradually dominates other processes. It is surprising that, compared with the analytic solution, a larger value of 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 -independence solution is successfully achieved, which matches the analytic solution with reasonable accuracy.
(a)
(b)
In Figures 4(a) and 4(b), the computed water surface profiles at by using the SLIC and C-SLIC schemes with different values of 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 -independence solution which matches the analytic solution almost exactly (see Figure 4(b)).
(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 Figure 2(b).
3.1. Governing Equations and Numerical Methods for the 2LSWEs
The 2LSWEs are generally written as (10) with vectors defined by [17]where denote the layer depths, with and , respectively, referring to the lower- and upper-layers (see Figure 2(b)); is the bed elevation above a reference level ; denote the layer surface elevations with respect to and have the relations and ; denote the layer discharges per unit width in the streamwise direction; denote the layer-averaged velocities; and are the layer densities and is the density ratio.
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 bywhere
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, , and, for both the C-FORCE and C-SLIC schemes, is set to be the approximately maximum local characteristic speed as [17, 18]where
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 . 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 . 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 with and is used as a reference solution. Here, is defined by (23) with calculated by (27); denotes the number of the mesh cells. In Figure 5, the interface profile of the reference solution at 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 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).
In Figures 6(a) and 6(b), the computed layer surface profiles at by using the FORCE and C-FORCE schemes with different values of 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 and . As we compare Figure 3 with Figure 6 and Figure 4 with Figure 7, it is found that the effects of 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 -independence solution is achieved by using the C-FORCE or C-SLIC scheme.
(a)
(b)
(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 -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.
Appendix
A View of the LLF Scheme from the Method of Characteristics
In this appendix, we provide a physical explanation of the LLF scheme (defined in (6)) from the method of characteristics. For brevity, the following discussion is presented for the LLF scheme for solving the linear convection problem (defined in (1a), (1b), and (1c)); for this problem, (2) reduces to (8).
Under a time-space domain depicted in Figure 8, one attains along the characteristics line , , where is the intersection point between the characteristics line and the line . can be linearly interpolated (of the first-order accuracy) from the known values and ; namely,and, thus,which is identical to (8). This suggests that when solving the linear convection problem, the LLF scheme is equivalent to the first-order accurate characteristics scheme. The interpolation algorithm in (A.1) also indicates that the LLF scheme is stable when in solving the linear convection problem.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported in part by National Natural Science Foundation of China (Grant no. 51409195), by the Fundamental Research Funds for the Central Universities (Grant no. 2042014kf0068), and by the China Postdoctoral Science Foundation (Grant no. 2014M550408).