Symbolic-Numeric Stability Investigation of Difference Schemes for the Euler and

Navier-Stokes Equations

V.G. Ganzha(*) and E.V. Vorozhtsov

Institute of Theoretical and Applied Mechanics, Russian Academy of Sciences, Novosibirsk 630090, Russia

(*) At the moment the University of Kassel, Kassel D 34127, Germany

Abstract

We present a number of the symbolic-numeric algorithms for the stability investigation of difference schemes approximating the 2D and 3D Euler and Navier-Stokes equations of the compressible fluids. All the algorithms implement in different ways the von Neumann stability analysis procedure. It is shown that the consideration of the curvilinear grid within the framework of the Fourier stability analysis method increases by a factor of 23/9 the number of nondimensional variables in the domain of the variation of which the stability region is then determined in comparison with the case of an uniform rectangular spatial grid for the difference discretization of the 3D thin-layer Navier-Stokes equations. The questions of the compression of intermediate algebraic expressions are discussed. A four-step stability analysis procedure for difference schemes on curvilinear grids is proposed. The influence of turbulence modeling on the sizes of the necessary stability region has been investigated. Computational examples are presented.

1 Introduction

The numerical modeling of various gas flows on the basis of the finite difference and finite volume discretizations of the Euler and compressible Navier-Stokes equations has now become a powerful tool of the solution of many real-life problems. The numerical stability of many current finite difference schemes has still not been investigated. An important reason is the complicated structure of the mathematical expressions that are required for an adequate formulation of such a stability investigation of difference schemes based on the Fourier method. In this connection, we have previously developed a number of the symbolic-numerical methods for the stability investigation of difference initialvalue problems. These are:

1^{o}. The method implementing the modified Routh algorithm for studying the location of the characteristic equation zeroes with respect to a unit disc in a complex plane [1].

2^{o}. The method using the Routh-Hurwitz determinants [2].

3^{o}. The method using the algebra of the resultants together with the Fourier method [3].

4^{o}. The probabilistic method for the stability analysis of both linear and nonlinear difference schemes [4].

5^{o}. The method based on a direct high-precision numerical solution of the characteristic equation [5].

6^{o}. The method using the algebra of the resultants and the catastrophe theory to study the manifold of the characteristic equation zeroes [6].

Until now, only the methods 1^{o}, 3^{o} and 6^{o} have been adjusted and used for the stability investigation of difference schemes approximating the 2D and 3D Euler and Navier-Stokes equations. One of the purposes of the present paper is to discuss the problems, which arise at the symbolic and numerical stages of these methods, and to present the ways to circumwent these problems. Our second objective is to discuss those additional difficulties, which arise in the case of the stability investigation of difference schemes on arbitrary curvilinear grids of quadrilateral cells for the Euler and Navier-Stokes equations.

2 Scheme Checker

The finite-difference or finite-volume schemes for real-life problems are often characterized by very big lengths of the corresponding difference equations. In such cases there arises the problem of the check-up of the correctness of the difference schemes, which are input into the program for the symbolic/numeric stability investigation or for the local approximation study of the difference scheme.

We have developed a four-stage symbolic analysis module ”SCHEME CHEC-KER” for this purpose. It is applicable to the finite-difference or finite-volume schemes approximating the Euler or Navier-Stokes equations of a compressible fluid. The scheme checker should be used before proceeding to the stability or approximation study and it prevents from making the occasional errors in the input difference equations.

The Euler equations for the two-dimensional nonstationary flow of an inviscid non-heat-conducting compressible fluid may be written as

∂u/∂t + ∂F(u)/∂x + ∂G(u)/∂y = 0,  (1)

where

u = ( \begin{matrix}ρ \\ ρu \\ ρv \\ ρE \end{matrix} ) , F(u) = ( \begin{matrix}ρu \\ p + ρu^{2} \\ ρuv \\ pu + ρuE \end{matrix} ) , G(u) = ( \begin{matrix}ρv \\ ρuv \\ p + ρv^{2} \\ pv + ρvE \end{matrix} ) .  (2)

Here x and y are the Cartesian spatial coordinates, ρ is the fluid density, u, v are the fluid velocity components along the x- and y-axis, respectively; E = ε + (u^{2} + v^{2})/2, ε is the specific internal energy, and p is the pressure. The system (1), (2) is completed in the following by the ideal gas equation of state

p = (γ − 1)ρε,  (3)

where γ is the ratio of the gas specific heats, usually γ = const > 1.

The proposed four stages of our scheme checker can be used optionally; the most complete check-up is performed when all the four stages are executed.

Stage 1. It is well known that the Euler equations (1)-(2) admit the solution u(x, y, t) = const. From the hydromechanical viewpoint such a solution describes an uniform (homogeneous) flow. Therefore, it is natural to require that a difference scheme approximating the system (1), (2) also possesses this property. The check-up of this property can easily be done with the aid of a symbolic computation. Let us indeed assume that u_{ml}^{n} = u_{0} for any node m, l of the scheme stencil. Then we substitute the constant components of the vector u_{0} instead of the grid values u_{m,l}^{n} into the difference scheme to be analyzed. If we obtain that u_{j,k}^{n+1} = u_{0}, then the scheme under consideration possesses the needed property. For example, it is easy to show that the MacCormack scheme [7] approximating the system (1)-(3) satisfies the conditions for the exact reproduction of a homogeneous flow when it is calculated by this scheme.

The Navier-Stokes equations governing the flow of a compressible, viscous, heat conducting fluid may be written in the form

∂u/∂t + ∂F(u)/∂x + ∂G(u)/∂y = ∂F_{v}(u)/∂x + ∂G_{v}(u)/∂y,  (4)

where F_{v}(u) and G_{v}(u) are the vectors of the viscous fluxes, whose detailed expressions can be found in [8]. It follows from (4) that the Navier-Stokes equations also allow the constant solution u(x, y, t) = u_{0}. Therefore, the corresponding difference equations can also be checked at Stage 1.

This stage can easily be implemented with the aid of the REDUCE computer algebra system [9]. Let us present the corresponding REDUCE program for the case of a check-up of the difference scheme

\frac{u_{j}^{n+1} − u_{j}^{n}}{τ} + a\frac{u_{j+1}^{n} − u_{j−1}^{n}}{2h} = ν \frac{u_{j+1}^{n} − 2u_{j}^{n} + u_{j−1}^{n}}{h^{2}}

approximating the advection-diffusion equation ∂u/∂t + a∂u/∂x = ν∂^{2}u/∂x^{2}, a = const, ν = const > 0.

