Two approaches are used to extend the essentially non-oscillatory (ENO) schemes to treat conservation laws with stiff source terms. One approach is the application of the Strang time-splitting method. Here the basic ENO scheme and the Harten modification using subcell resolution (SR), ENO/SR scheme, are extended this way. The other approach is a direct method and a modification of the ENO/SR. Here the technique of ENO reconstruction with subcell resolution is used to locate the discontinuity within a cell and the time evolution is then accomplished by solving the differential equation along characteristics locally and advancing in the characteristic direction. This scheme is denoted ENO/SRCD (subcell resolution - characteristic direction). All the schemes are tested on the equation of LeVeque and Yee (NASA-TM-100075, 1988) modeling reacting flow problems. Numerical results show that these schemes handle this intriguing model problem very well, especially with ENO/SRCD which produces perfect resolution at the discontinuity

Chang, Shih-Hung

English

NASA Technical Reports Server

ON THE APPLICATION OF E N 0  SCHEME WITH SUBCELL ! / ] I  
RESOLUTION TO CONSERVATION LAWS WITH STIFF SOURCE: TERMS 
Shih-Hung Chang* 
Department of Mathematics, Cleveland State University 
Cleveland, Ohio 44115 
ABSTRACT 
Two approaches are used to extend EN0 schemes to treat conservation 'laws with stiff source 
terms. One approach is the application of Strang's time-splitting method. Here the basic EN0 
scheme and Harten's modification using subcell resolution, ENO/SR scheme, are extended this 
way. The other approach is a direct method and a modification of ENOISR. Here the technique 
of EN0 reconstruction with subcell resolution is used to locate the discontinuity within a cell and 
the time evolution is then accomplished by solving the differential equation along characteristics 
locally and advancing in the characteristic direction. This scheme is denoted ENOISRCD. All 
the schemes are tested on the equation of LeVeque and Yee (NASA TM 100075, 1988) model- 
ing reacting flow problems. Numerical results show that these schemes handle this intriguing 
model problem very well, especially with ENO/SRCD which produces perfect resolution at the 
discontinuity. 
1. INTRODUCTION 
Recently, in studying numerical methods for reacting flow problems, LeVeque and Yee (ref.5) 
considered certain fundamental questions concerning the extension of current finite difference 
techniques developed for non-reacting flows to reacting flows. Namely, can one: (i) develop stable 
methods, (ii) obtain "high resolution" results with sharp discontinuities and second order accuracy 
in smooth regions, and (iii) obtain the correct jumps at the correct locations? They introduced 
and studied the following one-dimensional scalar conservation law with parameter-dependent 
source term 
where p is a parameter. This equation becomes stiff when the parameter p is large. Although 
this linear advection equation with a source term represents only a simple model of reacting flow 
problems, by studying its numerical solutions one would encounter some of the difficulties sure 
to occur in solving more realistic models. 
*Research was partially supported under Cleveland State University RCAC Grant, Ohio 
OBOR Grant in the Research Challenge Program, and NASA Space Act Agreement C99066G 
while the author was in residence at the Institute for Computational Mechanics in Propulsion 
(ICOMP), NASA Lewis Research Center, Cleveland, Ohio 44135 
https://ntrs.nasa.gov/search.jsp?R=19910011761 2020-03-19T19:12:43+00:00Z
In their study, two different approaches were used to construct second order accurate nu- 
merical methods. One approaeh was to use a modification of MacCormack's predictor-corrector 
method for conservation laws (ref.6), together with two TVD-like versions with appropriate lim- 
iters (ref.8,9). The other approach was based on the second order accurate Strang's time-splitting 
method (ref.7). Their numerical tests revealed that stable and second order schemes can be de- 
vised by using either of these approaches. However, in studying the ability of these methods in 
dealing with propagating discontinuities, it was reported that for a certain reasonably fixed mesh 
and for the very stiff case, all the methods produced solutions that look reasonable and yet are 
completely wrong, because the discontinuities are in the wrong locations. Their investigation 
pointed out that the main difficulty is the smearing of the discontinuity in the spatial direction, 
which in turn introduced a nonequilibrium state into the calculation. To avoid this difficulty, 
it will be necessary to increase the resolution near the discontinuity, at  least for the purpose of 
evaluating $J (u) . 
The purpose of this paper is to show that numerical methods can be devised to overcome 
the above mentioned difficulties. We will demonstrate this by extending EN0 schemes to treat 
conservation laws with source terms. We will construct numerical schemes which, when applied 
to Eq.(l), will produce stable solutions with excellent resolutions at the correct locations of 
discontinuities. We choose to describe the extensions for the following equation 
where the source term $J(u) arises from the chemistry of the reacting species. It can be handled 
similarly for a < 0. 
The basic EN0 scheme (ref.4) depends on an essentially non-oscillatory reconstruction pro- 
cedure. Harten (ref.3) recently modified this procedure using subcell resolution. The notion of 
subcell resolution is based on the observation that unlike point values, cell averages of a discontin- 
uous piecewise-smooth function contain information about the exact location of the discontinuity 
within the cell. Using this observation in his study of conservation laws, Harten (ref.3) obtained 
a modification of the basic EN0 scheme, which he denoted ENOISR, achieving significant im- 
provement in the resolution of contact discontinuities. Basically, when good approximations to 
the locations of discontinuities inside the cells are obtained, it is then possible to have good 
reconstruction of the solution at each time step. Here we will also demonstrate that when the 
information on the location of the discontinuity is used in treating the source term, significant 
improvement in numerical results can be obtained. 
One approach that we use to extend EN0 and ENO/SR is the application of Strang's time- 
splitting method(ref.7), in which one alternates between solving the conservation law without the 
source term and the ordinary differential equation modeling the chemistry. The other approach is 
a direct method and a modification of ENOISR. Here we use the technique of EN0 reconstruction 
with subcell resolution to locate the discontinuity within a cell and then accomplish the time 
evolution by solving the differential equation along characteristics locally and advancing in the 
characteristic direction. Since the subcell resolution and the characteristic direction are essential 
in the design of this scheme, it is denoted ENO/SRCD. En ref.1, ENO/SRCD was originally 
implemented using the time-splitting method. 
In section 2, we will first describe the construction of ENOISRGD and then the extensions 
of EN0 and ENO/SR schemes for Eq.(3). In section 3, we present the numerical results obtained 
from using these schemes on the model problem of LeVeque and Yee (ref.5). A conclusion will be 
given in sectiori 4. 
2. CONSTRUCTION OF THE SCHEMES 
We observe that along the characteristic x = xo + a t ,  the solution to (3) evolves according 
to the ODE 
d 
-U(XO + a t , t )  = $(U(XO + a t , t ) ) ,  dt (4) 
with initial data u(xo, 0). In the scheme ENOISRCD, this equation will be solved approximately 
from the time step t, to t,+'. At t,, suppose that we have obtained the numerical solution 
vn  = {v?}, where vy represents an approximation to iiy, the cell average of u over [xj-r ,  xj+l]  
at t,. Then, to  obtain vn+', we use the following steps: 
1. Obtain a reconstruction, R(x; vn)  , of the solution from the given values vn.  
2. Locate the discontinuity, if any, within each cell using the subcell resolution technique and 
modify the reconstruction R(x; vn)  to obtain k ( z ;  vn) . 
3. Advance a ( z ;  vn) via the ODE (4) along the characteristics to the tn+' level and then take 
cell averages to complete vn+'. 
In ENO/SRCD, steps 1 and 2 will follow the EN0 reconstruction procedure with subcell 
resolution of Harten (ref.3). The reconstructed solution function R (x; vn)  here is a piecewise 
quadratic polynomial obtained by using the primitive function approach. For the sake of com- 
pleteness, we will describe in straightforward terms the procedures used. For more details and 
general discussions on reconstruction with subcell resolution, see ref.3. 
Step 1. EN0 Reconstruction 
Over each cell [xj-i, xj+t],  choose i such that j - 2 5 i < j and minimizes the following: 
Let Rj(x; vn)  denote the reconstructed quadratic polynomial over this cell. Then 
where 
c j  = (v:+~ - 2v:+' + v:)/(Ax)~, 
Step 2. Modification by Subcell Resolution 
To detect a discontinuity in a cell [zj-t ,zj+:], we define 
Zj+ ?j 
R j - l ( ~ ;  vn)  dx + (x; vn)  d ~ ]  - V: . 
In our numerical experiment, the following criterion is used. If 
we consider that there is a discontinuity at Bj  in this cell satisfying 
The location Bj can then be approximated by using any standard root-finding method. We simply 
use the bisection method in our experiment. 
Now, if there is a discontinuity inside the cell at Bj, a modified reconstruction $(x; vn )  is 
used , where 
Otherwise, we use 
$(x;vn) = Rj(x;vn). 
Step 3. Time Evolution and Cell Averaging 
Consider the case a At < Ax and that there exists a discontinuity at Bj  inside the cell 
[xj-4, x~+ .L]  2 with 
Bj 5 x,++ - a At, 
as shown in Figure 1. At the t, level, R j - l ( ~ ;  vn) is used as the solution to the left of Oj and 
Rj+1(x; vn)  to the right of Bj. Since the solution to Eq.(3) evolves according to the ODE (4) along 
characteristics, we can obtain approximate solution values at the tn+1 level by solving Eq.(4) and 
advancing in the characteristic direction. If w ( ~ , t , + ~ )  denotes the solution of Eq.(4) at tn+1 over 
[ x ~ - ~ , x ~ + L ] ,  obtained by using the initial values Rj - l (~ ;  vn) on [zj-+ -aAt,Bj) and Rj+1(x; vn)  
on (Bj , xi+ t - a At], the numerical solution v?+' is then an approximation to 
Fig. 1. The case B j  5 xi++ - a At 
In our present implementation we use the following simple computation. Let xm and xp 
denote the midpoints in the intervals (xj- ; - a At, 8,) and (Oj ,  x j+ t  - a At) respectively (see 
Fig. 1). Then we compute 
wm = Rj-1 ( ~ m ;  un) + At $(Rj-l(xm; un)), 
wp = Rj+r (xp; vn )  + At $ (R j+ l (~p ;  vn)), 
u;+' = [w, (8, - xj-L + a At) + wp (zj+; - B j  - a a t ) ]  /Ax. 
In the above computation, w, and wp are obtained from the Euler's method. If the modified 
Euler method is desired, for example, one simply computes 
similarly for wi and wp, and then v;+' by (10). Other locations of 0, and the cases with regions 
without discontinuities can be treated similarly and easily. It is quite simple to modify the above 
scheme to obtain higher order versions. 
The extension of E N 0  and ENO/SR schemes to treat conservation laws with source terms 
will be done by using S trang's time-splitting method (ref.7). With respect to Eq. (3))  the numerical 
solution un+'  is computed from vn by 
where Sf (Ad) represents the nmerical solution operator for the conservation law without the 
source term 
ut + a u,  = 8, a > 0, 
over a time step At, and s+(?) represents the numerical solution operator for the ordinary 
differential equation 
.t = 9(4 (13) 
over At/2. Thus the basic EN0 and ENO/SR are used here as the operator Sf(&). The extended 
schemes will still be denoted EN0 and ENO/SR respectively. 
The following second order version of the EN0 scheme has been used in ref.2. In our com- 
putational experiment, it produces slightly better results than other second order versions. 
EN0  Scheme: 
For the operator Sf(At),  we use 
where 
with 
where the sj's used in the above computation come from (6) in step 1 and f R(vL ,  vR) denotes 
the flux at the origin in a Riemann problem with v~ to the left and v~ to the right. 
Now, let us describe the operator §$(At). Here we will follow steps 1 and 2 to find the 
discontinuity Oj, if any. Let us use the same notations introduced before and refer to Fig.1. Also, 
let z, and zp denote the midpoints in the intervals (xj- +, O j )  and (8, , xj+ $) respectively. Then, 
for approximating cell averages and for the case O j  5 xj++ - a At, we use 
(16) 
where w, and wp are the same as in (10). Again, other situations are handled similarly. 
The resulting algorithm then takes the following form: 
ENO/SR Scheme: 
The operator Sf(At) is now replaced by Barten9s second order EN0 scheme with subcell 
resolution (ref.3). It is given in the form of Eq.(14) with 
-EN0 
where f j+ +, will be the same as in (15) and the correction term ijj+: is computed as follows. If 
the discontinuity condition (8) is not satisfied, then 
else 
[ (Ax-a  At) (v; - aAts j /2 )  - bj-l(xj-+,xj+t - a ~ t ) ] / ~ t ,  
when Fj(xj+; - a A t )  F~(X,-~)  > 0, G j + l  2 = 
[ b j + l ( ~ j + ~  - a A t , x j + r )  2 - aAt(v? + (AX - aAt )  sj/2)]/At, 
otherwise, 
and the expression bj (yl , y2) is used to mean 
In the above computation, all the cj's, sj's, and aj's come from (6) in step 1. The operator 
S* (At) will be the same as in (16) and the final algorithm also takes the same form as in (17). 
3. COMPUTATIONAL RESULTS 
We use the same fixed mesh and initial data as in the model problem of LeVeque and Yee 
(ref.5) to test the ability of the above schemes in dealing with propagating discontinuities. Thus, 
we apply all the schemes to Eq.(l) together with the initial condition 
We take Ax = 0.02, At = 0.015, and the domain in a: to be from 0 to I. For comprison with 
ref.5, we also show the results at t = 0.3 and for the cases p = 1, 10, 100, and 1000. For all cases, 
ENO/SRCD produces a perfect resolution as shown in Figure 2. Figure 3 shows the computed 
results using EN0 and ENB/SR schemes for p = 1, 10, and 100. Here one can see that the 
results ffom ENO/SR are also almost perfect. For the very stifT case, p = 1000, both E N 0  and 
ENO/SR fail to produce stable solutions for the above mesh in our computational experiment. 
However, when we reduce the size of At to one half of the original, i.e., At = 0.0075, and march 
40 time steps, very good result is again obtained from E N 0  and perfect resolution is obtained 
from ENO/SR as shown in Figures 4a and 4b respectively. We understand that reducing At 
means the reduction of stiffness of the problem. The difficulty arises from the fact that in both 
-EN0 EN0 and ENO/SR schemes, the computation of the numerical flux f ,+ still produces 'large " 
error in the spatial direction. 
The computational results obtained here compare favorably to those in LeVeque and Yee 
(ref.5). 
Fig. 2. Numerical results at t = 0.3 using ENO/SRCD scheme with 
discontinuous initial data and p = 1,10,100,1000. 
. true solution, . . .: computed solution A X  = 0.02, At = 0.015, -----. 
Fig. 3. Numerical results at  t = 0.3 using EN0 (first column) and ENO/SR 
(second column) schemes with discontinuous initial data for 
,u = 1 (first row), ,u r=: 10 (second row), and p = 100 (third row). 
Az = 0.02, At = 0.015, ---------. . true solution, . . .: computed solution 
Fig. 4a. Numerical results at t = 0.3 using EN0 scheme with 
discontinuous initial data and p = 1000. 
. true solution, . .: computed solution Ax = 0.02, At = 0.0075, -. 
Fig. 4b. Numerical results at  t = 0.3 using ENO/SR scheme with 
discontinuous initial data and p -- 1000. 
. true solution, - .: computed solution Ax = 0.02, A t  = 0.0075, 
4. CONCLUSIONS 
We have extended the basic EN0 and Harten's ENO/SR schemes to treat conservation laws 
with source terms. Two approaches are used. One is to apply Strang9s time-splitting method, 
in which one alternates between solving the conservation law without the source term and the 
ordinary differential equation modeling the chemistry. The other is a modification of EN0 /SR 
and a direct method, which uses the technique of EN0 reconstruction with subcell resolution to 
locate the discontinuity within a cell and then accomplishes the time evolution by solving the 
differential equation along characteristics locally and advancing in the characteristic direction. 
We call this scheme ENO/SRCD. All the schemes are tested on the equation of LeVeque and Yee 
(ref.5) modeling reacting flow problems. The ENO/SRCD scheme produces perfect resolution at 
the propagating discontinuity. The extensions of basic EN0 and ENO/SR via time-splitting also 
perform very well, especially with ENO/SR showing almost perfect results, except for the very 
stiff case where some adjustment in the time step-size is needed. 
REFERENCES 
1. Chang, S.H., "On the Application of Subcell Resolution to Conservation Laws with Stiff 
Source Terms," NASA Technical Memorandum 102384, November 1989. 
2. Chang, S.H. and Liou, M.S., "A Numerical Study of EN0 and TVD Schemes for Shock 
Capturing," NASA Technical Memorandum 101355, September 1988. 
3. Harten, A., "EN0 Schemes with Subcell Resolution; J. Comp. Phys., 83(1989), pp. 148-184. 
4. Harten, A., Engquist, B., Osher, S., and Chakravarthy, S., "Uniformly High Order Accurate 
Essentially Non-Oscillatory Schemes 111," J. Comp. Phys., 71(1987), pp. 231-303. 
5. LeVeque, R.J. and Yee, H.C., "A Study of Numerical Methods for Hyperbolic Conservation 
Laws with Stiff Source Terms," NASA Technical Memorandum 100075, March 1988. 
6. MacCormack, R.W., "The Effect of Viscosity in Hypervelocity Impact Cratering," AIAA 
Paper 69-354, 1969. 
7. Strang, G., "On the Construction and Comparison of Difference Schemes," SIAM J. Numer. 
Anal., 5(1968), pp. 506-517. 
8. Yee, H.C., "Upwind and Symmetric Shock-Capturing Schemes," NASA Ames Technical 
Memorandum 89464, 1987. 
9. Yee, H.C. and Shinn, J.L., "Semi-Implicit and Fully Implicit Shock-Capturing Methods for 
Hyperbolic Conservation Laws with Stiff Source Terms," AIAA Paper 87-11 16, June 1987. 


On the application of ENO scheme with subcell resolution to conservation laws with stiff source terms

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server