Applied Mathematics
Diagram Design and Analysis of
Aerodynamics Problems
with Mathematica
Victor G. Ganzha
1, Evgeni V. Vorozhtsov
21 Institute of Informatics, Technical University of Munich, Munich 80290,
Arcis-str. 21, Germany e-mail:[email protected]
2 Institute of Theoretical and Applied Mechanics, Russian Academy of Sciences,
Novosibirsk 630090, Russia e-mail:[email protected]
Received: August 24, 2000
Summary.
The development of new programming paradigms like the object-oriented programming (OOP) and parallel programming opens up the new prospects for the development of reliable and ef-cient program packages on the basis of those general-purpose com-puter algebra systems (CASs), which support these paradigms. We propose to use a diagram representation of the objects and methods while developing new big programpackages for the solution of aerody-namics problems within the framework of Mathematica by using the OOP principles. We show at a number of computational examples that CAS Mathematica enables one to implement all the essential objects of such a goal-oriented package.Key words:
Euler equations, dierence schemes, curvilinear grids, stabilityMathematics Subject Classication (1991): 65M10, 68C20
1. Introduction
The increase in the eciency of using the computers at the numerical simulation of various uid ows is the central problem of computa-tional uid dynamics (CFD). There are several ways to increase this eciency.
The oldest way is the increase in the computational speed of the computer codes. One of the ways to achieve this goal is the paralleliza-tion of codes written in such languages of numerical computaparalleliza-tions as FORTRAN or C.
The present-day procedure of computer-aided aerodynamic design includes the use of a hierarchy of mathematical models for the sim-ulation of uid ows on computers 4]. These models range from the relatively simple full-potential equation to the complete compress-ible Navier{Stokes equations. In this connection, there emerges the second direction in the solution of the problem of increasing the ef-ciency of computer simulations in aerodynamic design: the use of the object-oriented programming (OOP) at the development of big program packages.
A rapid progress of general-purpose computer algebra systems (CASs) puts into foreground the idea of using these systems as a tool for the solution of CFD problems. It should be noted that there were already several attempts at developing the program packages for the construction and investigation of numerical methods for solving the partial dierential equations (PDEs) of the mathematical physics in the CAS environment. In 12], the REDUCE 3.3 based program package FIDE was presented. This package was designed to auto-mate the process of numerical solving of PDE systems by means of computer algebra. This package enables the user to study both ap-proximation and stability of nite dierence schemes and to generate a FORTRAN code, which implements the chosen nite dierence method. It should, however, be noted that the CAS REDUCE 3.3 does not support explicitly the object-oriented programming.
In 17], it was proposed to generate the needed FORTRAN codes for the solution of elliptic PDEs with the aid of a MACSYMA pro-gram. The ideas discussed in 17] were developed signicantly in 18], where the viewpoint has been expressed that the development of problem solving environments (PSEs) is needed. The SciNapse PSE is described in 18] for numerically solving partial dierential equa-tions. The main interface to the system is the \problem specication language". This language allows an initial-boundary value problem for a system of PDEs to be specied in terms of invariant dier-ential operators or in a particular coordinate system. The methods for discretizing the dierential equations and the boundary condi-tions can be specied in SciNapse by giving keywords for common methods or by giving detailed mathematical replacement rules, or some combination of both. The SciNapse system is implemented in Mathematica and the template are executable Mathematica code, so
they can easily be tested for correctness in Mathematica. SciNapse is a knowledge-based program synthesis system implemented with objects, transformation rules, and a reasoning system.
The purpose of the present paper is to discuss the problems associ-ated with the implementation of the OOP ideas and principles within the framework of CASMathematica 4.022] at the development of a program package for the solution of CFD problems.
Although Mathematica is written in an object-oriented version of C, this CAS does not provide any native support for object-oriented methods. We do note that the packageClasses.mby Roman Maeder
does implement a full-blown object-oriented extension to Mathemati-ca13]. This package does not make it possible to do any calculations that could not be done before, however, it makes it possible to com-pletely rearrange the way in which they are carried out.
We propose to use a diagram representation of the objects (PDEs, discretization methods, grid generation methods, visualization tech-niques) while developing new big program packages for the solution of multidimensional aerodynamics problems. Such a methodology is intended for its OOP implementation. We show at a number of ex-amples that all the important building blocks of a package designed for the numerical modeling of CFD problems can be implemented successfully withMathematica.
2. Hierarchies of models and methods
The governing equations of uid ows are nonlinear partial dieren-tial equations. The analytical solution of these equations is possible only for very particular cases of ow problems. This is the main rea-son for which the numerical methods for the solution of the governing equations of uid ows are the basic means for obtaining the solutions of practically important aerodynamics problems.
At present, all practically important CFD problems are solved numerically with the use of curvilinear spatial grids, which enable one to consider the geometrically complex curved boundaries in the most ecient way.
The process of the development of a computer program for the nu-merical solution of a specic aerodynamics problem on a curvilinear grid consists of the following three main stages:
1. The development of a program for the numerical generation of a curvilinear spatial grid
2. The development of a programimplementinga chosen numerical method for the solution of the Euler or Navier{Stokes equations on a given grid
3. The development of a program for the graphical visualization of the computational results, the computation of the needed integral functionals (the lift force coecients, the drag coecients, etc.).
Simplest uid ow models Full-potential equation Euler equations 2D Euler 3D Euler Navier{Stokes equations Parabolized Navier{Stokes equations Complete compressible Navier{Stokes equations ? ?
Fig.1. The hierarchy of uid models
Since a hierarchy of mathematicalmodels of uid ows can be used at dierent stages of computer-aided aerodynamic design, one can arrange the available mathematical models into a hierarchic diagram representation as shown in Fig. 1. For the sake of brevity, we do not show in Fig. 1 a number of further mathematical models, such as small-disturbance potential equation, the models of viscous-inviscid interaction, the Prandtl boundary layer equations, etc. 4].
At the stage of numerical discretization of the governing equations of uid ow, one can also use a number of the subclasses of numerical methods (Fig. 2). These numerical methods can be implemented on various classes of spatial grids 11] (Fig. 3).
The numerical methods can be implemented in many ways: with the use of the Runge{Kutta type discretizations, implicit methods, multigrid methods, preconditioners, etc.
Discretization methods Finite dierence methods Finite volume methods Finite element methods Approximation investigation Stability investigation ; ; @ @ R ? ? @ @ R ; ; ?
Fig.2. The subclasses of discretization methods Spatialgrids Logically rectangular grids Triangular grids Single-blo ck grids Multi-blo ck grids Solution-adaptive grids Structured grids Unstruct-ured grids ? @ @ @ R ? ) P P P P P q ; ; @ @ R
Fig.3. Classes of spatial grids
As soon as the numerical solution process is completed, the pas-sage to the visualization stage is performed (Fig. 4).
Since the basic goal of the overall mathematical process as de-picted in Figs. 1{4 is a rapid obtaining of the problem solution, the CAS like Mathematica plays an important role, because it has all the needed components for analytic and numeric computations and
Isolines Solutionproles functionalsIntegral Visualization ? H H H j
computer graphics. Therefore, one can implement each stage shown in Figs. 1{4 within the same CAS.
We will show in the subsequent sections how the dierent stages of the computer-aided solution of CFD problems can be implemented withMathematica 4.0.
3. Analytic tests
One of the important constituent parts of the development of a CFD program package is the debugging and validation of new comput-er codes. The validation problem is often solved by comparing the numerical solutions with some particular analytic solutions. The self-similar solutions can often be reduced to ordinary dierential equa-tions (ODEs). The analytic soluequa-tions of many ODE classes can be obtained successfully withMathematica 19{21,1].
In the case of relatively simple test problems involving the par-tial dierenpar-tial equations, the analytic solutions can also be obtained withMathematica 21]. The group methods implemented in CAS en-vironment enable one to obtain the particular solutions of much more complex equations, such as the Navier{Stokes equations 2].
4. Discretization
The construction of the discretizations of the governing partial dif-ferential equations mentioned in Fig. 1 can be implemented by the method of indeterminate coecients with the aid of CAS Mathemat-ica 5].
After a numerical discretization scheme has been constructed, it is necessary to investigate its approximation and stability properties. In the next two subsections, we show how these two tasks can be implemented with CAS Mathematica 4.0.
4.1. Approximation investigation
We consider the Euler equations for two-dimensional inviscid gas ow as the governing equations:
wherex and y are Cartesian coordinates and (2) w= 0 B B @ u v E 1 C C A f(w) = 0 B B @ u u2+p uv uH 1 C C A g(w) = 0 B B @ v vu v2+p vH 1 C C A: Here p, , u, v, E, and H denote the pressure, density, Cartesian velocity components, total energy and total enthalpy. For a perfect gas
(3) E= ( p
;1)+ 12(
u2+v2) H =E+ p
where is the ratio of specic heats.
q q q q q jk;1 j;1k jk j+ 1k jk+ 1 1 2 3 4
Fig.5. The centered scheme
Let us take an arbitrary cell of curvilinear grid (Fig. 5). Let the values ofwat the cell center be denoted bywjk, and let;jk andVjk
be the cell contour and the control volume bounded by contour;jk.
The Euler equations (1) can be written in integral form for the region
Vjk with boundary ;jk as (4) @t@ ZZ Vjkwdxdy+ I ;jk(f dy ;g dx) = 0:
In the result of discretization of (4) we can obtain the semi-discrete equation 9]
whereAjk is the cell area, and the operatorQrepresents an
approx-imation to the boundary integral dened by the second term of (4). For example, the ux balance for the x momentum component is represented in (5) as (6) @t@(Ajku) + 4 X k=1 (Qkuk+ykpk) = 0
where the ux velocity
(7) Qk =ykuk;xkvk
and the sum in (6) is over four sides of the cell, see Fig. 5. The values
xk andyk in (6) and (7) are the increments ofxandy along side k of the cell, with appropriate signs. Each quantity in (6) and (7) such as u2 or (u)2 is evaluated as the average of the values in the cells on the two sides of the face,
(8) (u)2= 12(u)jk+ (u)j+1k]:
The scheme (6){(8) reduces to a central dierence scheme on a Cartesian grid, therefore, it is second order accurate in space. As is known, the second-order schemes produce spurious oscillations of the numerical solution at the shock wave fronts and in their vicinity. In this connection it was proposed in 10] to introduce the articial dissipation terms in the nite volume scheme in order to damp the spurious oscillations. These terms are added to (5) as follows: (9) dwdt + (Qw=Ajk);Dw= 0
whereDis a dissipative operator the structure of which is similar for each of the four dependent variables:
(10) Dwjk=Dxwjk+Dywjk where (11) Dxwjk = (
d
j+1=2k ;d
j ;1=2k) Ajk (12) Dywjk= (d
jk+1=2 ;d
jk ;1=2) Ajk and(13)
d
j+1=2k= Aj+1=2k "(2) j+1=2k(wj +1k ;wjk) ;" (4) j+1=2k(wj +2k ;3wj +1k+ 3wjk ;wj ;1k)] Aj+1=2k = 12(Ajk+Aj +1k):The quantity in (13) is the time step of the numerical method. The coecients "(2) and "(4) are variable, they are adapted to the ow. Dene (14) jk =jpj +1k ;2pjk+pj ;1k j=(jpj +1k j+ 2jpjkj+jpj ;1k j): Then (15) "(2) j+1=2k= (2)max( j j +1k jj jkj) "(4) j+1=2k= max(0 (4) ;" (2) j+1=2k):
Typical values of the constants (2) and (4) are 10] (2) = 1=4,
(4) = 1=256. The values
d
jk1=2 in (12) are computed similarly to (13){(15) for example,Ajk+1=2 = (1=2)(Ajk+Ajk
+1).
The time stepping scheme approximating (9) is constructed as an explicit three{stage scheme of the Runge{Kutta type (see 9]): (16) w(0)=wn w(1)=w(0) ; 1 (Phw (0) ;Dw (0)) w(2)=w(0) ; 2 (Phw (1) ;Dw (0)) w(3)=w(0) ; (Phw (2) ;Dw (0)) where wn=w(xytn), tn =tn ;1 + n,n= 12:::, t0 = 0, n=
Phw = (1=Ajk)Qwjk, and 1, 2 are the nondimensional weight parameters.
A symbolic algorithm for the local approximation study of the nite volume approximation of @f(w)=@x and @g(w)=@y in accor-dance with (6) was presented in detail in 5]. This algorithm was implemented with Mathematica. In particular, the expansion of the grid functions entering the nite volume scheme may be performed by using the following function 5]:
Fx , y ] := Fc + Fx*(x - xc) + Fy*(y - yc)
+ 1/2*Fxx*(x - xc)2 + Fxy*(x - xc)*(y - yc) + 1/2*Fyy*(y-yc)2 + 1/6*Fxxx*(x - xc)3 + 1/2*Fxxy*(x - xc)*(y - yc)
Here Fx , y ] =f(w(xy)), (xc, yc) is the user specied center of the Taylor series expansion,
Fx =@f(w(xy))=@x Fxxy =@3f(w(xy))=@2x@y
etc. In the particular case of an uniform rectangular spatial grid the Taylor expansion gives the second order of approximation. But in the case of a general logically rectangular curvilinear grid, the proximation order reduces to the rst order. The corresponding ap-proximation error is presented in 5] for the apap-proximation of the
x-derivative we will not reproduce it here because of its bulky form. The expression for this error involves the derivatives @2f(w)=@x2,
@2f(w)=@x@y, and @2f(w)=@y2.
4.2. Stability investigation
Let us enumerate a number of items that explain the practical im-portance of the stability investigation of the numerical discretization schemes for PDEs. Knowledge of the stability region of a numerical method enables one
{
to avoid an increasing amplitude of error,{
to reduce the computer time while modeling nonstationary phe-nomena,{
to accelerate the convergence to a stationary solution when using the pseudounsteady method,{
to distinguish between physical and numerical oscillations at the onset of numerical instability (e.g., in turbulence modeling),{
to nd in a given class of numerical methods the specic nitedierence or nite volume schemes, which possess a maximum stability robustness 6].
If the stability condition of a specic numerical discretization scheme is known, it is usually incorporated into the computer code implementing the scheme. But many current nite dierence or nite volume schemes are so complicated that it is impossible to obtain the needed stability conditions by using the \manual" calculations. In such cases, the computer algebra systems prove to be the only means for obtaining the stability regions by combining the symbolic and numeric computations.
In what follows we briey outline the general scheme of a typical symbolic-numeric algorithm for stability investigation of nite dier-ence or nite volume schemes basing on the Fourier method 6]. We
demonstrate this procedure at the example of the stability investiga-tion of the nite volume scheme (6){(16).
If we \freeze" the components of the solution vector w in the dierence scheme coecients then we obtain the linearized dierence scheme
(17) wn+1=S(xyw
0)w
n
where w0 is the \frozen" solution vector. The dierence operator
S(xyw0) in (17) is called the step operator, and its structure deter-mines to a large extent the stability properties of dierence method. In what follows we investigate the stability of the dierence initial-value problem (18) wn+1=S(xyw 0)w n n= 012::: w0=
U
0(xy) ;1< xy <1where
U
0(xy) is the given initial condition att = 0. We apply the von Neumann stability analysis procedure 6] for obtaining the neces-sary stability condition of dierence problem (18). In this procedure, at rst the Fourier symbol G = F(S) of the step operator S is ob-tained. After that the characteristic equation(19) det(I;G) =
m
X
j=0
cjm;j = 0
is computed, where m = 4 in the case of the Euler system of equa-tions (1){(2), andI is themmidentity matrix. The von Neumann necessary stability conditions then read
(20) jj1 +O( ) = 1:::m:
These conditions are also sucient for stability of the dierence Cauchy problem (18) in the case of a normal or diagonalizable am-plication matrixG14].
Let us now describe some peculiarities of the symbolic stages of the symbolic-numeric algorithm, which we have implemented with Mathematica 3.0. A detailed presentation of this method for the case of the absence of articial dissipation terms may be found in 7].
Since the amplication matrixG is a function of the Jacobi ma-trices A and B for many nite dierence or nite volume schemes for Euler equations, there emerges the idea of working withAand B
matrices at the early stages of symbolic computations. In the process of such computations the expressions of the form
(21) c1AB+c2BA c1ABB +c2BAB
etc. arise. It is well known 6] that the gas dynamical matricesAand
B do not commute, therefore, the conventional commutative multi-plication is inapplicable.Mathematicatreats howeverAandB in (21) as commutative quantities and transforms them into (c1+c2)AB and (c1+c2)AB
2, which leads to incorrect results in our case.
In order to solve this problem we have proposed the idea of replac-ing the commutative multiplication operationTimesinMathematica
with the noncommutative operationDot. To illustrate the
implemen-tation of this idea in a Mathematicaprogram let us assume that cer-tain quantitiesA1andA2do not commute and consider the following
program:
c1 = HoldForm(A1+A2) (A1+A2) (A1+A2)] d1 = c1/.fPlus;>List, Times;>Timg
d2 = d1/. HoldFormx y ]];>HoldFormx, y] e = d2/. HoldForm;>Outer
f = e/. List;> Plus
pu3= f/. Timx , y , z ];>Dotx, y, z]
This program computes the product (A1+A2) (A1+A2) (A1+A2) in the correct form, which takes into account the noncommutativity of the quantitiesA1 andA2:
A1 . A1 . A1 + A1 . A1 . A2 + A1 . A2 . A1 +A1 . A2 . A2 + A2 . A1 . A1 + A2 . A1 . A2 + A2 . A2 . A1 + A2 . A2 . A2 The above program fragment is a key element of our Mathematica
program, which has enabled us to obtain the amplication matrixB
symbolically as a function of the matrices A and B, although these matrices themselves were not specied explicitly as matrices. Let us now describe the sequential stages of our symbolic algorithm whose output is the amplication matrixGof a given dierence scheme.
Step 1. Elimination of the intermediate values
u
(1) andu
(2)from the dierence equationu
n+1;
u
(3)= 0 of scheme (16). Let us intro-duce in theMathematica program the notations
uv1, n , i , j ] =
u
(1)ij ,
uv2, n , i , j ] =
u
(2)ij
uv1, n , i , j ] := uv0n, i, j] - alf1 dt (pu - d) uv0, n, i, j]
uv2, n , i , j ] := uv0n, i, j] - alf2 dt (pu uv1, n, i, j] - d uv0,n,i,j]) uv3, n , i , j ] := uv0n, i, j] - dt (pu uv2, n, i, j] - d uv0, n, i, j]) where pu stands for the matrix dierence operator Ph d is the
lin-earized operator of articial dissipation according to (10){(15) dt =
, alf1 = 1, alf2 = 2 ID is assumed to be the 4
4 identity matrix it will be specied explicitly as a matrix at a much later stage. We note that the quantity uv0n, i, j] =
u
nij is not declared here as a vec-tor of four components. We have denoted the result of the action of the dierence operator Ph on the grid functionu
nij simply as a usualproduct: Ph
u
nij = pu uv0n, i, j]. As we have shown in 5], suchdier-ence operators can be eciently implemented in Mathematica with the aid of transformation rules. The elimination of the values
u
(1) andu
(2) is performed in our program with the aid of the command sch = ID uv3, n, i, j]/.fuv0, n, i , j ] ;>uv0n, i, j]g//Expand and the left-hand side of the two-level dierence equationu
n+1;
u
(3)= 0 is obtained in the formID uv0n, i, j] - dt ID pu uv0n, i, j] + alf2 dt2 ID pu2 uv0n, i, j] - alf1 alf2 dt3 ID pu3 uv0n, i, j] - alf1 alf2 dt3 pu2 d
+ alf2 dt2 pu d - dt d
Rewriting this expression by using the original notation, we obtain that the step operatorS of the linearized nite volume scheme cor-responding to three-stage Runge{Kutta scheme (16) has the form (22) S =I; Ph+ 2( Ph) 2 ; 1 2( Ph) 3 ; 1 2 3P2 hD+ 2 2P hD; D
where Ph and D are the linear dierence operators obtained as a
result of the linearization of operatorsPh and Din (16).
Step 2. Substitution of the expression for the dierence operator
Ph into the dierence scheme in integral steps.
c1 = HoldForm(A1+A2) (A1+A2)]
d1 = c1/. fPlus ;>List, Times;>Timg d2 = d1/. HoldFormx , y ]];>HoldFormx, y] e = d2/. HoldForm;>Outer
f = e/. List;> Plus
pu2= f/. Timx , y ];>Dotx, y]
c1 = HoldForm(A1+A2) (A1+A2) (A1+A2)] d1 = c1/. fPlus ;>List, Times;>Timg d2 = d1/. HoldFormx , y ]];>HoldFormx, y]
e = d2/. HoldForm;>Outer f = e/. List;> Plus
pu3= f/. Timx , y , z ];>Dotx, y, z] sch1 = sch/.fpu ;>A1 + A2, pu 2 ;> pu2, pu 3 ;>pu3g A1 = t1 A A2 = t2 B sch12 = sch1/.fdt ID x;>dt x, dtk ID x ;>dtk xg sch2 = sch12/.falf1;>1/2, alf2;>1/2g sch3 = sch2/.f(x1 y1 ).(x2 y2 ).(x3 y3 );>x1.x2.x3 y1 y2 y3g sch4 = sch3/.f(x1 y1 ).(x2 y2);>x1.x2 y1 y2, (x1 y1).(x2 y2 ).(x3 y3 );>x1.x2.x3 y1 y2 y3g
As a result we obtain the right-hand side of (22) in the form involving various matrix products A.A, A.B, A.B.A, etc.
Step 3. Fourier transform of the step operator S(xyw0). For this purpose, the Fourier particular solution wnjk = W0
nei(j 1 +k 2 ), i=p ;1 W
0 = const, is substituted into the expression for the step operator the variables 1 and 2 represent the real wavenumbers.
Since S is a linear dierence operator, we can at rst compute the Fourier transform of each item entering S in order to save the computer memory. As a result of these symbolic computations, the amplication matrix G obtained by the Fourier transformation of operatorS (22) has the form
(23) G=I+ (d1A+d2B) + 2(d1A+d2B) 2+ 1 2(d1A+d2B) 3 + 1 2 3d 3(d1A+d2B) 2+ 2 2d 3(d1A+d2B) +d3I whereF(; Ph) =d 1A+d2B, F(; D) =d 3I, (24) d1= ; 2Ajk(yj+1k ;yjk)(cos 2 ;isin 2) + (yj+1k+1 ;yj +1k)(cos1+isin1) + (yjk+1 ;yj +1k+1)(cos2+isin2) + (yjk;yjk +1)(cos1 ;isin 1)] d2= ; 2Ajk(xjk;xj +1k)(cos2 ;isin 2) + (xj+1k ;xj +1k+1)(cos1+isin1) + (xj+1k+1 ;xjk +1)(cos2+isin2) + (xjk+1 ;xjk)(cos 1 ;isin 1)] d3="2x(2cos1 ;2);4" 4x(cos1 ;1) 2 +"2y(2cos2 ;2);4" 4y(cos2 ;1) 2:
Here A andB are the Jacobi matrices obtained from (2): (25) A= @f@w B(w) = @g@w :(w)
The coecients"2x"4x"2y"4y in (24) are obtained by freezing the coecients "(2) j1=2k," (4) j1=2k" (2) jk1=2" (4)
jk1=2 in the formulas of arti-cial dissipator: (26) "2x=" (2) j1=2k = const " 4x=" (4) j1=2k = const "2y =" (2) jk1=2 = const " 4y =" (4) jk1=2 = const:
Step 4. Symbolic computation of the entries of G. To save the needed memory at this stage, we at rst substitute the entries of the matricesAandBasa21=a21,b43=b43, etc., where we have denoted by ajk andbjkjk= 1:::4, the entries of the matricesA and B,
respectively. After that we nd the expressions for the rst row ofG
by replacingajk with their specic entries corresponding to the gas
dynamic Jacobi matricesA and B (the expressions for these entries may be found in 6]).
Thus, the rst row elements are the functions of the components of the solution vectorwin (2), the nondimensional parameters (26), and the dimensional quantitiesAjk,xjk+1
;xjk,yjk +1 ;yjk,xj +1k ;xjk, yj+1k
;yjk. In order to make the results of stability analysis inde-pendent both of a specic curvilinear grid and of the time step it is extremely important to express the entries of the amplication matrix
G in terms of the nondimensional similarity parameters 1::: M (M 1). It turns out that in the case of the Runge{Kutta nite vol-ume method dened by equations (10){(16) it is sucient to intro-duce along with the parameters (26) the following six nondimensional variables: (27) 1 = c(yjk +1 ;yjk) Ajk 3 = v(yjk +1 ;yjk) Ajk 5 = xj+1k ;xjk yjk+1 ;yjk 2 = u(yjk +1 ;yjk) Ajk 4 = yj+1k ;yjk yjk+1 ;yjk 6 = xjk+1 ;xjk yjk+1 ;yjk
wherecis the sound speed,c=p
p=for a perfect gas. The entries of matrixG can be eciently expressed in terms of the parameters
1::: 6 in aMathematica code by using the transformation rules. Let, for example,cp1= 1. Then we can introduce the notation cp1 into the entrygg1,1]]of Gwith the aid of transformation rule
gg1,1]] = gg1,1]]/. c ;>cp1*Ajk/(yj,k+1]] - yj,k]]), where Ajk = Ajk, yj,k]] = yjk, etc. The remaining parameters 2::: 6 are introduced in the entries of G in a similar way. The rst row of Gthus obtained is then stored in the le row1.m. After
that we compute in a similar way the next rows ofG. Such a strategy saves a lot of computer memory, so that only about 2 Mb are needed to compute one row of G. Once we have computed all the rows of matrixG, we can assemble them into a 44 matrix.
It turns out that G obtained in this way still depends explicitly on the time step . More precisely,G has the following structure:
(28) G= 0 B B B B B B B @ f11 f12 q f 13 q f 14 q 2 f21 q f22 1f 23 1 f 24 q f31 q f32 1f 33 1 f 34 q f41 q2 f42 q f43 q f44 1 1 C C C C C C C A whereq= (yjk+1
;yjk)=Ajkand thefij,ij= 1:::4 are functions of only. Similarly to 6] it is easy to show that there exists a one-parameter family of the diagonal matrices
Q= diag(qq+1q+1q+2) such that the matrix
(29) G0 =QGQ
;1
does already not depend on q. A simple choice for Q is obtained at = 0: Q = diag(1qqq2). The symbolic computation of G
0 in accordance with (29) can conveniently be implemented with Mathe-matica: G0=DiagonalMatrix f1qqq 2 g]GInverse DiagonalMatrix f1qqq 2 g]] As a result we obtain a matrixG0with nondimensional entries, which can be obtained formally from (28) by settingq= 1. In what follows we will omit the subscript 0 and writeGinstead of G0 for brevity.
In this way we have obtained the amplication matrix on a curvi-linear grid for Jameson's scheme (16), which takes 1898 lines of text, with 65 symbols in each line on the average.
At the numerical stages of our method, the zeroes of characteristic equation (19) are computed numerically with the aid of the Mathe-matica function Solve...]. It is well known 6] that these zeroes
are very sensitive to roundo errors when the machine arithmetic of oating-point numbers is used to compute the numerical values of coecients of (19). On the other hand, it is well known that the computer algebra systemMathematicaperforms exact arithmetic op-erations on rational numbers. In accordance with (24) the coecients
cj in (19) depend on cosm and sinm, m = 12. In order to avoid
the introduction of any roundo errors when computing the numeri-cal values of these functions we have determined these values as the following rational numbers:
(30) cosm = 8 < : ;1 m= 1;R 2(t m") 1 +R2(t m") m6= sinm = 8 < : 0 m= 2R(tm") 1 +R2(t m") m6= where R(tm") =RationalizeNtm 10;(e+1) ], "], m = 12 and tm = tan(m=2) " = 10;(e+1), e
0, is the user-specied accura-cy with which the built-inMathematica functionRationalize...]
converts a oating-point number tm into a rational number. It is
important that the calculation of cosm and sinm by (30) always
ensures the satisfaction of the relation cos2
m+ sin2
m= 1.
The values of nondimensional parameters (26) and (27) were al-so computed in (19) as the rational numbers. As a result, the co-ecients cj of equation (19) are exact for any complex numerical
method. Therefore, we can compute the zeroes of (19) exactly, be-cause the Mathematica function Solve...] implements the exact
solution formulas atm= 4.
For the verication of the above symbolic/numeric method we have used the analytic solution 6] for the necessary stability condi-tion of the linearized Runge{Kutta scheme, which can be obtained from (23) in the particular case whered3 = 0, i.e.,"2x="4x ="2y =
"4y = 0: (31) j 2 j+j 2 4 j+j 3 6 j+j 3 5 j + (j 1 j+j 1 4 j) 2+ ( j 1 6 j+j 1 5 j) 2]0:5 f( 1 2)
where (32) f( 1 2) = 1 1 2 1 ; 2+ ( 2( 2 ;4 1+ 8 2 1)) 0:5 2 2 0:5 :
In the particular case of a uniform rectangular spatial grid in the (xy) plane we have:
(33) xj+1k =xjk+h1 xjk+1 =xjk
yj+1k =yjk yjk+1 =yjk+h2
8jk
where h1 and h2 are the grid steps along the x- and y-axes, respec-tively. The substitution of (33) in (27) yields the formulas
1 = c h1 2 = u h1 3 = v h1 4 = 0 5 = h1 h2 6 = 0: Therefore, the stability condition (31) simplies at 1 = 2 = 1=2 to
(34) j 2 j+j 3 j 5+ 1 q 1 + 2 5 2: At the specied cell aspect ratio 5 the equality
(35) 1= (2 ;j 2 j;j 3 j 5)= q 1 + 2 5
determines a surface 1 ='( 2 3) of the stability region of scheme (16) on an uniform rectangular grid in the absence of articial dissi-pator. In particular, at the apex of pyramid (35), where 2= 3 = 0, we have 1 = 2=
p 1 + 2
5. In Table 1 we present the results of nu-merical computation of 1at the stability region boundary at a user-specied accuracy "= 10;2 for the computation of
1. In this table,
1num is the value of 1 on the stability region boundary obtained by the above presented symbolic-numerical method 1ex= 2=
p 1 + 2 5 1 = j 1ex ; 1num
j. It can be seen from the Table 1 that the absolute error 1 is smaller by an order of magnitude than the user-specied accuracy 10;2.
In the general case, where the articial dissipation coecients
"2x"2y,"4x,"4y are dierent from zero, we have at rst determined the stabilityregion boundary in dierent sections of a ten-dimensional Euclidean space of ("2x"2y,"4x "4y 1 2::: 6) points. This has required many runs on a computer. As a result of these runs, we have accumulated the information, which proved to be sucient for an an-alytic tting of all the numerical data obtained. At this tting we have considered the values 1 = 2 =
1
2 in (32), at which f(
1 2) = 2. In this case, we have obtained the von Neumann necessary stability
Table 1. Values of1 at point 2 =3 = 0 of the stability region
boundary of scheme (16) at dierent values of 5 1 n 5 1 =2 1 2 3 1num 1831 1024 5793 4096 7325 8192 971 1536 1:788086 1:41414307 0:894165 0:632162 1ex 1.788854 1.414214 0.894427 0.632456 1 0.000768 0.000093 0.000262 0.000294
condition of the three-stage Runge{Kutta nite volume method (16) on curvilinear grid in the form
(36) " 2x 0:5 2:7 + " 2y 0:5 2:7 + " 4x 1=8 8 + " 4y 1=8 8 + " j 1 j p (1 +j 4 j) 2+ ( j 5 j+j 6 j) 2 2;j 2 j(1 +j 4 j);j 3 j(j 5 j+j 6 j) # 4 1:
It is easy to see that in the particular case where "2x ="2y ="4x =
"4y = 0 formula (36) coincides with inequality (31). In accordance with (23) the amplication matrixGis a matrix polynomial of matrix
d1A+d2B. It is well known that the matricesd1Aandd2Bcan be di-agonalized simultaneously by the same similarity transformation 14]. By using the sucient stability criterion 14] for the diagonalizable amplication matrices it easy to prove that the necessary von Neu-mann stability condition (36) for the linearized Runge{Kutta scheme under consideration is also sucient for stability.
Substituting expressions (27) for the nondimensional similarity parameters 1::: 6 into (36) we can easily obtain the explicit ex-pression for the maximum time step allowed by the stability of the Runge{Kutta scheme under study. A simple analysis of this expres-sion shows that in the case of the presence of articial dissipator, that is in the case where at least one of the four coecients (26) is dierent from zero, the maximum time step size allowed by stability reduces in comparison with the case of the absence of articial dissipation terms.
5. Numerical grid generation
We will demonstrate the lling of one of the classes of spatial grids (Fig.3), logically rectangular single-block grids, with the aid of
Ma-thematica 4.0. Such curvilinear grids have gained a wide acceptance, see 11]. They are mapped with the aid of a non-degenerate transfor-mation
x=x() y=y()
from the plane of the Cartesian spatial coordinates (xy) onto a uni-form rectangular grid in the plane of curvilinear coordinates,.
It should be noted that both the local approximation study and the stability investigation by the von Neumann method are problem independent. The numerical grid generation procedures are, on the contrary, problem dependent and may require, for example, the im-plementation of multi-block grids in the cases of suciently complex geometries.
In this connection, we demonstrate the implementation of numer-ical grid generation at a number of specic aerodynamics problems. The computer algebra system Mathematica proved to be a very con-venient tool for a rapid generation of curvilinear grids. Owing to the powerful built-in computer graphics functions, the user can rapidly adjust the free grid parameters during singleMathematica session in an interactive way in order to obtain a curvilinear grid, which ensures the desired distributions of grid lines near important boundaries or surfaces.
In Fig. 6 we show the curvilinear grids generated with the aid of Mathematica program for problems of inviscid gas ow in a channel with a circular arc bump in lower wall. These grids were generat-ed with the id of the same Mathematica program implementing the multi-surface method 3] for numerical grid generation.
In Fig. 7 we show the curvilinear grids around the NACA 0012 airfoil. In the case of the grid of Fig. 7, (b) we have used the method of 15], in which the C-type grid was generated with the aid of the transformation
(37) x=B+Acosh()cos y =Asinh()sin
where A and B are constants. The coordinate lines = const wrap around the airfoil, whereas the lines = const are approximately normal to the lines= const. We now show a fragment of our Math-ematica program, in which the transformation (37) is implemented: Do eta = efcnss + ds, s, eta, a, thetai]]
xi, j]] = b + a*Cosheta]*Costhetai]] ] yi, j]] = a*Sinheta]*Sinthetai]] ]
0.5 1 1.5 2 2.5 3 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 0.2 0.4 0.6 0.8 1
(a) the grid of 50 15 nodes for
subsonic ow at the freestream Mach numberM
1= 0 :5
(b) curvilinear grid of 8040 nodes
for supersonic ow atM 1= 1
:65 Fig.6. Curvilinear grids for the problem of inviscid ow through a channel with
a circular arc bump in lower wall
2 4 6 8 -4 -2 2 4 2 4 6 8 -4 -2 2 4
(a) the grid obtained by multi-surface
method 3] (b) the grid obtained by the methodof 15]
Fig. 7. The curvilinear C-type grids of 6515 nodes around the NACA 0012
airfoil
6. Visualization
Owing to the availability of numerous powerful and versatile comput-er graphics functions, the Mathematica uscomput-er cn display graphically any one-dimensional plots, isolines, and surfaces obtained as a result of numerical simulation of a uid dynamics problem.
In Fig. 8 we present the computational results obtained on the curvilinear grid of Fig. 6, (b). As a criterion for the numerical solu-tion convergence to a stasolu-tionary limit we have checked the inequality
0.5 1 1.5 2 2.5 3 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 3 x 1.2 1.4 1.6 1.8 2 M
(a) predicted Mach number contours (b) predicted Mach number proles along lower wall (solid line) and up-per wall (dashed line)
0 100 200 300 400 n -5 -4 -3 -2 -1 Log (R )10 n (c) convergence history
Fig. 8. Inviscid supersonic ow through a channel with a circular arc bump in
lower wall at the freestream Mach numberM 1= 1
:65
Rn , where Rn = maxjkjn +1
jk ;njkj is the solution residual,
n = 012:::, and is a user-specied small positive number we have taken the value = 10;4.
Let us denote by max the maximum value of the time step ,
which is obtained from (36) if one replaces the symbol with the equality symbol. The actual aerodynamic computations were per-formed with = max, whereis a safety factor, 0< 1. In the case of Fig. 8, the value= 0:95 was specied.
In order to plot the isolines of a grid function given in the nodes of a curvilinear grid, we have written a Mathematica function that
0.2 0.4 0.6 0.8 1. 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 Cp 0 200 400 600 800n -4 -3 -2 -1 0 1 Log10Rn
(a) the distributions of the pressure coe cient Cp along the upper (|-)
and lower () airfoil surfaces
(b) convergence history 5 10 20 40 60 0 0.0001 0.0002 0.0003 0.0004 5 10
(c) the Mach number contours (d) the surface j(njk;n ;1
jk )=n ;1
jk j, n= 800
Fig.9. The computational results for the problem of transonic ow around the
airfoil NACA 0012 atM 1= 0
:80 and zero angle of attack
implements the algorithm described in 16]. In 8] we discuss the implementation of this algorithm with Mathematica and present a more rapid algorithm for plotting the isolines.
In Fig. 9 we show the numerical results obtained by scheme (16) on the Rizzi mesh of Fig. 7, (b). The maximum residualRndropped
by ve orders of magnitude during the rst 800 time steps: from the valueR1= 6:29
10 1 toR
800 = 3:90 10
results agree well with those obtained in 10] on a 6432 mesh. The stability condition (36) was used in our Mathematica code for this aerodynamics problem to compute the local time steps in accordance with the time stepping procedure proposed in 10].
It can be seen from Fig. 9, (d) that the main source of the numer-ical solution error is caused by the shock waves near the upper and lower airfoil surfaces. These errors decay with the increasing distance from the airfoil. This may be explained by the fact that the shock strength decreases as the distance from the airfoil increases. As a result, the numerical solution gradients within the zone of smeared shock wave become smaller, which leads to the diminution of the residualj(njk;n
;1
jk )= n;1
jk j.
7. Conclusions
We have shown above that the current general-purpose CASs may be used as the basis for the development of big program packages intended for a rapid solution of multi-dimensional uid dynamics problems. The incorporation of the ideas and principles of the object-oriented programming should result in more reliable and exible pro-gram packages. This work needs further big eorts and is now in progress. The purpose of the present paper was to outline a number of the directions toward the solution of the problem of the develop-ment of such a exible OOP based package within the framework of a general-purpose CAS.
A further objective would be the incorporation of some elements of the expert systems and articial intelligence to automatically choose the best methods for numerical grid generation and subsequent nu-merical solution of a specic ow problem.
References
1. G. Baumann (1999): Symmetry Analysis of Dierential Equations Using Mathematica, Springer-Verlag, New York.
2. V.V. Bublik (2000, in press): Group classication of the Navier{Stokes equations for compressible viscous heat-conducting gas, in: Computer Al-gebra in Scientic Computing / CASC'00, V.G.Ganzha, E.W.Mayr and E.V.Vorozhtsov (Eds.), Springer-Verlag, Berlin.
3. P.R. Eiseman (1979): A multi-surfce method of coordinate generation, J. Comput. Phys.33, 118{150.
4. C.A.J. Fletcher (1996):Computational Techniques for Fluid Dynamics, Vols. I, II, 3rd Edition, Springer-Verlag, Berlin, Heidelberg, New York.
5. V.G. Ganzha and E.V. Vorozhtsov (1996): Numerical Solutions for Partial Dierential Equations: Problem Solving Using Mathematica,CRC Press, Boca Raton, New York, London.
6. V.G. Ganzha and E.V. Vorozhtsov (1996): Computer-Aided Analysis of Dierence Schemes for Partial Dierential Equations, Wiley-Interscience, New York.
7. V.G. Ganzha and E.V. Vorozhtsov (1999): Application of computer algebra systems for stability analysi of dierence schemes on curvilinear grids, J. Sym-bolic Computation28, 401{433.
8. V.G. Ganzha, E.V. Vorozhtsov, and M. Wester (2000, in press): An assessment of the e ciency of computer algebra systems in the solution of scientic com-puting problems, in: Computer Algebra in Scientic Comcom-puting / CASC'00, V.G.Ganzha, E.W.Mayr and E.V.Vorozhtsov (Eds.), Springer-Verlag, Berlin. 9. A. Jameson and W. Schmidt (1985): Some recent developments in numerical methods for transonic ows, Computer Methods in Applied Mechanics and Engineering51, 467-493.
10. A. Jameson, W. Schmidt, and E. Turkel (1981): Numerical solution of the Euler equations by nite volume methods Using Runge Kutta time stepping schemes, AIAA Paper 81-1259.
11. P. Knupp and S. Steinberg (1994): Fundamentals of Grid Generation, CRC Press, Boca Raton, Ann Arbor.
12. R. Liska and L. Drska (1990): FIDE: A REDUCE package for automation of nite dierence method for solving PDE, in: Proc. Int. Symp. on Symbolic and Algebraic Computation, ISSAC'90, August 20{24, 1990, Tokyo, Japan, Eds. S. Watanabe and M. Nagata, ACM Press, New York, 169{176.
13. R. Maeder (1991): Programming in Mathematica, 2nd Edition, New York, Addison-Wesley.
14. R.D. Richtmyer and K.W. Morton (1967): Dierence Methods for Initial-Value Problems, 2nd Edition, Wiley-Interscience, New York.
15. A. Rizzi (1981): Computational mesh for transonic airfoils, in: Numerical Methods for the Computation of Inviscid Transonic Flows with Shock Waves, A GAMM Workshop, Eds. A. Rizzi and H. Viviand, Vieweg, Braunschweig, 222{263.
16. P.J. Roache (1976):Computational Fluid Dynamics, Hermosa, Albuquerque. 17. S. Steinberg and P. Roache (1993): Finite-dierence algorithm design and automatic code generation, in: Computer Algebra and Its Applications to Mechanics, CAAM'90. Proc. Int. Seminar, Novosibirsk, August 28{31, 1990 Irkutsk, September 1{3, 1990, Eds. V.G. Ganzha, V.M. Rudenko and E.V. Vo-rozhtsov, Nova Science Publishers, Inc., New York, 15{29.
18. S. Steinberg, R.L. Akers, E. Kant, C.J. Randall, and R.L. Young (1997): SciNapse: A problem-solving environment for partial dierential equations, IEEE Computational Science an Engineering 4, 32{42.
19. W. Strampp and V.G. Ganzha (1995):Dierentialgleichungen mit Mathemat-ica, Braunschweig, Vieweg.
20. W. Strampp, V.G. Ganzha, and E.V. Vorozhtsov (1997): Hohere Mathematik mit Mathematica. Band 3: Dierentialgleichungen und Numerik, Braun-schweig, Vieweg.
21. W. Strampp, V.G. Ganzha, and E.V. Vorozhtsov (1997): Hohere Mathematik mit Mathematica. Band 4: Funktionentheorie, Fourier- und Laplacetransfor-mationen, Braunschweig, Vieweg.
22. S. Wolfram (1999):The Mathematica Book, 4th Edition, Cambridge Univer-sity Press, Cambridge, UK New York.