operator u;
% a difference scheme with an input error
sc1:= (u(j,n+1) - u(j,n))/dt + a*(u(j+1,n) - u(j-1,n))/(2*h)
      - nu*(u(j+1,n) - 2*u(j,n) - u(j-1,n))/(h**2);
out res1;
for all j, k let u(j,k) = uu0;
 sc2:= sc1;
if sc2 neq 0 then
begin
write "Difference equation does not";
write "reproduce homogeneous flow";
pause;
end;
   for all j,k clear u(j,k);
sc1:= sc1;
shut res1;
end;

In this case we obtain that

       2*NU*UU0
SC2 := --------
           2

         H

This means that SC2≠ 0, therefore, the message

Difference equation does not
reproduce homogeneous flow

is printed out inside the file res1. The REDUCE command pause enables the user to interrupt or to continue the further symbolic processing of the difference equations.

Stage 2. In some cases it is possible to derive by hand the amplification matrix G as a function of the Jacobi matrices A and B, where

A = ∂F(u)/∂u, B = ∂G(u)/∂u.  (5)

Then we can compare the REDUCE - generated entries of G with the exact expressions for the entries of A and B. But the derivation of the analytic expression for the amplification matrix by hand for more complex difference schemes becomes very difficult and can naturally lead to errors. Thus this method has only a very limited scope of application.

Stage 3. If the amplification matrix G is a function of the Jacobi matrices (5), then it is useful to perform the similarity transformation of the amplification matrix G : G_{0} = T_{2}GT_{2}^{−1}. The form of the nonsingular matrices T_{2} and T_{2}^{−1} can be found, for example, in [7]. It is well known that the matrices G and G_{0} have the same eigenvalues. But the matrix G_{0} has a much simpler structure. Furthermore, it can be shown that the entries of G_{0} are nondimensional functions of the vector κ = (κ_{1}, …, κ_{M}) and vector ξ = (ξ_{1}, ξ_{2}), where M ≥ 1 and κ_{1}, …, κ_{M} are the nondimensional variables (for example, the Courant number) in the domain of the variation of which the stability region is to be determined; ξ_{j} = k_{j}h_{j}, j = 1, 2, h_{1} and h_{2} are the steps of an uniform grid along the xand y-axes, respectively; k_{1} and k_{2} are the wavenumbers entering the Fourier harmonic

u_{j,k}^{n} = u_{0}e^{i(jk_{1}h_{1}+kk_{2}h_{2}−ωnτ)},  (6)

τ is the time step of a difference scheme, ω is the wave frequency.

If the matrix G_{0} obtained as a result of the execution of Stage 3 does indeed not depend on h_{1}, h_{2} and τ, then the input difference equations are correct, and the matrix G_{0} can be used for further symbolic-numeric computations instead of G.

Stage 4. Truncation error analysis of the input difference equations on a computer. For example, in the case of the MacCormack scheme [7] it is known that the approximation order is equal to two both in space and in time.

Therefore, we can obtain symbolically the first differential approximation of scheme written as the input information for the REDUCE program and then compare the obtained approximation order with the order 2. Steinberg et al. [10] have used such an automated truncation error analysis for the symbolic testing of a FORTRAN code generated symbolically with the aid of the MACSYMA system. As it was noted in [10], the symbolic truncation error analysis, along with numerical tests, reduces the number of possible errors in the generated FORTRAN code, so that the probability of any remaining errors is very small.

It should be noted that the task of symbolic computation of a differential approximation, especially in the cases of difference schemes with intermediate steps, like the MacCormack scheme, may lead to laborious computational work and to large memory and CPU time expenses. Therefore, it appears that Stage 4 is generally realizable only on sufficiently powerful computers.

3 Symbolic Computation of the Amplification Matrix

Let us consider the difference scheme

C_{1}u^{n+1} + C_{2}u^{n} = 0,  (7)

where C_{1} and C_{2} are the matrix difference operators. In the Fourier method of the stability analysis the solution of the form (6) is substituted into the difference scheme (7). As a result one obtains from (7) the system

H_{1}(λ, h, τ, ξ)u^{n+1} + H_{2}(h, τ, ξ)u^{n} = 0  (8)

with λ = exp(−iωτ), h = (H_{1}, …, h_{L}), L is the number of spatial variables, 1 ≤ L ≤ 3. If the matrix H_{1} is nonsingular, equation (8) can be solved with respect to u^{n+1}:

u^{n+1} = Gu^{n},  (9)

where G = −H_{1}^{−1}H_{2} is the amplification matrix of scheme (7). Now let λ_{1}, …, λ_{m} (m = 4 in the case of the Euler equations (1), (2)) be the eigenvalues of G. The necessary von Neumann stability condition has the form [11]

|λ_{j}| ≤ 1 + O(τ), j = 1, …, m.  (10)

In the symbolic-numerical methods, which are described in the present paper, the amplification matrix G is found with the aid of the symbolic computations. After that the characteristic polynomial

f(λ, κ, ξ) = det(G − λI) = ∑_{j=0}^{m} c_{j}(κ, ξ)λ^{m−j}  (11)

is computed also in a symbolic form. The estimates of the length of the expanded polynomial (11), which were presented in [12], show that the symbolic computation of (11) makes increased demands for the size of the computer memory. A method was presented in [12] for the solution of the problem of the memory exhaustion in the process of the computation of the coefficients of (11). In this method, the symbolic computation of the entries of the amplification matrix G is at first optimized with the aid of the SCOPE package, which is now built in the REDUCE system. The coefficients of the characteristic equation (11) are then computed numerically with double precision by using the Newton identities. This method has enabled us to avoid the memory exhaustion even in the case of difference schemes for the 3D Navier-Stokes equations.

In the symbolic-numerical methods there exists a question: At what point to start the numerical computation? It appears that the method proposed in [12] gives a right answer to this question. In addition, the application of the SCOPE package reduces significantly the overall number of the arithmetic operations needed to compute the entries of the amplification matrix G. This is very important for the reduction of the roundoff errors accumulation at the numerical stages of the symbolic-numeric methods for the stability analysis. In this connection, the REDUCE system is preferable for the implementation of the above analysis methods. No other computer algebra system has an optimization program like SCOPE.

4 The modified Routh Algorithm

After we have obtained the characteristic equation coefficients we can investigate the stability of a difference scheme in different ways. One method, which was implemented in [1,2],reduces the von Neumann analysis to a generalized Routh-Hurwitz problem. Let

