- Email: [email protected]
Computation of travelling wave solutions of scalar conservation laws with a stiff source term
Computers & Fluids 32 (2003) 1161–1178 www.elsevier.com/locate/compfluid
Computation of travelling wave solutions of scalar conservation laws with a stiff source term Vicente Martınez b
a,*
, Antonio Marquina
b
a Departamento de Matem aticas, Universitat Jaume I, Campus de Riu Sec, 12071 Castell o, Spain Departamento de Matem atica Aplicada, Universitat de Val encia, Dr. Moliner, 50 46100 Burjassot (Valencia), Spain
Accepted 24 July 2002
Abstract In this paper we propose a nonoscillatory numerical technique to compute the travelling wave solution of scalar conservation laws with a stiff source term. This procedure is based on the dynamical behavior described by the associated stationary ODE and it reduces/avoids numerical errors usually encountered with these problems, i.e., spurious oscillations and incorrect wave propagation speed. We combine this treatment with either the first order Lax–Friedrichs scheme or the second order Nessyahu–Tadmor scheme. We have tested several model problems by LeVeque and Yee for which the stiffness coefficient can be increased. We have also tested a problem with a nonlinear flux and a discontinuous source term. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Travelling waves solution; Nonlinear conservation laws; Stiff source term
1. Introduction The models of nonequilibrium gas dynamics, where chemical reactions occur, contain physical magnitudes varying in very different scales of time and space. Numerical approximations to the solution of these problems may develop spurious oscillations near discontinuities and/or incorrect wave propagation speed (see [4,17]). This question has been addressed by various authors from different points of view (see [9,17,22] for the one-dimensional case and [3,5,21,24] for systems of equations).
*
Corresponding author. E-mail addresses: [email protected] (V. Martınez), [email protected] (A. Marquina).
0045-7930/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0045-7930(02)00079-8
1162
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
We shall restrict our analysis to one-dimensional scalar conservation laws with source terms (also called balance laws). In the case of models of nonequilibrium gas dynamics, the source term is stiff. We consider bounded solutions to the nonhomogeneous initial-boundary value problem ut þ f ðuÞx ¼ wðuÞ; ð1Þ uðx; 0Þ ¼ u0 ðxÞ; for ðx; tÞ 2 ½a; b 0; T with T > 0, f , w 2 C 1 ð½a; bÞ and the initial data u0 ðxÞ are supposed to be a piecewise-smooth function. We consider absorbing boundary conditions, i.e., homogeneous Neumann boundary conditions. The flux function f describes the fluid flow and the source term w represents the chemistry of the reacting species usually present in this kind of models. We assume the following hypothesis in order to get weak entropy solutions (see [8,15]): ðH1Þ 9M0 > 0 such that juj P M0 ) u wðuÞ 6 0: In this paper, we construct a numerical procedure to approximate a travelling wave solution of Eq. (1) with stiff source term wðuÞ. The procedure is based on a travelling wave approximation of the viscous equation with the same dynamical behavior as the original one. Travelling wave analysis is an analytical tool used by various authors in different problems and applications (see [3,6,7,10,11,18,23]). The paper is organized as follows: Section 2 contains the travelling viscous wave analysis needed to devise the procedure to compute the contribution of the source term for a finite difference scheme in conservation form; Section 3 describes the procedure; Section 4 includes some numerical tests that show the power of resolution of our algorithm; The concluding remarks appear in Section 5.
2. Travelling wave analysis In order to study inviscid solutions of Eq. (1), a classical procedure (see [23] for homogeneous conservation laws) is to consider the scalar viscous conservation law with a source term ut þ f ðuÞx ¼ wðuÞ þ euxx
ð2Þ
under the same initial-boundary conditions as in (1) with e > 0 small. Our goal in this section is to approximate the solutions of (2) by means of a homogeneous conservation law, where the flux function takes into account the contribution of the source term. In nonequilibrium gas dynamics, in which the source term is stiff, the change between different stages is very fast and connected with the stiffness. Our analysis concerns a model with two constant stages: u and uþ . If p is the discontinuity point in the solution x ¼ xðtÞ of (1) in the ðx; tÞplane, and s ¼ x0 ðtp Þ, then we assume that (2) has solutions, called travelling wave solutions of the form: uðx; tÞ ¼ uðnÞ;
n¼
x st ; e
satisfying the boundary conditions
ð3Þ
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
lim uðnÞ ¼ u ;
lim uðnÞ ¼ uþ
n! 1
n!þ1
1163
ð4Þ
and lim u0 ðnÞ ¼ 0;
lim u0 ðnÞ ¼ 0:
n! 1
n!þ1
ð5Þ
The functions uðxÞ, in each time step, approach a jump discontinuity in the limit as e goes to zero (see [14]). Furthermore, these solutions are usually the ‘‘physically relevant’’ ones (see [23]). To obtain existence of travelling wave solutions in Eq. (2), we assume that and f 0 ðuÞ is strictly increasing:
ðH2Þ
f 2 C 2 ðRÞ
ðH3Þ
wðuÞ has isolated simple zeros:
Now, we are ready to recall the following result: Theorem 2.1. If (H1), (H2), (H3) are satisfied and u P uþ , then there is a unique travelling wave solution of Eqs. (2) and (3). Proof. See [6,7,18].
If u is a solution of (2) and (3), we have uðnÞt þ f ðuðnÞÞx ¼ wðuðnÞÞ þ euðnÞxx and therefore, s 1 1 00 0 0 0 þ f ðuðnÞÞu ðnÞ ¼ wðuðnÞÞ þ e 2 u ðnÞ ; u ðnÞ e e e we multiply by e and we obtain the second order ODE: u00 ¼ ðf 0 ðuÞ sÞu0 ewðuÞ; that we can rewrite as u00 ¼ ðf 0 ðuÞu0 ewðuÞÞ su0 :
ð6Þ
On the other hand, we consider the following homogeneous problem: ut þ F ðu; x; tÞx ¼ 0; where
x st F ðu; x; tÞ ¼ f ðuÞ eW e
and WðnÞ ¼
Z
n
wðuðsÞÞ ds:
1
It is easy to check that a travelling wave of the viscous equation ut þ F ðu; x; tÞx ¼ euxx
ð7Þ
1164
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
satisfies the second order ordinary differential equation (6). Therefore, (1) and (7) have the same travelling wave solutions. Thus, to approximate travelling wave solutions to our original problem we will solve the homogeneous problem (6) by using numerical methods in conservation form. We shall explain this in the next section.
3. Numerical procedure: adaptive grid points for quadrature In this section we are going to describe how to compute numerical fluxes and how to integrate the source term. We use a uniform computational grid for the spatial variable x 2 ða; bÞ with stepsize h > 0, defined as xj ¼ a þ jh, j ¼ 0; 1; . . . ; N. Our time discretization is defined as tn ¼ nk; n ¼ 0; 1; 2 . . . M, where k ¼ Dt is the time step. We denote by unj ¼ uðxj ; tn Þ the approximation of u at the grid point ðxj ; tn Þ. The average of uð ; tn Þ on ½xj 12 ; xjþ12 is defined by Z 1 xjþ12 n uj ¼ uðx; tn Þ dx: h x 1 j
2
In each computational cell ½xj 12 ; xjþ12 ½tn ; tnþ1 , we consider k ¼ k=h and the numerical flux function Z 1 tnþ1 n F ðuðxjþ12 ; tÞ; xjþ12 ; tÞ dt; Fjþ 1 2 k tn hence n n unþ1 ¼ unj kðFjþ 1 Fj 1 Þ: j 2
ð8Þ
2
Therefore, in this model, the finite difference scheme for (8) in conservation form becomes n n n unþ1 ¼ unj kðfjþ 1 fj 1 Þ þ kCjþ1 ; j 2
2
ð9Þ
2
n where the numerical flux function fjþ 1 is given by 2
n fjþ 1 2
¼ f~ðunj mþ1 ; . . . ; unjþm Þ;
as a function of 2m variables, which is consistent with (1), i.e., f~ðu; . . . ; uÞ ¼ f ðuÞ, Z rn jþ1 2 Cnjþ1 ¼ wðuðs; tn ÞÞ ds 2
ð10Þ
rn
j 1 2
and rnjþ1 ¼
xjþ12 stn
2
e
:
For the sake of simplicity let us assume that s > 0 and that the discontinuity points pn :¼ ðxpn ; tn Þ, n ¼ 0; 1; . . . ; M, satisfy xpn 2 ða; bÞ, n ¼ 0; 1; . . . ; M, i.e., unj ¼ u
8xj < a; tn < tM
and unj ¼ uþ
8xj > b; tn < tM ;
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1165
Fig. 1. Range of influence of the Lax–Friedrichs first order scheme.
then, our procedure will be restricted to the computational domain X ¼ fðx; tÞ : a < x < b; 0 6 t 6 tM g: For a point x 2 ða; bÞ, the set Rðx Þ ¼ fðx; tÞ 2 X : jx x j 6 ctg is called the range of influence of x , for some constant value c (see [16]). For instance, if we consider a three-point scheme (see Fig. 1), we have Rðx Þ ¼ fðxj ; tn Þ 2 X : jxj x j 6 nh; j ¼ 0; 1; . . . ; N ; n ¼ 0; 1; . . . ; Mg: Let us notice that the auxiliary mesh frnjþ1 gj;n is only used to compute Cnjþ1 by controlling the 2 2 contribution of this term to the scheme (9). Since the source term is included in integral (10), a correct numerical treatment of Cnjþ1 is fundamental to avoid spurious oscillations. 2 In order to avoid the transfer of wrong values to the numerical scheme, the auxiliary mesh must satisfy rnjþ1 62 Rðxp0 Þ;
j ¼ 0; 1; . . . ; N 1; n ¼ 0; 1; . . . ; M:
2
ð11Þ
The introduction of an auxiliary mesh is the most interesting feature of this procedure. To show that it becomes simple we will describe the structure of this mesh: n o Lemma 3.1. The time dependent auxiliary mesh rnjþ1 has space step h=e. 2
j;n
The proof is trivial. Lemma 3.2. In each time step the auxiliary mesh is shifted by sk=e. The proof (see Fig. 2) follows from rnjþ1 rnþ1 ¼ jþ1 2
2
sðtnþ1 tn Þ ; e
j ¼ 0; 1; . . . ; N 1; n ¼ 0; 1; . . . ; M 1:
Lemma 3.3. If k ¼ k=h ¼ l=sq, with l and q as integer numbers, then the auxiliary mesh repeats its values with period q and carries out a displacement of l cells.
1166
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 2. Evolution in time of the auxiliary mesh with s > 0.
Proof. For the sake of convenience, we assume M q > 0 and N l > 0. Then, for j ¼ 0; 1; . . . ; N l and n ¼ 0; 1; . . . ; M q, we compute xjþ12þl stnþq xjþ12 þ stn xjþ12 þ lh stn sqk xjþ12 þ stn lh sqk n ¼ 0: rnþq
r ¼ ¼ 1 ¼ 1 jþ jþ2þl 2 e e e In summary, if we use k ¼ l=sq and the Lax–Friedrichs first order scheme (see Fig. 3), e may be chosen according to the following restriction: rnjþ1 62 ½xp0 d; xp0 þ Mh; 2
n ¼ 0; 1; . . . ; q and j ¼ 0; 1; . . . ; N l; for some d > 0; ð12Þ
that is simpler than (11). Let us notice that the auxiliary mesh has space step h=e, so it allows us to space its points as often as necessary, taking a suitable e to avoid numerical instabilities caused by (10). We propose the following numerical recipe to avoid spurious oscillations and to obtain a correct wave propagation speed: Step 1: Given s, choose h; k; N and M so that k satisfies Lemma 3.3 and a Courant–Friedrichs– Lewy (CFL) restriction.
Fig. 3. Evolution in time of the discontinuity points with s > 0.
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1167
Step 2: Given M and xp0 , choose e such that (12) is verified by frnjþ1 gj;n . 2 Step 3: Approximate Cnjþ1 by numerical quadrature using frnjþ1 gj;n . 2 2 Step 4: Evolve in time using Eq. (9) with a finite difference scheme in conservation form. Remark 3.1. Since Z x 1 Z jþ 2 wðuðx; tn ÞÞ dx ¼ xj 1
xjþ1 2
xj 1
2
2
Z rn x st jþ1 n 2 w u wðuðnÞÞ dn; dx ¼ e e rn 1 j 2
a suitable choice of e can handle stiff problems. Our numerical procedure is based on the integration of Eq. (7) instead of Eq. (1), i.e., on the absorption of the source term into the flux via integration. This procedure evokes the work of Koren [13]. However, this paper introduces the e parameter, which should be chosen small enough in order to control the stiffness of the source term (usually e h for the models that we have indicated in precedent sections). To end the section, let us notice that the procedure needs the wave speed, so if it is unknown an approximation should be found. 4. Numerical experiments To study the behavior of our numerical procedure, we consider the first order Lax–Friedrichs scheme (LAXF1) given, for Eq. (1), by n n n ¼ unj kðfjþ unþ1 1 fj 1 Þ þ kwðuj Þ; j 2
2
1 1 n ðf ðunj Þ þ f ðunjþ1 ÞÞ ðDiþ12 Þ; fjþ 1 ¼ 2 2 2k where Diþ12 ¼ uiþ1 ui : If we use LAXF1 together with the numerical procedure of Section 3 and the trapezoidal rule to approximate Cnjþ1 , we refer to it as LAXF1R. We will also study the effects that the previous 2 numerical procedure produces when a high order method is used. Here we will also consider the Nessyahu–Tadmor method (NT2), which can be written (see [19]) in the following predictor– corrector form 8 1 < unþ2 ¼ un 1kðf n Þ0 ; j j j 2 n o 1 nþ12 nþ1 1 n n 2 : u 1 ¼ ðuj þ ujþ1 Þ þ 1ððunj Þ0 ðunjþ1 Þ0 Þ k f ðunþ Þ
f ðu Þ þ kwðunj Þ; j jþ1 jþ2
2
8
ðunj Þ0 and ðfjn Þ0 can be approached by u0j ¼ minmodða Dujþ12 ; 12ðujþ1 uj Þ; a Duj 12 Þ; fj0 ¼ minmodða Dfjþ12 ; 12ðfjþ1 fj Þ; a Dfj 12 Þ;
1168
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
where 16a < 4
and minmodðr; sÞ ¼ 12 minðjrj; jsjÞðsignðrÞ þ signðsÞÞ:
In the same way as previously, when we use NT2 together with the recipe of Section 3 and the trapezoidal rule to approximate Cnjþ1 , we refer to it as NT2R. 2
4.1. Model problem 1 A simple model with source term representing the chemistry in gas dynamics (see [17]) is given by ut þ ux ¼ luðu 1Þðu 1=2Þ;
ðx; tÞ 2 R ð0; T Þ
and T > 0;
ð13Þ
which is a stiff problem for large values of the parameter l, called stiffness coefficient. In this model, the presence of different time scales also depends on the parameters l and h. LeVeque and Yee (see [17]) have considered stable second order accurate numerical methods, such as the MacCormack or splitting methods, but all of them are subject to incorrect propagation speeds of discontinuities (a similar behavior was reported in [4]). LeVeque and Yee have identified the quantity kl as the critical parameter affecting the stability of the problem. They have observed numerical difficulties when kl is much larger than 1. Berm udez and Vazquez (see [1]) have addressed this problem by using upwind discretizations instead of centred discretizations of the source term, but this technique is insufficient for large parameter values of l. The exact solution of (13) with the initial data 1; x < 0:3; ð14Þ u0 ðxÞ ¼ 0; x > 0:3; is the function uðx; tÞ ¼ u0 ðx tÞ. 4.1.1. Dynamical study of the travelling wave solutions From (6), we have u00 ¼ ewð/Þ ¼ eluðu 1Þðu 0:5Þ; if we consider y1 ¼ u and y2 ¼ u0 , we obtain the ordinary differential system 0 y1 y2 : ¼ y2 ely1 ðy1 1Þðy1 0:5Þ Solving it, we have 1 2 y 2 2
14elðy14 2y13 þ y12 Þ ¼ E;
ð15Þ
where E is a constant. The phase portrait of the travelling wave approximate solutions can be obtained from (15) and is shown in Fig. 4 for el ¼ 2. From (4) and (5) it is deduced that the trajectory of the phase portrait, which is the suitable travelling wave solution associated with problems (13) and (14), is the one that connects ðu ; 0Þ and ðuþ ; 0Þ. In this case, it is a separatrix connecting the equilibria (1; 0) and (0; 0), which we denote by uP 1 and it is given by
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1169
Fig. 4. Phase portrait of the travelling wave solutions for problem 1 and el ¼ 2.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi el 4 ðy 2y13 þ y12 Þ: y2 ¼ 2 1 This function reaches the minimum in rffiffiffiffiffi 1 el ; : 2 32 In Fig. 11, we can see uP 1 for several values of el. Keeping in mind that uP 1 approaches the jump discontinuity, as e ! 0, the movement between (1; 0) and (0; 0) has to be smooth to avoid spurious oscillations introduced by the source term in the numerical scheme for intermediate points between (1; 0) and (0; 0). Therefore, we can conclude that the parameter el should be small. 4.1.2. Numerical results We have used k ¼ 34, xp0 ¼ 0:3, tM ¼ 0:3, k ¼ 0:015 (M ¼ 20) and h ¼ 0:02 (N ¼ 100). Then we will be able to use ð 1; 1Þ as reference interval, since all the discontinuity points are in it from t ¼ t0 to t ¼ tM . We will also be able to use q ¼ 4, since k ¼ 34, so we only have to consider the first four time steps to check that the points of the auxiliary mesh (9) do not fall inside the critical interval ½0:3 d; 0:7 (according to (12)). In LAXF1R we take e ¼ 0:007, so that all the nearest values of the auxiliary mesh fall outside the critical interval (see Table 1). In Fig. 5 the results of LAXF1 agree with the references given in the introduction, but if we use LAXF1R with e ¼ 0:007, even for l ¼ 500; 000:0, the numerical oscillations will be avoided. However, the previous value of e is not small enough when tM ¼ 0:33 (M ¼ 22) because at this time the value 0.71428571 of the auxiliary mesh (see Table 1) falls inside the critical interval ½0:3 d; 0:74, and it is necessary to take LAXF1R with e ¼ 0:005 (see Fig. 6).
1170
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Table 1 Nearest points of the auxiliary mesh to the critical interval given in (12) for problem 1 Time step n¼0 n¼1 n¼2 n¼3 n¼4
tn
rnj 1
rnjþ1
0.000 0.015 0.030 0.045 0.060
0.00000000 )2.14285714 )1.42857143 )0.71428571 0.00000000
2.85714286 0.71428571 1.42857143 2.14285714 2.85714286
2
2
Fig. 5. Solution profiles with LAXF1 and LAXF1R to problem 1, cfl ¼ 0:75 and tM ¼ 0:3.
Fig. 6. Solution profiles to problem 1 with LAXF1R, cfl ¼ 0:75, l ¼ 500; 000:0 and tM ¼ 0:33.
Regarding the method NT2, Fig. 7 shows that erroneous values are obtained with the method NT2 when l ¼ 800, in contrast to LAXF1, where l ¼ 30 already leads to large errors. But the critical interval is in this case ½0:3 d; 1:1, since the scheme NT2 works with k < 0:5. As a consequence, the critical interval is wider. Therefore, to avoid the numerical oscillations, we have chosen NT2R with e ¼ 0:0035. In Fig. 8, the unstable behavior of the solution caused by the erroneous values introduced by the scheme NT2R when e ¼ 0:004 can be seen.
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1171
Fig. 7. Solution profiles with NT2 and NT2R to problem 1, cfl ¼ 0:375 and tM ¼ 0:3.
Fig. 8. Solution profiles with NT2R to problem 1, cfl ¼ 0:375, l ¼ 500; 000:0 and tM ¼ 0:3.
Another interesting feature in a numerical scheme is the correct propagation wave speeds. To compare wave speeds for various values of kl in time tM , we use in this study the same average speed as the one use by LeVeque and Yee in [17], which is given by ! N N X h X M 0 c¼ u uj : tM j¼1 j j¼1 Fig. 9 shows this average as a function of kl for tM ¼ 0:3, N ¼ 100, LAXF1R with e ¼ 0:007 and cfl ¼ 0:75, LAXF1 with cfl ¼ 0:75, NT2 with cfl ¼ 0:375, and NT2R with e ¼ 0:0035 and cfl ¼ 0:375. We can conclude that the best approximation is obtained when we use the recipe of Section 3. 4.2. Model problem 2 Another problem typical of a stiff nonlinear source term (see [9,17]) is given by ut þ ux ¼ luðu 1Þ;
ðx; tÞ 2 R ð0; T Þ
and T > 0;
ð16Þ
whose exact solution with the initial data given in (14) is again the function uðx; tÞ ¼ u0 ðx tÞ.
1172
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 9. Comparison of average wave speeds for problem 1.
4.2.1. Dynamical study of the travelling wave solutions We obtain the phase portrait of the travelling wave solutions of problem 2 in the same way as in problem 1. In this case, we need to solve the following ordinary differential system 0 y2 y1 ; ¼ y2 ely1 ðy1 1Þ whose solutions satisfy 1 2 y 2 2
16elð2y13 3y12 Þ ¼ E;
ð17Þ
where E is a constant. In a similar way as in problem 1, we obtain the phase portrait of the travelling wave approximate solutions by giving values to (17), which can be seen in Fig. 10 for el ¼ 4. In this figure, we can see the suitable trajectory of the travelling wave solution associated with problem 2 defined by (14), (16) and (17) which is a separatrix connecting (1; 0) and ( 0:5; 0). This trajectory is given by rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi el ð2y13 3y12 þ 1Þ y2 ¼ 3 and it is denoted by uP 2 . In this case, uP 2 reaches the minimum in rffiffiffiffiffi el 0; : 3 We can see uP 2 , for several values of el in Fig. 11. 4.2.2. Numerical results The behavior of our recipe of Section 3 in this problem is similar to the behavior in problem 1. For this reason, in order to analyze the effectiveness of our proposal, it is enough to consider the
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1173
Fig. 10. Phase plane of the travelling wave solutions for problem 2 and el ¼ 4.
Fig. 11. Comparison of trajectories in problems 1 and 2 for different values of el.
LAXF1R scheme. The results can be checked in Fig. 12, where we have used the same data as in the previous problem. Regarding propagation wave speeds and the average speed c, similar results are also obtained that we do not reproduce here. 4.3. Model problem 3 Next, we will analyze the effectiveness of our procedure on a problem with a nonlinear flux. The following test problem was studied by Klingenstein (see [12]) and it is given by BurgersÕ equation. In [12], BurgersÕ equation is coupled to a chemical kinetics equation, in order to analyze a combustion model.
1174
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 12. Solution profiles with LAXF1 and LAXF1R (e ¼ 0:007) to problem 2.
In this section, we consider a simple model with BurgersÕ equation and a discontinuous source term: ut þ 12u2 x ¼ lðu aðuÞÞ; ðx; tÞ 2 R ð0; T Þ and T > 0; where aðuÞ ¼
1; 0;
u P 0:25; u < 0:25:
Again, this problem is stiff for large values of the parameter l. And the exact solution with the initial data given by 1; x < 0; u0 ðxÞ ¼ 0; x > 0; is the function uðx; tÞ ¼ u0 ðx 12tÞ. 4.3.1. Dynamical study of the travelling wave solutions The phase portrait of the travelling wave solutions of problem 3 is obtained in two parts. First, in the interval ð 1; 0:25Þ and in the same way as in previous problems, we have the following ordinary differential system 0 y2 y1 ¼ : ð18Þ y1 y2 þ ely1 12 y2 y2 Next, in the interval ð0:25; þ1Þ, we have 0 y2 y1 : ¼ y1 y2 þ elðy1 1Þ 12 y2 y2
ð19Þ
We obtain the phase portrait of the travelling wave approximate solutions by using qualitative analysis of the ordinary differential equations (18) and (19), and numerical techniques (see, for
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1175
Fig. 13. Phase portrait of the travelling wave solutions for problem 3 and el ¼ 163 .
example [2,20]). In Fig. 13, the phase portrait for el ¼ 163 can be seen. In this case, the suitable trajectory of the travelling wave solution associated with problem 3 is a piecewise function formed by the separatrix that connects (1; 0) and P , and the one that connects Q and (0; 0). So, this travelling wave solution has a discontinuous derivative (see Fig. 14). Fortunately, when e ¼ 0 Eqs. (18) and (19) are the same, then for e small enough, it is possible to obtain a smooth trajectory connecting (1; 0) and (0; 0).
Fig. 14. Travelling wave solution for problem 3 and el ¼ 163 .
1176
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 15. Phase portrait of the travelling wave solutions for problem 3 and e ¼ 0.
In Fig. 15, the phase portrait of the travelling wave solutions of problem 3 when e ¼ 0 can be seen. Fig. 16 shows the suitable trajectory of the travelling wave solution associated with problem 3 for this case. In summary, this dynamical analysis allows us to obtain a suitable trajectory of the travelling wave solution to approach the jump discontinuity. Taking e small enough, we have a smooth movement between (1; 0) and (0; 0). So, we avoid the spurious oscillations introduced by the source term in the numerical scheme for intermediate points between (1; 0) and (0; 0).
Fig. 16. Travelling wave solution for problem 3 and e ¼ 0.
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1177
Fig. 17. Solution profiles with LAXF1 and LAXF1R (e ¼ 0:003) to problem 3, tM ¼ 0:32.
4.3.2. Numerical results The behavior of our recipe of Section 3 in this problem is similar to the behavior in the previous problems, including the propagation wave speeds and the average speed c. To analyze the effectiveness of our proposal in this problem, we consider again the LAXF1R scheme. We have used N ¼ 100 and M ¼ 20, so in order to verify (12) we need k ¼ 0:016 (k ¼ 45) and e ¼ 0:003. It is important to highlight that the value of e has been obtained prior to the numerical calculation and that once this calculation has been carried out, the results (see Fig. 17) are in agreement with the theoretical analysis. 5. Concluding remarks We have analyzed the use of a numerical procedure on some models with a stiff source term. In order to control the contribution of the source term, we have introduced an auxiliary mesh, which avoids the transfer of wrong values to the numerical scheme. The method has been tested on some problems with discontinuous solutions and it has been shown to yield correct wave speed, also reducing spurious oscillations near discontinuities. Obviously, when the source term is nonstiff, this procedure is not worth using and standard schemes for nonstiff source terms should be used. The dynamical analysis of the problems allows us to deduce that problem 2 is more sensitive than problems 1 and 3, with respect to the propagation of erroneous values of the source term when the parameter kl is increased. Problems 1 and 3 verify hypothesis (H1), while problem 2 does not verify it. Then, problem 2 does not have a travelling wave solution that connects u with uþ . To approximate a solution to this problem we need a smaller e. For this reason, the minimum of uP 2 should be smaller than the minimum of uP 1 in problem 1. To approximate a solution to problem 2, we need a trajectory very close to the phase portrait axis 0X, (1; 0) and (0; 0). Acknowledgements Vicente Martınez acknowledges the continuous support from DGES(MEC), project PB970397, and from the Universitat Jaume I Fundaci o, Caixa Castell o, project P1B-97-23. Antonio Marquina was supported in part by a grant from DGICYT, PB97-1402.
1178
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
References [1] Berm udez A, Vazquez ME. Upwind methods for hyperbolic conservation laws with source terms. Comput Fluids 1994;23(8):1049–71. [2] Boyce WE, DiPrima RC. Elementary differential equations and boundary value problems. New York: John Wiley & Sons; 1992. [3] Burman E, Sainsaulieu L. Numerical analysis of two operator splitting methods for a hyperbolic system of conservation laws with stiff relaxation terms. Comput Meth Appl Mech Eng 1995;128:291–314. [4] Colella P, Majda A, Roytburd V. Theoretical and numerical structure for reacting shock waves. SIAM J Sci Stat Comput 1986;7:1059–80. [5] Engquist B, Sj€ ogreen B. Robust difference approximations of stiff inviscid detonation waves. CAM Report No. 9103, UCLA March 1991. [6] Fan H. Existence and uniqueness of traveling waves and error estimates for Godunov schemes of conservation laws. Math Comput 1998;67(221):87–109. [7] Fan H, Hale JK. Large time behavior in inhomogeneous conservation laws. Arch Ration Mech Anal 1993; 125(3):201–16. [8] Gosse L. A priori error estimate for a well-balanced scheme designed for inhomogeneous scalar conservation laws. CR Acad Sci Paris Ser I Math 1998;327(5):467–72. [9] Griffiths DF, Stuart AM, Yee HC. Numerical wave propagation in an advection equation with a nonlinear source term. SIAM J Numer Anal 1992;29(5):1244–60. [10] Harabetian E. A numerical method for viscous perturbations of hyperbolic conservation laws. SIAM J Numer Anal 1990;27(4):870–84. [11] Harabetian E. A subcell resolution method for viscous systems of conservation laws. J Comput Phys 1992;103: 350–8. [12] Klingenstein P. Hyperbolic conservation with source terms: errors of the shock location. Seminar for applied mathematics. Research Report No. 94-07, 1994. [13] Koren B. A robust upwind discretization method for advection, diffusion and source terms. In: Vreugdenhil CB, Koren B, editors. Numerical methods for advection–diffusion problems. Notes on numerical fluids mechanics, vol. 45. Braunschweig: Vieweg Verlag; 1993. p. 117–38. [14] Kreiss HO, Lorenz J. Initial-boundary value problems and the Navier–Stokes equations. New York: Academic Press; 1989. [15] Kruzkov SN. First order quasilinear equations in several independant space variables. Math USSR Sbornik 1970;10:217–43. [16] LeVeque RJ. Numerical Methods for Conservation Laws. Basel: Birkh€ auser; 1990. [17] LeVeque RJ, Yee HC. A study of numerical methods for hyperbolic conservation laws with stiff source terms. J Comput Phys 1990;86:187–210. [18] Mascia C. Travelling wave solutions for a balance law. Proc Roy Soc Edinburgh, Sect A 1997;127(3):567–93. [19] Nessyahu H, Tadmor E. Non-oscillatory central differencing for hyperbolic conservation laws. J Comput Phys 1990;87:408–63. [20] Parker TS, Chua LO. Practical numerical algorithms for chaotic systems. New York: Springer-Verlag; 1989. [21] Pember RB. Numerical methods for hyperbolic conservation laws with stiff relaxation. II. Higher-order Godunov methods. SIAM J Sci Comput 1993;14(4):824–59. [22] Sj€ ogreen B, Engquist B. Numerical approximation of hyperbolic conservation laws with stiff terms. In: Proceed 3rd Inter Conf Hyperb Problems, Uppsala, vols. I&II; 1990. p. 848–60. [23] Smoller J. Shock waves and reaction–diffusion equations. New York: Springer-Verlag; 1994. [24] Yee HC, Shinn J. Semi-implicit and fully implicit shock-capturing methods for hyperbolic conservation laws with stiff source terms. AIAA, Paper 87-1116, June 1987.
Computation of travelling wave solutions of scalar conservation laws with a stiff source term Vicente Martınez b
a,*
, Antonio Marquina
b
a Departamento de Matem aticas, Universitat Jaume I, Campus de Riu Sec, 12071 Castell o, Spain Departamento de Matem atica Aplicada, Universitat de Val encia, Dr. Moliner, 50 46100 Burjassot (Valencia), Spain
Accepted 24 July 2002
Abstract In this paper we propose a nonoscillatory numerical technique to compute the travelling wave solution of scalar conservation laws with a stiff source term. This procedure is based on the dynamical behavior described by the associated stationary ODE and it reduces/avoids numerical errors usually encountered with these problems, i.e., spurious oscillations and incorrect wave propagation speed. We combine this treatment with either the first order Lax–Friedrichs scheme or the second order Nessyahu–Tadmor scheme. We have tested several model problems by LeVeque and Yee for which the stiffness coefficient can be increased. We have also tested a problem with a nonlinear flux and a discontinuous source term. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Travelling waves solution; Nonlinear conservation laws; Stiff source term
1. Introduction The models of nonequilibrium gas dynamics, where chemical reactions occur, contain physical magnitudes varying in very different scales of time and space. Numerical approximations to the solution of these problems may develop spurious oscillations near discontinuities and/or incorrect wave propagation speed (see [4,17]). This question has been addressed by various authors from different points of view (see [9,17,22] for the one-dimensional case and [3,5,21,24] for systems of equations).
*
Corresponding author. E-mail addresses: [email protected] (V. Martınez), [email protected] (A. Marquina).
0045-7930/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0045-7930(02)00079-8
1162
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
We shall restrict our analysis to one-dimensional scalar conservation laws with source terms (also called balance laws). In the case of models of nonequilibrium gas dynamics, the source term is stiff. We consider bounded solutions to the nonhomogeneous initial-boundary value problem ut þ f ðuÞx ¼ wðuÞ; ð1Þ uðx; 0Þ ¼ u0 ðxÞ; for ðx; tÞ 2 ½a; b 0; T with T > 0, f , w 2 C 1 ð½a; bÞ and the initial data u0 ðxÞ are supposed to be a piecewise-smooth function. We consider absorbing boundary conditions, i.e., homogeneous Neumann boundary conditions. The flux function f describes the fluid flow and the source term w represents the chemistry of the reacting species usually present in this kind of models. We assume the following hypothesis in order to get weak entropy solutions (see [8,15]): ðH1Þ 9M0 > 0 such that juj P M0 ) u wðuÞ 6 0: In this paper, we construct a numerical procedure to approximate a travelling wave solution of Eq. (1) with stiff source term wðuÞ. The procedure is based on a travelling wave approximation of the viscous equation with the same dynamical behavior as the original one. Travelling wave analysis is an analytical tool used by various authors in different problems and applications (see [3,6,7,10,11,18,23]). The paper is organized as follows: Section 2 contains the travelling viscous wave analysis needed to devise the procedure to compute the contribution of the source term for a finite difference scheme in conservation form; Section 3 describes the procedure; Section 4 includes some numerical tests that show the power of resolution of our algorithm; The concluding remarks appear in Section 5.
2. Travelling wave analysis In order to study inviscid solutions of Eq. (1), a classical procedure (see [23] for homogeneous conservation laws) is to consider the scalar viscous conservation law with a source term ut þ f ðuÞx ¼ wðuÞ þ euxx
ð2Þ
under the same initial-boundary conditions as in (1) with e > 0 small. Our goal in this section is to approximate the solutions of (2) by means of a homogeneous conservation law, where the flux function takes into account the contribution of the source term. In nonequilibrium gas dynamics, in which the source term is stiff, the change between different stages is very fast and connected with the stiffness. Our analysis concerns a model with two constant stages: u and uþ . If p is the discontinuity point in the solution x ¼ xðtÞ of (1) in the ðx; tÞplane, and s ¼ x0 ðtp Þ, then we assume that (2) has solutions, called travelling wave solutions of the form: uðx; tÞ ¼ uðnÞ;
n¼
x st ; e
satisfying the boundary conditions
ð3Þ
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
lim uðnÞ ¼ u ;
lim uðnÞ ¼ uþ
n! 1
n!þ1
1163
ð4Þ
and lim u0 ðnÞ ¼ 0;
lim u0 ðnÞ ¼ 0:
n! 1
n!þ1
ð5Þ
The functions uðxÞ, in each time step, approach a jump discontinuity in the limit as e goes to zero (see [14]). Furthermore, these solutions are usually the ‘‘physically relevant’’ ones (see [23]). To obtain existence of travelling wave solutions in Eq. (2), we assume that and f 0 ðuÞ is strictly increasing:
ðH2Þ
f 2 C 2 ðRÞ
ðH3Þ
wðuÞ has isolated simple zeros:
Now, we are ready to recall the following result: Theorem 2.1. If (H1), (H2), (H3) are satisfied and u P uþ , then there is a unique travelling wave solution of Eqs. (2) and (3). Proof. See [6,7,18].
If u is a solution of (2) and (3), we have uðnÞt þ f ðuðnÞÞx ¼ wðuðnÞÞ þ euðnÞxx and therefore, s 1 1 00 0 0 0 þ f ðuðnÞÞu ðnÞ ¼ wðuðnÞÞ þ e 2 u ðnÞ ; u ðnÞ e e e we multiply by e and we obtain the second order ODE: u00 ¼ ðf 0 ðuÞ sÞu0 ewðuÞ; that we can rewrite as u00 ¼ ðf 0 ðuÞu0 ewðuÞÞ su0 :
ð6Þ
On the other hand, we consider the following homogeneous problem: ut þ F ðu; x; tÞx ¼ 0; where
x st F ðu; x; tÞ ¼ f ðuÞ eW e
and WðnÞ ¼
Z
n
wðuðsÞÞ ds:
1
It is easy to check that a travelling wave of the viscous equation ut þ F ðu; x; tÞx ¼ euxx
ð7Þ
1164
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
satisfies the second order ordinary differential equation (6). Therefore, (1) and (7) have the same travelling wave solutions. Thus, to approximate travelling wave solutions to our original problem we will solve the homogeneous problem (6) by using numerical methods in conservation form. We shall explain this in the next section.
3. Numerical procedure: adaptive grid points for quadrature In this section we are going to describe how to compute numerical fluxes and how to integrate the source term. We use a uniform computational grid for the spatial variable x 2 ða; bÞ with stepsize h > 0, defined as xj ¼ a þ jh, j ¼ 0; 1; . . . ; N. Our time discretization is defined as tn ¼ nk; n ¼ 0; 1; 2 . . . M, where k ¼ Dt is the time step. We denote by unj ¼ uðxj ; tn Þ the approximation of u at the grid point ðxj ; tn Þ. The average of uð ; tn Þ on ½xj 12 ; xjþ12 is defined by Z 1 xjþ12 n uj ¼ uðx; tn Þ dx: h x 1 j
2
In each computational cell ½xj 12 ; xjþ12 ½tn ; tnþ1 , we consider k ¼ k=h and the numerical flux function Z 1 tnþ1 n F ðuðxjþ12 ; tÞ; xjþ12 ; tÞ dt; Fjþ 1 2 k tn hence n n unþ1 ¼ unj kðFjþ 1 Fj 1 Þ: j 2
ð8Þ
2
Therefore, in this model, the finite difference scheme for (8) in conservation form becomes n n n unþ1 ¼ unj kðfjþ 1 fj 1 Þ þ kCjþ1 ; j 2
2
ð9Þ
2
n where the numerical flux function fjþ 1 is given by 2
n fjþ 1 2
¼ f~ðunj mþ1 ; . . . ; unjþm Þ;
as a function of 2m variables, which is consistent with (1), i.e., f~ðu; . . . ; uÞ ¼ f ðuÞ, Z rn jþ1 2 Cnjþ1 ¼ wðuðs; tn ÞÞ ds 2
ð10Þ
rn
j 1 2
and rnjþ1 ¼
xjþ12 stn
2
e
:
For the sake of simplicity let us assume that s > 0 and that the discontinuity points pn :¼ ðxpn ; tn Þ, n ¼ 0; 1; . . . ; M, satisfy xpn 2 ða; bÞ, n ¼ 0; 1; . . . ; M, i.e., unj ¼ u
8xj < a; tn < tM
and unj ¼ uþ
8xj > b; tn < tM ;
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1165
Fig. 1. Range of influence of the Lax–Friedrichs first order scheme.
then, our procedure will be restricted to the computational domain X ¼ fðx; tÞ : a < x < b; 0 6 t 6 tM g: For a point x 2 ða; bÞ, the set Rðx Þ ¼ fðx; tÞ 2 X : jx x j 6 ctg is called the range of influence of x , for some constant value c (see [16]). For instance, if we consider a three-point scheme (see Fig. 1), we have Rðx Þ ¼ fðxj ; tn Þ 2 X : jxj x j 6 nh; j ¼ 0; 1; . . . ; N ; n ¼ 0; 1; . . . ; Mg: Let us notice that the auxiliary mesh frnjþ1 gj;n is only used to compute Cnjþ1 by controlling the 2 2 contribution of this term to the scheme (9). Since the source term is included in integral (10), a correct numerical treatment of Cnjþ1 is fundamental to avoid spurious oscillations. 2 In order to avoid the transfer of wrong values to the numerical scheme, the auxiliary mesh must satisfy rnjþ1 62 Rðxp0 Þ;
j ¼ 0; 1; . . . ; N 1; n ¼ 0; 1; . . . ; M:
2
ð11Þ
The introduction of an auxiliary mesh is the most interesting feature of this procedure. To show that it becomes simple we will describe the structure of this mesh: n o Lemma 3.1. The time dependent auxiliary mesh rnjþ1 has space step h=e. 2
j;n
The proof is trivial. Lemma 3.2. In each time step the auxiliary mesh is shifted by sk=e. The proof (see Fig. 2) follows from rnjþ1 rnþ1 ¼ jþ1 2
2
sðtnþ1 tn Þ ; e
j ¼ 0; 1; . . . ; N 1; n ¼ 0; 1; . . . ; M 1:
Lemma 3.3. If k ¼ k=h ¼ l=sq, with l and q as integer numbers, then the auxiliary mesh repeats its values with period q and carries out a displacement of l cells.
1166
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 2. Evolution in time of the auxiliary mesh with s > 0.
Proof. For the sake of convenience, we assume M q > 0 and N l > 0. Then, for j ¼ 0; 1; . . . ; N l and n ¼ 0; 1; . . . ; M q, we compute xjþ12þl stnþq xjþ12 þ stn xjþ12 þ lh stn sqk xjþ12 þ stn lh sqk n ¼ 0: rnþq
r ¼ ¼ 1 ¼ 1 jþ jþ2þl 2 e e e In summary, if we use k ¼ l=sq and the Lax–Friedrichs first order scheme (see Fig. 3), e may be chosen according to the following restriction: rnjþ1 62 ½xp0 d; xp0 þ Mh; 2
n ¼ 0; 1; . . . ; q and j ¼ 0; 1; . . . ; N l; for some d > 0; ð12Þ
that is simpler than (11). Let us notice that the auxiliary mesh has space step h=e, so it allows us to space its points as often as necessary, taking a suitable e to avoid numerical instabilities caused by (10). We propose the following numerical recipe to avoid spurious oscillations and to obtain a correct wave propagation speed: Step 1: Given s, choose h; k; N and M so that k satisfies Lemma 3.3 and a Courant–Friedrichs– Lewy (CFL) restriction.
Fig. 3. Evolution in time of the discontinuity points with s > 0.
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1167
Step 2: Given M and xp0 , choose e such that (12) is verified by frnjþ1 gj;n . 2 Step 3: Approximate Cnjþ1 by numerical quadrature using frnjþ1 gj;n . 2 2 Step 4: Evolve in time using Eq. (9) with a finite difference scheme in conservation form. Remark 3.1. Since Z x 1 Z jþ 2 wðuðx; tn ÞÞ dx ¼ xj 1
xjþ1 2
xj 1
2
2
Z rn x st jþ1 n 2 w u wðuðnÞÞ dn; dx ¼ e e rn 1 j 2
a suitable choice of e can handle stiff problems. Our numerical procedure is based on the integration of Eq. (7) instead of Eq. (1), i.e., on the absorption of the source term into the flux via integration. This procedure evokes the work of Koren [13]. However, this paper introduces the e parameter, which should be chosen small enough in order to control the stiffness of the source term (usually e h for the models that we have indicated in precedent sections). To end the section, let us notice that the procedure needs the wave speed, so if it is unknown an approximation should be found. 4. Numerical experiments To study the behavior of our numerical procedure, we consider the first order Lax–Friedrichs scheme (LAXF1) given, for Eq. (1), by n n n ¼ unj kðfjþ unþ1 1 fj 1 Þ þ kwðuj Þ; j 2
2
1 1 n ðf ðunj Þ þ f ðunjþ1 ÞÞ ðDiþ12 Þ; fjþ 1 ¼ 2 2 2k where Diþ12 ¼ uiþ1 ui : If we use LAXF1 together with the numerical procedure of Section 3 and the trapezoidal rule to approximate Cnjþ1 , we refer to it as LAXF1R. We will also study the effects that the previous 2 numerical procedure produces when a high order method is used. Here we will also consider the Nessyahu–Tadmor method (NT2), which can be written (see [19]) in the following predictor– corrector form 8 1 < unþ2 ¼ un 1kðf n Þ0 ; j j j 2 n o 1 nþ12 nþ1 1 n n 2 : u 1 ¼ ðuj þ ujþ1 Þ þ 1ððunj Þ0 ðunjþ1 Þ0 Þ k f ðunþ Þ
f ðu Þ þ kwðunj Þ; j jþ1 jþ2
2
8
ðunj Þ0 and ðfjn Þ0 can be approached by u0j ¼ minmodða Dujþ12 ; 12ðujþ1 uj Þ; a Duj 12 Þ; fj0 ¼ minmodða Dfjþ12 ; 12ðfjþ1 fj Þ; a Dfj 12 Þ;
1168
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
where 16a < 4
and minmodðr; sÞ ¼ 12 minðjrj; jsjÞðsignðrÞ þ signðsÞÞ:
In the same way as previously, when we use NT2 together with the recipe of Section 3 and the trapezoidal rule to approximate Cnjþ1 , we refer to it as NT2R. 2
4.1. Model problem 1 A simple model with source term representing the chemistry in gas dynamics (see [17]) is given by ut þ ux ¼ luðu 1Þðu 1=2Þ;
ðx; tÞ 2 R ð0; T Þ
and T > 0;
ð13Þ
which is a stiff problem for large values of the parameter l, called stiffness coefficient. In this model, the presence of different time scales also depends on the parameters l and h. LeVeque and Yee (see [17]) have considered stable second order accurate numerical methods, such as the MacCormack or splitting methods, but all of them are subject to incorrect propagation speeds of discontinuities (a similar behavior was reported in [4]). LeVeque and Yee have identified the quantity kl as the critical parameter affecting the stability of the problem. They have observed numerical difficulties when kl is much larger than 1. Berm udez and Vazquez (see [1]) have addressed this problem by using upwind discretizations instead of centred discretizations of the source term, but this technique is insufficient for large parameter values of l. The exact solution of (13) with the initial data 1; x < 0:3; ð14Þ u0 ðxÞ ¼ 0; x > 0:3; is the function uðx; tÞ ¼ u0 ðx tÞ. 4.1.1. Dynamical study of the travelling wave solutions From (6), we have u00 ¼ ewð/Þ ¼ eluðu 1Þðu 0:5Þ; if we consider y1 ¼ u and y2 ¼ u0 , we obtain the ordinary differential system 0 y1 y2 : ¼ y2 ely1 ðy1 1Þðy1 0:5Þ Solving it, we have 1 2 y 2 2
14elðy14 2y13 þ y12 Þ ¼ E;
ð15Þ
where E is a constant. The phase portrait of the travelling wave approximate solutions can be obtained from (15) and is shown in Fig. 4 for el ¼ 2. From (4) and (5) it is deduced that the trajectory of the phase portrait, which is the suitable travelling wave solution associated with problems (13) and (14), is the one that connects ðu ; 0Þ and ðuþ ; 0Þ. In this case, it is a separatrix connecting the equilibria (1; 0) and (0; 0), which we denote by uP 1 and it is given by
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1169
Fig. 4. Phase portrait of the travelling wave solutions for problem 1 and el ¼ 2.
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi el 4 ðy 2y13 þ y12 Þ: y2 ¼ 2 1 This function reaches the minimum in rffiffiffiffiffi 1 el ; : 2 32 In Fig. 11, we can see uP 1 for several values of el. Keeping in mind that uP 1 approaches the jump discontinuity, as e ! 0, the movement between (1; 0) and (0; 0) has to be smooth to avoid spurious oscillations introduced by the source term in the numerical scheme for intermediate points between (1; 0) and (0; 0). Therefore, we can conclude that the parameter el should be small. 4.1.2. Numerical results We have used k ¼ 34, xp0 ¼ 0:3, tM ¼ 0:3, k ¼ 0:015 (M ¼ 20) and h ¼ 0:02 (N ¼ 100). Then we will be able to use ð 1; 1Þ as reference interval, since all the discontinuity points are in it from t ¼ t0 to t ¼ tM . We will also be able to use q ¼ 4, since k ¼ 34, so we only have to consider the first four time steps to check that the points of the auxiliary mesh (9) do not fall inside the critical interval ½0:3 d; 0:7 (according to (12)). In LAXF1R we take e ¼ 0:007, so that all the nearest values of the auxiliary mesh fall outside the critical interval (see Table 1). In Fig. 5 the results of LAXF1 agree with the references given in the introduction, but if we use LAXF1R with e ¼ 0:007, even for l ¼ 500; 000:0, the numerical oscillations will be avoided. However, the previous value of e is not small enough when tM ¼ 0:33 (M ¼ 22) because at this time the value 0.71428571 of the auxiliary mesh (see Table 1) falls inside the critical interval ½0:3 d; 0:74, and it is necessary to take LAXF1R with e ¼ 0:005 (see Fig. 6).
1170
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Table 1 Nearest points of the auxiliary mesh to the critical interval given in (12) for problem 1 Time step n¼0 n¼1 n¼2 n¼3 n¼4
tn
rnj 1
rnjþ1
0.000 0.015 0.030 0.045 0.060
0.00000000 )2.14285714 )1.42857143 )0.71428571 0.00000000
2.85714286 0.71428571 1.42857143 2.14285714 2.85714286
2
2
Fig. 5. Solution profiles with LAXF1 and LAXF1R to problem 1, cfl ¼ 0:75 and tM ¼ 0:3.
Fig. 6. Solution profiles to problem 1 with LAXF1R, cfl ¼ 0:75, l ¼ 500; 000:0 and tM ¼ 0:33.
Regarding the method NT2, Fig. 7 shows that erroneous values are obtained with the method NT2 when l ¼ 800, in contrast to LAXF1, where l ¼ 30 already leads to large errors. But the critical interval is in this case ½0:3 d; 1:1, since the scheme NT2 works with k < 0:5. As a consequence, the critical interval is wider. Therefore, to avoid the numerical oscillations, we have chosen NT2R with e ¼ 0:0035. In Fig. 8, the unstable behavior of the solution caused by the erroneous values introduced by the scheme NT2R when e ¼ 0:004 can be seen.
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1171
Fig. 7. Solution profiles with NT2 and NT2R to problem 1, cfl ¼ 0:375 and tM ¼ 0:3.
Fig. 8. Solution profiles with NT2R to problem 1, cfl ¼ 0:375, l ¼ 500; 000:0 and tM ¼ 0:3.
Another interesting feature in a numerical scheme is the correct propagation wave speeds. To compare wave speeds for various values of kl in time tM , we use in this study the same average speed as the one use by LeVeque and Yee in [17], which is given by ! N N X h X M 0 c¼ u uj : tM j¼1 j j¼1 Fig. 9 shows this average as a function of kl for tM ¼ 0:3, N ¼ 100, LAXF1R with e ¼ 0:007 and cfl ¼ 0:75, LAXF1 with cfl ¼ 0:75, NT2 with cfl ¼ 0:375, and NT2R with e ¼ 0:0035 and cfl ¼ 0:375. We can conclude that the best approximation is obtained when we use the recipe of Section 3. 4.2. Model problem 2 Another problem typical of a stiff nonlinear source term (see [9,17]) is given by ut þ ux ¼ luðu 1Þ;
ðx; tÞ 2 R ð0; T Þ
and T > 0;
ð16Þ
whose exact solution with the initial data given in (14) is again the function uðx; tÞ ¼ u0 ðx tÞ.
1172
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 9. Comparison of average wave speeds for problem 1.
4.2.1. Dynamical study of the travelling wave solutions We obtain the phase portrait of the travelling wave solutions of problem 2 in the same way as in problem 1. In this case, we need to solve the following ordinary differential system 0 y2 y1 ; ¼ y2 ely1 ðy1 1Þ whose solutions satisfy 1 2 y 2 2
16elð2y13 3y12 Þ ¼ E;
ð17Þ
where E is a constant. In a similar way as in problem 1, we obtain the phase portrait of the travelling wave approximate solutions by giving values to (17), which can be seen in Fig. 10 for el ¼ 4. In this figure, we can see the suitable trajectory of the travelling wave solution associated with problem 2 defined by (14), (16) and (17) which is a separatrix connecting (1; 0) and ( 0:5; 0). This trajectory is given by rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi el ð2y13 3y12 þ 1Þ y2 ¼ 3 and it is denoted by uP 2 . In this case, uP 2 reaches the minimum in rffiffiffiffiffi el 0; : 3 We can see uP 2 , for several values of el in Fig. 11. 4.2.2. Numerical results The behavior of our recipe of Section 3 in this problem is similar to the behavior in problem 1. For this reason, in order to analyze the effectiveness of our proposal, it is enough to consider the
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1173
Fig. 10. Phase plane of the travelling wave solutions for problem 2 and el ¼ 4.
Fig. 11. Comparison of trajectories in problems 1 and 2 for different values of el.
LAXF1R scheme. The results can be checked in Fig. 12, where we have used the same data as in the previous problem. Regarding propagation wave speeds and the average speed c, similar results are also obtained that we do not reproduce here. 4.3. Model problem 3 Next, we will analyze the effectiveness of our procedure on a problem with a nonlinear flux. The following test problem was studied by Klingenstein (see [12]) and it is given by BurgersÕ equation. In [12], BurgersÕ equation is coupled to a chemical kinetics equation, in order to analyze a combustion model.
1174
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 12. Solution profiles with LAXF1 and LAXF1R (e ¼ 0:007) to problem 2.
In this section, we consider a simple model with BurgersÕ equation and a discontinuous source term: ut þ 12u2 x ¼ lðu aðuÞÞ; ðx; tÞ 2 R ð0; T Þ and T > 0; where aðuÞ ¼
1; 0;
u P 0:25; u < 0:25:
Again, this problem is stiff for large values of the parameter l. And the exact solution with the initial data given by 1; x < 0; u0 ðxÞ ¼ 0; x > 0; is the function uðx; tÞ ¼ u0 ðx 12tÞ. 4.3.1. Dynamical study of the travelling wave solutions The phase portrait of the travelling wave solutions of problem 3 is obtained in two parts. First, in the interval ð 1; 0:25Þ and in the same way as in previous problems, we have the following ordinary differential system 0 y2 y1 ¼ : ð18Þ y1 y2 þ ely1 12 y2 y2 Next, in the interval ð0:25; þ1Þ, we have 0 y2 y1 : ¼ y1 y2 þ elðy1 1Þ 12 y2 y2
ð19Þ
We obtain the phase portrait of the travelling wave approximate solutions by using qualitative analysis of the ordinary differential equations (18) and (19), and numerical techniques (see, for
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1175
Fig. 13. Phase portrait of the travelling wave solutions for problem 3 and el ¼ 163 .
example [2,20]). In Fig. 13, the phase portrait for el ¼ 163 can be seen. In this case, the suitable trajectory of the travelling wave solution associated with problem 3 is a piecewise function formed by the separatrix that connects (1; 0) and P , and the one that connects Q and (0; 0). So, this travelling wave solution has a discontinuous derivative (see Fig. 14). Fortunately, when e ¼ 0 Eqs. (18) and (19) are the same, then for e small enough, it is possible to obtain a smooth trajectory connecting (1; 0) and (0; 0).
Fig. 14. Travelling wave solution for problem 3 and el ¼ 163 .
1176
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
Fig. 15. Phase portrait of the travelling wave solutions for problem 3 and e ¼ 0.
In Fig. 15, the phase portrait of the travelling wave solutions of problem 3 when e ¼ 0 can be seen. Fig. 16 shows the suitable trajectory of the travelling wave solution associated with problem 3 for this case. In summary, this dynamical analysis allows us to obtain a suitable trajectory of the travelling wave solution to approach the jump discontinuity. Taking e small enough, we have a smooth movement between (1; 0) and (0; 0). So, we avoid the spurious oscillations introduced by the source term in the numerical scheme for intermediate points between (1; 0) and (0; 0).
Fig. 16. Travelling wave solution for problem 3 and e ¼ 0.
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
1177
Fig. 17. Solution profiles with LAXF1 and LAXF1R (e ¼ 0:003) to problem 3, tM ¼ 0:32.
4.3.2. Numerical results The behavior of our recipe of Section 3 in this problem is similar to the behavior in the previous problems, including the propagation wave speeds and the average speed c. To analyze the effectiveness of our proposal in this problem, we consider again the LAXF1R scheme. We have used N ¼ 100 and M ¼ 20, so in order to verify (12) we need k ¼ 0:016 (k ¼ 45) and e ¼ 0:003. It is important to highlight that the value of e has been obtained prior to the numerical calculation and that once this calculation has been carried out, the results (see Fig. 17) are in agreement with the theoretical analysis. 5. Concluding remarks We have analyzed the use of a numerical procedure on some models with a stiff source term. In order to control the contribution of the source term, we have introduced an auxiliary mesh, which avoids the transfer of wrong values to the numerical scheme. The method has been tested on some problems with discontinuous solutions and it has been shown to yield correct wave speed, also reducing spurious oscillations near discontinuities. Obviously, when the source term is nonstiff, this procedure is not worth using and standard schemes for nonstiff source terms should be used. The dynamical analysis of the problems allows us to deduce that problem 2 is more sensitive than problems 1 and 3, with respect to the propagation of erroneous values of the source term when the parameter kl is increased. Problems 1 and 3 verify hypothesis (H1), while problem 2 does not verify it. Then, problem 2 does not have a travelling wave solution that connects u with uþ . To approximate a solution to this problem we need a smaller e. For this reason, the minimum of uP 2 should be smaller than the minimum of uP 1 in problem 1. To approximate a solution to problem 2, we need a trajectory very close to the phase portrait axis 0X, (1; 0) and (0; 0). Acknowledgements Vicente Martınez acknowledges the continuous support from DGES(MEC), project PB970397, and from the Universitat Jaume I Fundaci o, Caixa Castell o, project P1B-97-23. Antonio Marquina was supported in part by a grant from DGICYT, PB97-1402.
1178
V. Martınez, A. Marquina / Computers & Fluids 32 (2003) 1161–1178
References [1] Berm udez A, Vazquez ME. Upwind methods for hyperbolic conservation laws with source terms. Comput Fluids 1994;23(8):1049–71. [2] Boyce WE, DiPrima RC. Elementary differential equations and boundary value problems. New York: John Wiley & Sons; 1992. [3] Burman E, Sainsaulieu L. Numerical analysis of two operator splitting methods for a hyperbolic system of conservation laws with stiff relaxation terms. Comput Meth Appl Mech Eng 1995;128:291–314. [4] Colella P, Majda A, Roytburd V. Theoretical and numerical structure for reacting shock waves. SIAM J Sci Stat Comput 1986;7:1059–80. [5] Engquist B, Sj€ ogreen B. Robust difference approximations of stiff inviscid detonation waves. CAM Report No. 9103, UCLA March 1991. [6] Fan H. Existence and uniqueness of traveling waves and error estimates for Godunov schemes of conservation laws. Math Comput 1998;67(221):87–109. [7] Fan H, Hale JK. Large time behavior in inhomogeneous conservation laws. Arch Ration Mech Anal 1993; 125(3):201–16. [8] Gosse L. A priori error estimate for a well-balanced scheme designed for inhomogeneous scalar conservation laws. CR Acad Sci Paris Ser I Math 1998;327(5):467–72. [9] Griffiths DF, Stuart AM, Yee HC. Numerical wave propagation in an advection equation with a nonlinear source term. SIAM J Numer Anal 1992;29(5):1244–60. [10] Harabetian E. A numerical method for viscous perturbations of hyperbolic conservation laws. SIAM J Numer Anal 1990;27(4):870–84. [11] Harabetian E. A subcell resolution method for viscous systems of conservation laws. J Comput Phys 1992;103: 350–8. [12] Klingenstein P. Hyperbolic conservation with source terms: errors of the shock location. Seminar for applied mathematics. Research Report No. 94-07, 1994. [13] Koren B. A robust upwind discretization method for advection, diffusion and source terms. In: Vreugdenhil CB, Koren B, editors. Numerical methods for advection–diffusion problems. Notes on numerical fluids mechanics, vol. 45. Braunschweig: Vieweg Verlag; 1993. p. 117–38. [14] Kreiss HO, Lorenz J. Initial-boundary value problems and the Navier–Stokes equations. New York: Academic Press; 1989. [15] Kruzkov SN. First order quasilinear equations in several independant space variables. Math USSR Sbornik 1970;10:217–43. [16] LeVeque RJ. Numerical Methods for Conservation Laws. Basel: Birkh€ auser; 1990. [17] LeVeque RJ, Yee HC. A study of numerical methods for hyperbolic conservation laws with stiff source terms. J Comput Phys 1990;86:187–210. [18] Mascia C. Travelling wave solutions for a balance law. Proc Roy Soc Edinburgh, Sect A 1997;127(3):567–93. [19] Nessyahu H, Tadmor E. Non-oscillatory central differencing for hyperbolic conservation laws. J Comput Phys 1990;87:408–63. [20] Parker TS, Chua LO. Practical numerical algorithms for chaotic systems. New York: Springer-Verlag; 1989. [21] Pember RB. Numerical methods for hyperbolic conservation laws with stiff relaxation. II. Higher-order Godunov methods. SIAM J Sci Comput 1993;14(4):824–59. [22] Sj€ ogreen B, Engquist B. Numerical approximation of hyperbolic conservation laws with stiff terms. In: Proceed 3rd Inter Conf Hyperb Problems, Uppsala, vols. I&II; 1990. p. 848–60. [23] Smoller J. Shock waves and reaction–diffusion equations. New York: Springer-Verlag; 1994. [24] Yee HC, Shinn J. Semi-implicit and fully implicit shock-capturing methods for hyperbolic conservation laws with stiff source terms. AIAA, Paper 87-1116, June 1987.