1 Introduction

In the past several decades, there has been great progress and success of numerical high resolution treatments of hyperbolic conservation laws (see Refs. [1]-[5] and references therein). Since solutions to nonlinear conservation laws are discontinuous in general, high-order schemes have to deal with the conflict of accuracy and local oscillations near discontinuities. It leads to the family of high resolution schemes that have high-order accuracy in smooth regions of solutions and are free from oscillations near discontinuities with as few mesh points as possi-ble indicating the waves. The monotonic upstream-centered scheme for the conservation law (MUSCL) is a successful representative in resolving such conflictions^[3,6]. The most important ingredient, in addition to abiding by the conservation property with the finite volume frame-work, is the data reconstruction or slope reconstruction that is usually obtained through the numerical solutions at grid points with monotone algorithms, for instance, slope limiters, flux limiters or essentially non-oscillatory (ENO) and weighted essentially non-oscillatory (WENO) algorithms^[7-8]. Such reconstruction techniques are easily implemented and hence widely used in practice.

We refine a little bit the method of high resolution schemes. In smooth regions of solutions, the evolution of slopes (derivatives) should obey the governing equations, and the Lax-Wendroff (L-W) approach^[9] or the Cauchy-Kovalevskaya (C-K) theorem^[10] should be used directly. Just near discontinuities, the slopes are obtained by using the technology of data reconstruction widely adopted in high resolution schemes, and only there limiters are adopted to refrain from local oscillations. We propose a slope relaxation equation to adaptively connect these two ingredients in the two different types of regions. This procedure is different from that for MUSCL-type schemes which reconstruct slopes from the solutions at grid points uniformly in the domain under consideration, but the idea has something in common with the discontinuous Galerkin (DG) method for which slopes are essentially derived from different equations.

Now we explain how to implement our idea. Consider nonlinear conservation laws,

(1)

where the variable t is the time coordinate, u_t is the time derivative of u, x=(x₁, x₂, ···, x_d)^T is the space coordinate, d is the number of space dimensions, is the gradient operator, u=(u₁, u₂,···, u_m)^T is the vector of conservative quantities, and f=(f₁, f₂,···, f_d) is the associated flux. In smooth regions of u, the so-called derivative evolution system is obtained by differentiating (1) directly, i.e.,

(2)

where is the gradient of u, and is the Jacobian matrix of f_i. Based on the L-W approach, the solutions σ_i then act as the slopes of u to be adopted for constructing high-order schemes for (1) in the smooth region. In Ref. [11], this procedure is used to devise a scheme directly.

The direct application of (2) raises some issues near discontinuities. The solution (u, σ₁,···, σ_d) contains δ-shocks, i.e., the vector (σ₁, σ₂,···, σ_d) becomes the singular Dirac measure, as pointed out in Ref. [11], due to the discontinuities of the solution. From the numerical point of view, the slope σ may be very large. Thus, it is improper to directly apply (2) for the slope reconstruction in discontinuous regions, and it would be impossible to derive the total variation diminishing (TVD) or total variation bounded (TVB) property of the resulting schemes. Hence it is reasonable to construct the slopes using solutions u at grid points near discontinuities.

Denote by u_xi the derivative of the solution u with respect to x_i. u_xi is obtained by using the L-W or C-K approach in smooth regions and with the reconstruction technology by interpolating from the numerical solution near discontinuities. Then, we construct a relaxation model for computing the slopes, i.e.,

(3)

where μ_i > 0 is the relaxation parameter, adaptively chosen as O(1) in smooth regions and very small (say, O(Δx³), where Δx is the mesh size) near discontinuities, u_xi is discretized by some certain monotone algorithms, such as the limiters or the WENO method. From the viewpoint of numerical computation, the significance of the relaxation term is at least two-fold as follows:

(i) In smooth regions of the solution of (1), the difference between σ_i and u_xi is at most of the order O(Δx),

(4)

Then, as the parameter μ_i is chosen as O(1), the difference between (3) and (2) is O(Δx).

Therefore, we can use (3) to replace (2) for the computing of the slopes in smooth regions of solutions.

