An implicit approximate factorization (AF) algorithm is constructed which has the following characteistics. In 2-D: The scheme is unconditionally stable, has a 3 x 3 stencil and at steady state has a fourth order spatial accuracy. The temporal evolution is time accurate either to first or second order through choice of parameter. In 3-D: The scheme has almost the same properties as in 2-D except that it is now only conditionally stable, with the stability condition (the CFL number) being dependent on the cell aspect ratios, delta y/delta x and delta z/delta x. The stencil is still compact and fourth order accuracy at steady state is maintained. Numerical experiments on a 2-D shock-reflection problem show the expected improvement over lower order schemes, not only in accuracy (measured by the L sub 2 error) but also in the dispersion. It is also shown how the same technique is immediately extendable to Runge-Kutta type schemes resulting in improved stability in addition to the enhanced accuracy

Abarbanel, Saul

Kumar, Ajay

English

NASA Technical Reports Server

NASA Contractor Report 181625 
ICASE REPORT NO. @-l& 
ICASE 
COMPACT HIGH ORDER SCHEMES FOR THE EULER EQUATIONS 
(EASA-CR-181625) COXPACT E l G B  CFiCEK SCIIEHES N88-1S138 
ZCE IBE BULBP ECDll'XiCYS P i c a 1  5eEcrt ( N A S A )  
17 P CSCL 09B 
Unclas 
63/61 01282E5 
Saul Abarbanel 
Ajay Kumar 
Contract No. NAS1-18107 
February 1988 
INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AtU) EMGINEERIWG 
NASA Langley Rerearch Center, Hampton, Virginia 23665 
Operated by the Universities Space Research Association 
https://ntrs.nasa.gov/search.jsp?R=19880009754 2020-03-20T08:01:21+00:00Z
COMPACT HIGH ORDER SCHEMES FOR THE EULER EQUATIONS 
Saul Abarbanel 
Tel-Aviv University 
and 
Ajay Kumar 
NASA Langley Research Center 
ABSTRACT 
An implicit approximate factorization (AI?) algorithm is constructed which has the following 
characteristics. 
0 In 2-D: The scheme is unconditionally stable, has a 3 x 3 stencil and at steady state has a 
fourth order spatial accuracy. The temporal evolution is time accurate either to 1st or 2nd 
order through choice of parameter. 
0 In 3-D: The scheme has almost the same properties as in 2-D except that it is now only 
conditionally stable, with the stability condition (the CFL number) being dependent on the 
%ell aspect ratios," Ay/Az and Aa/Az.  The stencil is still compact and fourth order accuracy 
at steady state is maintained. 
Numerical experiments on a 2-D shock-reflection problem show the expected improvement over 
It is also shown how the same technique is immediately extendable to Runge-Kutta type schemes 
lower order schemes, not only in accuracy (measured by the LZ error) but also in the dispersion. 
resulting in improved stability in addition to the enhanced accuracy. 
'This research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS1-18107 while the first author was in residence at the Institute for Computer Applications in Science and 
Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. Partial support was also provided 
under USAF' Grant No. AF'OSR87-0218. 
1 
1 INTRODUCTION 
It can be shown [l] that numerical approximations to the linearized Euler equations of gas dynamics 
give rise to dispersive errors which in the 2-D supersonic case depend on a similarity parameter 
n = (Ay/Az)d- (under the assumption v << u everywhere, where u and v are the z 
and y components of the velocity vector). The difference between the dispersion relations of any 
numerical algorithm and that of the original partial differential equations can be plotted as curves 
in the Fourier plane with n as a parameter. 
In particular the results in [l] indicate that for central-difference schemes, the dispersive error 
are contributed mostly by the third power of the errors in the Fourier variables 0 and 4. It  is, 
therefore, natural to think of fourth order spatially accurate algorithms as having better dispersive 
properties. By utilizing the structure of the Euler equations one can obtain, on a Cartesian grid, a 
fourth order approximation which instead of using a 5 x 5 stencil (and 5 x 5 x 5 in 3-D) relies on 
a compact support of 3 x 3 (and 3 x 3 x 3 in 3-D). The advantages of the combination of fourth 
order accuracy together with compact support are quite obvious in terms of total computer work 
and memory. 
In Section 2, we describe the construction of an approximate factorization (AF) central differ- 
ence scheme and examine its theoretical (linear) stability properties. In Section 3, we derive the 
corresponding 4-step Runge-Kutta scheme and show how the Jameson-Schmidt-Turkel algorithms 
[2] may be easily modified to that form which has, in addition to the higher accuracy, markedly en- 
hanced stability limits. In Section 4, we describe some numerical experiments using the AF-version. 
Section 5 summarizes our findings. 
2 DERIVATION OF THE APPROXIMATE-FACTORIZATION 
SCHEME 
2.1 The Two-Dimensional Case 
Consider a general hyperbolic conservation law in 2-D: 
G + fi +G= 0. 
A - A  
A 
In the case of Euler equations, for example, the vectors ;f, f ( u ) ,  9 ( u )  are given by 
where p, u, v, E ,  and p are respectively the density, velocity components in the z and y directions, 
the total energy per unit volume and the pressure. In addition there is the equation of state relating 
algebraically (in the case of gas-dynamics) the internal energy to the pressure and density. One 
may also write (1) in non-conservation form as 
< +A +B .',= 0 
A A 
where A and B are the Jacobians of f and 9 with respect to . 
1 
(3) 
Consider first a forward-time second order central finite difference approximation to (1) 
where we have dropped the suparrow indicating a vector and used the conventional differencing 
notation 
2 
Uj,k’ t(i”,k = u(jAz, kAy, nAt). ( 5 )  
We have a Cartesian grid with nodes a t  z = jAz ,y  = kAy and t = nAt. We shall now introduce 
the usual shift operations notation: 
Dtu = un+l- j ,  k uj”,k 
6Zu = u j + i , k  - uj-i,k, ‘u = u j , k + i  - u j ,k - l .  
Eq. (4) then may be written as: 
(7) 
6 
B3 Dtu = - A P d z f  - -P&g 
where A = At/Az and lR = Ay/Az. If we take a Taylor series expansion about the grid point 
(z, y, t) = (jAz, kAy, nAt), we can construct the modified equation corresponding to (7): 
or 
Thus if we want to approximate (1) to a higher than second order (and in particular to fourth 
order in space-see comment in Introduction concerning dispersive errors) we must modify (7) by 
subtracting out the terms on the right hand side of (8). At this point we realize that using a 
straightforward approach to approximating fzzz and guy,, will lead to bigger stencils. However, 
using the original partial differential equation, (l), we have 
fizz = -ut22 - gyzz 
gvvv = - w y y  - fzvv. 
(9) 
(10) 
and similarly 
Since fizz and guuU need be approximated only to second order (because of their coefficients, Az2/3 
and Ay2/3) the required stencil for the terms on the right hand side is only 3 x 3. (This observation 
was previously made, in another context by Jones, et al. [3].) So, replacing fizz and guuu in (8) 
by -6:DtulAz’At - 6:p,6,,g/Az2Ay and -6; Dtu/Ay2At - 6y2pZ6, f/Ay2Az, respectively; and 
also replacing ut: by Dtut/At -, Dt(- fi - gu)/At + -Dt [(PA f/Az) + (py6,,g/Ay)] we obtain 
the following approximation to (1) which is spatially fourth order accurate and temporally second 
order accurate: 
Az2 6z2 Ay2 6y2 
f + y v g ]  
2 
Using (3) and the definition of Dt we rearrange the above into the delta form 
where I is the identity matrix. This is an implicit unfactored scheme. Since we are interested in 
marching to steady state we approximatefactor the left hand side to obtain: 
We introduced into (12) the parameter u. If u = 1 then we retain the temporal second order 
accuracy while u = 2 gives first order in time. Note that (12) involves the inversion of block- 
tridiagonal matrices. The right hand side represents the steady state operator to fourth-order. The 
whole scheme involves only a 3 stencil. 
2.2 The 3-D Case 
The starting point, corresponding to (1) is 
A A + f z  + s^y + h,= 0. 
Following the same steps as in the 2-D case, we get the modified equation 
where using (13), we have 
fizz = -Utzz - Qyzz - hzzz 
Qvvv = -utvu - hzvv - f zvv  
hzzz = - w a r  - f zrz  - Q y r r .  
Repeating all the steps leading to (12) we obtain its three dimensional analog: 
In (18) the matrix C is the Jacobian a 5; /a 
manner analogous to (6). Q is the cell aspect ratio in the z direction, Q = A z / A x .  
and the shift operators pr,iir are defined in a 
3 
2.3 
We consider the linearized @e., frozen coefficients) version of (12) .  We carry out a Von-Neumann 
analysis in the usual manner by casting ( 1 2 )  in Fourier space. A typical Fourier component of uzk 
Stability Analysis for the 2-D Case 
where 
with 
-T < +,e < T and -00 < t , m  < 00. 
The mapping of the various shift operations is as follows: 
1 
We also define the amplification matrix G by fin+’ = Gir”. With these notations, the frozen 
coefficient version of (12) is mapped into the Fourier space as follows: 
( I  - :AS2 + i u A A e d z )  ( I  - 5Bq2 2 + 
= -2ix ( I  - f.2) + ;qJ1-.  ( I  - f e y  -
In the general case the matrices A and B do not commute, thus rendering the analysis of (26) 
almost intractable. It is instructive, however, to consider the scalar case. Since the aspect ratio 
R = Ay/Az is arbitrary, and since the original partial differential equation ( l ) ,  in the scalar linear 
case, could be transformed to the wave equation 
uz + UT + up = 0 
with 3 = At, z = x and p = Ay/B,  we can without loss of generality (in this special linear scalar 
case) rewrite (26) as: 
Equation (27) can be rearranged to solve for the amplification factor G : 
a + i ( i  - l)b 
a + i?jb G =  
4 
where 
and b is the Fourier map of the steady-state operator, Le.: 
(30) b = 2 [x(dG(l- p2) 2 + Z d G ( 1 -  at2)] .
We see immediately from (28)  that for u = 1 (i.e., the case of second order temporal accuracy) we 
have a Crank-Nicholson type scheme with IlGll = 1. For u > 1 (Le., first order accuracy in time), 
we have IlGll < 1. We have thus demonstrated the unconditional stability of (27) for all values of 
the cell aspect ratio. 
2.4 Stability for the 3-D Case 
The three dimensional analog to (27) is 
[ 4 - 7  1 - Z q 2  - Z c 2 )  + &/=(I - I e 2  - - 2 € 2 )  + $c$c-F(l - $2 2 - 3$2)] 2 
3 Rz 3 =-2 iX  ( 1-(( 
(31) 
where 
The stability of three dimensional amplification factor G, as defined by (31), is difficult to 
analyze and we resorted to numerical evaluations of [GI2, using up to 8 x log Fourier modes. The 
numerical study of (31) was carried out on the Cray 2.  We found that as in the case of forward 
Euler approximate factorization second order scheme [4], the amplification factor was conditionally 
stable. For example for R = Q = 1, the stability limit is X 5 .43. These stability properties are 
obtained with the aid of artificial viscosity (AV) term which is of sixth order but still resides on 
the compact stencil. The AV term is added to the right hand side (explicit term) of (18) and is of 
the form 
is the Fourier dual variable defined analogously to ( and q .  
(32) 
1 
P,.xg6,26,26,2 
where pa,, is of order unity. We found 1.5 < pa,, < 2 to be most efficacious. Without the artificial 
viscosity term, the allowed value of X is about one order of magnitude less (e.g., for R = Q = l,X 
without using AV is about .035). 
3 COMPACT FOUR STEP FOURTH ORDER RUNGE-KUTTA 
ALGORITHM 
The basic four step explicit Runge-Kutta scheme, as proposed in [2], takes the form: 
uo = u" 
(33) 
where R is the finite difference (or finite volume) representation of the steady operator. For example, 
in the 2-D second order spatial accuracy case 
where f and g are the flux vectors encountered in section 2. If we want fourth order spatial accuracy, 
then it follows directly that the residual R is modified to read 
and we still retain the compact support. 
Similarly, in the three dimensional case we have 
Thus (33), with R given either by (35) for the 2-D case or by (36) for the 3-D case, retains all the 
features of the second order scheme but gains us the fourth order accuracy. In addition one can 
easily verify by simple analysis that for a given Cartesian grid and flow conditions the new fourth 
order formulation enhances the stability condition. In the 2-D case we have, using (35) rather than 
= 1.36. (37) 
(AtI(4) a (~ t )4 th  order 
( a t ) ( ~ )  (At)2mi order 
(34) 
-- 
In the 3-D case the gain is even more favorable, 
Thus the algorithm efficiency gains are two fold. First, for a given acceptable error level the fourth 
order accuracy allows a coarser grid, i.e., fewer node points. Second, not only At is increased due 
to the larger cell size but in addition it gains due to (37) (or (38) in the 3-D case). 
6 
4 NUMERICAL RESULTS FOR 2-D CASE 
The two-dimensional fourth order compact scheme is applied to a shock reflection problem sketched 
in Figure 1. It shows a 5" shock at Mach 1.95 reflecting from a flat plate. Results are presented 
both without any explicit artificial viscosity (AV) and with a fourth order AV term of the type 
-+(6:s,") added to the right hand side. Addition of this AV term reduced the accuracy of the 
compact scheme to third order (note that in 3-D the sixth order AV terms precludes this reduction 
in accuracy). Calculations are also made with a second order implicit Euler AF' scheme without and 
with the preceding AV term. Figures 2a and 2b show the results of the compact scheme whereas 
Figures 3a and 3b show the corresponding results from the second order scheme. It is seen that for 
t = 0, both fourth and second order schemes produce very oscillatory results but with e = 0.36, the 
results from the compact scheme (Figure 2b) improve dramatically and, in fact, are much better 
than the corresponding results obtained from the second order scheme (Figure 3b). The shocks 
captured by the compact scheme are sharper and the convergence is also seen to improve. 
Results from a study of grid aspect ratio effect on the compact scheme are also presented here 
in terms of the similarity parameter IC of reference 1. Three values of IC are considered, namely 1.67, 
1.01, and 0.42. Figure 2b corresponds to IC = 1.67 and figures 4 and 5 correspond to IC = 1.01 and 
0.42, respectively. In all these cases, e is set equal to 0.36. It is clear from the figures displaying 
the effect of IC that the best results are obtained for n near unity. In reference 1, a linear theory 
(e.g., for weak shock) predicts the same results. It is interesting that we find numerically that this 
is also the case in the present nonlinear problem. 
SUMMARY 
1. The steady state solution of the Euler equations of gas dynamics may be achieved to fourth 
order accuracy using a compact grid stencil of 3 x 3 and 3 x 3 x 3 in the 2-D and 3-D cases 
respectively. We presented two examples of such algorithms: one implicit (Euler approximate 
factorizations scheme) and one explicit (Four-stage Runge-Kutta). 
2. Numerical experiments were carried out for the 2-D shock reflection problem, using the im- 
plicit algorithm. Comparisons are made with a corresponding second order scheme. The 
results show that the compact higher order scheme offers marked improvement in both accu- 
racy and convergence rate. 
3. In connection with this work we would like to make the following remarks. It is known that 
the finite difference scheme cannot obtain high order accuracy for conservation-form equations 
computed on a non-uniform grid. This observation coupled with the ease of obtaining fourth 
order compact schemes even in 3-D for the Euler equations revives an old debate: Can one use 
uniform grid and apply conveniently boundary conditions to arbitrary shapes. The potential 
gain in reduced number of computational nodes and enhanced convergence rate appears large 
enough to study this question again. 
7 
References 
[l] R. Shapiro, “Prediction of dispersive errors in numerical solutions of the Euler equations,” sub- 
mitted for presentation at the 2nd International Conference on Hyperbolic Problems, Aachen, 
West Germany, March 1988. 
[2] A. Jameson, W. Schmidt, and E. Turkel, “Numerical solutions of the Euler equations by finite- 
volume methods using Runge-Kutta time-stepping schemes,” AIAA Paper No. 81-1259, June 
1981. 
[3] D. J. Jones, J. C. South, and E. B. Klunker, “On the numerical solution of elliptic partial 
differential equations by the method of lines,” J. Comp. Physics, Vol. 9, No. 3, June 1982, pp. 
496-527. 
[4] R. M. Beam and R. F. Warming, “An implicit factored scheme for the compressible Navier- 
Stokes equations,” AIAA J., Vol. 16, 1978, pp. 393-401. 
y=.5 
y=.32 
y=o 
x=o x=l 
Figure 1: Shock reflection problem 
8 
E- 
O b  
.80 
= 0.25 7 
1.05 1.30 
P 
111111111111) 
1 .s 1.80 
1.05 
t' 
" 
2 -  .25 
0 
.80 
X "  
UII 
0.85 
1 e 3 0  1.55 1.80 
P 
I E -5.50 
-6.75 
h 
Figure 2a: Pressure distribution and residual plot for compact scheme (e  = 0, IC = 1.67). 
9 
'wF x = 0.25 I . W F  x = 0.85 
L 
Figure 2b: Pressure distribution and residual plot for compact scheme (e = 0.36,~ = 1.67). 
10 
1 y = 0.14 
1.55 
L 
1 i i e 1.30 
E I 
L 
0 
.80 
-3 -00 F 
x =  
u 
0.85 
u 
P 
4 . 2 5   
(Y -5.50 
-6 -75 
Figure 3a: Pressure distribution and residual plot for second order scheme (e = 0, n = 1.67). 
11 
x =  
> s .25 
1 . q -  
I- 
k Y = 0.14 
t 
e 1.3Ok 
t 
.80 E .25 .50 .7!5 1 .oo 
0 
X 
.so F 
E 
x = 0.85 
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1  
-80 1 .os 1 .#) 1.55 1 .BO 
P 
-5.50 
L \ t \ 
-6.75 E 
L 
Figure 3b: Pressure distribution and residual plot for second order. scheme (e = 0 . 3 6 , ~  = 1.67). 
12 
i 
I i ,
1 
i 
i 
I  
i 
I 
I 
I 
I 
= 0.25 1 
L 
0 .BO  1.05 1 . 3 0  1 .s 1.80 
P 
y = 0.14 
I 
1 
L l  
1 
L 
.BO I I I I I I  
0 .2 
1 1 h 
.50 
X 
11 
.7 
d 
I .m 
- 3 . O F  
-Y -5  \ 
LT - 6 . O F  
-7.+ 
E l  I 
L 
Figure 4: Pressure distribution and residual plot for compact scheme ( E  = 0.36, K = 1.01). 
13 
'"E x = 0.25 I 
E 
1 1 1 1 1 1 1 1 1 1 1 l 1 1 1 1 1 1 1 1  
.80 1 .a 1.30 1.55 1.80 
P 
1 . e o F  
y = 0.14 E 
1.55 h r 
c I 
1 E x = 0.85 
t 
0 0  1.80 
1 .os 1.30 1.55 
P 
.80 
L 
Figure 5: Pressure distribution and residual plot for compact scheme ( E  = 0.36, n = 0.42). 
14 
Report Documentation Page 
1. Report No. 
NASA CR-181625 
ICASE Report  No. 88-13 
2. Government Accession No. 
7. AuthorM 
Saul  Abarbanel and Ajay Kumar 
7. Key Words (Suggested by Authork)) 
Eule r  equa t ions ,  compact schemes, 
h igh  order  accuracy,  d i s p e r s i o n  
9. Performing Organization Name and Address 
I n s t i t u t e  f o r  Computer Appl ica t ions  i n  Science 
Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665-5225 
and Engineering 
18. Distribution Statement 
61 - Computer Programming and 
64 - Numerical Analysis  
Unc la s s i f i ed  - unlimited 
Software 
2. Sponsoring Agency Name and Address 
National  Aeronaut ics  and Space Adminis t ra t ion 
Langley Research Center  
Hampton, VA 23665-5225 
9. Security Classif. (of this report) 20. Security Classif. (of this page) 21. No. of pages 
Unc 1 ass i f i ed Unc las s i f i ed  16 
3. Recipient's Catalog No. 
22. Price 
A02 
5. Report Date 
February 1988 
6. Performing Organization Code 
8. Performing Organization Report No. 
88-13 
10. Work Unit No. 
505-90-2 1-0 1 
11. Contract or Grant No. 
NAS1-18107 
13. Type of Report and Period Covered 
Contrac tor  Report 
14. Sponsoring hgency Code 
5. Supplementa Notes 
Langley Technica l  Monitor: 
Richard W. Barnwell 
Submitted t o  J. of S c i e n t i f i c  
Computing 
F i n a l  ReDort 
6. Abstract 
An implicit approximate factorization (Al?) algorithm is constructed which has the following 
characteristics. 
0 In 2-D: The scheme is unconditionally stable, has a 3 x 3 stencil and at steady state has a 
fourth order spatial accuracy. The temporal evolution is time accurate either to 1st or 2nd 
order through choice of parameter. 
0 In 3-D: The scheme has almost the same properties as in 2-D except that it is now only 
conditionally stable, with the stability condition (the CFL number) being dependent on the 
"cell aspect ratios," Ay/Az and Aa/Az.  The stencil is still compact and fourth order accuracy 
a t  steady state is maintained. 
NASA FORM 1626 OCT 86 NASA-Langley, 1988 


Compact high order schemes for the Euler equations

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server