numerical integrations of systems of conservation laws of ... - CiteSeerX
c ?? Society for Industrial and Applied Mathematics
SIAM J. APPL. MATH. Vol. ??, No. ??, pp. 1{16, ?? ??
??
NUMERICAL INTEGRATIONS OF SYSTEMS OF CONSERVATION LAWS OF MIXED TYPE SHI JINy
Abstract. The systems of conservation laws have been used to model dynamical phase transitions in, for example, the propagating phase boundaries in solids and the van der Waals uid. When integrating such mixed hyperbolic-elliptic systems the Lax-Friedrichs scheme is known to give the correct solutions selected by a viscosity-capillarity criterion except a spike at the phase boundary which does not go away even with a re ned mesh [15]. We identify the source of this spike as an inconsistency between the Lax-Friedrichs discretization and the viscosity-capillarity equations, and show a simple change of variable that can eliminate this spike. We then implement a high resolution scheme for the mixed type problems that select the same viscosity-capillarity solutions as the Lax-Friedrichs scheme with higher resolutions. Furthermore, a exibility in the ( rst order) scheme is used to obtain solutions for a wide range of the viscosity-capillarity equations. Key words. Systems of conservation laws of mixed type, viscosity-capillarity admissibility criterion, the Lax-Friedrichs scheme, shock capturing schemes. AMS(MOS) subject classi cations. 35L65, 65M06 1. Introduction. Systems of conservation laws of mixed type have been used to describe certain dynamic phase transitions. The best known mixed type system is the van der Waals gas at suciently low temperature [16]. Mixed type system is also used for the propagating phase boundaries in solids, where the stress-strain relation is not monotone [4]. In this paper we are concerned with the equations of isothermal motion for a compressible uid in Lagrangian coordinates wt ? u x = 0 ; (1.1) ut + p(w)x = 0 ; where u is the velocity, w is the speci c volume, and p is the pressure. Here the equation of state is a non-monotone function: p0 (w) < 0 if w < m or w > M ; the hyperbolic region; (1.2) p0 (w) > 0 if m < w < M ; the elliptic region:
p
Notice that the eigenvalues of (1.1) are ?p0(w), which become complex in the elliptic region and result in the classical Hardmard instability. We say that the solution is in phase 1, 2 or 3 at time t if p(w) lies in the interval (0; m); (m; M) or (M; 1) respectively. In the hyperbolic region it is expected that shocks may form [6]. A phase boundary is a shock that connects the solutions in di erent phases. When selecting the admissible weak solution to Eqs.(1.1), it is well known that the viscosity criterion is too restrictive to be used since it rules out propagating phase boundaries near the equilibrium co-existence line [13]. Instead, Slemrod [13] and Truskinovsky [17] independently proposed to use the viscosity-capillarity criterion. y School of Mathematics and Center for Dynamical System and Nonlinear Studies, Georgia Institute of Technology, Atlanta, GA 30332, USA ([email protected]). This work was supported by NSF grant no DMS-9404157. 1
2
shi jin
The viscosity-capillarity regularization for (1.1) is: wt ? ux = 0 ; (1.3) ut + p(w)x = uxx ? !2 wxxx : Here > 0 denotes the viscosity and !2 denotes the coecient of interfacial capillarity. Both and ! are positive constants. A ouf and Ca isch conducted numerical studies on the Riemann problem solutions and stability for (1.3) [2]. More recently Cockburn and Gau developed numerical schemes for (1.3) that work for all ! > 0 [3]. In addition, Abeyaratne and Knowles have proposed to use a kinetic relation to select the unique solution [1]. We are concerned with the numerical integrations of the mixed type problem (1.1) and (1.2). In this pursuit of particular interest is the work of Slemrod [14], who showed that the Lax-Friedrichs scheme for (1.1) gives the solution of (1.3) with ! = 41 . To see this, letting ! = 41 and introducing the change of variable (1.4) v = u ? 21 wx ; then Eqs.(1.3) can be rewritten as wt ? vx = 12 wxx ; (1.5) vt + p(w)x = 21 vxx : The Lax-Friedrichs discretization on (1.1), as shown in [14], is a second order discretization of (1.5) for some under an appropriate Courant-Friedrichs-Lewy (CFL) condition. Thus the Lax-Friedrichs scheme gives the weak solution of (1.1) selected by the viscosity-capillarity criterion (1.3) with ! = 41 . This analysis was later con rmed by the numerical experiments done by Slemrod and Flaherty [15]. Unfortunately the Lax-Friedrichs scheme always produces a spurious spike near the phase boundary, the nature of which is still unclear. To the author's knowledge no method based on a direct discretization of system (1.1) has been developed that can eliminate the spike. In recent years considerable interests have been shown to study the mixed type problem numerically. Despite its overdissipation and the unpleasant spike, the LaxFriedrichs scheme is still the most frequently used method for mixed type problems [7]. Attempts have been made to design modern high resolution shock capturing methods such as the total-variation-diminishing (TVD) methods and the total-variationbounded (TVB) methods. The direct implementation of these methods to mixed type problems faces many diculties, since these methods are usually based on the Riemann solver and local eld-by- eld characteristic decompositions. It has been shown by Shearer that the Riemann problem of (1.1) does not have a unique solution [10]. Also the local eld-by- eld decomposition is not available in the elliptic region. Shu found a hyperbolic ux splitting, and implemented the essentially-nonoscillatory methods (ENO) to each of the split uxes respectively. Local eld-by- eld decomposition is applied in the hyperbolic region while in the elliptic region the reconstruction is on the primitive variables. Pitman and Ni applied the high order Godunov methods to a relaxation regularization of the mixed type equations and the results are nonoscillatory [9]. Their method may select solutions di erent from those selected by the viscosity-capillarity solution and their method has to resolve the small relaxation time. In this paper we intend to answer some basic questions regarding to numerical discretization of the mixed type system (1.1). Speci cally, we would like to address
numerical solutions of mixed type systems
3
the following issues: 1) Identify the real source of the numerical spike produced by the Lax-Friedrichs scheme and show how to eliminate it. 2) Implement a second order shock capturing scheme for mixed type system (1.1) that pick up the solution selected by the viscosity-capillarity criterion with higher accuracy and sharper shocks and phase boundaries. 3) Attempt to extend these methods to handle a wide range of ! 2 (0; 14 ], which expands the applicability of our methods for mixed type problems. In an e ort to understand the relation between the solution of (1.1) and the capillarity coecient !, Cockburn and Gau have conducted a numerical convergence study for the viscosity-capillarity equations (1.3) as ! 0 with di erent ! [3]. They constructed a new system equivalent to (1.3) under (1.4) with enough free parameters so as to work for all ! > 0. The numerical discretization is based on this new system, and the change of variable (1.4) is performed to convert to the solutions of (1.3). Interesting enough is that their results do not have the spike, and their method works e ectively for all ! > 0. Their method, however, is based on a discretization for a viscosity-capillarity equations, thus is at most rst order to (1.1). Our work is partly motivated by their study, but we emphasize more on the practical side. We aim at developing modern, robust and high order numerical schemes based on a direct discretization of the original system (1.1) (rather than the viscosity-capillarity equations) that automatically pick up the correct viscosity-capillarity solutions. First we point out that the numerical spike produced by the Lax-Friedrichs scheme is caused by an inconsistency of the scheme with the viscosity-capillarity equations (1.3), in the sense that the transformation (1.4) that leads from (1.3) to (1.5), or vice versa, is not inherently possessed in the Lax-Friedrichs discretization. A numerical change of variable (1.6) u = v + 21 wx in the last step of the integration should be applied to obtain satisfactory results. Such a change of variable involves a discretization of the right hand side of (1.6) which must be consistent to the discretization of (1.1). If this is performed suitably then the numerical solution, while maintaining other features of the Lax-Friedrichs scheme, does not have the spike. In our e ort to implement a second order shock capturing method for (1.1), we use the relaxation schemes developed by the author and Xin [5]. The relaxation schemes are natural candidates for such problems since they do not need the Riemann solver, nor the local eld-by- eld decomposition. At the relaxation level, the Riemann solver and the eld-by- eld decomposition are trivially obtained due to the linear convection of the relaxation system and an operator splitting that separates the linear convection from the nonlinear sti source terms. The limiting scheme as the relaxation time approaches zero is a TVD scheme with reconstruction on primitive variables. Due to these special features the relaxation schemes have many advantages when applied to the mixed type problems. Here by the relaxation schemes we refer to the \ = 0" version of the relaxation scheme (called the relaxed schemes in [5]). Here is the relaxation time. This zerorelaxation scheme can also be derived by applying van Leer's MUSCL scheme [18] to a ux-splitting procedure, similar to Shu's work in [11]. For the rst order scheme, by
4
shi jin
performing the modi ed equation analysis, we nd that it selects the same viscositycapillarity solution as the Lax-Friedrichs scheme. We also nd that with higher order time integrator such as the second order Runge-Kutta method such a selection of the viscosity-capillarity solution can be made with a more relaxed CFL condition. Since our second order scheme reduces to the rst order scheme at discontinuities it may select the same shocks and phase boundaries as the rst order scheme (or the Lax-Friedrichs scheme), but with sharper resolutions. When forming the relaxation scheme in [5] we chose a positive diagonal constant matrix A that satis es Liu's subcharacteristic condition [8]. When passing to the zero relaxation limit A enters the numerical viscosity coecients. For problems sensitive to the parameter !, the exibility of A allows us to obtain numerical results for a wide range of ! 2 (0; 41 ] with the rst order scheme, thus expand the applicability of the scheme to a wider class of mixed type problems. However, we did not observe the same convergence for the second order scheme. We now outline this paper. In section 2 we review the modi ed equation analysis done in [15] on the Lax-Friedrichs scheme for ! = 14 . We then identify the real source of the numerical spike produced by the Lax-Friedrichs scheme and show how to eliminate it. In section 3 we perform the same analysis to a rst order scheme (derived as the zero relaxation limit of the rst order relaxation scheme or by a ux splitting procedure) and show that it gives the same viscosity-capillarity solution as the Lax-Friedrichs scheme. We nd that with higher order time integrator the CFL restriction can be greatly relaxed. In section 4 we implement a second order scheme (derived as the zero relaxation limit of the second order relaxation scheme or by a
ux splitting procedure) to the mixed type problem. Numerical experiments show that, for problems not sensitive to !, it picks up the same solution as the rst order relaxation scheme or the Lax-Friedrichs scheme but with sharper resolutions in shocks and phase boundaries. Finally in section 5 we show how the ( rst order) scheme can be used to select the viscosity-capillarity solution for all ! 2 (0; 14 ]. In our numerical discretizations, xj + 21 is the grid point, xj is the cell center, x = xj + 21 ? xj ? 12 is the cell width, tn = nt (n 0), the time step t = tn ? tn?1. Here we assume a uniform grid without loss of generality. As usual, for a quantity Q, the cell-edge value is Qj + 12 = Q(xj + 21 ), and the cell average value is 1 Z xj+ 21 Q(x) dx : (1.7) Qj = x xj? 1 We also introduce the CFL number (1.8)
2
t : = x
2. The Lax-Friedrichs Scheme. 2.1. The Modi ed Equation Analysis. We begin with a brief review of the
modi ed analysis done by Slemrod [15] on the Lax-Friedrichs scheme. Consider the Lax-Friedrichs discretization of (1.1): Wjn+1 ? Wjn = 21 (Vjn+1 ? Vjn?1) + 12 (Wjn+1 ? 2Wjn + Wjn?1 ) ; (2:1a) n n n n n n +1 n 1 1 Vj ? Vj = ? 2 (p(Wj +1 ) ? p(Wj ?1)) + 2 (Vj +1 ? 2Vj + Vj ?1) :(2:1b) If (2.1) has a smooth solution W and V , one can obtain the following modi ed equation
numerical solutions of mixed type systems
5
of (2.1) by the Taylor expansion: Wt ? Vx = 21 ?1 x Wxx + 21 t [p0(W)Wx ]x + O(t2) = 21 t [ (; W)Vx]x + O(t2) ; (2.2) Vt + p(W)x = 21 ?1 x Vxx + 21 t [p0(W)Vx ]x + O(t2 ) = 21 t [ (; W)Wx]x + O(t2) ; where (2.3) (; W) = ?2 + p0(W) ; and ! constant as x; t ! 0. This shows that (2.1) is rst order to (1.1) and second order to Wt ? Vx = 21 [ (; W)Wx ]x ; (2.4) Vt + p(W)x = 21 [ (; W)Vx ]x : It is desirable that (; W) > 0. When p0 > 0 (the elliptic region) this is true automatically. If p0 < 0 (the hyperbolic region) one requires (2.5) max f?p0(W) ; W such that p0 (W) < 0g 2 < 1 W which is just the CFL condition. Rewrite (2.4) as Wt ? Vx = 21 ?2 t Wxx + 12 t [p0(W)Wx ]x ; (2.6) Vt + p(W)x = 21 ?2 t Vxx + 12 t [p0(W)Vx ]x : If one chooses (2.7) ?2 t = t for 0 < < 21 or (2.8) = t 1?2 ; then the rst terms in the right hand side of (2.6) dominate and the Lax-Friedrichs scheme becomes a rst order approximation to (1.5) with = t . Since (1.5) is equivalent to (1.3) upon the change of variable (1.4), this analysis shows that the LaxFriedrichs scheme (2.1) indeed picks up the viscosity-capillarity solution with ! = 41 . Applying the Lax-Friedrichs scheme (2.1) to the Riemann problem W(x; 0) = WL = 0:7; U(x; 0) = UL = ?0:327; for x < 0 ; W(x; 0) = WR = ?0:7; U(x; 0) = UR = ?0:327; for x > 0 ; where p(w) has the van der Waals like form (2.9) p(w) = w ? w3 : p p The elliptic region is w 2 (? 33 ; 33 ). The exact solution to this problem, with the viscosity-capillarity criterion, contains a left and right moving shocks, with a steady phase boundary in between at x = 0 [15], as shown by the solid lines in Figure 1. 1 ; = 1 and the solutions Numerically we take ?5 < x < 5; x = 0:02; t = 350 7 at t = 2 are depicted by the circles in Figure 1. In the W coordinate the results seem nicely to depict the exact solution. The V coordinate also seems to provide a reasonable representation of the exact solution except at x = 0 where a spike is formed.
6
shi jin
This is exactly the numerical experiment carried out in [15]. The appearance of the spike was interpreted in [15] from the approximation in the neighborhood of x = 0 given by (2.6a) (2.10) Vx ? 21 ?2 t Wxx : As W approaches a step, V approaches a spike. 2.2. The Real Reason for the Spike. By carefully examining Figure 1, we nd that as W approaches a step, V approaches a spike which looks like a Dirac delta function. At the same time, Wx also approaches a delta like function with a negative sign. If one adds V with a suitably scaled Wx one might get a nice cancellation of the spike. This observation motivates our identi cation of the source of the spike. Notice that the Lax-Friedrichs scheme is only an approximation to (1.5), which is not yet the viscosity-capillarity equations (1.3) unless the change of variable (1.4) is carried out. It is crucial to notice that (1.5) clearly bears the feature to cancel the delta like function in V with the negative delta like function in Wx . Thus we identify the source of the spike as the inconsistency of the Lax-Friedrichs scheme to the viscositycapillarity equations (1.3). Such an inconsistency may be eliminated by the numerical change of variable (2.11) U = V + 12 Wx : 2.3. The Elimination of the Spike. To eliminate the numerical spike one has to perform the change of variable (2.11) at the end of the integration and depict U as the velocity. Notice that (2.11) involves the spatial derivative Wx which should be discretized properly. For example, our numerical experiment shows that a naive centered di erence will not give a satisfactory result. Rewriting the right hand side of (2.1a) as (Vj + 21 ? Vj ? 12 ) where (2.12) Vj + 21 = 21 (Vj +1 + Vj ) + 12 ?1(Wj +1 ? Wj ) : This is exactly a second order discretization to the right hand side of (2.11) at x = xj + 12 if = ?1 x. A discretization to (2.11) should be consistent to (2.12) [3]. Thus we take (2.13a) Uj + 12 = 21 (Vj +1 + Vj ) + 21 Wj +1x? Wj ; with (2.13b) = ?1x = ?2 t : Our Lax-Friedrichs algorithm then consists of the following steps: by
1) Given initial data W(0; x) and U(0; x) nd the average quantity Wj0 and Vj0
(2.14)
1 Z xj+ 21 W(0; x) dx ; Wj0 = x xj? 12
1 Z Vj0 = x
xj+ 12 xj? 21
U(0; x) dx :
2) For n = 1; 2; ; N where Nt = T, applying the Lax-Friedrichs scheme (2.1) to obtain Wjn and Vjn .
7
numerical solutions of mixed type systems
0.8
1
0.6
0.5
0.4
0
0.2
W
V
1.5
-0.5
0
-1
-0.2
-1.5 -5
-4
-3
-2
-1
0 x
1
2
3
4
-0.4 -5
5
-4
-3
-2
-1
0 x
1
2
3
4
5
Fig. 1. The Lax-Friedrichs solutions ('o') of (2.9) at t = 2 with 500 cells and the CFL number = 0:202. The solid line is the exact solution. The numerical solution of V has a spike near the phase boundary x = 0.
0.05
0
-0.05
-0.1
U
-0.15
-0.2
-0.25
-0.3
-0.35 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
Fig. 2. Same as Figure 1 with an extra step (2.13) at the output time. The U solutions
are plotted. The numerical spike has been eliminated (note the scale change).
3) Use (2.13) to nd UjN+ 21 and output WjN and UjN+ 21 . In Figure 2 we use the same data as in Figure 1 with the extra step 3) in the new algorithm given above. In the U-coordinate the spike is completely eliminated and the numerical results are quite satisfactory.
3. An Improved First Order Scheme. In an e ort to develop a higher order scheme for (1.1), we begin with the relaxation system introduced by the author and Xin [5]. Letting U = (w; u)T and F(U) = (?u; p(w))T , and introducing a new variable V 2 R2, we couple U and V by the following second order hyperbolic system (3.1)
@t U + @ x V = 0 ; @t V + A@x U = ? 1 (V ? F(U)) :
8
shi jin
In (3.1), is a small positive parameter called the relaxation time, (3.2) A = diag(a21 ; a22) for a1 > 0; a2 > 0 is a positive diagonal constant matrix satisfying Liu's subcharacteristic condition (3.3) A F0(U)2 where F 0(U) is the Jacobian matrix of F. Note that (3.1) is always hyperbolic and well-posed for any > 0, even if the original system contains elliptic region. Standard numerical methods for hyperbolic systems may be applied to the linear convection part, coupled with some sti ODE solver for the source term. Numerical discretizations of (3.1) are called the relaxation schemes in [5], which does not involve any Riemann solver and is a good approximation to the original system when is small. In this article will not play any role, since we will only use the zero-relaxation limit of the relaxation schemes, but A is an important parameter in the method applied here. The zero-relaxation limit of the rst order relaxation scheme takes the form [5] Wjn+1 ? Wjn = 21 (Vjn+1 ? Vjn?1) + 21 a1 (Wjn+1 ? 2Wjn + Wjn?1) ; (3.4) Vjn+1 ? Vjn = ? 12 (p(Wjn+1 ) ? p(Wjn?1)) + 21 a2 (Vjn+1 ? 2Vjn + Vjn?1) : A simple choice of the matrix A is [5] (3.5)
p
a1 = a2 = a = max f ?p0 (W); W such that p0 (W) < 0g ; W
i.e., they are chosen to be the maximal (real) eigenvalue of (1.1). Remark: Note that (3.4) can also be derived from a ux-splitting (3.6a) F = F + + F? where p p (3.6b) F+ = 12 (F(U) + AU) ; F? = 21 (F(U) ? AU) ; by applying upwind schemes to F+ and F? respectively. This ux splitting bears some similarity with the one used in [11]. We apply the modi ed equation analysis to scheme (3.4). With the a given in (3.5), the Taylor expansion on (3.4) gives Wt ? Vx = 12 ax Wxx + 12 t [p0(W)Wx ]x + O(t2) = 12 t [ (a; ; W)Wx]x + O(t2) ; (3.7) Vt + p(W)x = 21 ax Vxx + 12 t [p0(W)Vx ]x + O(t2) = 12 t [ (a; ; W)Vx]x + O(t2) ; where (3.8) (a; ; W) = a?1 + p0(W) : Note (?1 ; ; W) = (; W). Again if one takes (3.9) = at1? for 0 < < 21 ;
numerical solutions of mixed type systems
9
then scheme (3.4) approximates
Wt ? Vx = 21 t Wxx ; Vt + p(W)x = 21 t Vxx ; with a rst order accuracy, similar to the Lax-Friedrichs scheme. Thus one should expect that the scheme (3.4) selects the same solutions as the Lax-Friedrichs scheme. (3.10)
In (3.7), 21 axWxx and 12 axVxx in the error terms re ect the spatially rst order nature of scheme (3.4), while 12 t[p0(W)Vx ]x and 12 t[p0(W)Vx ]x re ect the temporally rst order nature of scheme (3.4). In order to obtain a numerical dissipation with constant coecient (to be consistent with (1.5)), it makes sense to use a higher order temporal integrator so the nonlinear part of the numerical dissipations does not appear. We use the second order TVD Runge-Kutta method [12] for this purpose. Denote the numerical ux of (3.4) by Vj + 21 = 12 (Vj + Vj +1 ) + 12 a(Wj +1 ? Wj ) ; (3.11) pj + 21 = 12 (p(Wj ) + p(Wj +1 )) ? 21 a(Vj +1 ? Vj ) ; then the second order Runge-Kutta method gives Wj(1) = Wjn + (Vjn+ 12 ? Vjn? 21 ) ; Vj(1) = Vjn ? (pnj+ 12 ? pnj? 21 ) ; (1) (2) (1) (1) (1) (3.12) Wj(2) = Wj(1) + (Vj(1) + 21 ? Vj ? 12 ) ; Vj = Vj ? (pj + 21 ? pj ? 12 ) ; Wjn+1 = 12 (Wjn + Wj(2)) ; Vjn+1 = 21 (Vjn + Vj(2)) : The Taylor expansion on (3.11), (3.12) gives Wt ? Vx = 21 ax Wxx + O(t2) ; (3.13) Vt + p(W)x = 21 ax Vxx + O(t2) : This is consistent to (1.5) when ax = a?1 t = t for 0 < < 1 or (3.14) = at1? : This is a much more relaxed CFL restriction compared to (3.9). From now on we will always use the second order temporal discretization (3.12). Remark. The CFL restrictions (2.7), (3.9) and (3.14) are for the purpose of asymptotic analysis. In our numerical experiments carried out below we do not nd the necessity for such restrictions. It seems that any CFL number that guarantees the numerical stability will select the same solution, although the solutions with smaller CFL number are more smeared out. Now our algorithm for the rst order relaxation scheme consists of the following three steps: 1) Initialization by (2.14). 2) For n = 1; 2; ; N where Nt = T, apply (3.11) and (3.12). 3) Find UjN+ 21 by
(3.15)
UjN+ 12 = 21 (VjN + VjN+1 ) + 12 a(WjN+1 ? WjN ) :
shi jin
1.5
0.8
1
0.6
0.5
0.4
0
0.2
V
W
10
-0.5
0
-1
-0.2
-1.5 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
-4
-3
-2
-1
-0.4 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
0.05
0
-0.05
U
-0.1
-0.15
-0.2
-0.25
-0.3
-0.35 -5
0 x
1
2
3
4
5
Fig. 3.
The numerical solution ('o') of (2.9) by the rst order scheme (3.11) and (3.12) at t = 2 with 500 cells and the CFL number = 1. The solid lines are the exact solution. Note the solution U has no spike.
We now test this scheme p by the example p spatial p of section 2 with the same grid, and choose a = maxW ?p0 (W) = 2, the CFL number = a = 2 xt = 1. Although V coordinate has a spike near x = 0, both W coordinate and U-coordinate, as depicted in Figure 3, give satisfactory results. The resolutions here are sharper than Figure 1 and 2 because here the CFL number = 1 is much larger. 4. A Second Order Scheme. We now implement a second order scheme, which is the zero-relaxation limit of the second order relaxation scheme derived in [5]. Instead of the piecewise constant approximation in the rst order scheme, the second order scheme uses van Leer's MUSCL type piecewise linear approximation [18]. For temporal discretization we still use the second order Runge-Kutta method (3.12). Spatially the numerical uxes are given by [5] Vj + 21 = 21 (Vj + Vj +1) + 21 a (Wj +1 ? Wj ) + 41 x ( j+ + j?+1 ) ; (4.1) pj + 12 = 21 (p(Wj ) + p(Wj +1)) ? 12 a (Vj +1 ? Vj ) ? 14 x ( j+ + j?+1 ) ;
where and are the slopes of aW V and aV p(W) respectively. For a general
numerical solutions of mixed type systems
11
quantity Q, de ne (4.2)
jL = Qj ?xQj ?1 ; jR = Qj +1x? Qj ; jO = 12 (jL + jR ) ;
then the limited slope of Q is given by 8 > if jL jR 0 ; O jL jR ;2 otherwise : min 1; 2 j : jO jO
Another slope limiter, called the minmod limiter, is de ned by ! !! jL jL R (4.4) j = j max 0; min 1; 2 R ; min R ; 2 : j j Note that if there is a discontinuity both limited slopes (4.3) and (4.4) become zero and the second order scheme reduces to the rst order scheme. Thus one might expect that the second order scheme (4.1) will select the same shocks and phase boundaries as the rst order scheme or the Lax-Friedrichs scheme, i.e., the viscosity-capillarity solutions with ! = 41 . This is the case when the solution of the problem is not sensitive to !, as demonstrated by the Slemrod-Flaherty test problem to be given below. However, if the solution is indeed a function of !, then the test performed in the next section indicates that the second order scheme may not converge to the correct solution. Note that all these arguments are not theoretically justi ed, thus are solely our numerical observations. Remark: One can also derive (4.1) by applying the MUSCL to the split uxes F+ and F? in (3.6b) respectively. Here the slope limiters are applied to the primitive variables thus no local eldby- eld decomposition is needed. Actually such a decomposition is not available at all in the elliptic region. In summary, our second order method consists of the initialization (2.14), the numerical uxed de ned by (4.1) combined with the second order Runge-Kutta method (3.12) for n = 1; 2; ; N. Finally we set (4.5) UjN+ 21 = 21 (VjN + VjN+1 ) + 12 a (WjN+1 ? WjN ) + 41 x ( j+;N + j?+1;N ) :
We test this scheme (with slope limiter (4.3)) for the problem given in section 2 with the same data as the rst order scheme (3.11)-(3.12) except now the CFL number = 21 . The solutions are depicted in Figure 4. The W-coordinate gives a sharp, non-oscillatory solution. The V -coordinate gives a sharp solution which is slightly oscillatory near the phase boundary with a spike at the phase boundary x = 0, and the U-coordinate is similar to V without the spike. We have also performed re ned computation and the slightly oscillatory behavior of V and U eventually die out as one re nes the grids. If one uses the minmod slope limiter (4.4) then all the solutions are monotone (except the spike in V ). We also test all these schemes for the van der Waals equation of state 0:9 ; (4.6) p(W) = W ?10:25 ? W 2 where W 2 (0:494279; 1:405065) is the elliptic region. We solve the following Riemann
12
shi jin
1.5
0.3
0.2
1
0.1 0.5
V
W
0 0
-0.1 -0.5 -0.2
-1
-1.5 -5
-0.3
-4
-3
-2
-1
0 x
1
2
3
4
5
-4
-3
-2
-1
-0.4 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
0.05
0
-0.05
U
-0.1
-0.15
-0.2
-0.25
-0.3
-0.35 -5
0 x
1
2
3
4
5
Fig. 4. The numerical solutions ('o') by the second order scheme (4.1) at t = 2 with slope limiter (4.3), 500 cells and the CFL number = 0:5. The solid lines are the exact solution.
problem:
(VL ; WL) = (1; 0:54); (VR ; WR ) = (1; 1:8517) : This example was used by Shu in [11] to test his ux-split ENO method, where he obtained an oscillatory behavior in W near the phase boundary as well as an undershoot. In our test we use the same grid as in [11], choosing a to be the maximal real eigenvalue of the system (a = 1:264911) and the CFL number = 0:5. The results by the rst order scheme are monotone but very dissipative (we omit the gures). The results by the second order scheme with the slope limiter (4.3) are depicted in Figure 5, where the 'o' lines stand for the numerical results with 200 cells, and the solid lines are the results by 2000 cells of the same method (served as the converged exact solution). The W-coordinate shows a sharp, non-oscillatory solution, while the V and U coordinates are sharp but oscillatory near the phase boundary. The re ned computation (the solid line) shows similar oscillatory behavior (with smaller magnitude) in V and U as the coarse one. The oscillatory behavior of V and U remains pretty much the same even if the minmod limiter (4.4) is used. 5. The Numerical Schemes for ! 2 (0; 41 ]. So far we have restricted our
13
2
1.08
1.8
1.06
1.6
1.04
1.4
1.02
1.2
1
V
W
numerical solutions of mixed type systems
0.98
1
0.96
0.8
0.94
0.6
0.92
0.4 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
-4
-3
-2
-1
-5
-4
-3
-2
-1
0 x
1
2
3
4
5
1.08
1.06
U
1.04
1.02
1
0.98
0.96 -5
0 x
1
2
3
4
5
Fig. 5. Numerical solutions of (4.6) by the second order scheme (4.1) at t = 4 with slope limiter (4.3) and the CFL number = 0:5. The 'o' lines are the solutions with 200 cells, the solid lines are the solutions with 2000 cells.
study to the case ! = 41 . For certain mixed type problems di erent ! corresponds to di erent materials, it is often necessary to compute the solution of the mixed type problem for ! other than 41 . This issue was addressed in [3], where Cockburn and Gau constructed a system equivalent to (1.3) with enough free parameters to capture all ! > 0. Based on a direct discretization of the original system (1.1), our strategy, however, is to make the scheme itself capture certain range of ! automatically. Thanks the exibility of the relaxation schemes this goal can be achieved for ! 2 (0; 41 ] rather straightforwardly. In the original derivation of the relaxation schemes in [5], we allow a1 and a2 to be di erent in (3.1), as long as p 0 (5.1) a1 ; a2 a = max f ?p (W) ; W such that p0(W) < 0g : W If the second order Runge-Kutta method (3.12) is used then, corresponding to (3.13), we now have Wt ? Vx = 21 a1 xWxx ; (5.2) Vt + p(W)x = 12 a2 xVxx ;
14
shi jin
after ignoring the O(t2) terms. De ne (5.3) U = V + 21 a1xWx ; then one can rewrite (5.2) as Wt ? Ux = 0 ; (5.4) Ut + p(W)x = 12 (a1 + a2)xUxx ? 14 a1 a2x2Wxxx : Let = 21 (a1 + a2)x and !2 = 41 a1 a2x2, then (5.5) ! = (a a+1aa2 )2 ; 1 2 and (5.4) recovers (1.3). Clearly for a1 ; a2 > 0 the range of ! now is (0; 41 ]. If one chooses p 0 (5.6) a1 = a = max f ?p (W) ; W such that p0(W) < 0g ; W and let (5.7) then (5.1) is satis ed and (5.8)
Solving (5.8) for r 1 gives (5.9)
a2 = a 1 r
for r 1 ;
! = (1 +r r)2 : 1 [ 1 ? 2! + p1 ? 4! ] : r = 2!
Therefore we have the following algorithm: 1) Given any ! 2 (0; 14 ], nd r 1 from (5.9), and then a1 and a2 from (5.6) and (5.7) respectively. 2) Initialization by (2.14). 3) Use the scheme (a combination of (3.4) and (4.1)) Vj + 12 = 12 (Vj + Vj +1 ) + 12 a1 (Wj +1 ? Wj ) + 14 x ( j+ + j?+1 ) ; (5.10) pj + 21 = 12 (p(Wj ) + p(Wj +1 )) ? 12 a2 (Vj +1 ? Vj ) ? 41 x ( j+ + j?+1 ) ; spatially and the Runge-Kutta method (5.8) temporally for n = 1; 2; ; N. Here ( 0 for the rst order relaxation scheme ; (5.11) = 1 for the second order relaxation scheme : and are the limited slopes of a1 W V and a2V p(W) respectively with the slope limiters de ned by (4.3) or (4.4). 4) Find UjN+ 12 by (5.12)
UjN+ 21 = 12 (VjN + VjN+1) + 21 a1(WjN+1 ? WjN ) + 14 x ( j+;N + j?+1;N ) :
By the above analysis, this algorithm will guarantee that the rst order scheme produces the correct viscosity-capillarity solution for any ! 2 (0; 41 ]. However, we
numerical solutions of mixed type systems
15
do not observe the same convergence for the second order scheme in the numerical example given below. We test these schemes for a piecewise-linear equation of state 8 ? 20w if w < 0:1; > < (5.13) p(w) = > ? 3 + 10w if 0:1 w < 0:2; : ? 5w if w 0:2: The elliptic region is the domain w 2 (0:1; 0:2). The Riemann problem for the above equation of state (stress-strain relation) was solved by Abeyaratne and Knowles [1]. In our numerical experiments we take the initial data (5.14) (VL ; WL ) = (0; 0:05); (VR ; WR ) = (0:5; 0:07) : This problem is challenging because the solution not only consists of both shocks and phase boundaries but also depends sensitively on !. One needs far more grid points to obtain a reasonable resolution than the ones we studied in previous sections. In Figure 6 we depict the numerical results of the rst order and the second order (with slope limiter (4.3)) schemes for ! = 1=4 and ! = 0:0625 respectively at t = 0:05. In each run p we take 2000 spatial cells between [?0:5; 0:5], a1 = 20, the CFL number is 1 for the rst order scheme and 1=2 for the second order scheme. For the rst order scheme we observe the numerical convergence, and the solution is more dissipative with smaller !. The second order scheme, with a higher resolution though, fails to converge near the phase boundaries. We also used the minmod slope limiter but the results essentially remain the same. We think that the slope limiter should be modi ed so it is more sensitive to the discontinuities, or in (3.1) should be turned on to give more viscosity and capillarity. This remains unsolved and is a subject of future research. Acknowledgement. I am grateful to Bernardo Cockburn and Huiing Gau for showing me their recent work [3] and giving me the exact solution code for the test problem in section 5. In addition, I would like to thank John Lowengrub, Marshall Slemrod and Zhouping Xin for some helpful discussions on mixed type problems.
[1] [2] [3] [4] [5] [6] [7] [8]
REFERENCES R. Abeyaratne and J.K. Knowles, Kinetic relations and the propagation of phase boundaries in solids, Arch. Rat'l Mech. Anal., 114 (1991), pp. 119{154. M. Affouf and R. Caflisch, A numerical study of Riemann problem solutions and stability for a system of viscous conservation laws of mixed type, SIAM J. Appl. Math., 51 (1991), pp. 605-634. B. Cockburn and H. Gau, Model numerical schemes for phase transitions in solids, SIAM J. Sci. Comp., submitted. R. James, Dynamic phase changes with non-monotone stresses, Arch. Rat'l Mech. Anal., 66 (1983), pp. 59{125. S. Jin and Z.P. Xin, The relaxation schemes for systems of hyperbolic conservation laws in arbitrary space dimensions,, Comm. Pure Appl. Math., to appear. P.D. Lax, Hyperbolic systems of conservation laws and the mathematical theory of shock waves, 11 (1973). B.L. Keyfitz and M. Shearer, eds., Nonlinear Evolution Equations that Change Type, The IMA Volumes in Mathematics and its Applications, Springer-Verlag, 1990. T.-P. Liu, Hyperbolic conservation laws with relaxation, Comm. Math. Phys., 108 (1987), pp. 153{175.
16
shi jin
0.4
0.35
0.3
0.3
0.25
0.25
W
W
0.35
0.4 (a)
0.2
0.15
0.1
0.1
0.05
0.05
-0.4
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0 -0.5
0.5
0.4
0.35
(c)
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0.5
0.3
0.25
0.25
0.2
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0.5
0.2
0.15
0.15
0.1
0.1
0.05
0.05
-0.4
(d)
0.35
0.3
0 -0.5
-0.4
0.4
W
W
0.2
0.15
0 -0.5
(b)
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0.5
0 -0.5
-0.4
Fig. 6. Numerical solutions of w for the piecewise-linear equation of state (5.13) with initial data (5.14). In each run 2000 spatial cells are used. At t = 0:05, the exact solution is plotted by the solid lines, and the numerical solutions by the dashed lines. In (a) (! = 1=4) and (b) (! = 0:0625) the rst order scheme is applied with the CFL number = 1. In (c) (! = 1=4) and (d) (! = 0:0625) the second order scheme with slope limiter (4.3) is applied with the CFL number = 1=2.
[9] E.B. Pitman and Y. Ni, Visco-elastic relaxation with a van der Waals type stress, Int. J. Eng. Sci., 32 (1994), pp. 327{338. [10] M. Shearer, Non-uniqueness of admissible solutions of conservation laws of mixed type, Arch. Rat'l Mech. Anal., 93 (1986), pp. 45{59. [11] C.-W. Shu, A numerical method for systems of conservation laws of mixed type admitting hyperbolic ux splitting, J. Comp. Phys., 100 (1992), pp. 424{429. [12] C.-W. Shu and S. Osher, Ecient implementation of essentially non-oscillatory shock capturing methods, J. Comp. Phys., 77 (1988), pp. 439{471. [13] M. Slemrod, Admissibility criteria for propagating phase boundaries in a van der Waals uid, Arch. Rat'l Mech. Anal., 81 (1983), pp. 301{315. [14] M. Slemrod, Lax-Friedrichs and the viscosity-capillarity criterion, in Physical Mathematics and Nonlinear Partial Di erential Equations, S.M. Rankin and J.H. Lightbourne, eds., 1985. [15] M. Slemrod and J.E. Flaherty, Numerical integration of a Riemann problem for a van der Waals uid, in Phase Transformations, C.A. Elias and G. John, eds., 1986. [16] A. Sommerfeld, Thermodynamics and Statistical Mechanics, Academic Press, New York and
numerical solutions of mixed type systems
17
London, 1964. [17] L. Truskinovski, Equilibrium phase interface, Dokl. Akad. Nauk. SSSH, 256 (1982), pp. 306{ 310. [18] B. van Leer, Toward the ultimate conservative di erence schemes V. A second-order sequal to Godunov's method, J. Comp. Phys., 32 (1979), pp. 101{136.
SIAM J. APPL. MATH. Vol. ??, No. ??, pp. 1{16, ?? ??
??
NUMERICAL INTEGRATIONS OF SYSTEMS OF CONSERVATION LAWS OF MIXED TYPE SHI JINy
Abstract. The systems of conservation laws have been used to model dynamical phase transitions in, for example, the propagating phase boundaries in solids and the van der Waals uid. When integrating such mixed hyperbolic-elliptic systems the Lax-Friedrichs scheme is known to give the correct solutions selected by a viscosity-capillarity criterion except a spike at the phase boundary which does not go away even with a re ned mesh [15]. We identify the source of this spike as an inconsistency between the Lax-Friedrichs discretization and the viscosity-capillarity equations, and show a simple change of variable that can eliminate this spike. We then implement a high resolution scheme for the mixed type problems that select the same viscosity-capillarity solutions as the Lax-Friedrichs scheme with higher resolutions. Furthermore, a exibility in the ( rst order) scheme is used to obtain solutions for a wide range of the viscosity-capillarity equations. Key words. Systems of conservation laws of mixed type, viscosity-capillarity admissibility criterion, the Lax-Friedrichs scheme, shock capturing schemes. AMS(MOS) subject classi cations. 35L65, 65M06 1. Introduction. Systems of conservation laws of mixed type have been used to describe certain dynamic phase transitions. The best known mixed type system is the van der Waals gas at suciently low temperature [16]. Mixed type system is also used for the propagating phase boundaries in solids, where the stress-strain relation is not monotone [4]. In this paper we are concerned with the equations of isothermal motion for a compressible uid in Lagrangian coordinates wt ? u x = 0 ; (1.1) ut + p(w)x = 0 ; where u is the velocity, w is the speci c volume, and p is the pressure. Here the equation of state is a non-monotone function: p0 (w) < 0 if w < m or w > M ; the hyperbolic region; (1.2) p0 (w) > 0 if m < w < M ; the elliptic region:
p
Notice that the eigenvalues of (1.1) are ?p0(w), which become complex in the elliptic region and result in the classical Hardmard instability. We say that the solution is in phase 1, 2 or 3 at time t if p(w) lies in the interval (0; m); (m; M) or (M; 1) respectively. In the hyperbolic region it is expected that shocks may form [6]. A phase boundary is a shock that connects the solutions in di erent phases. When selecting the admissible weak solution to Eqs.(1.1), it is well known that the viscosity criterion is too restrictive to be used since it rules out propagating phase boundaries near the equilibrium co-existence line [13]. Instead, Slemrod [13] and Truskinovsky [17] independently proposed to use the viscosity-capillarity criterion. y School of Mathematics and Center for Dynamical System and Nonlinear Studies, Georgia Institute of Technology, Atlanta, GA 30332, USA ([email protected]). This work was supported by NSF grant no DMS-9404157. 1
2
shi jin
The viscosity-capillarity regularization for (1.1) is: wt ? ux = 0 ; (1.3) ut + p(w)x = uxx ? !2 wxxx : Here > 0 denotes the viscosity and !2 denotes the coecient of interfacial capillarity. Both and ! are positive constants. A ouf and Ca isch conducted numerical studies on the Riemann problem solutions and stability for (1.3) [2]. More recently Cockburn and Gau developed numerical schemes for (1.3) that work for all ! > 0 [3]. In addition, Abeyaratne and Knowles have proposed to use a kinetic relation to select the unique solution [1]. We are concerned with the numerical integrations of the mixed type problem (1.1) and (1.2). In this pursuit of particular interest is the work of Slemrod [14], who showed that the Lax-Friedrichs scheme for (1.1) gives the solution of (1.3) with ! = 41 . To see this, letting ! = 41 and introducing the change of variable (1.4) v = u ? 21 wx ; then Eqs.(1.3) can be rewritten as wt ? vx = 12 wxx ; (1.5) vt + p(w)x = 21 vxx : The Lax-Friedrichs discretization on (1.1), as shown in [14], is a second order discretization of (1.5) for some under an appropriate Courant-Friedrichs-Lewy (CFL) condition. Thus the Lax-Friedrichs scheme gives the weak solution of (1.1) selected by the viscosity-capillarity criterion (1.3) with ! = 41 . This analysis was later con rmed by the numerical experiments done by Slemrod and Flaherty [15]. Unfortunately the Lax-Friedrichs scheme always produces a spurious spike near the phase boundary, the nature of which is still unclear. To the author's knowledge no method based on a direct discretization of system (1.1) has been developed that can eliminate the spike. In recent years considerable interests have been shown to study the mixed type problem numerically. Despite its overdissipation and the unpleasant spike, the LaxFriedrichs scheme is still the most frequently used method for mixed type problems [7]. Attempts have been made to design modern high resolution shock capturing methods such as the total-variation-diminishing (TVD) methods and the total-variationbounded (TVB) methods. The direct implementation of these methods to mixed type problems faces many diculties, since these methods are usually based on the Riemann solver and local eld-by- eld characteristic decompositions. It has been shown by Shearer that the Riemann problem of (1.1) does not have a unique solution [10]. Also the local eld-by- eld decomposition is not available in the elliptic region. Shu found a hyperbolic ux splitting, and implemented the essentially-nonoscillatory methods (ENO) to each of the split uxes respectively. Local eld-by- eld decomposition is applied in the hyperbolic region while in the elliptic region the reconstruction is on the primitive variables. Pitman and Ni applied the high order Godunov methods to a relaxation regularization of the mixed type equations and the results are nonoscillatory [9]. Their method may select solutions di erent from those selected by the viscosity-capillarity solution and their method has to resolve the small relaxation time. In this paper we intend to answer some basic questions regarding to numerical discretization of the mixed type system (1.1). Speci cally, we would like to address
numerical solutions of mixed type systems
3
the following issues: 1) Identify the real source of the numerical spike produced by the Lax-Friedrichs scheme and show how to eliminate it. 2) Implement a second order shock capturing scheme for mixed type system (1.1) that pick up the solution selected by the viscosity-capillarity criterion with higher accuracy and sharper shocks and phase boundaries. 3) Attempt to extend these methods to handle a wide range of ! 2 (0; 14 ], which expands the applicability of our methods for mixed type problems. In an e ort to understand the relation between the solution of (1.1) and the capillarity coecient !, Cockburn and Gau have conducted a numerical convergence study for the viscosity-capillarity equations (1.3) as ! 0 with di erent ! [3]. They constructed a new system equivalent to (1.3) under (1.4) with enough free parameters so as to work for all ! > 0. The numerical discretization is based on this new system, and the change of variable (1.4) is performed to convert to the solutions of (1.3). Interesting enough is that their results do not have the spike, and their method works e ectively for all ! > 0. Their method, however, is based on a discretization for a viscosity-capillarity equations, thus is at most rst order to (1.1). Our work is partly motivated by their study, but we emphasize more on the practical side. We aim at developing modern, robust and high order numerical schemes based on a direct discretization of the original system (1.1) (rather than the viscosity-capillarity equations) that automatically pick up the correct viscosity-capillarity solutions. First we point out that the numerical spike produced by the Lax-Friedrichs scheme is caused by an inconsistency of the scheme with the viscosity-capillarity equations (1.3), in the sense that the transformation (1.4) that leads from (1.3) to (1.5), or vice versa, is not inherently possessed in the Lax-Friedrichs discretization. A numerical change of variable (1.6) u = v + 21 wx in the last step of the integration should be applied to obtain satisfactory results. Such a change of variable involves a discretization of the right hand side of (1.6) which must be consistent to the discretization of (1.1). If this is performed suitably then the numerical solution, while maintaining other features of the Lax-Friedrichs scheme, does not have the spike. In our e ort to implement a second order shock capturing method for (1.1), we use the relaxation schemes developed by the author and Xin [5]. The relaxation schemes are natural candidates for such problems since they do not need the Riemann solver, nor the local eld-by- eld decomposition. At the relaxation level, the Riemann solver and the eld-by- eld decomposition are trivially obtained due to the linear convection of the relaxation system and an operator splitting that separates the linear convection from the nonlinear sti source terms. The limiting scheme as the relaxation time approaches zero is a TVD scheme with reconstruction on primitive variables. Due to these special features the relaxation schemes have many advantages when applied to the mixed type problems. Here by the relaxation schemes we refer to the \ = 0" version of the relaxation scheme (called the relaxed schemes in [5]). Here is the relaxation time. This zerorelaxation scheme can also be derived by applying van Leer's MUSCL scheme [18] to a ux-splitting procedure, similar to Shu's work in [11]. For the rst order scheme, by
4
shi jin
performing the modi ed equation analysis, we nd that it selects the same viscositycapillarity solution as the Lax-Friedrichs scheme. We also nd that with higher order time integrator such as the second order Runge-Kutta method such a selection of the viscosity-capillarity solution can be made with a more relaxed CFL condition. Since our second order scheme reduces to the rst order scheme at discontinuities it may select the same shocks and phase boundaries as the rst order scheme (or the Lax-Friedrichs scheme), but with sharper resolutions. When forming the relaxation scheme in [5] we chose a positive diagonal constant matrix A that satis es Liu's subcharacteristic condition [8]. When passing to the zero relaxation limit A enters the numerical viscosity coecients. For problems sensitive to the parameter !, the exibility of A allows us to obtain numerical results for a wide range of ! 2 (0; 41 ] with the rst order scheme, thus expand the applicability of the scheme to a wider class of mixed type problems. However, we did not observe the same convergence for the second order scheme. We now outline this paper. In section 2 we review the modi ed equation analysis done in [15] on the Lax-Friedrichs scheme for ! = 14 . We then identify the real source of the numerical spike produced by the Lax-Friedrichs scheme and show how to eliminate it. In section 3 we perform the same analysis to a rst order scheme (derived as the zero relaxation limit of the rst order relaxation scheme or by a ux splitting procedure) and show that it gives the same viscosity-capillarity solution as the Lax-Friedrichs scheme. We nd that with higher order time integrator the CFL restriction can be greatly relaxed. In section 4 we implement a second order scheme (derived as the zero relaxation limit of the second order relaxation scheme or by a
ux splitting procedure) to the mixed type problem. Numerical experiments show that, for problems not sensitive to !, it picks up the same solution as the rst order relaxation scheme or the Lax-Friedrichs scheme but with sharper resolutions in shocks and phase boundaries. Finally in section 5 we show how the ( rst order) scheme can be used to select the viscosity-capillarity solution for all ! 2 (0; 14 ]. In our numerical discretizations, xj + 21 is the grid point, xj is the cell center, x = xj + 21 ? xj ? 12 is the cell width, tn = nt (n 0), the time step t = tn ? tn?1. Here we assume a uniform grid without loss of generality. As usual, for a quantity Q, the cell-edge value is Qj + 12 = Q(xj + 21 ), and the cell average value is 1 Z xj+ 21 Q(x) dx : (1.7) Qj = x xj? 1 We also introduce the CFL number (1.8)
2
t : = x
2. The Lax-Friedrichs Scheme. 2.1. The Modi ed Equation Analysis. We begin with a brief review of the
modi ed analysis done by Slemrod [15] on the Lax-Friedrichs scheme. Consider the Lax-Friedrichs discretization of (1.1): Wjn+1 ? Wjn = 21 (Vjn+1 ? Vjn?1) + 12 (Wjn+1 ? 2Wjn + Wjn?1 ) ; (2:1a) n n n n n n +1 n 1 1 Vj ? Vj = ? 2 (p(Wj +1 ) ? p(Wj ?1)) + 2 (Vj +1 ? 2Vj + Vj ?1) :(2:1b) If (2.1) has a smooth solution W and V , one can obtain the following modi ed equation
numerical solutions of mixed type systems
5
of (2.1) by the Taylor expansion: Wt ? Vx = 21 ?1 x Wxx + 21 t [p0(W)Wx ]x + O(t2) = 21 t [ (; W)Vx]x + O(t2) ; (2.2) Vt + p(W)x = 21 ?1 x Vxx + 21 t [p0(W)Vx ]x + O(t2 ) = 21 t [ (; W)Wx]x + O(t2) ; where (2.3) (; W) = ?2 + p0(W) ; and ! constant as x; t ! 0. This shows that (2.1) is rst order to (1.1) and second order to Wt ? Vx = 21 [ (; W)Wx ]x ; (2.4) Vt + p(W)x = 21 [ (; W)Vx ]x : It is desirable that (; W) > 0. When p0 > 0 (the elliptic region) this is true automatically. If p0 < 0 (the hyperbolic region) one requires (2.5) max f?p0(W) ; W such that p0 (W) < 0g 2 < 1 W which is just the CFL condition. Rewrite (2.4) as Wt ? Vx = 21 ?2 t Wxx + 12 t [p0(W)Wx ]x ; (2.6) Vt + p(W)x = 21 ?2 t Vxx + 12 t [p0(W)Vx ]x : If one chooses (2.7) ?2 t = t for 0 < < 21 or (2.8) = t 1?2 ; then the rst terms in the right hand side of (2.6) dominate and the Lax-Friedrichs scheme becomes a rst order approximation to (1.5) with = t . Since (1.5) is equivalent to (1.3) upon the change of variable (1.4), this analysis shows that the LaxFriedrichs scheme (2.1) indeed picks up the viscosity-capillarity solution with ! = 41 . Applying the Lax-Friedrichs scheme (2.1) to the Riemann problem W(x; 0) = WL = 0:7; U(x; 0) = UL = ?0:327; for x < 0 ; W(x; 0) = WR = ?0:7; U(x; 0) = UR = ?0:327; for x > 0 ; where p(w) has the van der Waals like form (2.9) p(w) = w ? w3 : p p The elliptic region is w 2 (? 33 ; 33 ). The exact solution to this problem, with the viscosity-capillarity criterion, contains a left and right moving shocks, with a steady phase boundary in between at x = 0 [15], as shown by the solid lines in Figure 1. 1 ; = 1 and the solutions Numerically we take ?5 < x < 5; x = 0:02; t = 350 7 at t = 2 are depicted by the circles in Figure 1. In the W coordinate the results seem nicely to depict the exact solution. The V coordinate also seems to provide a reasonable representation of the exact solution except at x = 0 where a spike is formed.
6
shi jin
This is exactly the numerical experiment carried out in [15]. The appearance of the spike was interpreted in [15] from the approximation in the neighborhood of x = 0 given by (2.6a) (2.10) Vx ? 21 ?2 t Wxx : As W approaches a step, V approaches a spike. 2.2. The Real Reason for the Spike. By carefully examining Figure 1, we nd that as W approaches a step, V approaches a spike which looks like a Dirac delta function. At the same time, Wx also approaches a delta like function with a negative sign. If one adds V with a suitably scaled Wx one might get a nice cancellation of the spike. This observation motivates our identi cation of the source of the spike. Notice that the Lax-Friedrichs scheme is only an approximation to (1.5), which is not yet the viscosity-capillarity equations (1.3) unless the change of variable (1.4) is carried out. It is crucial to notice that (1.5) clearly bears the feature to cancel the delta like function in V with the negative delta like function in Wx . Thus we identify the source of the spike as the inconsistency of the Lax-Friedrichs scheme to the viscositycapillarity equations (1.3). Such an inconsistency may be eliminated by the numerical change of variable (2.11) U = V + 12 Wx : 2.3. The Elimination of the Spike. To eliminate the numerical spike one has to perform the change of variable (2.11) at the end of the integration and depict U as the velocity. Notice that (2.11) involves the spatial derivative Wx which should be discretized properly. For example, our numerical experiment shows that a naive centered di erence will not give a satisfactory result. Rewriting the right hand side of (2.1a) as (Vj + 21 ? Vj ? 12 ) where (2.12) Vj + 21 = 21 (Vj +1 + Vj ) + 12 ?1(Wj +1 ? Wj ) : This is exactly a second order discretization to the right hand side of (2.11) at x = xj + 12 if = ?1 x. A discretization to (2.11) should be consistent to (2.12) [3]. Thus we take (2.13a) Uj + 12 = 21 (Vj +1 + Vj ) + 21 Wj +1x? Wj ; with (2.13b) = ?1x = ?2 t : Our Lax-Friedrichs algorithm then consists of the following steps: by
1) Given initial data W(0; x) and U(0; x) nd the average quantity Wj0 and Vj0
(2.14)
1 Z xj+ 21 W(0; x) dx ; Wj0 = x xj? 12
1 Z Vj0 = x
xj+ 12 xj? 21
U(0; x) dx :
2) For n = 1; 2; ; N where Nt = T, applying the Lax-Friedrichs scheme (2.1) to obtain Wjn and Vjn .
7
numerical solutions of mixed type systems
0.8
1
0.6
0.5
0.4
0
0.2
W
V
1.5
-0.5
0
-1
-0.2
-1.5 -5
-4
-3
-2
-1
0 x
1
2
3
4
-0.4 -5
5
-4
-3
-2
-1
0 x
1
2
3
4
5
Fig. 1. The Lax-Friedrichs solutions ('o') of (2.9) at t = 2 with 500 cells and the CFL number = 0:202. The solid line is the exact solution. The numerical solution of V has a spike near the phase boundary x = 0.
0.05
0
-0.05
-0.1
U
-0.15
-0.2
-0.25
-0.3
-0.35 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
Fig. 2. Same as Figure 1 with an extra step (2.13) at the output time. The U solutions
are plotted. The numerical spike has been eliminated (note the scale change).
3) Use (2.13) to nd UjN+ 21 and output WjN and UjN+ 21 . In Figure 2 we use the same data as in Figure 1 with the extra step 3) in the new algorithm given above. In the U-coordinate the spike is completely eliminated and the numerical results are quite satisfactory.
3. An Improved First Order Scheme. In an e ort to develop a higher order scheme for (1.1), we begin with the relaxation system introduced by the author and Xin [5]. Letting U = (w; u)T and F(U) = (?u; p(w))T , and introducing a new variable V 2 R2, we couple U and V by the following second order hyperbolic system (3.1)
@t U + @ x V = 0 ; @t V + A@x U = ? 1 (V ? F(U)) :
8
shi jin
In (3.1), is a small positive parameter called the relaxation time, (3.2) A = diag(a21 ; a22) for a1 > 0; a2 > 0 is a positive diagonal constant matrix satisfying Liu's subcharacteristic condition (3.3) A F0(U)2 where F 0(U) is the Jacobian matrix of F. Note that (3.1) is always hyperbolic and well-posed for any > 0, even if the original system contains elliptic region. Standard numerical methods for hyperbolic systems may be applied to the linear convection part, coupled with some sti ODE solver for the source term. Numerical discretizations of (3.1) are called the relaxation schemes in [5], which does not involve any Riemann solver and is a good approximation to the original system when is small. In this article will not play any role, since we will only use the zero-relaxation limit of the relaxation schemes, but A is an important parameter in the method applied here. The zero-relaxation limit of the rst order relaxation scheme takes the form [5] Wjn+1 ? Wjn = 21 (Vjn+1 ? Vjn?1) + 21 a1 (Wjn+1 ? 2Wjn + Wjn?1) ; (3.4) Vjn+1 ? Vjn = ? 12 (p(Wjn+1 ) ? p(Wjn?1)) + 21 a2 (Vjn+1 ? 2Vjn + Vjn?1) : A simple choice of the matrix A is [5] (3.5)
p
a1 = a2 = a = max f ?p0 (W); W such that p0 (W) < 0g ; W
i.e., they are chosen to be the maximal (real) eigenvalue of (1.1). Remark: Note that (3.4) can also be derived from a ux-splitting (3.6a) F = F + + F? where p p (3.6b) F+ = 12 (F(U) + AU) ; F? = 21 (F(U) ? AU) ; by applying upwind schemes to F+ and F? respectively. This ux splitting bears some similarity with the one used in [11]. We apply the modi ed equation analysis to scheme (3.4). With the a given in (3.5), the Taylor expansion on (3.4) gives Wt ? Vx = 12 ax Wxx + 12 t [p0(W)Wx ]x + O(t2) = 12 t [ (a; ; W)Wx]x + O(t2) ; (3.7) Vt + p(W)x = 21 ax Vxx + 12 t [p0(W)Vx ]x + O(t2) = 12 t [ (a; ; W)Vx]x + O(t2) ; where (3.8) (a; ; W) = a?1 + p0(W) : Note (?1 ; ; W) = (; W). Again if one takes (3.9) = at1? for 0 < < 21 ;
numerical solutions of mixed type systems
9
then scheme (3.4) approximates
Wt ? Vx = 21 t Wxx ; Vt + p(W)x = 21 t Vxx ; with a rst order accuracy, similar to the Lax-Friedrichs scheme. Thus one should expect that the scheme (3.4) selects the same solutions as the Lax-Friedrichs scheme. (3.10)
In (3.7), 21 axWxx and 12 axVxx in the error terms re ect the spatially rst order nature of scheme (3.4), while 12 t[p0(W)Vx ]x and 12 t[p0(W)Vx ]x re ect the temporally rst order nature of scheme (3.4). In order to obtain a numerical dissipation with constant coecient (to be consistent with (1.5)), it makes sense to use a higher order temporal integrator so the nonlinear part of the numerical dissipations does not appear. We use the second order TVD Runge-Kutta method [12] for this purpose. Denote the numerical ux of (3.4) by Vj + 21 = 12 (Vj + Vj +1 ) + 12 a(Wj +1 ? Wj ) ; (3.11) pj + 21 = 12 (p(Wj ) + p(Wj +1 )) ? 21 a(Vj +1 ? Vj ) ; then the second order Runge-Kutta method gives Wj(1) = Wjn + (Vjn+ 12 ? Vjn? 21 ) ; Vj(1) = Vjn ? (pnj+ 12 ? pnj? 21 ) ; (1) (2) (1) (1) (1) (3.12) Wj(2) = Wj(1) + (Vj(1) + 21 ? Vj ? 12 ) ; Vj = Vj ? (pj + 21 ? pj ? 12 ) ; Wjn+1 = 12 (Wjn + Wj(2)) ; Vjn+1 = 21 (Vjn + Vj(2)) : The Taylor expansion on (3.11), (3.12) gives Wt ? Vx = 21 ax Wxx + O(t2) ; (3.13) Vt + p(W)x = 21 ax Vxx + O(t2) : This is consistent to (1.5) when ax = a?1 t = t for 0 < < 1 or (3.14) = at1? : This is a much more relaxed CFL restriction compared to (3.9). From now on we will always use the second order temporal discretization (3.12). Remark. The CFL restrictions (2.7), (3.9) and (3.14) are for the purpose of asymptotic analysis. In our numerical experiments carried out below we do not nd the necessity for such restrictions. It seems that any CFL number that guarantees the numerical stability will select the same solution, although the solutions with smaller CFL number are more smeared out. Now our algorithm for the rst order relaxation scheme consists of the following three steps: 1) Initialization by (2.14). 2) For n = 1; 2; ; N where Nt = T, apply (3.11) and (3.12). 3) Find UjN+ 21 by
(3.15)
UjN+ 12 = 21 (VjN + VjN+1 ) + 12 a(WjN+1 ? WjN ) :
shi jin
1.5
0.8
1
0.6
0.5
0.4
0
0.2
V
W
10
-0.5
0
-1
-0.2
-1.5 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
-4
-3
-2
-1
-0.4 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
0.05
0
-0.05
U
-0.1
-0.15
-0.2
-0.25
-0.3
-0.35 -5
0 x
1
2
3
4
5
Fig. 3.
The numerical solution ('o') of (2.9) by the rst order scheme (3.11) and (3.12) at t = 2 with 500 cells and the CFL number = 1. The solid lines are the exact solution. Note the solution U has no spike.
We now test this scheme p by the example p spatial p of section 2 with the same grid, and choose a = maxW ?p0 (W) = 2, the CFL number = a = 2 xt = 1. Although V coordinate has a spike near x = 0, both W coordinate and U-coordinate, as depicted in Figure 3, give satisfactory results. The resolutions here are sharper than Figure 1 and 2 because here the CFL number = 1 is much larger. 4. A Second Order Scheme. We now implement a second order scheme, which is the zero-relaxation limit of the second order relaxation scheme derived in [5]. Instead of the piecewise constant approximation in the rst order scheme, the second order scheme uses van Leer's MUSCL type piecewise linear approximation [18]. For temporal discretization we still use the second order Runge-Kutta method (3.12). Spatially the numerical uxes are given by [5] Vj + 21 = 21 (Vj + Vj +1) + 21 a (Wj +1 ? Wj ) + 41 x ( j+ + j?+1 ) ; (4.1) pj + 12 = 21 (p(Wj ) + p(Wj +1)) ? 12 a (Vj +1 ? Vj ) ? 14 x ( j+ + j?+1 ) ;
where and are the slopes of aW V and aV p(W) respectively. For a general
numerical solutions of mixed type systems
11
quantity Q, de ne (4.2)
jL = Qj ?xQj ?1 ; jR = Qj +1x? Qj ; jO = 12 (jL + jR ) ;
then the limited slope of Q is given by 8 > if jL jR 0 ; O jL jR ;2 otherwise : min 1; 2 j : jO jO
Another slope limiter, called the minmod limiter, is de ned by ! !! jL jL R (4.4) j = j max 0; min 1; 2 R ; min R ; 2 : j j Note that if there is a discontinuity both limited slopes (4.3) and (4.4) become zero and the second order scheme reduces to the rst order scheme. Thus one might expect that the second order scheme (4.1) will select the same shocks and phase boundaries as the rst order scheme or the Lax-Friedrichs scheme, i.e., the viscosity-capillarity solutions with ! = 41 . This is the case when the solution of the problem is not sensitive to !, as demonstrated by the Slemrod-Flaherty test problem to be given below. However, if the solution is indeed a function of !, then the test performed in the next section indicates that the second order scheme may not converge to the correct solution. Note that all these arguments are not theoretically justi ed, thus are solely our numerical observations. Remark: One can also derive (4.1) by applying the MUSCL to the split uxes F+ and F? in (3.6b) respectively. Here the slope limiters are applied to the primitive variables thus no local eldby- eld decomposition is needed. Actually such a decomposition is not available at all in the elliptic region. In summary, our second order method consists of the initialization (2.14), the numerical uxed de ned by (4.1) combined with the second order Runge-Kutta method (3.12) for n = 1; 2; ; N. Finally we set (4.5) UjN+ 21 = 21 (VjN + VjN+1 ) + 12 a (WjN+1 ? WjN ) + 41 x ( j+;N + j?+1;N ) :
We test this scheme (with slope limiter (4.3)) for the problem given in section 2 with the same data as the rst order scheme (3.11)-(3.12) except now the CFL number = 21 . The solutions are depicted in Figure 4. The W-coordinate gives a sharp, non-oscillatory solution. The V -coordinate gives a sharp solution which is slightly oscillatory near the phase boundary with a spike at the phase boundary x = 0, and the U-coordinate is similar to V without the spike. We have also performed re ned computation and the slightly oscillatory behavior of V and U eventually die out as one re nes the grids. If one uses the minmod slope limiter (4.4) then all the solutions are monotone (except the spike in V ). We also test all these schemes for the van der Waals equation of state 0:9 ; (4.6) p(W) = W ?10:25 ? W 2 where W 2 (0:494279; 1:405065) is the elliptic region. We solve the following Riemann
12
shi jin
1.5
0.3
0.2
1
0.1 0.5
V
W
0 0
-0.1 -0.5 -0.2
-1
-1.5 -5
-0.3
-4
-3
-2
-1
0 x
1
2
3
4
5
-4
-3
-2
-1
-0.4 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
0.05
0
-0.05
U
-0.1
-0.15
-0.2
-0.25
-0.3
-0.35 -5
0 x
1
2
3
4
5
Fig. 4. The numerical solutions ('o') by the second order scheme (4.1) at t = 2 with slope limiter (4.3), 500 cells and the CFL number = 0:5. The solid lines are the exact solution.
problem:
(VL ; WL) = (1; 0:54); (VR ; WR ) = (1; 1:8517) : This example was used by Shu in [11] to test his ux-split ENO method, where he obtained an oscillatory behavior in W near the phase boundary as well as an undershoot. In our test we use the same grid as in [11], choosing a to be the maximal real eigenvalue of the system (a = 1:264911) and the CFL number = 0:5. The results by the rst order scheme are monotone but very dissipative (we omit the gures). The results by the second order scheme with the slope limiter (4.3) are depicted in Figure 5, where the 'o' lines stand for the numerical results with 200 cells, and the solid lines are the results by 2000 cells of the same method (served as the converged exact solution). The W-coordinate shows a sharp, non-oscillatory solution, while the V and U coordinates are sharp but oscillatory near the phase boundary. The re ned computation (the solid line) shows similar oscillatory behavior (with smaller magnitude) in V and U as the coarse one. The oscillatory behavior of V and U remains pretty much the same even if the minmod limiter (4.4) is used. 5. The Numerical Schemes for ! 2 (0; 41 ]. So far we have restricted our
13
2
1.08
1.8
1.06
1.6
1.04
1.4
1.02
1.2
1
V
W
numerical solutions of mixed type systems
0.98
1
0.96
0.8
0.94
0.6
0.92
0.4 -5
-4
-3
-2
-1
0 x
1
2
3
4
5
-4
-3
-2
-1
-5
-4
-3
-2
-1
0 x
1
2
3
4
5
1.08
1.06
U
1.04
1.02
1
0.98
0.96 -5
0 x
1
2
3
4
5
Fig. 5. Numerical solutions of (4.6) by the second order scheme (4.1) at t = 4 with slope limiter (4.3) and the CFL number = 0:5. The 'o' lines are the solutions with 200 cells, the solid lines are the solutions with 2000 cells.
study to the case ! = 41 . For certain mixed type problems di erent ! corresponds to di erent materials, it is often necessary to compute the solution of the mixed type problem for ! other than 41 . This issue was addressed in [3], where Cockburn and Gau constructed a system equivalent to (1.3) with enough free parameters to capture all ! > 0. Based on a direct discretization of the original system (1.1), our strategy, however, is to make the scheme itself capture certain range of ! automatically. Thanks the exibility of the relaxation schemes this goal can be achieved for ! 2 (0; 41 ] rather straightforwardly. In the original derivation of the relaxation schemes in [5], we allow a1 and a2 to be di erent in (3.1), as long as p 0 (5.1) a1 ; a2 a = max f ?p (W) ; W such that p0(W) < 0g : W If the second order Runge-Kutta method (3.12) is used then, corresponding to (3.13), we now have Wt ? Vx = 21 a1 xWxx ; (5.2) Vt + p(W)x = 12 a2 xVxx ;
14
shi jin
after ignoring the O(t2) terms. De ne (5.3) U = V + 21 a1xWx ; then one can rewrite (5.2) as Wt ? Ux = 0 ; (5.4) Ut + p(W)x = 12 (a1 + a2)xUxx ? 14 a1 a2x2Wxxx : Let = 21 (a1 + a2)x and !2 = 41 a1 a2x2, then (5.5) ! = (a a+1aa2 )2 ; 1 2 and (5.4) recovers (1.3). Clearly for a1 ; a2 > 0 the range of ! now is (0; 41 ]. If one chooses p 0 (5.6) a1 = a = max f ?p (W) ; W such that p0(W) < 0g ; W and let (5.7) then (5.1) is satis ed and (5.8)
Solving (5.8) for r 1 gives (5.9)
a2 = a 1 r
for r 1 ;
! = (1 +r r)2 : 1 [ 1 ? 2! + p1 ? 4! ] : r = 2!
Therefore we have the following algorithm: 1) Given any ! 2 (0; 14 ], nd r 1 from (5.9), and then a1 and a2 from (5.6) and (5.7) respectively. 2) Initialization by (2.14). 3) Use the scheme (a combination of (3.4) and (4.1)) Vj + 12 = 12 (Vj + Vj +1 ) + 12 a1 (Wj +1 ? Wj ) + 14 x ( j+ + j?+1 ) ; (5.10) pj + 21 = 12 (p(Wj ) + p(Wj +1 )) ? 12 a2 (Vj +1 ? Vj ) ? 41 x ( j+ + j?+1 ) ; spatially and the Runge-Kutta method (5.8) temporally for n = 1; 2; ; N. Here ( 0 for the rst order relaxation scheme ; (5.11) = 1 for the second order relaxation scheme : and are the limited slopes of a1 W V and a2V p(W) respectively with the slope limiters de ned by (4.3) or (4.4). 4) Find UjN+ 12 by (5.12)
UjN+ 21 = 12 (VjN + VjN+1) + 21 a1(WjN+1 ? WjN ) + 14 x ( j+;N + j?+1;N ) :
By the above analysis, this algorithm will guarantee that the rst order scheme produces the correct viscosity-capillarity solution for any ! 2 (0; 41 ]. However, we
numerical solutions of mixed type systems
15
do not observe the same convergence for the second order scheme in the numerical example given below. We test these schemes for a piecewise-linear equation of state 8 ? 20w if w < 0:1; > < (5.13) p(w) = > ? 3 + 10w if 0:1 w < 0:2; : ? 5w if w 0:2: The elliptic region is the domain w 2 (0:1; 0:2). The Riemann problem for the above equation of state (stress-strain relation) was solved by Abeyaratne and Knowles [1]. In our numerical experiments we take the initial data (5.14) (VL ; WL ) = (0; 0:05); (VR ; WR ) = (0:5; 0:07) : This problem is challenging because the solution not only consists of both shocks and phase boundaries but also depends sensitively on !. One needs far more grid points to obtain a reasonable resolution than the ones we studied in previous sections. In Figure 6 we depict the numerical results of the rst order and the second order (with slope limiter (4.3)) schemes for ! = 1=4 and ! = 0:0625 respectively at t = 0:05. In each run p we take 2000 spatial cells between [?0:5; 0:5], a1 = 20, the CFL number is 1 for the rst order scheme and 1=2 for the second order scheme. For the rst order scheme we observe the numerical convergence, and the solution is more dissipative with smaller !. The second order scheme, with a higher resolution though, fails to converge near the phase boundaries. We also used the minmod slope limiter but the results essentially remain the same. We think that the slope limiter should be modi ed so it is more sensitive to the discontinuities, or in (3.1) should be turned on to give more viscosity and capillarity. This remains unsolved and is a subject of future research. Acknowledgement. I am grateful to Bernardo Cockburn and Huiing Gau for showing me their recent work [3] and giving me the exact solution code for the test problem in section 5. In addition, I would like to thank John Lowengrub, Marshall Slemrod and Zhouping Xin for some helpful discussions on mixed type problems.
[1] [2] [3] [4] [5] [6] [7] [8]
REFERENCES R. Abeyaratne and J.K. Knowles, Kinetic relations and the propagation of phase boundaries in solids, Arch. Rat'l Mech. Anal., 114 (1991), pp. 119{154. M. Affouf and R. Caflisch, A numerical study of Riemann problem solutions and stability for a system of viscous conservation laws of mixed type, SIAM J. Appl. Math., 51 (1991), pp. 605-634. B. Cockburn and H. Gau, Model numerical schemes for phase transitions in solids, SIAM J. Sci. Comp., submitted. R. James, Dynamic phase changes with non-monotone stresses, Arch. Rat'l Mech. Anal., 66 (1983), pp. 59{125. S. Jin and Z.P. Xin, The relaxation schemes for systems of hyperbolic conservation laws in arbitrary space dimensions,, Comm. Pure Appl. Math., to appear. P.D. Lax, Hyperbolic systems of conservation laws and the mathematical theory of shock waves, 11 (1973). B.L. Keyfitz and M. Shearer, eds., Nonlinear Evolution Equations that Change Type, The IMA Volumes in Mathematics and its Applications, Springer-Verlag, 1990. T.-P. Liu, Hyperbolic conservation laws with relaxation, Comm. Math. Phys., 108 (1987), pp. 153{175.
16
shi jin
0.4
0.35
0.3
0.3
0.25
0.25
W
W
0.35
0.4 (a)
0.2
0.15
0.1
0.1
0.05
0.05
-0.4
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0 -0.5
0.5
0.4
0.35
(c)
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0.5
0.3
0.25
0.25
0.2
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0.5
0.2
0.15
0.15
0.1
0.1
0.05
0.05
-0.4
(d)
0.35
0.3
0 -0.5
-0.4
0.4
W
W
0.2
0.15
0 -0.5
(b)
-0.3
-0.2
-0.1
0 x
0.1
0.2
0.3
0.4
0.5
0 -0.5
-0.4
Fig. 6. Numerical solutions of w for the piecewise-linear equation of state (5.13) with initial data (5.14). In each run 2000 spatial cells are used. At t = 0:05, the exact solution is plotted by the solid lines, and the numerical solutions by the dashed lines. In (a) (! = 1=4) and (b) (! = 0:0625) the rst order scheme is applied with the CFL number = 1. In (c) (! = 1=4) and (d) (! = 0:0625) the second order scheme with slope limiter (4.3) is applied with the CFL number = 1=2.
[9] E.B. Pitman and Y. Ni, Visco-elastic relaxation with a van der Waals type stress, Int. J. Eng. Sci., 32 (1994), pp. 327{338. [10] M. Shearer, Non-uniqueness of admissible solutions of conservation laws of mixed type, Arch. Rat'l Mech. Anal., 93 (1986), pp. 45{59. [11] C.-W. Shu, A numerical method for systems of conservation laws of mixed type admitting hyperbolic ux splitting, J. Comp. Phys., 100 (1992), pp. 424{429. [12] C.-W. Shu and S. Osher, Ecient implementation of essentially non-oscillatory shock capturing methods, J. Comp. Phys., 77 (1988), pp. 439{471. [13] M. Slemrod, Admissibility criteria for propagating phase boundaries in a van der Waals uid, Arch. Rat'l Mech. Anal., 81 (1983), pp. 301{315. [14] M. Slemrod, Lax-Friedrichs and the viscosity-capillarity criterion, in Physical Mathematics and Nonlinear Partial Di erential Equations, S.M. Rankin and J.H. Lightbourne, eds., 1985. [15] M. Slemrod and J.E. Flaherty, Numerical integration of a Riemann problem for a van der Waals uid, in Phase Transformations, C.A. Elias and G. John, eds., 1986. [16] A. Sommerfeld, Thermodynamics and Statistical Mechanics, Academic Press, New York and
numerical solutions of mixed type systems
17
London, 1964. [17] L. Truskinovski, Equilibrium phase interface, Dokl. Akad. Nauk. SSSH, 256 (1982), pp. 306{ 310. [18] B. van Leer, Toward the ultimate conservative di erence schemes V. A second-order sequal to Godunov's method, J. Comp. Phys., 32 (1979), pp. 101{136.