(ii) In the regions containing discontinuities, the relaxation parameter μ_i is chosen to be very small. Then, the right side of (3) is dominated, and σ_i is approximated by the approximation of u_xi obtained from monotone algorithms. Moreover, when (3) is used to compute the slopes, the relaxation takes effect of suppressing local oscillations.

Thus, (1) and (3) are used to devise any high resolution schemes for (1). The following ex-ample illustrates the necessity through such a procedure. Let us consider the Burgers7 equation with a continuous initial condition,

(5)

(6)

The second-order staggered central scheme in Ref. [3] is used here. Tables 1 and 2 show errors at t=0.12 between the exact solution and numerical solutions reconstructed by the solution of (3) and the minmod limiter, respectively. From the tables, it is observed that the scheme could keep the second-order accuracy for both L^∞ and L¹ errors as the solutions to (3) are used for reconstruction. Moreover, the solutions reconstructed by (3) (see Table 1) are more accurate than those reconstructed by the minmod limiter (see Table 2) when the same mesh is used.

To validate this procedure devising schemes for (1) rigorously, we use the framework of the staggered central scheme^[3] to discretize (1) and (3) just for simplicity. Other approaches can be used similarly. It is shown that the resulting scheme is TVB stable and satisfies a generalized entropy inequality that is less accurate than the second-order one from the classical discrete entropy inequality. Hence, it is entropy-satisfied. Here, we point out in passing that very few results are known about the entropy convergence of high resolution schemes. There are examples for the central Nessyahu-Tadmor (NT) scheme^[3], MUSCL-type schemes^[12-13], and DG schemes^[14]. All results for them are obtained with some assumptions.

This paper is organized into six sections. Besides the introduction here, the two-value scheme is formulated in Section 2, and the stability properties are shown in Sections 3 and 4. Numerical results are performed in Section 5 to validate the new scheme. Finally, we present a conclusive discussion in Section 6.

2 Basic formulation for one-dimensional conservation laws

We start with the initial value problem for one-dimensional conservation laws,

(7)

(8)

coupled with the relaxed derivative evolution system,

(9)

(10)

where the interval I=[a, b] is finite, u=(u₁, u₂, ···, u_m)^T is an m-vector of conservative quantities, f=(f₁, f₂, ···, f_m)^T is the vector of flux functions, σ=u_x is the derivative of u, and μ > 0 is a relaxation parameter. Denote by the corresponding Jacobian matrix.

In the design of second-or higher-order schemes for (7)-(8), one of the key points is how to approximate the derivatives of u in the reconstruction stage at each time step. Usually, they are constructed by the cell averages of u with certain monotonicity algorithms (limiters or ENO and WENO methods) so as to obtain non-oscillatory schemes^{[3-5, 15-16]}.

However, if the initial values of u only have finite discontinuities, so does u. In accordance with the L-W approach and the C-K theorem, it is natural to derive these derivatives in the smooth regions of the solution u directly from the equation,

(11)

which is obtained from the governing equation (7). On the other hand, through adaptively choosing the relaxation parameter μ, the slope σ in the discontinuous region would be computed by limiters or other monotonicity algorithms which are regarded as the approximations of u_x in (9).

Now we show details of the two-value scheme, developed by (7)-(9). The framework of the staggered central scheme, the NT scheme, is adopted here^[3]. Certainly, any other type of schemes (e.g., the Godunov-type) can be also used to solve (7) and (9). The reason of making such a choice is that no Riemann problems are solved, and field-by-field characteristic decompositions are avoided so that we can concentrate on the new features of the resulting scheme. The corresponding stability properties can also be easily established.

The idea of the staggered central scheme is to use staggered control volumes for avoiding solving Riemann problems. Assume that initially u₀(x) is defined on a finite interval I=[a, b]. Take an equally spaced grid over I with cells The center of I_j is , where Δx=(b — a)/(J + 1) is the spatial mesh size, and .

Given the initial data {u_jⁿ} and {σ_jⁿ} for (u, σ) at the time level t=t_n, a piecewise linear approximation u(x, t_n) for u and a piecewise constant approximation σ(x, t_n) for σ are reconstructed at t=t_n, i.e.,

(12)

Let t_n+1=t_n + Δt_n. Use u(x, t_n) and σ(x, t_n) as initial values, and take (x_j, x_j+1) × [t_n, t_n+1) as the control volume for the integrations of (7) and (9). We arrive at