g(w, κ, ξ) = (w − 1)^{m}f((w + 1)/(w − 1), κ, ξ)  (12)

and let w_{j}, j = 1, …, m, be the zeros of the polynomial (12). Then the condition Rew_{j} ≤ 0 corresponds to the condition |λ_{j}| ≤ 1, j = 1, …, m. We have used in [1] for the numerical solution of the generalized Routh-Hurwitz problem for polynomial (12) the modified Routh algorithm. The modification reduces to a scaling of certain 2 × 2 matrices, which can be introduced within the framework of the Routh table. It turned out that this approach is much more accurate than the use of the double precision arithmetic. The latter is inapplicable in cases of complex difference schemes for 2D and especially for 3D fluid dynamics problems.

With the aid of the modified Routh algorithm the necessary stability condition (NSC) of scheme [7] was found in the form

|κ_{2}| + |κ_{3}|κ_{4} + κ_{1}(1 + κ_{4}^{2})^{0.5} ≤ 1,  (13)

which coincides with the sufficient stability condition found in [13] with the aid of a number of test calculations by the MacCormack scheme. The nondimensional variables κ_{1}, …, κ_{4} entering formula (13) are defined as follows:

κ_{1} = cτ/h_{1}, κ_{2} = uτ/h_{1}, κ_{3} = vτ/h_{1}, κ_{4} = h_{1}/h_{2}.  (14)

It should be noted that the variable κ_{5} = γ, where the constant γ enters the gas equation of state (3), should also be included in a general case in the set of the variables κ_{1}, …, κ_{M}. It is fortunate that γ does not enter explicitly the NSC (13) of the MacCormack scheme. When we initiated the calculations for the purpose of determining the NSC for the scheme of [7] we have introduced γ in the set κ_{1}, …, κ_{5}. But the subsequent runs at different values of γ showed that in the case of the MacCormack scheme the NSC does not depend on γ. This enabled us to reduce by one the dimension M of the Euclidean space E^{M} in which the points of the boundary of the necessary stability region were determined. We do not know whether the ratio of the specific heats γ will always be absent in the necessary stability condition of any difference scheme approximating the Euler equations (1)-(3).

More complex difference schemes often contain some nondimensional weight parameters α_{1}, …, α_{N}. These are, for example, the coefficients of the added artificial dissipation terms as in [14]. Then these parameters should generally be included in the set of the variables κ_{1}, …, κ_{M}:

κ_{1} = cτ/h_{1}, κ_{2} = uτ/h_{1}, κ_{3} = vτ/h_{1}, κ_{4} = h_{1}/h_{2},
κ_{4+ν} = α_{ν}, ν = 1, …, N.

In such cases much more computational work is required to obtain the NSC in the space E^{4+N}.

5 The Use of the Resultants

It was proposed in [15] to use the algebra of the resultants together with the catastrophe theory methods for the stability investigation of difference initialvalue problems. In the subsequent work [3] we have proposed a symbolicnumerical method based on the resultants only. Since Rew_{j} ≤ 0, the polynomial (12) should have at least one purely imaginary zero on the boundary Γ of the stability region. Let us set w = iσ, where σ is a real number, and consider the polynomial φ(σ, κ, ξ) = g(iσ, κ, ξ). The necessary stability region boundary Γ is determined by those values of κ and ξ for which φ(σ, κ, ξ) has a real zero σ. These zeroes are determined by the system

Reφ(σ, κ, ξ) = 0, Imφ(σ, κ, ξ) = 0.  (15)

It is known that this system has the solution if and only if the resultant R(κ, ξ) of the polynomials Reφ and Imφ is equal to zero:

R(κ, ξ) = res(Reφ, Imφ) = 0.  (16)

Let

Re φ = b_{0}σ^{p_{1}} + b_{1}σ^{p_{1}−1} + ⋯ + b_{p_{1}},
Im φ = c_{0}σ^{p_{2}} + c_{1}σ^{p_{2}−1} + ⋯ + c_{p_{2}},

where p_{1} ≤ m, p_{2} ≤ m. The Sylvester matrix S of these polynomials is a (p_{1} + p_{2}) × (p_{1} + p_{2}) - matrix of the form

S = ( \begin{matrix}b_{0} & b_{1} & … & b_{p_{1}} & 0 & … & … & 0 \\ 0 & b_{0} & b_{1} & … & b_{p_{1}} & … & … & 0 \\ … & … & … & … & … & … & … & … \\ 0 & 0 & … & 0 & b_{0} & … & … & b_{p_{1}} \\ c_{0} & c_{1} & … & … & c_{p_{2}} & 0 & … & 0 \\ 0 & c_{0} & … & … & … & c_{p_{2}} & … & 0 \\ … & … & … & … & … & … & … & … \\ 0 & … & … & … & c_{0} & … & … & c_{p_{2}} \end{matrix} )\} p_{2} rows \} p_{1} rows  (17)

Then the resultant R is the determinant of the matrix S: R = det S. The coefficients of polynomial (11) are usually periodic functions of the spectral variables ξ_{1}, …, ξ_{L} with periods T_{1}, …, T_{L} (for example, T_{1} = 2π). Consider a parallelepiped

Π = {0 ≤ ξ_{l} ≤ T_{l}, l = 1, …, L}  (18)

in an L-dimensional Euclidean space of the ξ points. Now we can construct as in [12] the following binary function to characterize the stability and instability points:

b(κ) = \{ \begin{matrix}1, & sign R(κ, ξ) = const ∀ξ ∈ Π \\ 0, & otherwise. \end{matrix}  (19)

The κ-points at which b(κ) = 1 are the points where the difference scheme is stable in the sense of the satisfaction of the von Neumann conditions (10).

There are a number of the ways of using the resultant within the framework of some symbolic-numerical method. For example, it was proposed in [12] to use the SCOPE package only for the compression of the expressions for the entries of the matrix (17). The further computation of the determinant R = det(S) was proposed in [12] to be carried out by some numerically stable library facility.

We have also tried another strategy, in which the analytic expression for the resultant R was computed symbolically with the aid of the REDUCE system and SCOPE package. In the case of the three-dimensional MacCormack scheme [7] this has led to a very large generated FORTRAN subroutine containing more than 2,500 lines. It is clear that a numerical computation by this subroutine would lead to an unacceptable level of the roundoff errors even in the case of using the double precision arithmetic.

We have also implemented the third algorithm, in which the analytic expressions for the entries of S were compressed with the aid of SCOPE as in [12]. But the numerical value of the resultant R was further computed with the aid of our own FORTRAN subroutine. The substantial part of this subroutine is the use of the program balanc from [16] for the equilibration of the values of the entries of the Sylvester’s matrix S. This enabled us to obtain the acceptable results on the determination of the boundary Γ even without a recourse to the double precision arithmetic.

Table 1

\begin{matrix}κ_{5} & κ_{6} & κ_{1c} & κ_{1e} \\ 1.0 & 1.0 & 0.5763 & 0.5774 \\ 2.0 & 1.0 & 0.4022 & 0.4082 \\ 2.0 & 2.0 & 0.3260 & 0.3333 \end{matrix}

We applied this symbolic-numerical method to the stability analysis of the MacCormack scheme [7] for the three-dimensional Euler equations. The necessary stability condition for this scheme was found to have the form [12]

|κ_{2}| + |κ_{3}|κ_{5} + |κ_{4}|κ_{6} + κ_{1}(1 + κ_{5}^{2} + κ_{6}^{2})^{0.5} ≤ 1,  (20)

where

κ_{1} = cτ/h_{1}, κ_{2} = uτ/h_{1}, κ_{3} = vτ/h_{1}, κ_{4} = wτ/h_{1},
κ_{5} = h_{1}/h_{2}, κ_{6} = h_{1}/h_{3}, κ_{7} = γ,  (21)

and u, v, w are the velocity components along the x- , y- and z-axes, respectively.

In Table 1 we present the values of the quantity κ_{1} at κ_{2} = κ_{3} = κ_{4} = 0 and at different values of the quantities κ_{5} and κ_{6}; κ_{1e} = (1 + κ_{5}^{2} + κ_{6}^{2})^{−0.5} in Table 1 and the values κ_{1c} have been obtained by the symbolic-numeric method. The data of Table 1 confirm the result (20). The overall length of the automatically generated FORTRAN subroutines for the computation of the entries of the amplification matrix G amounted to 700 lines, with 60 symbols in each line.

6 Euler Equations and Curvilinear Grids

The major part of applied fluid dynamics problems are now solved numerically on curvilinear spatial grids. Therefore, it is important to have the methods for the stability analysis of difference schemes written down on such grids.

Let us write the nonstationary Euler equations in a fixed-in-time system of curvilinear coordinates ξ, η and the time t:

\frac{∂uJ}{∂t} + \frac{∂F̃}{∂ξ} + \frac{∂G̃}{∂η} = 0,  (22)

where

F̃ = Jξ_{x}F(u) + Jξ_{y}G(u), G̃ = Jη_{x}F(u) + Jη_{y}G(u),  (23)

and the flux vectors F and G are defined by formulas (2). The Jacobian J of the transformation

x = x(ξ, η), y = y(ξ, η)  (24)

is computed by the formula

J = x_{ξ}y_{η} − y_{ξ}x_{η}.  (25)

In order to understand the structure of the step operator of a difference scheme approximating the Euler system of equations (22) let us take a threestage Runge-Kutta scheme of the form [14]

u^{(0)} = u^{n};
u^{(1)} = u^{(0)} − α_{1}τP_{h}u^{(0)};
u^{(2)} = u^{(0)} − α_{2}τP_{h}u^{(1)};
u^{(3)} = u^{(0)} − τP_{h}u^{(2)};
u^{n+1} = u^{(3)},

(26)

where α_{1} and α_{2} are the weight parameters, and P_{h}u is the difference operator approximating the spatial differential operator

P u = \frac{1}{J} (\frac{∂F̃}{∂ξ} + \frac{∂G̃}{∂η})  (27)

with the aid of central differences. In accordance with [14], an efficient scheme of the form (26) is obtained at α_{1} = α_{2} = 0.6. The scheme (26) with α_{1} = α_{2} = 0.6 was used in [17] for the computation of turbulent gas flows with shock waves with the use of the finite volume method.

Let us carry out a Fourier stability analysis of the linearized scheme obtained from scheme (26). The linearized operator P̄u has the form

P̄_{h}u_{j,k} = A_{1}(u_{j,k}, ξ_{j}, η_{k})\frac{u_{j+1,k} − u_{j−1,k}}{2} + A_{2}(u_{j,k}, ξ_{j}, η_{k})\frac{u_{j,k+1} − u_{j,k−1}}{2} ,

(28)

where

A_{1} = \frac{1}{J} \frac{∂F̃}{∂u} = ξ_{x}A + ξ_{y}B, A_{2} = \frac{1}{J} \frac{∂G̃}{∂u} = η_{x}A + η_{y}B,  (29)

and the Jacobi matrices A and B are defined by formulas (5). At the frozen coefficients in (28), the scheme (26) becomes a linear difference scheme. Eliminating the intermediate values u^{(1)} and u^{(2)} we obtain the following two-level difference scheme [18]:

u^{n+1} = S(ξ, η, u_{0})u^{n},  (30)

where the step operator S(ξ, η, u_{0}) has the form

S(ξ, η, u_{0}) = I − τP_{h} + α_{2}(τP_{h})^{2} − α_{1}α_{2}(τP_{h})^{3}.  (31)

Let us denote by F(−τP) the Fourier symbol of the operator P_{h}, and let

C = F(−τP).

Then it is easy to compute with regard for (28) that

C = −iτ(sin ξ_{1}A_{1} + sin ξ_{2}A_{2}).  (32)

The amplification matrix G of scheme (28) then has in accordance with (31) the form

G = I + C + α_{2}C^{2} + α_{1}α_{2}C^{3}.  (33)

Let μ_{j}, j = 1, 2, 3, 4, be the eigenvalues of the matrix C. Then the eigenvalues λ_{j}, j = 1, 2, 3, 4, of the matrix G (33) have the form

λ_{j} = 1 + μ_{j} + α_{2}μ_{j}^{2} + α_{1}α_{2}μ_{j}^{3}, j = 1, 2, 3, 4.

The matrix C represents a linear combination of the Jacobi matrices A and B:

C = d_{1}A + d_{2}B,  (34)

where with regard for (29)

d_{1} = −iτξ_{x} sin ξ_{1} − iτη_{x} sin ξ_{2}, d_{2} = −iτξ_{y} sin ξ_{1} − iτη_{y} sin ξ_{2}.  (35)

The eigenvalues μ_{j} of the matrix (34) are known to have the form [19]

\begin{matrix}μ_{1} = μ_{2} = d_{1}u + d_{2}v, \\ μ_{3} = μ_{1} + c(d_{1}^{2} + d_{2}^{2})^{0.5}, \\ μ_{4} = μ_{1} − c(d_{1}^{2} + d_{2}^{2})^{0.5}, \end{matrix}  (36)

where c is the sound velocity in the gas. It follows from (36) that all the eigenvalues μ_{j}, j = 1, …, 4, are purely imaginary. Let us set as in [18] μ_{j} = iσ_{j}, j = 1, …, 4, where σ_{j} is a real number. Then the von Neumann stability condition |λ_{j}| ≤ 1 implies

α_{1}^{2}α_{2}^{2}σ^{4} + σ^{2}(α_{2}^{2} − 2α_{1}α_{2}) + 1 − 2α_{2} ≤ 0.

With regard for (36) we can then obtain the necessary stability condition of scheme (26) in the form [18]

|d_{1}u| + |d_{2}v| + c(|d_{1}|^{2} + |d_{2}|^{2})^{0.5} ≤
\frac{1}{α_{1}} [\frac{2α_{1} − α_{2} + \sqrt{α_{2}(α_{2} − 4α_{1} + 8α_{1}^{2})}}{2α_{2}}]^{0.5}.  (37)

Basing on the formulas (36) and (37) we now want to determine the form and the number of the nondimensional variables κ_{1}, …, κ_{M} which enter the stability condition (37). It is easy to see from the expressions (35) and (36) that we can introduce the following nondimensional variables:

κ_{1,1} = τξ_{x}c, κ_{1,2} = τξ_{y}c, κ_{1,3} = τη_{x}c, κ_{2,1} = τξ_{x}u,
κ_{2,2} = τη_{x}u, κ_{3,1} = τξ_{y}v, κ_{3,2} = τη_{y}v, κ_{4} = η_{y}/ξ_{x}.  (38)

We have 8 variables in the total. In terms of these variables, we can obtain from (37) the necessary stability condition of scheme (26) on a general curvilinear grid in the closed form:

|κ_{2,1}| + |κ_{2,2}| + |κ_{3,1}| + |κ_{3,2}| + c((|κ_{1,1}| + |κ_{1,3}|)^{2} + (|κ_{1,2}| +
|κ_{1,1}κ_{4}|)^{2})^{0.5} ≤

\frac{1}{α_{1}} [\frac{2α_{1} − α_{2} + \sqrt{α_{2}(α_{2} − 4α_{1} + 8α_{1}^{2})}}{2α_{2}}]^{0.5}.

For the subsequent analysis it is convenient to express the derivatives ξ_{x}, ξ_{y}, η_{x}, η_{y} in terms of the derivatives x_{ξ}, x_{η}, y_{ξ}, y_{η} [20]:

ξ_{x} = \frac{1}{J} y_{η}, ξ_{y} = −\frac{1}{J} x_{η}, η_{x} = −\frac{1}{J} y_{ξ}, η_{y} = \frac{1}{J} x_{ξ}.  (39)

The derivatives x_{ξ}, x_{η}, y_{ξ}, y_{η} are usually approximated by central differences [20]:

(x_{ξ})_{j,k} = (x_{j+1,k} − x_{j−1,k})/2, (y_{η})_{j,k} = (y_{j,k+1} − y_{j,k−1})/2,  (40)

etc. A more explicit expression for uτξ_{x}, by using (25) and (39), reads

uτ(ξ_{x})_{j,k} = \frac{uτ(y_{η})_{j,k}}{(x_{ξ}y_{η} − x_{η}y_{ξ})_{j,k}} .  (41)

For the particular case of a uniform rectangular mesh in the (x, y) plane, we have

\begin{matrix}x_{j+1,k} = x_{j,k} + h_{1}, & x_{j,k+1} = x_{j,k}, \\ y_{j+1,k} = y_{j,k}, & y_{j,k+1} = y_{j,k} + h_{2} & ∀j, k. \end{matrix}  (42)

Substituting (40) into (41) and using (42), we obtain a simplified version of (41):

uτ(ξ_{x})_{j,k} = uτh_{2}/(h_{1}h_{2}) = uτ/h_{1}.  (43)

The expression in the right hand side of (43) coincides with the value of κ_{2} in (14). By considering the remaining variables in the set (38) in a similar way we can see that (38), for the uniform rectangular grid of (42), reduces to:

κ_{1,1} = cτ/h_{1}, κ_{1,2} = κ_{1,3} = 0, κ_{2,1} = uτ/h_{1},
κ_{2,2} = 0, κ_{3,1} = 0, κ_{3,2} = vτ/h_{2}, κ_{4} = h_{1}/h_{2}.  (44)

It is easy to see that the set of the variables (44) coincides with the four variables of (14), as was to be expected.

Thus, we can see that the number of the nondimensional variables κ_{1}, …, κ_{M} in the case of the consideration of a fixed curvilinear spatial grid has increased by a factor of 2 in comparison with the case of an uniform rectangular grid.

It follows from the above consideration that at least for some classes of difference schemes on curvilinear grids for the Euler equations the characteristic equation of a scheme can be written in the form (11), where

κ = (κ_{1,1}, κ_{1,2}, κ_{1,3}, κ_{2,1}, κ_{2,2}, κ_{3,1}, κ_{3,2}, κ_{4}).

7 Navier-Stokes Equations and Curvilinear Grids

Let us assume that the Navier-Stokes equations are to be solved in a physical domain D with (x, y, z) points. We further assume that there exist functions

x = x(ξ, η, ζ), y = y(ξ, η, ζ), z = z(ξ, η, ζ),  (45)

which enable us to uniquely map D onto a domain in the space with points (ξ, η, ζ), where Oξηζ is a rectangular Cartesian coordinate system. Because of the memory/speed limitations of computers currently available, simplified versions of the full system of the compressible Navier-Stokes equations are used for the numerical simulation of flows with high Reynolds number. Following [21] let us assume that ξ, η, and ζ represent the streamwise, normal, and spanwise coordinates, respectively, for the flow around an aircraft wing. Since the dominant viscous effects for turbulent flows with high Reynolds number arise from viscous diffusion normal to the body surface, a thin-layer assumption may be employed by retaining only the viscous diffusion terms in the η-direction. Then the simplified 3D Navier-Stokes equations have the form

\frac{∂u⃗}{∂t} + \frac{∂F⃗ (u⃗)}{∂ξ} + \frac{∂G⃗(u⃗)}{∂η} + \frac{∂H⃗(u⃗)}{∂ζ} = \frac{∂G⃗_{v}(u⃗)}{∂η},  (46)

where

u⃗ = J [ \begin{matrix}ρ \\ ρu \\ ρv \\ ρw \\ ρE \end{matrix} ] , F⃗ = J [ \begin{matrix}ρU \\ ρuU + ξ_{x}p \\ ρvU + ξ_{y}p \\ ρwU + ξ_{z}p \\ (ρE + p)U \end{matrix} ] ,

G⃗ = J [ \begin{matrix}ρV \\ ρuV + η_{x}p \\ ρvV + η_{y}p \\ ρwV + η_{z}p \\ (ρE + p)V \end{matrix} ] , H⃗ = J^{−1} [ \begin{matrix}ρW \\ ρuW + ζ_{x}p \\ ρvW + ζ_{y}p \\ ρwW + ζ_{z}p \\ (ρE + p)W \end{matrix} ] ,  (47)

and

U = ξ_{x}u + ξ_{y}v + ξ_{z}w, V = η_{x}u + η_{y}v + η_{z}w, W = ζ_{x}u + ζ_{y}v + ζ_{z}w, (48)

J is the Jacobian of the transformation (45),

J = x_{ξ}(y_{η}z_{ζ} − y_{ζ}z_{η}) − y_{ξ}(x_{η}z_{ζ} − x_{ζ}z_{η}) + z_{ξ}(x_{η}y_{ζ} − y_{η}x_{ζ}).  (49)

In Eqs. (46)-(48), u, v and w are the fluid velocity components along the x-,y-and z-axes, respectively. The expression for the viscous flux G_{v} is as follows:

G⃗_{v} = \sqrt{γ} με̄J \frac{M_{∞}}{Re_{∞}} [ \begin{matrix}0 \\ ϕ_{1G}u_{η} + η_{x}ϕ_{2G} \\ ϕ_{1G}v_{η} + η_{y}ϕ_{2G} \\ ϕ_{1G}w_{η} + η_{z}ϕ_{2G} \\ ϕ_{1G}[(\frac{q^{2}}{2} )_{η} + (\frac{γ}{γ−1})\frac{ε̃}{σε̄}T_{η}] + V ϕ_{2G} \end{matrix} ] ,  (50)

where

q^{2} = u^{2} + v^{2} + w^{2}, ϕ_{1G} = η_{x}^{2} + η_{y}^{2} + η_{z}^{2},

ϕ_{2G} = \frac{1}{3}(η_{x}u_{η} + η_{y}v_{η} + η_{z}w_{η}).

T is the fluid temperature, M_{∞} is the freestream Mach number, Re_{∞} is the freestream Reynolds number. The transformation formulas read as follows:

\begin{matrix}ξ_{x} = (y_{η}z_{ζ} − y_{ζ}z_{η})/J & η_{x} = (y_{ζ}z_{ξ} − y_{ξ}z_{ζ})/J & ζ_{x} = (y_{ξ}z_{η} − y_{η}z_{ξ})/J \\ ξ_{y} = (x_{ζ}z_{η} − x_{η}z_{ζ})/J & η_{y} = (x_{ξ}z_{ζ} − x_{ζ}z_{ξ})/J & ζ_{y} = (x_{η}z_{ξ} − x_{ξ}z_{η})/J \\ ξ_{z} = (x_{η}y_{ζ} − x_{ζ}y_{η})/J & η_{z} = (x_{ζ}y_{ξ} − x_{ξ}y_{ζ})/J & ζ_{z} = (x_{ξ}y_{η} − x_{η}y_{ξ})/J \end{matrix} (51)

In the sets of equations (46)-(47), (50), distances have been put in a dimensionless form as in [21]. The turbulence modeling is here performed by using the algebraic turbulence model of Baldwin and Lomax [22]. This model is based on the concepts of an eddy viscosity and eddy conductivity. In the momentum equations, the molecular viscosity μ is replaced by the effective viscosity μ_{e}, μ_{e} = μ + μ_{t} = μ(1 + μ_{t}/μ) = με̄, where ε̄ = 1 + μ_{t}/μ and μ_{t} denotes the turbulent viscosity. Similarly, in the energy equation, the molecular conductivity k is replaced by the effective conductivity k_{e},

k_{e} = k + k_{t} = \frac{c_{p}}{σ} μ + \frac{c_{p}}{σ_{t}} μ_{t} = \frac{c_{p}μ}{σ}(1 + \frac{σ}{σ_{t}} \frac{μ_{t}}{μ} ) = \frac{c_{p}μ}{σ}ε̃ = kε̃,

where ε̃ = 1 + (σμ_{t})/(σ_{t}μ). Here σ is the laminar Prandtl number and σ_{t} the turbulent Prandtl number; σ = 0.72, σ_{t} = 0.9.

It was shown in [5] that in the particular case of a Cartesian uniform rectangular grid, the necessary stability region should be determined in the ninedimensional Euclidean space E^{9} of the variables

κ_{1} = \frac{cτ}{h_{1}} , κ_{2} = \frac{uτ}{h_{1}} , κ_{3} = \frac{vτ}{h_{1}} , κ_{4} = \frac{wτ}{h_{1}} , κ_{5} = \frac{h_{1}}{h_{2}} ,

κ_{6} = \frac{h_{1}}{h_{3}} , κ_{7} = γ, κ_{8} = \sqrt{γ} \frac{M_{∞}}{Re_{∞}} \frac{τμε̄}{h_{1}^{2}ρ}, κ_{9} = \frac{ε̃}{ε̄σ} ,  (52)

where h_{1}, h_{2} and h_{3} are the steps of an uniform grid along the axes x, y, and z, respectively; c is the sound speed.

In the case of the general curvilinear coordinates ξ, η and ζ it was shown in [23] that the coefficients of the characteristic polynomial (11) depend on the following dimensionless variables:

\begin{matrix}κ_{1,1} = τξ_{x}c, & κ_{1,2} = τξ_{y}c, & κ_{1,3} = τξ_{z}c, & κ_{1,4} = τη_{x}c, \\ κ_{1,5} = τη_{y}c, & κ_{1,6} = τη_{z}c, & κ_{1,7} = τζ_{x}c, & κ_{1,8} = τζ_{y}c, \\ κ_{1,9} = τζ_{z}c, & κ_{2,1} = τξ_{x}u, & κ_{2,2} = τη_{x}u, & κ_{2,3} = τζ_{x}u, \\ κ_{3,1} = τξ_{y}v, & κ_{3,2} = τη_{y}v, & κ_{3,3} = τζ_{y}v, & κ_{4,1} = τξ_{z}w, \\ κ_{4,2} = τη_{z}w, & κ_{4,3} = τζ_{z}w, & κ_{5} = η_{y}/ξ_{x}, & κ_{6} = ζ_{z}/ξ_{x}, \\ κ_{7} = γ, & κ_{8} = \sqrt{γ} \frac{M_{∞}}{Re_{∞}} \frac{τμε̄}{ρ}(ξ_{x})^{2}, & κ_{9} = \frac{ε̃}{ε̄σ} . \end{matrix} (53)

In the particular case of a uniform rectangular mesh in the (x, y, z) space, we have

\begin{matrix}x_{m+1,q,r} = x_{m,q,r} + h_{1}, x_{m,q+1,r} = x_{m,q,r}, x_{m,q,r+1} = x_{m,q,r}, \\ y_{m+1,q,r} = y_{m,q,r}, y_{m,q+1,r} = y_{m,q,r} + h_{2}, y_{m,q,r+1} = y_{m,q,r}, \\ z_{m+1,q,r} = z_{m,q,r}, z_{m,q+1,r} = z_{m,q,r}, z_{m,q,r+1} = z_{m,q,r}. \end{matrix}  (54)

We can substitute into (53) the expressions (51) for the derivatives ξ_{x}, ξ_{y}, …, ζ_{z}. Then we can approximate the derivatives x_{ξ}, x_{η}, …, z_{ζ} by the central differences similar to (40). Then it is easy to see with regard for (54) that in the particular case of an uniform rectangular grid the set (53) reduces to:

\begin{matrix}κ_{1,1} = cτ/h_{1}, & κ_{1,2} = κ_{1,3} = κ_{1,4} = 0, & κ_{1,6} = κ_{1,7} = κ_{1,8} = 0, \\ κ_{1,5} = cτ/h_{2}, & κ_{1,9} = cτ/h_{3}, & κ_{2,1} = uτ/h_{1}, \\ κ_{2,2} = κ_{2,3} = 0, & κ_{3,1} = κ_{3,3} = 0, & κ_{3,2} = vτ/h_{2}, \\ κ_{3,3} = κ_{4,1} = κ_{4,2} = 0, & κ_{4,3} = wτ/h_{3}, & κ_{5} = h_{1}/h_{2}, \\ κ_{6} = h_{2}/h_{3}, κ_{7} = γ, & κ_{8} = \sqrt{γ} \frac{M_{∞}}{Re_{∞}} \frac{τμε̄}{ρh_{1}^{2}}, & κ_{9} = \frac{ε̃}{ε̄σ} . \end{matrix} (55)

Since cτ/h_{3} = (cτ/h_{1})(h_{1}/h_{3}), the set of the variables (55) indeed coincides with the nine variables of (52), as was to be expected. For a general curvilinear grid, the set (53) encompasses 23 variables. A similar set (52) for the case of an uniform rectangular grid contains only 9 variables. Thus the consideration of a curvilinear grid increases by a factor of 23/9 the dimension of the Euclidean space E^{M} in which the stability region boundary is to be determined.

It was shown in [23] that in the case of a four-stage Jameson’s scheme approximating the system (46)-(50) the expanded characteristic polynomial (11) would require about 4 · 5! × 6^{135} ≈ 5.392 · 10^{107} bytes of computer memory. Therefore, the use of a package for the compression of the expressions for the entries of the amplification matrix G is necessary. The only such package is SCOPE, which is now available in the REDUCE library. After the entries of the matrix G are obtained, the coefficients of the polynomial (11) can be computed numerically with the use of the Newton identities as in [5].

It is important that the analytic form of the coefficients c_{m}(κ, ξ) of(11) remains the same for any curvilinear grid. In this connection we can propose the following procedure for the stability analysis on a general curvilinear grid.

Step 1. Stability analysis of a given difference scheme for a uniform rectangular spatial grid. This step is needed for the check-up of the necessary stability condition which is obtained on a general curvilinear grid. As shown above, the general necessary stability condition has to coincide with that obtained for a uniform rectangular grid, provided we substitute (54) into (53).

Step 2. Symbolic computation of the analytic expressions for the coefficients of the characteristic polynomial (11), and the execution of the substitutions for introducing the dimensionless variables κ_{1}, …, κ_{23} in accordance with (55).

Step 3. Computation of the coordinates κ_{1}, …, κ_{23} of the points of the necessary stability region boundary in the space of the 23 dimensionless variables (55).

Step 4. Determination of the coefficients of analytic formulas which fit well the numerical data obtained in the third step. As a result we obtain the approximate necessary stability condition in the form of closed analytic expressions,

g_{ν}(κ⃗) ≥ 0, ν = 1, …, N  (56)

(N ≥ 1) involving the 23 variables κ_{1}, …, κ_{23} from the set (53).

The numerical data are processed at Step 4 in the least-squares sense. As an example we have considered in [23] the numerical and analytical data obtained in [5] for a four-stage Jameson’s scheme in the particular case of a uniform rectangular grid. The corresponding approximate necessary stability condition in the case of a viscous, heat conducting gas was obtained in [23] in the form

\frac{(|κ_{2}| + |κ_{3}|κ_{5} + |κ_{4}|κ_{6} + (κ_{1} − 0.58344κ_{9}^{0.08517})\sqrt{1 + κ_{5}^{2} + κ_{6}^{2}})^{2}}{8} +

(\frac{κ_{5}^{2}κ_{8}}{0.5222425})^{2} ≤ 1.  (57)

This formula can be used for the check-up of the general necessary stability condition (56) to be obtained at Step 4. It follows from (57) that the consideration of the turbulent viscosity (κ_{8} > 0) can generally impose additional restriction on the time step, because the time step τ enters the numerator of the formula for κ_{8} (see (52)). On the other hand, the consideration of the effective heat conductivity (0 < κ_{9} < 1/σ ≈ 1.388889) relaxes slightly the restriction on the time step τ.

8 Concluding Remarks

Because of a large number of the variables κ_{1}, …, κ_{23} in the space of the variation of which the necessary stability region boundary is determined in the case of the Navier-Stokes equations solved on curvilinear grids it is clear that the main computational work is required at steps 3 and 4 of the stability analysis procedure described in Section 7. In this connection it is reasonable to develop simpler procedures, which would give the accuracy of the stability predictions sufficient for engineering applications.

Another possible direction of the future research is the extension of the stability investigation method outlined in Section 7 for the case of moving curvilinear spatial grids. The motion of grids is calculated in the following situations: (i) when the solution-adaptive grids are used; (ii) when the grid lines are adjusted to the motion of some geometric boundaries (as, for example, in the wing flutter problems, or in the problems on the blood flow through blood vessels).

References

1. V.G. Ganzha, E.V. Vorozhtsov, and R. Liska. Stability analysis of difference initial-value problems with the aid of the REDUCE system. In: Computer Algebra and Its Applications to Mechanics. Novosibirsk, Aug. 28-31; Irkutsk, Sept. 1-3, 1990. Proc. Int. Seminar / Eds. V.G.Ganzha, V.M.Rudenko and E.V. Vorozhtsov. Nova Science Publishers, New York, 1993, p. 31-41.

2. V.G. Ganzha and E.V. Vorozhtsov. The stability analysis of difference schemes by numerical solution of the generalised Routh-Hurwitz problem. Computer Physics Communications 43, 209-216, 1987.

3. V.G. Ganzha, E.V. Vorozhtsov, and C. Zenger. Stabilitätsuntersuchung von Differenzenverfahren mit Hilfe der Resultantenalgebra unter Verwendung des REDUCE-Systems. Preprint of the Institut für Informatik Technische Universität München Tum-I9120, Munich, 1991, p. 1-23.

4. V.G. Ganzha and E.V. Vorozhtsov. A probabilistic symbolic-numerical method for the stability analysis of difference schemes for PDEs. Proceedings of the 1993 Internat. Sympos. on Symbolic and Algebraic Computation, July 6-8, 1993, Kiev, Ukraine/ Ed. M.Bronstein. ACM Press, New York, 9–13, 1993.

5. V.G. Ganzha, E.V. Vorozhtsov, J. Boers, and J.A. van Hulzen. Symbolicnumeric stability investigations of Jameson’s schemes for the thin-layer Navier-Stokes equations. -In: Proc. ISSAC’94, Oxford, July 20-22, 1994 / Eds. J.von zur Gathen and M.Giesbrecht. ACM Press, New York, 1994, p. 234-241.

6. E.V. Vorozhtsov, B.Yu. Scobelev, and V.G. Ganzha. Symbolic-numerical method for the stability analysis of difference schemes on the basis of the catastrophe theory. J. Comput. Phys. 116, 26-38(1995).

7. R.W. MacCormack, The effect of viscosity in hypervelocity impact cratering, AIAA Paper, No. 354, 1969.

8. R. Peyret and T.D. Taylor. Computational Methods for Fluid Flow. Springer-Verlag, New York, Heidelberg, Berlin, 1983.

9. A.C. Hearn. REDUCE 3.5, User’s Manual. Rand Corporation, Santa Monica, Cal., 1993.

10. S. Steinberg and P.J. Roache. Finite-difference algorithm design and automated code generation. Computer Algebra and Its Applications to Mechanics (CAAM’ 90), Novosibirsk-Irkutsk, August 28 - September 3, 1990. Proceedings of an International Seminar/Eds V.G. Ganzha, V.M. Rudenko and E.V. Vorozhtsov. Nova Science Publishers, Inc., New York, 15–29, 1993.

11. R.D. Richtmyer and K.W. Morton. Difference Methods for Initial-Value Problems, Second Edition. Interscience Publishers, John Wiley and Sons, New York, 1967.

12. V.G. Ganzha, E.V. Vorozhtsov, and J.A. van Hulzen. A new symbolicnumeric approach to stability analysis of difference schemes. - In: Proceedings ISSAC’92, July 27-29, 1992, Berkeley, California/ Ed. P.S. Wang. ACM Press, New York, 1992, p. 9-15.

13. R.W. MacCormack. Numerical solution of the interaction of a shock wave with a laminar boundary layer. Proceedings of the Second International Conference on Numerical Methods in Fluid Dynamics /Ed. M. Holt. Lecture Notes in Physics, vol. 8. Springer-Verlag, New York, 151–163, 1971.

14. A. Jameson and W. Schmidt. Some recent developments in numerical methods for transonic flows. Computer Methods in Applied Mechanics and Engineering, 51, 467-493, 1985.

15. B.Yu. Scobelev and E.V. Vorozhtsov. Application of the catastrophe theory methods to the stability analysis of difference schemes. The Second Japan-Soviet Union Joint Symp. on Computational Fluid Dynamics. Tsukuba, Aug. 27–31, 1990. Preprints. Published by Japan Society of Computational Fluid Dynamics. The University of Tsukuba, Tsukuba, 53–60, 1990.

16. J.H. Wilkinson and C. Reinsch. Handbook for automatic computation. Linear algebra. Springer-Verlag, Heidelberg, New York, Berlin, 1971.

17. F. Grasso and C.G. Speziale. Supersonic flow computations by two-equati-on turbulence modeling. In: AIAA 9th Computational Fluid Dynamics Conference, June 13-15, 1989, Buffalo, N.Y. A Collection of Technical Papers. Published by AIAA, Washington, DC 20024, 1989, p. 232-239.

18. V.G. Ganzha and E.V. Vorozhtsov. On the stability of Jameson schemes. In: Bridging Mind and Model: Papers in Applied Mathematics / Ed. P.J. Costa. St. Thomas Technology Press, St. Paul, MN, USA, 1993, p. 231-300.

19. R.F. Warming, R.M. Beam, and B.J. Hyett. Diagonalization and simultaneous symmetrization of the gas-dynamic matrices. Mathematics of Computation 29, 1037–1045, 1975.

20. J.F. Thompson, Z.U.A. Warsi, and C.W. Mastin. Numerical Grid Generation. Foundations and Applications. North-Holland, New York, 1985.

21. V.N. Vatsa. Accurate numerical solutions for transonic viscous flow over finite wings. J. Aircraft 24, No. 6, p. 377-385(1987).

22. B. Baldwin and H. Lomax. Thin layer approximation and algebraic model for separated turbulent flows, AIAA Paper No. 78-257, 1978.

23. E.V. Vorozhtsov. Symbolic-numeric stability investigation of difference schemes for compressible 3D Navier-Stokes equations. In: Proc. Internat. Workshop on Computer Algebra in Science and Engineering - Algorithms, Systems, Applications. Bielefeld, August 28-31, 1994 /Ed. F. Hehl. World Science Publishers, Singapore, 1995 (in press).