- Email: [email protected]
Symbolic-numerical computation of the stability regions for Jameson's schemes
I
Mathematics
ELSEVIER
and Computers
Symbolic-numerical
in Simulation
MATHEMATICS AND COMPUTERS IN SIMULATJON
3
42 (1996) 607615
computation of the stability regions for Jameson’s schemes
V.G. Ganzha *, E.V. Vorozhtsov Institute of Theoretical
and Applied Mechanics,
Novosibirsk
630090,
Russia
Abstract We describe a symbolic-numerical method for the Fourier stability analyses of difference initial-value problems approximating the initial-value problems for hyperbolic or parabolic PDEs. The Fourier method is reduced to the algebra of the resultants. We further use the REDUCE computer algebra system for the symbolic computation of the resultant and for the generation of a FORTRAN function to compute the value of the resultant. Basing on this FORTRAN function we further construct a binary function to characterize the stability and instability points. Using this function we generate a bilevel digital picture, and the stability region boundaries are then detected in this picture with the aid of the efficient algorithm proposed previously by Pavlidis (1982). The above symbolic-numerical method has enabled us to obtain for the first time the stability regions of the considered Jameson’s schemes as applied to the two-dimensional advection-diffusion equation. Analytical formulas are proposed for the approximation of the boundaries of the obtained stability regions. It is shown that these formulas underestimate insignificantly the actual sizes of the stability regions. Therefore, these formulas can efficiently be used in practical computations by Jameson’s schemes. Keywords:
Partial differential
equations;
Difference
schemes; Stability; Symbolic-numerical
computation
1. Introduction The numerical simulation is at present a powerful tool for studying various processes in physics and fluid dynamics. Finite difference schemes and the method of finite volume are now widely accepted in the computational fluid dynamics. The convergence of the solutions of difference equations to the exact solution of the original initial- and boundary-value problems for the partial differential equations in the limit when the steps of the mesh in the spatial variables xl, . . . , XL (L 3 1) and the time 1 tend to zero necessitates the stability and the approximation [ 11. However, for many widespread schemes of the finite difference method and of the finite volume method the stability has not been investigated until now even within the framework of the linear theory. The main reason for this is the complexity of mathematical expressions arising in the analysis of the above types of numerical discretizations of the PDEs by the Fourier method. * Corresponding author. [email protected]. 0378-4754/96/$15.00
Present
address:
GH-Universitat,
D-34127
0 1996 Elsevier Science B.V. All rights reserved
PII SO378-4754(96)00037-7
Kassel,
Germany.
Fax: +561-804-3597;
e-mail:
608
VG. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
The practical solution of a problem on the stability investigation of the numerical approximations to complex multidimensional problems is impossible without the application of computers. A review of the state-of-the-art of the use of computers for the stability analyses of difference initial-value problems and of difference initial- and boundary-value problems may be found in [2,3]. A significant progress in the solution of problems arising in the computer implementation of the Fourier method has been achieved in the work [4]. Owing to the use of symbolic manipulations on a computer with the aid of the computer algebra system REDUCE [5], the packages GENTRAN [6] and SCOPE [7], as well as the Newton identities [4], it has become possible to alleviate substantially the problem on the intermediate expression swell arising in the symbolic realization on computer of the Fourier method of the stability analysis. This enabled the authors of [4] to obtain for the first time necessary stability conditions of the MacCormack schemes [8] for the two- and three-dimensional nonstationary problems of gas dynamics. We analyze in the present work the Von Neumann stability [l] of a number of schemes from the families of Jameson’s schemes [9-l l] as applied to the two-dimensional advectiondiffusion equation. It should be noted that the schemes of Jameson’s families have gained a relatively wide acceptance for the numerical solution of applied problems, see, in particular, the works [ 12-171. This is explained by the fact that in the process of the construction of Jameson’s schemes the following useful design principles have been implemented [lo]: (1) The conservation laws of gas dynamics are satisfied in the discrete conservative form. (2) Nonphysical solutions are eliminated by introducing the appropriate dissipative terms in the discrete approximation. (3) The final steady state does not depend on the time marching procedure. (4) The total enthalpy is constant in the stationary solution. (5) The uniform flow is the exact solution of the difference equations on an arbitrary grid. (6) The schemes are applicable both for stationary and nonstationary problems.
2. Reduction of the Fourier stability analysis to the algebra of the resultants A detailed presentation of this procedure may be found in [3]. Therefore, we restrict ourselves here to a short description of it. Consider the difference initial-value problem of the form u n+l=Hu”,
n=0,1,2
u” = uu(x),
)...,
(1) (2)
wheren = (xl, . . . , XL), L 2 1, xl, . . . , XL, are the spatial variables; u” = U(X, t”), t* = nr, T is the time step, H is a matrix difference operator, uu(x) is the initial-value vector function entering the formulation of an initial-value problem for the original PDE system. In accordance with the Fourier method the solutions of the form u(x, t) = U exp(i(kx - or)} are substituted into scheme (1). Here k = (kl , . . . , kL) is a real wave vector, o is the wave frequency, U is a constant vector and i = 2/-i. As a result one obtains from (1) the formula
V.G. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
609
where B is the amplification matrix of the difference scheme (1). Let
Cm(A)=
&j(K.&P-j,
m3 1,
(3)
j=O
be acharacteristic polynomial of the matrix B, where 6 = (61, . . . , to), ej = kj hj , j = 1, . . . , L, h 1, . . . , hL arethestepsofauniformmeshalongtheaxesxl,. . .,XL.In(3)K = (~1,. . . , KM), M>l,andKl,. . . ,KM are nondimensional complexes in the space of the variation of which the stability region is determined. Now make as in [2- 41 the Mobius transformation h = (w + l)/(w - 1) in (3). Let g(w)
= (w - lYC,((w
+ l)/(w
- 1)).
(4)
The necessary Von Neumann stability condition, as formulated in [ 11, has the form ]hj 1< 1 + O(t), j = 1 . . 1m, where hj, j = 1, . . . , m, are the eigenvalues of the matrix B. Denote by Wj the zeroes of the m, corresponds with the condition j = l,..., polynomial (4), j = 1, . . . , m.ThentheconditionRewj
(5)
Usually the coefficients of polynomial (3) are periodic functions of the spectral variables 61, . . . ,& TL . Consider a parallelepiped periodsTt,..., 17:
(0 < 61 < c, 1 = 1, . . . , l}
with
(6)
in an L-dimensional Euclidean space E ’ of the ,$ points. Now we can construct the following binary function to characterize the stability and instability points (see also [4]): sign R(k, 6) = const Vt E fl, The K-points at which f(k)
= 1 are the points where the difference
(7) scheme is stable.
3. Symbolic stages It follows from (7) that we need the values of the resultant R to determine whether the difference scheme (1) is stable or unstable at a given K-point. In the present version we have implemented the same strategy for symbolically computing the resultant as in [3]. Therefore, we only enumerate the basic steps of the symbolic procedure following [3]. Step 1: Symbolic computation of the amplification matrix B of a difference scheme. Step 2: Computation of the characteristic polynomial (3). Step 3: Computation of the polynomial (4). Step 4: Computation of the symbolic expression for the resultant R(K, 0 in (5) with the aid of a standard function RESULTANT of the REDUCE system [5].
610
KG. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
5: Generation of the FORTRAN function result (K., <) with the aid of the REDUCE system. This function enables us to compute the numerical value of the resultant at a given point (K, [) at the numerical stage of our symbolic-numerical method, which is briefly described in Section 4.
Step
4. An outline of Pavlidis’ algorithm We can generate on the basis of the binary function (7) an M-dimensional bilevel digital picture . . . iM}, where the indices (il, . . . , iM) refer to the centers of cells of a uniform rectangular computing mesh in the given hypercube P of K-points, where we wish to determine the stability region of difference scheme. In what follows we briefly describe a modified Pavlidis’ algorithm for the segmentation of two-dimensional (i.e. planar) bilevel pictures [ 181. Step I: Tracing of the external contours. The start cell of each such contour for the subsequent tracing is found by scanning the digital image along the horizontal rows of cells (pixels) from the left to the right and from the bottom to the top of the image. It is convenient to determine the direction for further contour tracing in terms of the observer who moves along the stability region boundary in such a way that the stable points lie to the right of the observer. Step 2: The obtained boundary cells of external contours are used as the start cells for the scanning inside the external contours to detect the possible unstable holes. This process is repeated recursively until no new contours are detected. Step 3: Let us denote by AKI and AK~ the steps of a uniform rectangular mesh in the (~1, ~2) plane. Then the quantization effects lead to an error of the order O(AKI) + ~(AKz) in determining the locations of the stability region boundaries. In this connection we have augmented Pavlidis’ algorithm by a refinement procedure, in which the bisection process is applied along a normal to the contour at each point of the contour. As a result it is possible to compute the positions of the points of the stability region boundaries at a user-specified accuracy E (for example, E = 10p3). Ifi,
5. Stability analysis of Jameson’s schemes for the advection4iffusion
equation
It was shown in [l l] that Jameson’s schemes [9,10] are optimal with respect to stability in cases when they are applied for the numerical solution of the Euler equations for an inviscid fluid. However, in a number of works (see, in particular, [ 12-171) these schemes were applied also for the numerical integration of the Navier-Stokes equations. In this connection the questions of the stability investigation of Jameson’s schemes as applied to the Navier-Stokes equations are of present interest. We investigate below the stability of two Jameson’s schemes applied to a simpler model equation, the two-dimensional advection+liffusion equation of the form au/at
+ Alh/&q
+ B&.4/&x2 = v(~~u/~x~
+ a2u/L-h,2).
(8)
Here A and B are the constant components of the advection velocity vector along the XI- and x2-axes, respectively; ‘u is the diffusion coefficient, v = const. > 0. In accordance with the lines method for the solution of PDEs we will consider instead of the original PDE the ordinary differential equation du/dt
+ Pu = 0,
(9)
KG. Ganzha, E.K Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
where the operator P is the operator of spatial differencing.
PU =A(ui+l,j Let
US
- ui-l.j)/(2hl)
u(Ui+l,j
- 2Uij +
denote by F(-t
F(-rP)
611
In the case of Eq. (8) it has the form
+ B(ui,j+l - ui,j-l)/(2h2)
ui-l,j)/(h:)
-
v(Ui,j+l
-
2Uij +
Ui,j-l)/(hi).
(10)
P) the Fourier symbol of the operator P in (10). Then we have that Bt sin& + sin62
= -i
(11)
h2
It follows from (11) that the stability regions of Jameson’s schemes for Eq. (8) may be considered in the three-dimensional space of the variables K1 =
At/h],
Bt/hz,
K2 =
K3 =
vs/(h:).
In what follows we will specify the cell aspect ratio three-stage Jameson’s scheme of the form [ 193 u(o) =
Un
u(l) =
ucd) _ t PU(O),
u(2) = $p)
+ $(1)
u(3) = &CO) + $(2) U(n+l) =
~4
=
h 1/ h2. We consider at first a relatively simple
- $spu(l),
(12)
- $rp&),
u(3)
Let 0 = F(--5 P). Then we obtain from (12) the characteristic
equation
where Q(a) = 1 + (J + $r2 + ;03.
(13)
In Fig. 1 we show some sections of the f boundary of the stability region of the three-stage Jameson’s scheme (12), (lo), which were obtained with the aid of the above presented symbolic-numerical method. It can be shown analytically on the basis of (13) that at ~4 = 1 the value of the coordinate ~3 at the point of the intersection of the f boundary with the OKg-axis is equal to K3 =
lr7()1/[4(1 +
K;)]
2
0.31409317,
(14)
- 1 S -2.5127453. In all the computational examples where au = (-4 + fl)‘/3 - (4 + m)‘j3 presented below and obtained by our symbolic-numerical method the absolute error E of the computation of the (~1, ~2, ~3) coordinates on f was specified as & = 10-3. At point ~1 = ~2 = 0 the value of ~3 on f computed by our method was ~3 = 0.3 148. Thus the actual error in determining ~3 on the Kg-axis does not exceed the given E = 1Om3. It follows from Fig. 1 that the height of the f surface over the plane ~3 = 0 depends substantially on the magnitude of the cell aspect ratio ~4 = hl/ h2. This impedes the obtaining of an analytic approximation of the f boundary (for example, by the least-squares method). Therefore, it is desirable to find other
612
KG. Ganzha, E.V Vorozhtsov/Mathematics
Fig. 1. 1: section = 2, K2 = 0.
K4
=
1, K2
=
0; 2: the projection
and Computers in Simulation 42 (1996) 607615
of the section
~4
1, KI
=
=
K2
on
the
plane ~2 = 0; 3: section
K4
nondimensional variables K; , K;, K; with respect to which the boundary f is less sensitive to the variation of the cell aspect ratio ~4. In connection with foregoing we introduce the variable K;
=
“t(l/h:
+
instead of the variable Ki
=
l/h;)
~3.
=
K3(1
+ K;)
(15)
We-obtain from (14) that the relationship
$7ol
(16)
takes place on r at ~1 = ~2 = 0 independently of the magnitude of the cell aspect ratio ~4 = hl/h2. If Fig. 2 we present a number of sections of the f boundary of the stability region of the three-stage Jameson’s scheme (12), (10) which pass through the OK;-axis. It may be seen that the maximum height of the r surface over the plane K; = 0 now remains the same under the variation of the cell aspect ratio within wide limits. The behavior of the r boundary at ~4 E [O.l, 101 shows that the possible approximation of the stability region of the scheme under consideration can be represented by analytic formula
-2
Fig. 2. 1: section ~4 = 0.1, ~2 = 0, 2: section ~4 = 10, ~2 = 0; 6: the arc of the ellipse ~f/3
-1
~4 +
0
1, ~2 = 0; 3: section Ki/b2 = 1.
=
2
1
~4
=
2, k-2 = 0; 4: section K4 = 4, K2 = 0,5: section
VG. Ganzha, E. V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
613
(17) where b = $l~ul = 0.6281863. It follows from Fig. 2 that Eq. (17) describes the f boundary more accurately with increasing ~4. In the plane K; = 0 it yields for any ~4 the exact result 1~1 I + 1~2 I < &. For a quantitative estimation of the error of approximation (17) let us introduce in the plane ~2 = 0 the polar coordinates K$
K1 =?-COSP,
=rsin(p.
(18)
Then the step t can be calculated in the section t < R(@/[(A/W2
~2
0
=
+ u~(Iz,~ + h,2)210.5,
(19)
where r = Z?(q) is the equation of the r boundary in the section K:! = 0 in the polar coordinates. Let iV Then points (K,i, K&), i = 1, . . . , N, be calculated in the section ~2 = 0. Let tan p = K&/Kli. Ri E R(qi) =
(K:i
+
(Kji)2)o.5.
(20)
Substituting (18) into (17) we obtain in the section
(
r=l/
2
where a = A.
a2
b2
=
(21)
’
1
In particular, at (p = vi = arctan(w&/~ti)
r E r(q$) = Ri/[(Kli/a)2
0 the equation of the ellipse in polar coordinates:
0.5
sin2 q
cos+-
~2
•i-
we can obtain from (21) the relationship
(K;i/b)2]o’5
with the aid of (20). We will calculate the mean-square error 6R
=
1 N G x[(Ri 1=l
0.5
. 100%
- ri)/Ri]*
(22)
1
in the section ~2 = 0. Eq. (22) yields, in particular, at ~4 = 1 the value 6R = 20.6975%. This means Eq. (17) underestimates the maximal value of t allowed by the stability by about 20% on the average taking into account (19). Thus that Eq. (17) can be used in engineering computations at least in the range 0.1 < ~4 < 10. The convenience of the practical use of Eq. (17) consists of the fact that owing to the fact that ~1, ~2 and ~4 enter with equal powers the left-hand side of (17) one can easily find from (17) the explicit analytic expression for the maximal step T allowed by the stability. Now consider the optimal five-stage Jameson’s scheme [9-l l]
(23) Jm)
= u(o) -
Un+l =&n)
y1 T
hdm-‘),
614
VG. Ganzha, E.K Vorozhtsov/Mathematics
_ ..
IJ
”
Fig. 3. 1: section ~4
=
1, KI
=
~4 ~2
= 0.1, ~2 = 0; 2: section Kj = 1, K2 = 0; 3: section ~4 = 4, K2 = 0; 4: the projection of the section on the plane ~2 = 0; 5: the arc of the ellipse (~1 /3.3)2 + (~;/0.62)~ = 1; 6: the arc of the ellipse
+ (~;/0.62)~
(2~1 /3.3)2
and Computers in Simulation 42 (1996) 607-615
= 1.
where YI
=l,
y2=;.
y3=g
y4=
&
y5 = ;.
(24)
In Fig. 3 we show several sections of the r boundary of the stability region of Jameson’s scheme (23), (24), where K; is again determined by Eq. (15). In this case the exact value of K; on r at ~1 = ~2 = 0 is K; = 0.6477989. In the numerical computation by our method we have obtained at the point f fl OK; the value K; = 0.6481 at 0.1 f ~4 < 10. The r curve for the case ~4 = 10 differs very little from the curve f for the case ~4 = 4 in the section ~2 = 0. Therefore, the curve for the case ~4 = 10 is not plotted in Fig. 3. It follows from Fig. 3 that the volume of the stability region increases with the reduction of the value ~4; the same phenomenon takes place also in the case of the three-stage scheme (of Fig. 2). By analogy with the scheme (12) let us take in (21) a = 4, b = Iiool = 0.6477989. In this case we obtain that at ~4 > 4 this ellipse will intersect the numerically determined curves of the r boundary in the section ~2 = 0. This means that some points of the ellipse (21) will be outside the stability regions of corresponding schemes. In this connection one can complicate the problem on determining the constants a and b in (21). Let us introduce the functional 0.5
1 N z
C[(Ri i=l
- ri)/Ri
(25)
J2
I
We formulate for this functional the following constrained minimization problem: find
min F(a, b)
under the constraints Ri - pi > 0,
i = l,...,N,
O
0 < b < 0.6477989.
(26)
We have found numerically the solution of the problem (26) with the error 61 = 0.1 in determining the value of a and with the error ~2 = 0.01 in determining the value of b. The answer is a = 3.3, b = 0.62. Thus
KG. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
615
in the range 0.1 < ~4 < 10 one can use in the practical engineering computations the following approximate stability condition of the five-stage Jameson’s scheme: (27) In particular, in the section ~2 = 0 the mean-square error S R computed by Eq. (22) is equal to 11.95 15%. The local error of Eq. (27) in the neighbourhood of the point r fl OK; is (0.6478 - 0.62)/0.6478 2 4.3%. References [I] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial-value Problems, 2nd edition (Interscience, New York, 1967). [2] S.I. Mazurik and E.V. Vorozhtsov, Symbolic-numerical computations in the stability analyses of difference schemes, in: S. Watanabe and M. Nagata, eds., Proc. Internat. Symp. on Symbolic and Algebraic Computation (Tokyo, Japan, 20-24 August 1990) (Addison-Wesley, Reading, MA, 1990) 177-184. [3] V.G. Ganzha, E.V. Vorozhtsov and Chr. Zenger, Stabilitatsuntersuchung von Differenzenverfahren mit Hilfe der Resultantenalgebra unter Verwendung des REDUCE-Systems, Preprint, Institut fur Informatik, Technische Universitat Miinchen TUM- 19 120, Munich, 199 1. [4] V.G. Ganzha, E.V. Vorozhtsov and J.A. van Hulzen, A new symbolic-numeric approach to stability analysis of difference schemes, in: P.S. Wang, ed., Proc. Intemat. Symp. on Symbolic and Algebraic Computation (ISSAC’92) (Berkeley, California, 27-29 July 1992) (ACM Press, New York, 1992) 9-15. [5] A.C. Heam, REDUCE Users Manual, Version 3.3, The Rand Corporation, Santa Monica, 1987. [6] B.L. Gates, GENTRAN User’s Manual, REDUCE Version, The RAND Corporation, Santa Monica, 1987. [7] J.A. van Hulzen, B.J.A. Hulshof, B.L. Gates and M.C. van Heerwaarden, A code optimization package for REDUCE, in: G.H. Gonnet, ed., Proc. ISSAC’ (ACM Press, New York, 1989) 163-170. [8] R.W. MacCormack, The effect of viscosity in hypervelocity impact cratering, AIAA Paper, No. 354, 1969. [9] A. Jameson, W. Schmidt and E. Turkel, Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes, AIAA Paper, No. 1259, 198 1. [lo] A. Jameson and W. Schmidt, Some recent developments in numerical methods for transonic flows, Comput. Methods Appl. Mech. Engrg. 5 1 (l-3) (1985) 467493. [l l] J. Pike and P.L. Roe, Accelerated convergence of Jameson’s finite-volume Euler scheme using van der Houwen integrators, Comput. Fluids 13 (2) (1985) 223-236. [ 121 V.N. Vatsa, Accurate numerical solutions for transonic flow over finite wings, J. Aircraft 24 (6) (1987) 377-385. [ 131 L. Martinelli, Validation of a multigrid method for the Reynolds-averaged equations, AIAA paper No. 88 04 14, 1988. [ 141 D.G. Holmes and S.D. Connell, Solution of the 2D Navier-Stokes equations on unstructured adaptive grids, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-15 June 1989); A Collection of Technical Papers, (AIAA, Washington, DC 20024, 1989) 25-39. [ 151 L.N. Long, M.M.S. Khan and H.T. Sharp, A massively parallel three-dimensional Euler/Navier-Stokes method, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-15 June 1989); A Collection of Technical Papers (AIAA, Washington, DC 20024, 1989) 89-102. [ 161 R.K. Agarwal, Development of a Navier-Stokes code on a connection machine, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-15 June 1989); A Collection of Technical Papers (AIAA, Washington, DC 20024, 1989) 103-108. [17] F. Grass0 and C.G. Speziale, Supersonic flow computations by twwquation turbulence modeling, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-l 5 June 1989); A Collection of Technical Papers (AIAA, Washington, DC 20024, 1989) 232-239. [ 181 T. Pavlidis, Algorithms for Graphics and Image Processing (Springer, Berlin, 1982). [ 191 C.W. Shu and S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes, J. Comput. Phys. 77 (2) (1988) 43947 I.
Mathematics
ELSEVIER
and Computers
Symbolic-numerical
in Simulation
MATHEMATICS AND COMPUTERS IN SIMULATJON
3
42 (1996) 607615
computation of the stability regions for Jameson’s schemes
V.G. Ganzha *, E.V. Vorozhtsov Institute of Theoretical
and Applied Mechanics,
Novosibirsk
630090,
Russia
Abstract We describe a symbolic-numerical method for the Fourier stability analyses of difference initial-value problems approximating the initial-value problems for hyperbolic or parabolic PDEs. The Fourier method is reduced to the algebra of the resultants. We further use the REDUCE computer algebra system for the symbolic computation of the resultant and for the generation of a FORTRAN function to compute the value of the resultant. Basing on this FORTRAN function we further construct a binary function to characterize the stability and instability points. Using this function we generate a bilevel digital picture, and the stability region boundaries are then detected in this picture with the aid of the efficient algorithm proposed previously by Pavlidis (1982). The above symbolic-numerical method has enabled us to obtain for the first time the stability regions of the considered Jameson’s schemes as applied to the two-dimensional advection-diffusion equation. Analytical formulas are proposed for the approximation of the boundaries of the obtained stability regions. It is shown that these formulas underestimate insignificantly the actual sizes of the stability regions. Therefore, these formulas can efficiently be used in practical computations by Jameson’s schemes. Keywords:
Partial differential
equations;
Difference
schemes; Stability; Symbolic-numerical
computation
1. Introduction The numerical simulation is at present a powerful tool for studying various processes in physics and fluid dynamics. Finite difference schemes and the method of finite volume are now widely accepted in the computational fluid dynamics. The convergence of the solutions of difference equations to the exact solution of the original initial- and boundary-value problems for the partial differential equations in the limit when the steps of the mesh in the spatial variables xl, . . . , XL (L 3 1) and the time 1 tend to zero necessitates the stability and the approximation [ 11. However, for many widespread schemes of the finite difference method and of the finite volume method the stability has not been investigated until now even within the framework of the linear theory. The main reason for this is the complexity of mathematical expressions arising in the analysis of the above types of numerical discretizations of the PDEs by the Fourier method. * Corresponding author. [email protected]. 0378-4754/96/$15.00
Present
address:
GH-Universitat,
D-34127
0 1996 Elsevier Science B.V. All rights reserved
PII SO378-4754(96)00037-7
Kassel,
Germany.
Fax: +561-804-3597;
e-mail:
608
VG. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
The practical solution of a problem on the stability investigation of the numerical approximations to complex multidimensional problems is impossible without the application of computers. A review of the state-of-the-art of the use of computers for the stability analyses of difference initial-value problems and of difference initial- and boundary-value problems may be found in [2,3]. A significant progress in the solution of problems arising in the computer implementation of the Fourier method has been achieved in the work [4]. Owing to the use of symbolic manipulations on a computer with the aid of the computer algebra system REDUCE [5], the packages GENTRAN [6] and SCOPE [7], as well as the Newton identities [4], it has become possible to alleviate substantially the problem on the intermediate expression swell arising in the symbolic realization on computer of the Fourier method of the stability analysis. This enabled the authors of [4] to obtain for the first time necessary stability conditions of the MacCormack schemes [8] for the two- and three-dimensional nonstationary problems of gas dynamics. We analyze in the present work the Von Neumann stability [l] of a number of schemes from the families of Jameson’s schemes [9-l l] as applied to the two-dimensional advectiondiffusion equation. It should be noted that the schemes of Jameson’s families have gained a relatively wide acceptance for the numerical solution of applied problems, see, in particular, the works [ 12-171. This is explained by the fact that in the process of the construction of Jameson’s schemes the following useful design principles have been implemented [lo]: (1) The conservation laws of gas dynamics are satisfied in the discrete conservative form. (2) Nonphysical solutions are eliminated by introducing the appropriate dissipative terms in the discrete approximation. (3) The final steady state does not depend on the time marching procedure. (4) The total enthalpy is constant in the stationary solution. (5) The uniform flow is the exact solution of the difference equations on an arbitrary grid. (6) The schemes are applicable both for stationary and nonstationary problems.
2. Reduction of the Fourier stability analysis to the algebra of the resultants A detailed presentation of this procedure may be found in [3]. Therefore, we restrict ourselves here to a short description of it. Consider the difference initial-value problem of the form u n+l=Hu”,
n=0,1,2
u” = uu(x),
)...,
(1) (2)
wheren = (xl, . . . , XL), L 2 1, xl, . . . , XL, are the spatial variables; u” = U(X, t”), t* = nr, T is the time step, H is a matrix difference operator, uu(x) is the initial-value vector function entering the formulation of an initial-value problem for the original PDE system. In accordance with the Fourier method the solutions of the form u(x, t) = U exp(i(kx - or)} are substituted into scheme (1). Here k = (kl , . . . , kL) is a real wave vector, o is the wave frequency, U is a constant vector and i = 2/-i. As a result one obtains from (1) the formula
V.G. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
609
where B is the amplification matrix of the difference scheme (1). Let
Cm(A)=
&j(K.&P-j,
m3 1,
(3)
j=O
be acharacteristic polynomial of the matrix B, where 6 = (61, . . . , to), ej = kj hj , j = 1, . . . , L, h 1, . . . , hL arethestepsofauniformmeshalongtheaxesxl,. . .,XL.In(3)K = (~1,. . . , KM), M>l,andKl,. . . ,KM are nondimensional complexes in the space of the variation of which the stability region is determined. Now make as in [2- 41 the Mobius transformation h = (w + l)/(w - 1) in (3). Let g(w)
= (w - lYC,((w
+ l)/(w
- 1)).
(4)
The necessary Von Neumann stability condition, as formulated in [ 11, has the form ]hj 1< 1 + O(t), j = 1 . . 1m, where hj, j = 1, . . . , m, are the eigenvalues of the matrix B. Denote by Wj the zeroes of the m, corresponds with the condition j = l,..., polynomial (4), j = 1, . . . , m.ThentheconditionRewj
(5)
Usually the coefficients of polynomial (3) are periodic functions of the spectral variables 61, . . . ,& TL . Consider a parallelepiped periodsTt,..., 17:
(0 < 61 < c, 1 = 1, . . . , l}
with
(6)
in an L-dimensional Euclidean space E ’ of the ,$ points. Now we can construct the following binary function to characterize the stability and instability points (see also [4]): sign R(k, 6) = const Vt E fl, The K-points at which f(k)
= 1 are the points where the difference
(7) scheme is stable.
3. Symbolic stages It follows from (7) that we need the values of the resultant R to determine whether the difference scheme (1) is stable or unstable at a given K-point. In the present version we have implemented the same strategy for symbolically computing the resultant as in [3]. Therefore, we only enumerate the basic steps of the symbolic procedure following [3]. Step 1: Symbolic computation of the amplification matrix B of a difference scheme. Step 2: Computation of the characteristic polynomial (3). Step 3: Computation of the polynomial (4). Step 4: Computation of the symbolic expression for the resultant R(K, 0 in (5) with the aid of a standard function RESULTANT of the REDUCE system [5].
610
KG. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
5: Generation of the FORTRAN function result (K., <) with the aid of the REDUCE system. This function enables us to compute the numerical value of the resultant at a given point (K, [) at the numerical stage of our symbolic-numerical method, which is briefly described in Section 4.
Step
4. An outline of Pavlidis’ algorithm We can generate on the basis of the binary function (7) an M-dimensional bilevel digital picture . . . iM}, where the indices (il, . . . , iM) refer to the centers of cells of a uniform rectangular computing mesh in the given hypercube P of K-points, where we wish to determine the stability region of difference scheme. In what follows we briefly describe a modified Pavlidis’ algorithm for the segmentation of two-dimensional (i.e. planar) bilevel pictures [ 181. Step I: Tracing of the external contours. The start cell of each such contour for the subsequent tracing is found by scanning the digital image along the horizontal rows of cells (pixels) from the left to the right and from the bottom to the top of the image. It is convenient to determine the direction for further contour tracing in terms of the observer who moves along the stability region boundary in such a way that the stable points lie to the right of the observer. Step 2: The obtained boundary cells of external contours are used as the start cells for the scanning inside the external contours to detect the possible unstable holes. This process is repeated recursively until no new contours are detected. Step 3: Let us denote by AKI and AK~ the steps of a uniform rectangular mesh in the (~1, ~2) plane. Then the quantization effects lead to an error of the order O(AKI) + ~(AKz) in determining the locations of the stability region boundaries. In this connection we have augmented Pavlidis’ algorithm by a refinement procedure, in which the bisection process is applied along a normal to the contour at each point of the contour. As a result it is possible to compute the positions of the points of the stability region boundaries at a user-specified accuracy E (for example, E = 10p3). Ifi,
5. Stability analysis of Jameson’s schemes for the advection4iffusion
equation
It was shown in [l l] that Jameson’s schemes [9,10] are optimal with respect to stability in cases when they are applied for the numerical solution of the Euler equations for an inviscid fluid. However, in a number of works (see, in particular, [ 12-171) these schemes were applied also for the numerical integration of the Navier-Stokes equations. In this connection the questions of the stability investigation of Jameson’s schemes as applied to the Navier-Stokes equations are of present interest. We investigate below the stability of two Jameson’s schemes applied to a simpler model equation, the two-dimensional advection+liffusion equation of the form au/at
+ Alh/&q
+ B&.4/&x2 = v(~~u/~x~
+ a2u/L-h,2).
(8)
Here A and B are the constant components of the advection velocity vector along the XI- and x2-axes, respectively; ‘u is the diffusion coefficient, v = const. > 0. In accordance with the lines method for the solution of PDEs we will consider instead of the original PDE the ordinary differential equation du/dt
+ Pu = 0,
(9)
KG. Ganzha, E.K Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
where the operator P is the operator of spatial differencing.
PU =A(ui+l,j Let
US
- ui-l.j)/(2hl)
u(Ui+l,j
- 2Uij +
denote by F(-t
F(-rP)
611
In the case of Eq. (8) it has the form
+ B(ui,j+l - ui,j-l)/(2h2)
ui-l,j)/(h:)
-
v(Ui,j+l
-
2Uij +
Ui,j-l)/(hi).
(10)
P) the Fourier symbol of the operator P in (10). Then we have that Bt sin& + sin62
= -i
(11)
h2
It follows from (11) that the stability regions of Jameson’s schemes for Eq. (8) may be considered in the three-dimensional space of the variables K1 =
At/h],
Bt/hz,
K2 =
K3 =
vs/(h:).
In what follows we will specify the cell aspect ratio three-stage Jameson’s scheme of the form [ 193 u(o) =
Un
u(l) =
ucd) _ t PU(O),
u(2) = $p)
+ $(1)
u(3) = &CO) + $(2) U(n+l) =
~4
=
h 1/ h2. We consider at first a relatively simple
- $spu(l),
(12)
- $rp&),
u(3)
Let 0 = F(--5 P). Then we obtain from (12) the characteristic
equation
where Q(a) = 1 + (J + $r2 + ;03.
(13)
In Fig. 1 we show some sections of the f boundary of the stability region of the three-stage Jameson’s scheme (12), (lo), which were obtained with the aid of the above presented symbolic-numerical method. It can be shown analytically on the basis of (13) that at ~4 = 1 the value of the coordinate ~3 at the point of the intersection of the f boundary with the OKg-axis is equal to K3 =
lr7()1/[4(1 +
K;)]
2
0.31409317,
(14)
- 1 S -2.5127453. In all the computational examples where au = (-4 + fl)‘/3 - (4 + m)‘j3 presented below and obtained by our symbolic-numerical method the absolute error E of the computation of the (~1, ~2, ~3) coordinates on f was specified as & = 10-3. At point ~1 = ~2 = 0 the value of ~3 on f computed by our method was ~3 = 0.3 148. Thus the actual error in determining ~3 on the Kg-axis does not exceed the given E = 1Om3. It follows from Fig. 1 that the height of the f surface over the plane ~3 = 0 depends substantially on the magnitude of the cell aspect ratio ~4 = hl/ h2. This impedes the obtaining of an analytic approximation of the f boundary (for example, by the least-squares method). Therefore, it is desirable to find other
612
KG. Ganzha, E.V Vorozhtsov/Mathematics
Fig. 1. 1: section = 2, K2 = 0.
K4
=
1, K2
=
0; 2: the projection
and Computers in Simulation 42 (1996) 607615
of the section
~4
1, KI
=
=
K2
on
the
plane ~2 = 0; 3: section
K4
nondimensional variables K; , K;, K; with respect to which the boundary f is less sensitive to the variation of the cell aspect ratio ~4. In connection with foregoing we introduce the variable K;
=
“t(l/h:
+
instead of the variable Ki
=
l/h;)
~3.
=
K3(1
+ K;)
(15)
We-obtain from (14) that the relationship
$7ol
(16)
takes place on r at ~1 = ~2 = 0 independently of the magnitude of the cell aspect ratio ~4 = hl/h2. If Fig. 2 we present a number of sections of the f boundary of the stability region of the three-stage Jameson’s scheme (12), (10) which pass through the OK;-axis. It may be seen that the maximum height of the r surface over the plane K; = 0 now remains the same under the variation of the cell aspect ratio within wide limits. The behavior of the r boundary at ~4 E [O.l, 101 shows that the possible approximation of the stability region of the scheme under consideration can be represented by analytic formula
-2
Fig. 2. 1: section ~4 = 0.1, ~2 = 0, 2: section ~4 = 10, ~2 = 0; 6: the arc of the ellipse ~f/3
-1
~4 +
0
1, ~2 = 0; 3: section Ki/b2 = 1.
=
2
1
~4
=
2, k-2 = 0; 4: section K4 = 4, K2 = 0,5: section
VG. Ganzha, E. V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
613
(17) where b = $l~ul = 0.6281863. It follows from Fig. 2 that Eq. (17) describes the f boundary more accurately with increasing ~4. In the plane K; = 0 it yields for any ~4 the exact result 1~1 I + 1~2 I < &. For a quantitative estimation of the error of approximation (17) let us introduce in the plane ~2 = 0 the polar coordinates K$
K1 =?-COSP,
=rsin(p.
(18)
Then the step t can be calculated in the section t < R(@/[(A/W2
~2
0
=
+ u~(Iz,~ + h,2)210.5,
(19)
where r = Z?(q) is the equation of the r boundary in the section K:! = 0 in the polar coordinates. Let iV Then points (K,i, K&), i = 1, . . . , N, be calculated in the section ~2 = 0. Let tan p = K&/Kli. Ri E R(qi) =
(K:i
+
(Kji)2)o.5.
(20)
Substituting (18) into (17) we obtain in the section
(
r=l/
2
where a = A.
a2
b2
=
(21)
’
1
In particular, at (p = vi = arctan(w&/~ti)
r E r(q$) = Ri/[(Kli/a)2
0 the equation of the ellipse in polar coordinates:
0.5
sin2 q
cos+-
~2
•i-
we can obtain from (21) the relationship
(K;i/b)2]o’5
with the aid of (20). We will calculate the mean-square error 6R
=
1 N G x[(Ri 1=l
0.5
. 100%
- ri)/Ri]*
(22)
1
in the section ~2 = 0. Eq. (22) yields, in particular, at ~4 = 1 the value 6R = 20.6975%. This means Eq. (17) underestimates the maximal value of t allowed by the stability by about 20% on the average taking into account (19). Thus that Eq. (17) can be used in engineering computations at least in the range 0.1 < ~4 < 10. The convenience of the practical use of Eq. (17) consists of the fact that owing to the fact that ~1, ~2 and ~4 enter with equal powers the left-hand side of (17) one can easily find from (17) the explicit analytic expression for the maximal step T allowed by the stability. Now consider the optimal five-stage Jameson’s scheme [9-l l]
(23) Jm)
= u(o) -
Un+l =&n)
y1 T
hdm-‘),
614
VG. Ganzha, E.K Vorozhtsov/Mathematics
_ ..
IJ
”
Fig. 3. 1: section ~4
=
1, KI
=
~4 ~2
= 0.1, ~2 = 0; 2: section Kj = 1, K2 = 0; 3: section ~4 = 4, K2 = 0; 4: the projection of the section on the plane ~2 = 0; 5: the arc of the ellipse (~1 /3.3)2 + (~;/0.62)~ = 1; 6: the arc of the ellipse
+ (~;/0.62)~
(2~1 /3.3)2
and Computers in Simulation 42 (1996) 607-615
= 1.
where YI
=l,
y2=;.
y3=g
y4=
&
y5 = ;.
(24)
In Fig. 3 we show several sections of the r boundary of the stability region of Jameson’s scheme (23), (24), where K; is again determined by Eq. (15). In this case the exact value of K; on r at ~1 = ~2 = 0 is K; = 0.6477989. In the numerical computation by our method we have obtained at the point f fl OK; the value K; = 0.6481 at 0.1 f ~4 < 10. The r curve for the case ~4 = 10 differs very little from the curve f for the case ~4 = 4 in the section ~2 = 0. Therefore, the curve for the case ~4 = 10 is not plotted in Fig. 3. It follows from Fig. 3 that the volume of the stability region increases with the reduction of the value ~4; the same phenomenon takes place also in the case of the three-stage scheme (of Fig. 2). By analogy with the scheme (12) let us take in (21) a = 4, b = Iiool = 0.6477989. In this case we obtain that at ~4 > 4 this ellipse will intersect the numerically determined curves of the r boundary in the section ~2 = 0. This means that some points of the ellipse (21) will be outside the stability regions of corresponding schemes. In this connection one can complicate the problem on determining the constants a and b in (21). Let us introduce the functional 0.5
1 N z
C[(Ri i=l
- ri)/Ri
(25)
J2
I
We formulate for this functional the following constrained minimization problem: find
min F(a, b)
under the constraints Ri - pi > 0,
i = l,...,N,
O
0 < b < 0.6477989.
(26)
We have found numerically the solution of the problem (26) with the error 61 = 0.1 in determining the value of a and with the error ~2 = 0.01 in determining the value of b. The answer is a = 3.3, b = 0.62. Thus
KG. Ganzha, E.V Vorozhtsov/Mathematics
and Computers in Simulation 42 (1996) 607-615
615
in the range 0.1 < ~4 < 10 one can use in the practical engineering computations the following approximate stability condition of the five-stage Jameson’s scheme: (27) In particular, in the section ~2 = 0 the mean-square error S R computed by Eq. (22) is equal to 11.95 15%. The local error of Eq. (27) in the neighbourhood of the point r fl OK; is (0.6478 - 0.62)/0.6478 2 4.3%. References [I] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial-value Problems, 2nd edition (Interscience, New York, 1967). [2] S.I. Mazurik and E.V. Vorozhtsov, Symbolic-numerical computations in the stability analyses of difference schemes, in: S. Watanabe and M. Nagata, eds., Proc. Internat. Symp. on Symbolic and Algebraic Computation (Tokyo, Japan, 20-24 August 1990) (Addison-Wesley, Reading, MA, 1990) 177-184. [3] V.G. Ganzha, E.V. Vorozhtsov and Chr. Zenger, Stabilitatsuntersuchung von Differenzenverfahren mit Hilfe der Resultantenalgebra unter Verwendung des REDUCE-Systems, Preprint, Institut fur Informatik, Technische Universitat Miinchen TUM- 19 120, Munich, 199 1. [4] V.G. Ganzha, E.V. Vorozhtsov and J.A. van Hulzen, A new symbolic-numeric approach to stability analysis of difference schemes, in: P.S. Wang, ed., Proc. Intemat. Symp. on Symbolic and Algebraic Computation (ISSAC’92) (Berkeley, California, 27-29 July 1992) (ACM Press, New York, 1992) 9-15. [5] A.C. Heam, REDUCE Users Manual, Version 3.3, The Rand Corporation, Santa Monica, 1987. [6] B.L. Gates, GENTRAN User’s Manual, REDUCE Version, The RAND Corporation, Santa Monica, 1987. [7] J.A. van Hulzen, B.J.A. Hulshof, B.L. Gates and M.C. van Heerwaarden, A code optimization package for REDUCE, in: G.H. Gonnet, ed., Proc. ISSAC’ (ACM Press, New York, 1989) 163-170. [8] R.W. MacCormack, The effect of viscosity in hypervelocity impact cratering, AIAA Paper, No. 354, 1969. [9] A. Jameson, W. Schmidt and E. Turkel, Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes, AIAA Paper, No. 1259, 198 1. [lo] A. Jameson and W. Schmidt, Some recent developments in numerical methods for transonic flows, Comput. Methods Appl. Mech. Engrg. 5 1 (l-3) (1985) 467493. [l l] J. Pike and P.L. Roe, Accelerated convergence of Jameson’s finite-volume Euler scheme using van der Houwen integrators, Comput. Fluids 13 (2) (1985) 223-236. [ 121 V.N. Vatsa, Accurate numerical solutions for transonic flow over finite wings, J. Aircraft 24 (6) (1987) 377-385. [ 131 L. Martinelli, Validation of a multigrid method for the Reynolds-averaged equations, AIAA paper No. 88 04 14, 1988. [ 141 D.G. Holmes and S.D. Connell, Solution of the 2D Navier-Stokes equations on unstructured adaptive grids, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-15 June 1989); A Collection of Technical Papers, (AIAA, Washington, DC 20024, 1989) 25-39. [ 151 L.N. Long, M.M.S. Khan and H.T. Sharp, A massively parallel three-dimensional Euler/Navier-Stokes method, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-15 June 1989); A Collection of Technical Papers (AIAA, Washington, DC 20024, 1989) 89-102. [ 161 R.K. Agarwal, Development of a Navier-Stokes code on a connection machine, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-15 June 1989); A Collection of Technical Papers (AIAA, Washington, DC 20024, 1989) 103-108. [17] F. Grass0 and C.G. Speziale, Supersonic flow computations by twwquation turbulence modeling, in: Proc. AIAA 9th Computational Fluid Dynamics Conf. (Buffalo, NY, 13-l 5 June 1989); A Collection of Technical Papers (AIAA, Washington, DC 20024, 1989) 232-239. [ 181 T. Pavlidis, Algorithms for Graphics and Image Processing (Springer, Berlin, 1982). [ 191 C.W. Shu and S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes, J. Comput. Phys. 77 (2) (1988) 43947 I.