(13)

(14)

where the ratio satisfies the Courant-Friedrichs-Lewy (CFL) constraint,

(15)

Here, ρ(f', (u_jⁿ)) is the spectral radius of f', (u_jⁿ). The numerical fluxes and g_jⁿ in (13) and (14) should approximate the time averages of f (u(x_j, t)) and f'(u(x_j, t))σ(x_j, t) over (t_n, t_n+1). Moreover, the two functions are continuous with t ∈ (t_n, t_n+1) under the constraint (15). Hence, the time averages of f (u(x_j, t)) can be integrated approximately by the midpoint rule at the expense of the local truncation error O(Δt_n²), and the time averages of the f' (u(x_j, t)) can be integrated by the rectangle approximation at the expense of O(Δt_n), if a second-order scheme is preferred. The two numerical fluxes and g_jⁿ are given by

(16)

(17)

where the L-W approach is used.

Define

(18)

Then, the numerical solution of (11). By (14), can be solved as

(19)

Here, u_x in (19) is discretized properly with a monotonicity algorithm for suppressing oscilla-tions. For example, we use the following analysis:

(20)

to approximate We would like to emphasize that the monotonicity algorithm is not needed in the smooth regions of solutions.

The relaxation parameter μ is chosen to be

so that the slope is dominated by in smooth regions and close to in regions containing discontinuities.

This is consistent with the methodology of high resolution schemes. Although we choose different ways to compute the slopes in smooth regions and regions containing discontinuities separately, the use of (9) uniformly combines two methods by the relaxation term so as to choose proper slopes in the two situations.

Thus, is derived at the new time level t=t_n+1. In the diagram, we obtain

(21)

We call the resulting scheme a two-value scheme.

To summarize, the introduction of (9) brings new features as follows:

(i) System (9) can effectively describe the propagation of gradients of conservation quan-tities. It is, therefore, natural to use σ computed by (9) at the data reconstruction stage in smooth regions for high resolution MUSCL-type schemes, instead of the reconstruction from the numerical solution at grid points for the whole computing region under consideration.

(ii) Thanks to the relaxation term in (9), the smooth and discontinuous parts are combined adaptively so that the monotonicity algorithm just works near discontinuities, and the accuracy reduction is avoided in smooth regions.

(iii) The resulting scheme is second-order accurate and TVB stable, and satisfies cell-entropy inequalities. Hence it is entropy stable.

Remark 1  Although the algorithm presented here is not exactly the same as the Hermite-Harten-Lax-van Leer-contact (HLLC) scheme in Ref. ^[11], the stability analysis in the next section still could be regarded as a supportive evidence that the algorithm in Ref. [11] is entropy-convergent.

3 Some properties of two-value scheme

In this section, we are going to prove that if the relaxation parameter μ is properly cho-sen, solutions computed by the two-value scheme (13)-(19) essentially satisfy the maximum principle, and the scheme is TVB stable for scalar conservation laws,

(22)

Here, u is a scalar variable. Then, it can be proved that these approximations converge to a weak solution to (22) as Δx→ 0 by use of the framework of compactness (see Ref. [17]). With the establishment of the entropy inequalities, we will find out that numerical solutions are just the approximations of the unique entropy solution.

In the previous section, the two-value scheme (13)-(19) is presented in detail. Specially, it can be shown that the numerical flux is obtained by a linear approximation of f(u(x_j, t)) and is obtained from the exact solution to some approximate conservation law (see Ref. [18]), i.e.,

(23)

if the flux function f ∈ C² is convex for scalar conservation laws.

Before we establish stability theorems of the proposed scheme, the smooth and discontinuous cells of solutions will be defined. Assume that , , and M be a constant depending on the second-order derivatives of u₀(x) (if they exist).

Definition 1  Fix the index j at the time level t=t_n+1. If the following inequalities hold:

(24)

(25)

where and then the interval [x_j, x_j+1] is a smooth cell for u(x, t_n+1). Otherwise, [x_j, x_j+1] is a discontinuous cell.

Recall that is defined in (19) as follows:

(26)

(27)

(28)

By Definition 1, in the smooth cells, is equivalent to within an O(Δx) error,

(29)

It can be proved that the sum of with respect to the index j is bounded by the total variation of u(x, t_n).

Lemma 1  Assume that the CFL condition

(30)

is met, where is a uniform constant, and is a positive constar^[18] Then, the following estimate holds:

(31)

Proof  According to the definition (27), it follows that

(32)

Based on Definition 1 and Lemma 1, the TVB property of the proposed scheme is proved with the following theorem.

Theorem 1  Let , and . The CFL condition (30) is met. Denote T=(N +1)Δt and L=b -a. Then, there exists a constant

If

(33)

holds, where with being a constant only dependent on M, M_B, T, and L, and the spatial step satisfies

(34)

then the numerical solutions computed by (13) -(19) essentially satisfy the maximum principle, i.e.,

(35)

and the total variations of the approximate solutions reconstructed by (12) satisfy

(36)

Remark 2  (i) The selection of μ_c does not affect the stability of the two-value scheme here. In the next section, we will present a constraint on μ_c to ensure that a (generalized) discrete entropy inequality holds.

(ii) Based on the property of solutions to the conservation laws, we know that there are only countably discontinuous cells for u at each time level.

Proof  Obviously, (35)-(36) hold for n=0. Denote

(37)

(38)

In the above notations, the indices are replaced by j only for simplicity as p is odd. Assume that M_B=B₀≤B₁≤···≤B_n≤M^B + CΔx with C being a constant and V₀≤V₁≤···≤V_n≤M_T. We can write σ_jⁿ as

(39)

If the cell I_j is a smooth interval, the relaxation parameter μ=μ_c. Thus can be estimated by

(40)

Otherwise, I_j is a discontinuous cell. Then, μ_j=μ_d. In the light of (33) and (30), we arrive at

(41)

if the mesh size satisfies

Based on the above analysis, is bounded by

(42)

Hence, , where . By analogy, B_n+1 is bounded by

If the constant C is picked by C=2TC₀f'/(2M_B), then we have

Finally, the essential maximum principle (35) holds.

So far, we have proved that the numerical solutions computed by the proposed scheme essentially satisfy the maximum principle. Next, we will show the TVB property of the two-value scheme (13)-(19).

For t=t_n+1 assume that , j=k, k+1··· k + l − 1 and , k + l, ··· k+l-1 for the moment (the other cases could be proved similarly), i.e., the cells {I_k, I_k+1, ···, I_k+l-1} are discontinuous cells, and the others are continuous cells. Then, the sum of slopes is bounded by

(43)

Here, L=b — a is the length of the interval I. Recall that the data u(x, t_n+1) are defined as

The total variation of u(x, t_n+1) can be estimated as

(44)

The last inequality holds due to (43), and is a constant only dependent on M, M_B, T, and L. Moreover, the middle term in (44) could be divided into three parts,

(45)

(24) and (30) are used for the first part of (45) to arrive at

(46)

Similarly, the second term in (45) is computed by

(47)

The last part in (45) can be controlled by

(48)

(49)

Based on (33) and (43), plugging (45)-(49) into (44) yields

where is a constant only dependent on M, M_B, T, and L. Thus, for t=t_N+1=T,

(50)

Here, . Thus, the proposed two-value scheme has been proved to be TVB stable.

4 Cell-entropy inequality for two-value scheme

So far, Theorem 1 demonstrates that solutions computed by the two-value scheme (13)-(19) essentially satisfy the maximum principle and are TVB stable. Hence, they converge to a weak solution to (22) in view of the L-W theorem. As long as the following cell entropy inequality is proved:

(51)

this weak solution must be the unique entropy solution to (22)^[17]. Here, S(u) is a convex entropy function, and is the numerical entropy flux which is consistent with the corresponding entropy flux,

This goal will be achieved by using the strategy in Ref. [3].

Denote and . In accordance with Ref. [3], define a grid function g,

and define some auxiliary functions,

Then, we have the following identity:

where the residual term is given by

and the numerical flux

is consistent with the entropy flux F (u). Therefore, (51) is valid if

(52)

is verified.

Theorem 2  Consider the scalar conservation law and its derivative system (7)-(9). Let , and f ∈ C² satisfies f" > 0. The CFL condition (30) is met. The relaxation parameters are selected as

(53)

Then, the approximate solutions computed by (13)-(19) satisfy the following approximate cell entropy inequality:

(i) If f (u)=au is linear, the numerical entropy inequality (51) holds with a remainder,

(54)

(ii) If f(u) is convex (and of course nonhnear), in (20) is redefined as follows:

(55)

where

(56)

Then, the cell entropy inequality (54) still holds for the numerical solutions.

Remark 3  It should be noted that the numerical entropy inequality in (54) is not the traditional numerical entropy inequality which has the extra residue O(Δx²Δt^α). Fortunately, the remainder will go to zero as Δx→0 and Δt→0 after the summation with respect to j and n, which implies that the entropy inequality still holds in the limit. It means that the numerical solutions computed by the new scheme converge to the unique entropy solution.

Proof  We will prove it only in terms of the quadratic entropy function . Indeed, thanks to Ref. [19], the entropy inequality (51) with one convex entropy can guarantee the satisfaction of entropy inequalities for all convex entropy functions.

Since both u_jⁿ(s) and u_jⁿ(r, s) lie between u_jⁿ and u_j+1ⁿ, we know, according to the definition of g(v), that

Then, (52) equals

(57)

Assume that . Then, we only need to prove

(58)

Note that g_jⁿ and g_j+1ⁿ can be rewritten as

where . Then, the first two terms in (58) can be written as

(59)

and the third term can be similarly computed as

(60)

Since the sum of in (59)-(60) satisfies

inserting (59) and (60) into (58) arrives at

(61)

We note that

(62)

If we choose μ_c=Δt^1+α, α < 1, and , (61) finally becomes

(63)

Here, the constant C is only associated with M_B, M, L, and T. Since σ_j^{n, l} and σ_j+1^{n, +1} agree in sign with (63) has the upper bound,

(64)

The last inequality holds due to the CFL condition (30).

If f is a linear function, f"(u)=0. By (20), (64) finally becomes

which means that the numerical entropy inequality (54) holds.

On the other hand, if f is a nonlinear function, and there is one extra term in the estimate (64). To obtain the numerical entropy inequality (54), a stronger monotone algorithm (55)-(56) in Ref. [20] is applied. This is in contrast with the linear case for which the monotone algorithm (20) is sufficient to guarantee (64). The modification is consistent with the assumption in Ref. [3] since we use the central discretization.

5 Numerical results

In this section, numerical results of some benchmark problems are presented to display the performance of the current method. We want to emphasize that the main purpose of this paper is to establish a derivative evolution model for MUSCL-type schemes and explain why it is reasonable by the stability analysis. We choose the staggered central scheme to discretize the proposed models just for simplicity. However, other frameworks are available, and the better resolution can be derived.

The numerical results derived by the two-value relaxed staggered central scheme (TVRX) are shown in this section compared with the traditional finite volume NT scheme and the generalized Riemann problem (GRP) method. We take the Courant number cfl=0.475 in practice. The time step is determined by the CFL condition. The relaxation parameter μ is selected as 1/Δx³ in the smooth regions and Δx³ in the discontinuous regions.

5.1 Scalar conservation laws in one-dimension

Example 1  Accuracy test

Here, we solve Burgers’ equation with a periodic boundary condition,

(65)

(66)

The exact solution is smooth up to t=2/π. Therefore, the output time T=0.6 and T=1.5 is chosen to present the solutions before and after the shock formation, respectively.

Since there is no singular point at T=0.6, (9) can be used to develop σ, which makes our scheme second-order accurate even at extreme points. We list the errors in Table 1, which proves that this two-value scheme keeps the second-order accuracy both in L₁ and L_∞ and norms.

At T=1.5, the shock appears. Here, (11) is used to keep a non-oscillatory scheme. Nu-merical results are shown at T=0.6 (see Fig. 1(a)) and T=1.5 (see Figs. 1(b) and 1(c)). It is shown that our scheme works well to resolve the shock wave.

Example 2  Rarefaction wave

Here, we solve the same equation (65) as that in Example 1 but with the Riemann initial data,

(67)

in order to show how our scheme works well on another typical wave, the rarefaction wave. Solutions at T=0.2 are displayed in Fig. 2. By the numerical results, we can find out that our scheme also displays good performance for the rarefaction wave. And if the relaxation parameter μ is properly selected, the two-value scheme is non-oscillatory.

5.2 Compressible Euler systems

In this subsection, numerical results of the compressible Euler systems for a polytropic gas are given, with γ=1.4,

(68)

where ρ, u, , and E are the density, velocity, pressure, and total energy, respectively.

5.2.1 One-dimensional cases

Example 3  Accuracy test

We consider the following initial value problem for the Euler equations subject to the initial data:

(69)

to test accuracy of the current two-value scheme (see Table 3) compared with the finite vol-ume NT scheme (see Table 4). The convergence rate coincides with the results for the scalar equations.

Example 4  Sod problem

We consider the shock tube problem^[21]. Initially, the gas is at rest with ρ=1, u=0.75, p=1 for 0≤x≤0.3 and ρ=0.125, u=0, p=0.1 for 0.3 < x < 1. Its solution is very mild and consists of a left rarefaction, a contact, and a right shock.

Numerical results are shown at the time T=0.2 in Fig. 3. 100 grid points are used in the computation. It is shown that our scheme works well in the smooth area, and performs similarly with the NT scheme at discontinuities. The GRP method works better than the above two schemes. This is understandable since the GRP method is designed based on its second-order generalized Riemann problem solver.

Example 5  Low density and internal energy problem

The initial condition is given by ρ=1, u=—2, p=0.4 for 0≤x≤0.5 and ρ=1, u=2, p=0.4 for 0.5 < x≤1^[22]. The solution to this problem consists of two strong rarefaction waves and a trivial stationary contact discontinuity.

Numerical results are shown at the time t=0.15 in Fig. 4. Through the numerical result, it is easy to see that our scheme basically performs well for the low density problem.

5.2.2 Two-dimensional cases

Actually, the two-dimensional simulation of the Euler equations is much more interesting and challenging. Here, we are going to present two classical test problems.

Example 6  Accuracy test

We solve the following initial value problem for the two-dimensional Euler equations to test the accuracy of our scheme:

(70)

The errors and convergence rates are shown in Table 5 compared with the NT scheme using the minmod limiter in Table 6.

Example 7  Double Mach reflection of strong shock

This problem was first studied by Woodward and Colella^[23] and later by many others. We use exactly the same initial conditions as those in Ref. [23]. The computational domain for this problem is [0, 4] × [0, 1], however, we only show the results on [0, 3] × [0, 1]. The reflecting wall lies at the bottom of the computational domain starting from , a right-moving Mach number 10 shock is positioned at y=0 and makes a 60° angle with the x-axis. For the bottom boundary, the exact post-shock condition is imposed for the part from x=0 to and a reflective boundary condition is used for the rest. At the top boundary of our computational domain, the flow values are set to describe the exact motion of the Mach number 10 shock. Inflow/outflow boundary conditions are used for the left/right boundaries. The numerical result with an 80 x 240 grid is displayed in Fig. 5. Details can be found in Ref. [23].

6 Conclusions

The aim of this paper is to establish a reliable derivative evolution model by introducing the slope relaxation for high resolution schemes of MUSCL-type. The schemes developed by the proposed model are high-order accurate in smooth regions of solutions and refrain from local oscillations near discontinuities. The relaxation term combines these two different ingredients together so that the L-W or C-K theorem is used by using the governing equation, and limiters are used near discontinuities for suppressing oscillations. The staggered central type discretiza-tion is used here for easy analysis. Thus the TVB property and cell entropy inequality are derived to show that the slope relaxation model is implementable. The numerical results show that the new scheme can deal with both the smooth solution and discontinuities.

The significance of this paper is to show the reliability of the slope relaxation model but not provide the specific discretization method. We can certainly use other discretizations such as using the HLLC solver or even the exact Riemann solver in Ref. [11]. Then, this framework is much more closely related to that in Ref. [11]. Hence, two-dimensional problems on unstructured meshes can be solved. However, the stability analysis is much more difficult. Moreover, the same idea can also be applied for even higher order schemes.

Acknowledgements The first author appreciates Prof.Haiyang HUANG and Prof. Shuanghu WANG for their kind guidance.