Comparison of Different Krylov Subspace Methods Embedded in an Implicit Finite Volume Scheme for the Computation of Viscous and Inviscid Flow Fields on Unstructured Grids

JOURNAL OF COMPUTATIONAL PHYSICS ARTICLE NO.

140, 311–345 (1998)

CP975862

Comparison of Different Krylov Subspace Methods Embedded in an Implicit Finite Volume Scheme for the Computation of Viscous and Inviscid Flow Fields on Unstructured Grids Andreas Meister Institut f¨ur Angewandte Mathematik, Universit¨at Hamburg, Bundesstrasse 55, 20146 Hamburg, Germany E-mail: [email protected] Received November 27, 1996

In the past few years a great variety of different Krylov subspace methods have been developed and investigated for several model equations. This paper is devoted to the comparison of current preconditioned Krylov subspace methods concerning several inviscid and viscous flow problems of interest in engineering applications. Therefore, the design of an implicit finite volume approximation of the Navier– Stokes equations on unstructured grids is described whereby a new combination of an isotropic triangulation with unisotropic subgrids is presented to achieve high resolution for high Reynolds number flows. For the first time, based on a specific selection of different inviscid and viscous flow fields, a reliable answer can be given to the fundamental question concerning the choice of iterative method depending on the underlying flow field in the area of the Euler and Navier–Stokes equations to get a stable and fast numerical scheme. °c 1998 Academic Press

1. INTRODUCTION

One of the most important assignments of the future in computational fluid dynamics is the highly accurate simulation of turbulent and chemically reacting flow fields as well as unsteady fluid flow about moving geometries arising from high capacity aircrafts with strongly pitching airfoils. The development of a modern numerical method which could be used as a basic tool for distinguished actual and subsequent real engineering applications is governed by several fundamental requirements like high resolution of inviscid and viscous flow fields, treatment of complicated geometries, acceptable run-time on modern computer architectures, and 311 0021-9991/98 $25.00 c 1998 by Academic Press Copyright ° All rights of reproduction in any form reserved.

312

ANDREAS MEISTER

extendability to chemically reacting flow fields, turbulent flows, and unsteady flows about moving geometries. In finite volume approximations the needs concerning the treatment of complicated geometries and adaptation as well as the needs for robust and accurate algorithms are perfectly matched. Finite volume methods are formulated on general control volumes and well-known successful approximative Riemann solver can easily be employed in the general framework [69]. During the past 15 years a wide variety of upwind schemes for the discretization of the convective part have been successfully developed. Most of the upwind schemes can be characterized as flux difference splitting biased or flux vector splitting biased. The first wellknown flux vector splitting scheme was proposed by Steger and Warming [71]. This method suffers from pertubations around the sonic point which arise from the fact that the scheme is not continuously differentiable at a sonic point. Later on van Leer [75] presented a further upwind scheme which avoids these errors and H¨anel et al. [33] modified the energy-flux in van Leer’s scheme to preserve total enthalpy in steady state solutions. All these methods are easily implemented but obviously too dissipative for a straightforward employment in a numerical method for viscous fluid flows. During the same period Roe [59] and Osher [55] proposed two different flux difference splitting schemes. Both methods yield an impressive improvement of the accuracy at the cost of increased computational effort. The numerical flux of Roe’s method is not continuously differentiable and therefore not particularly proper for an implicit formulation. In 1991 Liou and Steffen [47] presented a new flux splitting scheme called AUSM (advection upstream splitting method) which successfully combines the accuracy of the flux difference splitting schemes with the simplicity of the flux vector splitting methods. The main drawback of this scheme is the numerical overshoot at shock waves. Recently, Liou and Wada [78] proposed the AUSMDV scheme which overcomes the difficulties of the AUSM method and therefore represents a high accurate scheme in which the computational effort is less than in finite difference splitting biased approaches. In addition to the approximation of the fluxes over the cell interfaces, a time stepping scheme is necessary. Hereafter, the methods are categorized as explicit and implicit schemes. Explicit methods are restricted by a CFL number of 1 which is a restrictive upper bound for several applications. Implicit methods inherently satisfy the stability condition of Courant, Friedrich, and Lewy [14] for any time step, because the numerical domain of dependence always covers the physical one. At present, two different techniques are available. One is the dual-time stepping scheme based on ideas of Jameson [43] and recently proposed by Gaitonde [27] and Arnone et al. [2]. This technique introduces a dual time and reformulates the governing equations so that the time dependent solution of the original equations is equivalent to the steady state solution in pseudo time. Thereafter, different acceleration techniques like multigrid or local time stepping can be used. Nevertheless, the most popular treatment of the implicit time stepping is the linearization of the numerical flux function, which in general leads to a linear system of equations including a large, sparse, nonsymmetric, and indefinite matrix. Therefore, the performance of an implicit method decisively depends on the properties of the incorporated solver for the linear system of equations. For large and sparse matrices direct solvers are too expensive and thus iterative methods are required. More than 40 years ago, Hestenes and Stiefel [37] proposed an efficient method for solving large and sparse linear systems of equations, the conjugate gradient method. Unfortunately, the scheme is inherently restricted to Hermitian positive definite matrices, but its publication provided the decisive

DIFFERENT KRYLOV SUBSPACE METHODS

313

impulse for the development of a wide variety of Krylov subspace methods. In particular the last decade has seen tremendous progress in the development of Krylov subspace methods for arbitrary nonsingular matrices as, for example, the GMRES method proposed by Saad and Schultz [64]. In this method, a function has to be minimized which could be derived from a Petrov–Galerkin condition instead of a Galerkin condition as employed in the CG scheme. Three years later, Sonneveld [70] derived the CGS method from the BiCG scheme [21], whereby two advantages were achieved: first, the acceleration of the convergence rate and, second, the avoidance of adjoint matrix vector multiplications. Motivated by the observation that CGS tends to cause oscillations in the convergence behaviour, van der Vorst [73] presented a further variant of the BiCG algorithm, called BiCGSTAB. Unfortunately, the BiCG method suffers from possible breakdowns and passes this disadvantage to the CGS and BiCGSTAB algorithms. In order to avoid probable breakdown situations, Freund and Nachtigal [24] proposed a BiCG-like approach called QMR. This method combines desirable properties of the GMRES method and the BiCG scheme, so that the algorithm is characterized by a minimization property like GMRES and a short recurrence relation like BiCG. Hence, the scheme ensures smooth convergence behaviour and requires low storage per iteration. The BiCG method needs adjoint matrix vector multiplications and passes this requirement to the QMR algorithm. Therefore, it is not surprising that the same strategy which generates CGS and BiCGSTAB from BiCG was employed to improve QMR. Hence, the QMR-like approaches TFQMR and QMRCGSTAB were developed by Freund [29] and Chan et al. [13], respectively. A further Krylov subspace method which can be interpretated as a combination of GMRES (l) and Bi-CG was developed by Sleipjen and Fokkema [68] and is called BICGstab (`). During the past decade, detailed numerical investigations of Krylov subspace methods for different model equations have been presented [3, 11, 13, 19, 22, 29, 24, 30, 44, 45, 53, 56, 60, 61, 63, 65, 66, 70, 72–74, 80, 82]. Furthermore, a large number of publications prove that these methods are widespread in numerical methods for real engineering applications [1, 12, 18, 42, 49, 50, 52, 54, 58, 67, 76, 77]. However, a systematic investigation and comparison of Krylov subspace methods in combination with applicable preconditioning techniques is still missing. The goals of this paper are to describe the development and analysis of a modern numerical scheme, which satisfies the requirements mentioned above and to close the gap of information concerning the behaviour of current preconditioned Krylov subspace methods in the context of real engineering applications. The outline of this paper is as follows: In Section 2, the governing equations are briefly explained. In Section 3, we describe the discretization of the convective and diffusion terms as well as the employed time stepping scheme. Furthermore, the properties of the AUSMDV method used are outlined. To achieve high resolution, a second-order TVD approach concerning a new combination of isotropic triangulations with unisotropic subgrids is presented. In Section 4, we discuss several Krylov subspace methods and preconditioning techniques. In Section 5, we show the performance and applicability of the new numerical method for the simulation of distinguished flow fields. Thereafter, a special selection of test cases is considered to take into account of the great variety of physical and numerical flow phenomenons of interest. Here, the numerical scheme will be applied to transonic and supersonic inviscid test cases as well as viscous fluid flows, so, the Reynolds number will be varied from 500 up to 6 × 106 . Furthermore, the obtained results are compared with experimental data and numerical simulations of other authors and a critical investigation of different modern preconditioned Krylov

314

ANDREAS MEISTER

subspace methods is presented, which enables reliable statements concerning the stability and accuracy of these iterative solvers to be made. Finally, in Section 6 we finish with some concluding remarks. 2. GOVERNING EQUATIONS

In the present study, the flow is assumed to be governed by the two-dimensional timedependent compressible Navier–Stokes equations. We assume that G ⊂ R2 is a domain which is bounded by polygonal curve segments ∂G k , 1 ≤ k ≤ n, such that ∂G = ∪nk=1 ∂G k . Furthermore, we define D := {(x, t) ∈ R3 | x = (x 1 , x2 )T ∈ G ∪ ∂G, t ≥ 0} and denote the open set Ä ⊂ R4 as the state space. The quantities ρ, v1 , v2 , p, E, H : D → R describe the density, velocity in x1 - and x2 direction, pressure, total energy, and enthalpy of the fluid. The integral form of the Navier– Stokes equations is Z σ

∂t u dx +

Z X 2 σ j=1

∂x j f jc (u) dx =

1 Re∞

Z X 2 σ j=1

∂x j f jv (u) dx,

(2.1)

where u : D → Ä, u = (ρ, ρv1 , ρv2 , ρ E)T is the vector of the conserved variables and f jc , f jv : Ä → R4 , j = 1, 2, are the convective and viscous fluxes which are given by 

  0 ρv j  ρv v + δ p   τ 1j  1 j  1j  f jc (u) =   , f jv (u) =  τ2 j  ρv2 v j + δ2 j p   P2 ρHvj i=1 vi τi j +



µγ ∂ e Pr∞ x j

  , 

respectively. The quantity e denotes the internal energy, which is given by e = E − 12 (v12 + v22 ) and H is defined to be H = E + p/ρ. The pressure is defined by the equation of state p = (γ −1)ρ(E − 12 (v12 + v22 )), where γ denotes the ratio of specific heats. The temperature is given by T = γ (γ − 1)Ma2∞ e, where Ma∞ denotes the Mach number at infinity. The elements of the shear stress tensor are τi j = µ(∂x j vi + ∂xi v j ) + δi j λ(∂x1 v1 + ∂x2 v2 ), with the viscosity assumed to follow the Sutherland law µ = T 1.5 (1 + S)/(T + S), where S = 110K /T∞ and T∞ denotes the temperature at infinity measured in degrees Kelvin. Moreover, the relation between the thermal conductivity and the viscosity is defined by the Stokes’ hypothesis to be λ = − 23 µ and Re∞ and Pr∞ denote the Reynolds and Prandtl numbers at infinity, respectively. Interchanging the time derivative of the conserved variables with the integration, in combination with the use of the rotational invariance of the Euler equations, leads to the form of the Navier–Stokes equations, d dt

Z

Z σ

u dx +

∂σ

T−1 (n)f1c (T(n) u) ds =

2 Z 1 X f v (u)n j ds, Re∞ j=1 ∂σ j

(2.2)

where n = (n 1 , n 2 )T represents the outwards unit normal vector at ∂σ and T denotes the rotational matrix.

DIFFERENT KRYLOV SUBSPACE METHODS

315

3. FINITE VOLUME METHOD

¯ = G ∪ ∂G and the time part In order to solve Eqs. (2.2) numerically, the space part G of the set D have to be discretised. We start from an arbitrary conforming triangulation Th of the domain G which is called the primary mesh and which consists of triangles T . Furthermore, Nh denotes the set of all nodes of the triangulation Th . We denote the three edges of the triangle T by eT ,k , k = 1, 2, 3, and define + R0

E(i) = {eT ,k | k ∈ {1, 2, 3}, T ∈ Th , node xi ∈ eT ,k } and V (i) = {T ∈ Th | node xi is vertex of T }. For the generation of the triangulation we make use of two different strategies. The first scheme is based on an algorithm developed by Friedrich [25] for the generation of inner points in two-dimensional triangulations and provides mostly isotropic triangles. In the case of viscous flow fields, we prefer to use structured grids in boundary layers which are triangulated by the introduction of diagonals. As is shown in example 5.4, we finally combine these subgrids with the triangulation defined by the method described in [25]. We denote triangulations which are based on structured grids by Ths and triangulations which are generated by the use of the algorithm described in [25] by Thu . In general, the description of the control volumes can be accomplished in both cases by a unique definition. However, for practical use it makes sense to define each type of triangulation, independently. We define a discrete control volume σi as the open subset of R2 , including the node xi = (xi1 , xi2 )T and bounded by the straight lines, which are defined by the connection of the midpoint of the edge eT ,k ∈ E(i) with the inner point xs = (xs1 , xs2 )T of the corresponding triangle T (see Fig. 1). In the case when xi is at the boundary of the computational domain, the line defined by the connection of the midpoint of the boundary edge and the node itself is also a part of ∂σi . The definition of the boxes σi which correspond to the above-mentioned triangulations Thu and Ths differs only in the definition of the point xs . In the case of Thu we use xs =

X

αms xm

m∈{i, j,k}

FIG. 1. General form of a control volume of the triangulation Thu (left) and its boundary (right).

316

ANDREAS MEISTER

with αms =

1 2(|li | + |l j | + |lk |)

X

|lm¯ |, |li | = |x j − xk |

¯ m∈{i, j,k} ¯ =m m6

(see Fig. 1). For distorted triangles this definition exhibits the advantage that the deformation P of the control volume in comparison to those which we get by choosing xs = 13 3m=1 xm is much smaller. Considering Ths we use xs =

X

αms xm

αis = α sj = 12 ,

αks = 0,

m∈{i, j,k}

where xi and x j represent the two nodes which are connected by the introduced diagonal (see Fig. 2). Thereby, we have the possibility of using structured subgrids in order to get a suitable boundary layer resolution in the case of high Reynolds number flows. The union of all boxes σi is called the secondary mesh. Obviously, 4 := ∪i σi covers G. Let NR(i) denote the index set of all nodes neighbouring node xi , i.e. those nodes x j for which ∂σi ∩∂σ j 1 ds 6 = 0 is valid. Usually, for j ∈ N (i) the boundary between the control volume σi and σ j consists of two line segments which are denoted by likj , k = 1, 2, and nikj , k = 1, 2, represent the accompanying unit normal vectors. Numerical schemes are not able to compute the exact solution of the Navier–Stokes equations. Hence, we consider the space of piecewise polynomial functions t = {w(., t) : 4 → R4 |w j (., t)|σi ∈ Pm , ∀σi ∈ 4 and j = 1, . . . , 4} Hh,m

for all t ∈ R+ 0 and m ∈ N0 , where Pm represents the space of polynomials with degree less than m + 1. Utilizing our notion of control volumes and introducing the cell average operator by Z 1 (Mh u) (t)|σi := u(x, t) dx ∀σi ∈ 4 |σi | σi

FIG. 2. General form of a control volume of the triangulation Ths .

317

DIFFERENT KRYLOV SUBSPACE METHODS

and the inviscid and viscous operators by Z ¡ c ¢ 1 Lh u (t)|σi = T−1 (n(x))f1c (T(n(x))u(x, t)) ds ∀σi ∈ 4, |σi | ∂σi 2 Z ¡ v ¢ 1 X 1 v Lh u (t)|σi = f (u(x, t))n m (x) ds ∀σi ∈ 4 |σi | m=1 ∂σi Re∞ m

(3.1)

(3.2)

the Navier–Stokes equations (2.1) can be rewritten in the form ¡ ¢ ¡ ¢ d (Mh u) (t)|σi = − Lch u (t)|σi + Lvh u (t)|σi dt

∀σi ∈ 4.

(3.3)

Furthermore ui (t) = u(t)|σi := (Mh u) (t)|σi and the initialization of the piecewise constant functions is given by application of the cell average operator to the given initial function ui (0) := (Mh u)(0)|σi

∀σi ∈ 4.

3.1. Discretization of the Inviscid Terms With regard to the chosen discretization of the computational domain, the convective part of Eq. (3.3) can be written as ¡

2 Z 1 X X (t)|σi = T−1 (n(x))f1c (T(n(x))u(x, t)) ds |σi | j∈N (i) k=1 likj

¢

Lch u

∀σi ∈ 4.

Obviously, the evaluation of the line integrals leads to trouble when u is discontinuous across ∂σi . For the solution of these Riemann problems we introduce the concept of a numerical flux function, which is often called a Riemann solver. Concerning the fact that the implicit numerical scheme to be described should be extendable to turbulent and chemically reacting flows, as well as unsteady flow fields on moving grids, we require some basic properties of the numerical flux function as follows: • • • • • •

numerical efficiency and accuracy, extendability to unsteady flows on moving grids, extendability to chemically reacting and turbulent flows, high resolution for shock waves and contact discontinuities, continuous differentiability, conservation of enthalphy for steady flows.

Furthermore, several instances have been presented where various Riemann solvers show undesirable behaviour and therefore lead to unreliable results. Recently, Quirk [57] investigated some failings like expansion shocks, production of negative internal energies which directly leads to negative pressures, oscillations behind slowly moving shocks, the carbuncle phenomenon, kinked mach stems appearing at the reflection of a shock at a fixed wall and the so-called odd–even decoupling which leads to instability effects when a shock is aligned with the mesh. Thus, in the dilemma between requirements and failings it is very important to review the great variety of presently available numerical flux functions in order to find a sufficiently accurate and robust Riemann solver.

318

ANDREAS MEISTER

Since the early eighties various discretization schemes based on upwind techniques have been developed. These methods include the advantage that no added artificial dissipation is necessary to achieve stability as known from central difference approximations [43]. Therefore, we decided to focus on upwind schemes which are mostly subdivided into the concepts of finite difference splitting and finite vector splitting. Exact Riemann solvers such as the one proposed by Godunov [28] are too inefficient for the use in modern numerical schemes because they require an iterative procedure for each time step even in explicit schemes. The flux vector splitting methods are very robust and simple algorithms for the price of high numerical dissipation in boundary and shear layers. Prominent exponents of this class are the methods proposed by Steger and Warming [71] and van Leer [75] which could be straightforwardly implemented in an explicit as well as implicit framework, but they are too dissipative to be applied to viscous flow fields. Furthermore, these schemes are not able to preserve total enthalpy for steady flows. To overcome this difficulty, H¨anel et al. [33, 32] proposed a simple modification of van Leer’s scheme in the energy-flux that preserves the total enthalpy and improves the accuracy of the computed wall temperature distribution in viscous supersonic flow fields [33]. However, the massflux splitting was not changed and, therefore, the scheme also suffers from the disadvantage of large numerical dissipation on shear layers and contact discontinuities which occur in boundary layers. Schemes based on finite difference splitting such as the ones proposed by Roe [59] and Osher [55] lead to very accurate algorithms at the cost of increased complexity and computational effort, especially in an implicit method. Furthermore, the method of Roe computes the exact solution of a linearized Riemann problem that does not lead to a continuously differentiable numerical flux function, which is an unfavourable property for an implicit method. In 1991 Liou and Steffen [47] proposed a new flux splitting scheme called AUSM (advection upstream splitting method) which successfully combines the accuracy of the flux difference splitting schemes with the simplicity of the flux vector splitting methods. Thus, the method is able to capture a stationary contact discontinuity without numerical diffusion and avoids the carbuncle phenomenon. Furthermore, it is suitable for the use in an implicit algorithm, but unfortunately it suffers from some difficulties like oscillations in the pressure distribution behind a moving shock. Recently, Liou and Wada [78] presented an improved AUSM which they called AUSMDV. This method combines the so-called AUSMD, AUSMV, and H¨anel’s scheme, including an additional numerical dissipation to remove the glitches at the sonic point within rarefaction waves. “D” and “V” means that the corresponding scheme is flux difference or flux vector splitting biased, respectively. The method has desirable properties like no numerical dissipation on stationary as well as moving contact discontinuities, preserving enthalpy in steady flows, high resolution for shock waves and the scheme avoids pressure oscillations observed in the AUSM scheme. Furthermore, the resulting numerical flux function is continuously differentiable and can be extended to simulate chemically reacting [78, 34], as well as turbulent flow fields [52], and was recently modified for unsteady flow fields on moving grids [50]. Because of these favourable properties we use the AUSMDV scheme in an implicit version in the present numerical algorithm and, thus, we approximate the inviscid operator Lch by means of 2 ¢ ¢ 1 X X ¯¯ k ¯¯ c ¡ li j H uˆ i (t), uˆ j (t); nikj Lch u (t)|σi = |σi | j∈N (i) k=1

¡

(3.4)

DIFFERENT KRYLOV SUBSPACE METHODS

319

with the numerical flux function ¡ ¢ ¡ ¢ ¡ ¢ Hc uˆ i (t), uˆ j (t); nikj = T−1 nikj HAUSMDV uˆ i (t), uˆ j (t); nikj .

(3.5)

Herein, uˆ i (t) and uˆ j (t) represent the conserved variables at the midpoint xi jk of the line segment likj which are defined by uˆ i (t) = lim W −1 (w(x, t)), uˆ j (t) = lim W −1 (w(x, t)), x→xi jk x∈σi

x→xi jk x∈σ j

t where w(., t) ∈ Hh,m is a recovery function in the primitive variables (ρ, v1 , v2 , p) and W denotes a bijective mapping between the conserved and the primitive variables. A finite volume scheme always computes cell averages. Therefore, it seems to be natural t and use functions w(., t) ∈ to interprete those as piecewise constant functions u(., t) ∈ Hh,0 t Hh,0 defined by w(x, t) = W(u(., t))|σi , x ∈ σi for all σi ∈ 4. This procedure leads to a first-order scheme. The central idea in increasing the spatial accuracy lies in the construction of a recovery t , m > 0, such that on each control volume σi ∈ 4 the accuracy function w(., t) ∈ Hh,m condition w(x, t) − u(x, t) = O(h m+1 ) is valid. For the sake of simplicity, we choose polynomial recovery. In particular, we are interested in a second-order scheme in which linear t are sought. Thus, for each control volume σi we use a ansatz functions w(., t) ∈ Hh,1 polynomial

¡ ¢T wm∗ (., t)|σi = wm∗ (x1 , x2 , t) = C + ∇x wm (t)|σi · x1 − x1i0 , x2 − x2i0 ,

(3.6)

with C ∈ R, (x1i0 , x2i0 )T = xi0 ∈ σi and wm∗ ∈ {ρ, v1 , v2 , p}. In the context of CPU-time we have to find out the explicit form of the polynomial w∗ = (w1∗ , w2∗ , w3∗ , w4∗ )T by an efficient calculation from the known cell averages. Therefore, we insist on the piecewise constant function w ˜ m (t)|σi = (W(u(., t))|σi )m being the constant C in the formulation (3.6). Furthermore, we require that the linear ansatz function preserves the cell average which leads to Z 0=

σi

(wm∗ (., t)|σi − w ˜ m (t)|σi ) dx = ∇x wm (t)|σi

Z σi

(x − xi0 ) dx.

Because we consider nonvanishing gradients, it follows that xi0 has to be the barycenter of the box σi . Let T ∈ Th be an arbitrary triangle; then ∇x wT = (∇x w1T , . . . , ∇x w4T )T denotes the vector of the four gradients on T which can be uniquely defined by the interpolant of the cell averages of the three control volumes σi ∈ 4 with σi ∩ T 6= ∅. For the computation of each gradient, we introduce two different formulations depending on the type of triangulation (Thu or Ths ). First, we consider a control volume σi where the corresponding node xi satisfies the condition that there exists a triangle T ∈ Th with xi ∈ T . In this case we define the gradients by ∇x wm |σi =

1 X |T ∩ σi |∇x wmT |σi | T ∈T h

for m = 1, . . . , 4.

(3.7)

320

ANDREAS MEISTER

Otherwise, we use a central or, if necessary at the boundaries, a one-sided second-order difference for the computation of the gradients. It is well known that particularly in shock regions a straightforward implementation of this recovery procedure generates new extrema in the physical variables and therefore yields oscillations in the numerical approximation and often the numerical method fails. To avoid this behaviour, a function depending on the data is introduced to limit the gradient. At present, various limiters are available which satisfy the total variation diminishing (TVD) condition [46, 40]. We employ the slope limiter Lm σi described by Barth and Jespersen [5]. Consequently, the recovery polynomial is finally given in the form ¡ i jk ¡ i jk i jk ¢ i jk i0 i 0 ¢T ˜ m (t)|σi + L m . wm x1 , x2 , t = w σi ∇x wm (t)|σi · x 1 − x 1 , x 2 − x 2 3.2. Discretization of the Viscous Terms The numerical approximation of the viscous fluxes is performed in the form of a central difference. Analogously to the inviscid part we consider the formulation Z X 2 2 ¡ v ¢ 1 X X 1 f v (u(x, t))n ikj,` (t) ds Lh u (t)|σi = |σi | j∈N (i) k=1 Re∞ likj `=1 l of the dissipative part of Eq. (3.2). Let τikj be the triangle T which contains the line likj , then P the flux (1/Re∞ ) 2`=1 f`v (u)n ikj` over the edge likj of the control volume σi is evaluated in the barycenter of the triangle τikj . For each primitive variable φ ∈ {ρ, v1 , v2 , p} the gradients ∇x φ needed for the calculation of the viscous fluxes are given by the unique linear distribution of the variable φ on the triangle τikj in the form

(∇x φ)τikj

1 = det

µ

(φ2 − φ1 )(x3,2 − x1,2 ) − (φ3 − φ1 )(x2,2 − x1,2 ) (φ3 − φ1 )(x2,1 − x1,1 ) − (φ2 − φ1 )(x3,1 − x1,1 )

¶ ,

(3.8)

where det = (x2,1 − x1,1 )(x3,2 − x1,2 ) − (x3,1 − x1,1 )(x2,2 − x1,2 ) and the mean value φ P is determined by the arithmetic mean value (φ)τikj = 13 3m=1 φm , where φm represents the

mean averaged physical value at the box σm and xm = (xm,1 , xm,2 )T , m = 1, 2, 3, are the nodes of the triangle τikj . The combination of the described procedure leads directly to a formulation of the numerical viscous flux and, therefore, yields the approximation ¡

2 ¢ ¢ 1 X X ¯¯ k ¯¯ v ¡ li j H uτikj ,1 , uτikj ,2 , uτikj ,3 ; n ikj . Lvh u (t)|σi = |σi | j∈N (i) k=1

Note that a general second-order Gaussian quadrature rule for the approximation of the line integral requires the evaluation of the physical values at the midpoint of the line segment likj . Numerical investigations have shown that the described procedure leads to similar accuracy in the results and therefore we make use of the advantage that the chosen discretization avoids many additional calculations.

DIFFERENT KRYLOV SUBSPACE METHODS

321

3.3. Time Discretization For the description of the time stepping scheme, we integrate the weak form of the Navier–Stokes equations from time level t to time level t + 1t Z

t+1t t

Z

t+1t

∂t ui dt = t

Ã

¶ ! 2 Z µ X 1 v c ∂x f (u) − ∂x j f j (u) dx dt Re∞ j j j=1 σi

∀σi ∈ 4. (3.9)

With regard to the steady state computations it is not necessary to use a higher order time integration. Evaluating the right-hand side at time level t leads to an explicit scheme. However, these schemes are limited in CFL number by 1 which is a severe upper bound for simulations of viscous flows with high Reynolds numbers, as well as applications on adapted meshes where grid cells can be very small. One way to overcome this restriction is an implicit discretization where impressive acceleration, compared to an explicit time stepping scheme, can be achieved [49, 50]. Therefore, the right-hand side, and thus, the numerical flux function have to be evaluated at time level t + 1t. Hence, a linearization can be employed which leads to a linear system of equations in the form Aii 1ui +

X

Ai j 1u j = bi , i = 1, . . . , I,

j∈N (i)

where 1ui := ui (t + 1t) − ui (t) and Aii , Ai j ∈ R4×4 . Consequently, for each time step a linear system Ay = b

(3.10)

has to be solved, where A is a large, sparse, and nonsymmetric matrix. 4. ITERATIVE SOLVERS

The performance of an implicit method depends on the properties of the solver of the linear system of equations which has to be stable and fast. In the last few years a great variety of Krylov subspace methods were developed and used in different flow solvers. However, reliable statements concerning the performance and applicability of these methods in the context of current well-known flow problems are still missing. We will give a brief introduction to the Krylov subspace methods used in the following study. For a detailed description of these algorithm we refer to [4, 63, 51]. A projection method for solving Eq. (3.10) is an algorithm which seeks an approximate solution ym ∈ y0 + K m with regard to the condition (b − Aym ) ⊥ L m ,

(4.1)

where y0 ∈ Rn denotes an arbitrary vector and K m and L m represent m-dimensional subspaces of Rn . The orthogonality condition is given by the Euclidean inner product y ⊥ x ⇐⇒ (y, x)2 = 0. A Krylov subspace method belongs to the class of projection methods. Here K m represents the Krylov subspace given by K m = span{r0 , Ar0 , . . . , Am−1 r0 } with r0 = b − Ay0 .

322

ANDREAS MEISTER

At present, the leading Krylov subspace methods in current CFD applications are GMRES (generalized minimal residual), CGS (conjugate gradient squared), BiCGSTAB (biconjugate gradient stabilized), TFQMR (transpose-free quasi minimal residual), and QMRCGSTAB (quasi minimal residual conjugate gradient stabilized). The GMRES algorithm was developed by Saad and Schulz [65] in 1986. The method can be integrated in the above framework by choosing L m = {y ∈ Rn | ∃ x ∈ K m with y = Ax}. Thus, the residual satisfies the Petrov–Galerkin condition rm ⊥ span{Ar0 , A2 r0 , . . . , Am r0 } and can be written in the form rm = pm (A) r0 , where pm denotes a polynomial of degree m with pm (0) = I. With regard to the implementation of this method it is advantageous to consider GMRES as being based on a transformation of the system (3.10) into an equivalent minimization 2 problem. First, we define the function f : Rn → R+ 0 by f (y) = ||b − Ay||2 and choose an arbitrary initial vector y0 . Starting with m = 0 the residual r¯ m = miny=y0 +z,z∈K m f (y) is computed. Now m is increased until r¯ m is below a given accuracy. Then the optimal approximate solution ym = arg miny=y0 + z,z∈K m f (y) is computed. Both formulations are equivalent. A proof is given in [51]. Because of the expense of calculating the residual increase with the Krylov subspace dimension, it is efficient to limit this dimension. In the case when this limit is reached and the stopping criterion is still not satisfied the approximation ym is calculated and used as the initial value during a repetition. This technique is called “GMRES with restart.” Various numerical experiments indicate that a maximum allowable Krylov subspace dimension of 10 seems to be a good choice. Therefore, this upper bound is employed during all of the following computations. The CGS scheme is based on the BiCG method [21] which represents also a Krylov subspace method with a different Petrov–Galerkin condition L m = K mT := span{˜r0 , AT r˜ 0 , . . . , (AT )m−1 r˜ 0 }, where often r˜ 0 = r0 is taken. Similar to the above-mentioned formulation of the residual vector, BiCG constructs a sequence of residual vectors which can be written as rm = pm (A)r0 .

(4.2)

In comparison with the GMRES iteration, the BiCG scheme includes the advantage that the approximate solution can be computed by a three-term recurrence relation and thus the required amount of storage does not grow with m. On the other hand and in contrast to the GMRES algorithm, BiCG suffers from several breakdown possibilities for general nonsymmetric matrices and requires multiplications with AT . Note that in the case of symmetric matrices and r˜ 0 = r0 , BiCG reduces to the CG method. Sonneveld [70] observed that the second disadvantage can be avoided by a simple modification of the BiCG algorithm. Therefore, one has to use the polynomial pm arising from the BiCG scheme and replace (4.2) by rm = pm2 (A)r0 .

(4.3)

This procedure leads also to a three-term recurrence relation and avoids multiplications by AT .

DIFFERENT KRYLOV SUBSPACE METHODS

323

As was published by many authors [23, 29, 53, 72] for special systems of equations, the CGS algorithm tends to oscillations and sometimes diverges while other iterative solvers converge. In order to damp this oscillatory behaviour van der Vorst [73] suggested a variant of the CGS method called BiCGSTAB. Here Eq. (4.3) has to be replaced by rm = pm (A)qm (A)r0 , where pm and qm represent polynomials of degree m, respectively. Therefore, an additional degree of freedom is available, which is used for a concluding minimization of the residual. To overcome the breakdown possibilities Freund and Nachtigal [24] combined the BiCG scheme with the minimization strategy used in the GMRES method, which results in a BiCG-like approach called QMR. The method uses a look-ahead Lanczos algorithm to compute a basis of K m . In contrast to the GMRES iteration where the Arnoldi algorithm is used, the Lanczos algorithm incorporated in the BiCG scheme produces two sequences of vectors which are orthogonal to each other, whereby the calculated basis of the Krylov subspace is nonorthogonal, in general. A complete minimization of the function f employed in the GMRES method leads to trouble since the expense to compute the minimal rm grows with 4I m 2 , where I is the number of grid points (A ∈ R4I ×4I ). The algorithm proposed by Freund and Nachtigal uses a quasi minimization of f , whereby nonoptimal approximate solutions ym ∈ y0 + K m are achieved with the advantage of low cost. To derive the TFQMR algorithm [29] we consider the QMR method and replace the contained basis by a new one built by using vectors computed in the CGS iteration. Hence, the TFQMR algorithm inherits the advantage of the CGS method that the multiplication by AT is avoided. Analogously, Chan et al. [13] developed the QMRCGSTAB algorithm. Motivated by the procedure to define the TFQMR method from the CGS scheme they constructed a basis of K m by vectors used in BiCGSTAB. Some numerical experiments [13] show that the advantage of smoother convergence behaviour obtained by using BiCGSTAB instead of CGS also carries over from TFQMR to QMRCGSTAB. 4.1. Preconditioning One of the main acceleration techniques for iterative methods is the use of a preconditioning matrix. As with the above-presented iterative solvers a great variety of different preconditioners is also available such as scaling, incomplete formulations of the Cholesky, LU, or LQ factorization and some modifications like fill-in or threshold strategies. Furthermore, an iterative solver itself (i.e., SOR or SSOR), a multilevel or a multigrid method can be used as a preconditioner. A detailed presentation of various techniques is given in [62, 63]. In contrast to the iterative methods, a lot of these preconditioning techniques are too expensive or not even practical for the described framework. For example, the incomplete Cholesky factorization and a modified version are investigated as a preconditioner for the CG and CGS method for two model problems by van der Vorst [72], but obviously these preconditioning techniques are inherently restricted to linear systems of equations with a symmetric matrix. Multilevel preconditioners described by Bramble et al. in a series of publications [6–9] suffer from the same problem, namely that they are developed in

324

ANDREAS MEISTER

the background of elliptic problems and therefore restricted to symmetric positive definite matrices. A multigrid approach can certainly be employed as a preconditioner but it requires a significant additional effort in the development with regard to unstructured meshes, especially in three space dimensions. Therefore, a multigrid method is usually used as a solver, instead of a preconditioner. Note that the presented scheme can be used for the simulation of unsteady flow fields (usually by using a second-order backward difference formula (BDF(2)) or a Crank–Nicolson time discretization instead of a simple backward Euler approach) [51]. The standard multigrid approach [31, 10] is only applicable for the simulation of steady flow fields or unsteady fluid flows where the solution is sufficiently smooth [48] such as rarefaction waves but certainly not for moving shocks or contact discontinuities. For unsteady flow fields, like transonic flows about pitching airfoils, we therefore suggest an extension of the presented implicit scheme as proposed in [50, 51, 26] or the dual time stepping method [27, 2]. The incomplete LQ preconditioner suffers from the property that sometimes the algorithm produces a singular matrix Q [60]. A criterion to ensure the nonsingularity is unfortunately often in conflict with the sparsity requirement [62]. Concerning this difficulty the ILQ factorization is not very popular and therefore not implemented in the presented scheme. The most popular and successful preconditioner for large, sparse, nonsymmetric, and indefinite matrices is the incomplete LU factorization which is a pair of a lower left (L) and a upper right (U) matrix satisfying the following three conditions: 1. Uii presents the unit matrix for all i, 2. Li j = Ui j = Ai j , if Ai j is a null matrix, 3. (LU)i j = Ai j , if Ai j is not a null matrix. We employ a right preconditioning and therefore the linear system is transformed into AU−1 L−1 y˜ = b, y = U−1 L−1 y˜ . Considering these three requirements, one can easily see that the structure of L + U is identical with that of A and therefore the ILU preconditioner requires the same storages as the matrix A. Because the use of additional fill-in elements inside the ILU factorisation is mainly efficient for band matrices and does not extend to general sparse matrices, Saad [62] proposed a threshold strategy, but, as has been mentioned by him, this procedure may fail for many of the matrices that arise from real applications. Therefore, Saad extended this method with a pivoting which requires an additional computational effort. Thus, the only modification of the standard ILU preconditioning one can propose in the current framework is a frozen ILU technique, where the expensive calculation of the ILU factorization is enforced only after either a fixed or variable number of time steps [49], whereby the optimal number depends on the underlying flow problem. Numerical experiments indicate that the acceleration obtained by the frozen formulation of the ILU preconditioner does not depend on the chosen Krylov subspace method and, thus, during the comparison we focus our view on the nonfrozen ILU preconditioner. In order to decrease the amount of storage, especially in the case of complex twodimensional or three-dimensional flow fields, we furthermore investigate the behaviour of a simple scaling as a preconditioner which obviously can be used without any additional

325

DIFFERENT KRYLOV SUBSPACE METHODS

TABLE 1 Inviscid Test Cases Configuration

Ma∞

Angle of attack

References

Gridpoints

Bi-NACA0012 double airfoil Combustion chamber

0.55 2.0

6◦ 0◦

[15] [79]

13577 31202

amount of storage. This technique is given by a left preconditioning, S−1 Ay = S−1 b, where S = blockdiag(Aii ), i = 1, . . . , I . 5. NUMERICAL RESULTS

The following test cases show the performance and applicability of the presented numerical scheme for the simulation of distinguished flow fields and enables the investigation of the performance of the different preconditioned Krylov subspace methods. We consider a special selection of test cases to take account of the great variety of physically and numerically interesting flow phenomena. The first two test cases present a transonic flow about a double airfoil and a supersonic flow in a combustion chamber. Furthermore, we consider three viscous test cases. Here, the influence of the viscosity depends directly on the Reynolds number. The variation of this nondimensional number enables the investigation of viscous dominated, as well as convection dominated, flow problems. As in the case of inviscid flows (Table 1) the viscous test cases are summarized in Table 2. 5.1. Bi-NACA0012 Double Airfoil The first test case is concerned with the inviscid calculation of a transonic steady state flow about the Bi-NACA0012 double airfoil at a freestream Mach number of Ma∞ = 0.55 and an angle of attack of α = 6◦ . The test case was taken from [15] and exhibits the advantage that a flow about as well as between the two airfoils can be simultaneously investigated. The unstructured grid includes 13,577 nodes and 26,632 triangles. The density distribution obtained with the presented implicit method is shown in Fig. 4, and Fig. 3 shows the grid

TABLE 2 Viscous Test Cases Configur.

Ma∞

Angle of attack

NACA0012 airfoil Cylinder Flat plate

0.85

0◦

500

273 K

8.03 5.0

0◦

193.750 6 × 106

122.1 K 80 K

Re∞

Tˆ∞

Tˆw

294.4 K 288 K

References

Gridpoints

[20]

8742

[81, 41, 17] [16]

27120 7350

326

ANDREAS MEISTER

FIG. 3. Partial view of the triangulation.

used. Here we can see a good resolution of the dominant flow features which are the shock between the airfoils and the upper shock. 5.2. Combustion Chamber This test case shows the inviscid flow about a wedge in a channel. The configuration was taken from an experimental investigation of a supersonic combustion chamber [79]. The freestream Mach number is Ma∞ = 2.0 and a horizontal inflow is considered. Thus, the flow is characterized by strong shock–shock and shock–boundary interactions. Furthermore, an expansion fan with steep pressure, density, and velocity gradients takes place at the backward edge of the wedge. The presented density distribution (Fig. 6) was calculated on the adapted grid shown in Fig. 5. Figure 6 demonstrates that the shock–shock and the shock–boundary interactions are sharply resolved and the rarefaction waves are clearly captured. 5.3. NACA0012 Airfoil A steady laminar flow about a NACA0012 airfoil was chosen as the first case to test the described numerical method for viscous flow computations. Adiabatic walls with a reference temperature of Tˆw = 273 K are considered. The further freestream conditions are a Reynolds number of Re∞ = 500, a Mach number of Ma∞ = 0.85, and an angle of attack of α = 0◦ .

DIFFERENT KRYLOV SUBSPACE METHODS

FIG. 4. Isolines of the density.

FIG. 5. Triangulation of the channel.

FIG. 6. Density distribution.

327

328

ANDREAS MEISTER

FIG. 7. Partial view of the triangulation.

The obtained Mach number distribution is shown in Fig. 8, and Fig. 7 presents the adapted grid. For the validation of the numerical results we consider the C p -distribution at the surface of the airfoil. The curves in Fig. 9 show a good agreement between the numerical results of the described scheme in comparison with the ViB-Code [39], which is a cell-centered Jameson-type scheme on structured grids. 5.4. Shock–Shock Interaction This test case shows the numerical simulation of the flow over the inlet of a scramjet. The flow was measured by Wieting and Holden et al. [81, 41] in 1987 and 1988. The onflow Mach number was chosen to be Ma∞ = 8.03, so that a strong bow shock in front of the intake cowl occurs. Generated by the intake ramp an impinging shock hits the bow shock at an angle of 18.11◦ . Furthermore, the Reynolds number computed with the radius of the cylinder (r = 38.1 mm) is Re∞ = 193.750 and we consider isothermal walls, where the wall temperature was given to be Tˆw = 294.4 K and the reference temperature was chosen to be Tˆ∞ = 122.1 K. Edney [17] subdivided this shock–shock interaction into six groups, depending on the angle where the impinging shock hits the bow shock. Concerning this subdivision the investigated flow belongs to type IV. The physically interesting part of the flow is the pressure and heat flux distribution along the surface of the cylinder. The pressure, as well as the heat flux in the stagnation point, are about 10 times higher than those which arise in the same

DIFFERENT KRYLOV SUBSPACE METHODS

329

FIG. 8. Mach number distribution.

flow without an impinging shock. Figure 10 shows a partial view of the grid used and Fig. 11 presents the obtained Mach number distribution. Furthermore, the curves (Figs. 12 and 13) quantify the heat flux and pressure distribution along the surface of the cylinder, whereby the numerical results are divided by the heat flux q0 and pressure p0 at the stagnation point in a flow without an impinging shock, respectively. Both numerical results are in very good agreement with the experimental data.

FIG. 9. C p -Distribution at the surface of the airfoil.

330

ANDREAS MEISTER

FIG. 10. Partial view of the triangulation.

5.5. Flat Plate Boundary Layer In addition to the mentioned viscous flow fields at low and middle Reynolds numbers, we close the series of test cases with a laminar flow about a flat plate at a high Reynolds number of Re∞ = 6 × 106 , which was built with a reference length of one meter. The flat plate has a constant wall temperature of Tˆw = 288 K and a length of 35 cm. Furthermore we consider a reference temperature of Tˆ∞ = 80 K and the freestream Mach number and the angle of attack were chosen to be Ma∞ = 5.0 and α = 0◦ , respectively. For validation of the numerical results we compare the obtained skin friction coefficient (Fig. 14) and the Stanton number (Fig. 15) with simulations presented in [16]. Both figures show good agreement with the other computational results. 5.6. Convergence Histories and Discussion For the comparison of the convergence histories the weighted discrete L 2 norm of the density residual was taken. This convergence criterion is often used in numerical methods for the Euler and Navier–Stokes equations because it is equal to zero when the steady state solution is reached. Implicit methods are not restricted by the stability condition of Courant, Friedrich, and Lewy. Nevertheless, a time-stepping strategy is necessary since

DIFFERENT KRYLOV SUBSPACE METHODS

331

FIG. 11. Mach number distribution.

the incorporated linearization of the numerical flux function introduces an error of order O((1u)2 ). Hence, at the beginning of each time step we calculate the upper bound for an explicit scheme and multiply this number by a factor η which depends on the investigated flow field and is adapted during the computation to achieve fast convergence. For the first time step η was always chosen to be 1. With regard to the physical flow phenomena we increase η between a factor of 2 after each fiftieth time step (test case 5.2) and a factor of 2 after each third time step (test case 5.3). The maximum η lay between 250 (test case 5.2) and 105 (test case 5.5). The calculations were done on a Silicon Graphic Power Challenge XL. For the solution of the linear systems of equations we always took y0 = 0 as initial guess and the required tolerance ε was chosen to be ε = 10−6 . It has been proven that this tolerance is sufficiently small to guarantee a good resolution of the arising physical phenomena inside the considered test cases. Many linear systems of equations have to be solved until the numerical result represents the physical flow phenomena. Since the obtained linear systems can be extremely different from one time step to another, the following convergence curves show the convergence histories versus CPU-time, instead of convergence history for the Krylov subspace methods with reference to a special linear system which is generated during the calculation (see Figs. 16–37). Hence, for all Krylov subspace methods considered, excluding restarted

FIG. 12. Comparison of the heat flux distribution along the surface of the cylinder.

FIG. 13. Comparison of the pressure distribution along the surface of the cylinder.

FIG. 14. Distribution of the skin friction coefficient C f .

332

DIFFERENT KRYLOV SUBSPACE METHODS

333

FIG. 15. Distribution of the Stanton number St.

GMRES, a steeper convergence curve indicates a smaller product built by the time needed for one iteration step with the number of iterations needed to achieve the required tolerance ε. The beam diagrams on the right-hand side quantify the percentage comparison of the CPU-time needed. The results always correspond with the horizontal line in the bystanding figure of the convergence histories, whereby all percentage data are achieved by dividing each CPU-time by the CPU-time of the fastest algorithm. A missing beam implies that the corresponding iterative solver could not achieve the required density residual in an acceptable time. At first, we consider the presented diagrams concerning the simulation where an ILUpreconditioning of the Krylov subspace method was used. With regard to the CPU-time we can subdivide the five methods into three groups.

FIG. 16. Convergence histories for the ILU-preconditioned schemes. Test case 5.1.

FIG. 17. Percentage comparison of the CPU-time.

FIG. 18. Convergence histories for the schemes with scaling as a preconditioner. Test case 5.1.

FIG. 19. Percentage comparison of the CPU-time.

334

DIFFERENT KRYLOV SUBSPACE METHODS

335

FIG. 20. Convergence histories for the schemes without preconditioning. Test case 5.1.

The first group contains only the BiCGSTAB method. This scheme turned out to be the most efficient method in all calculations. CGS, QMRCGSTAB, and GMRES are summarized in the second group. These schemes required between 1.8% and 26.3% more CPUtime, compared with the BiCGSTAB method. For all investigated test cases the TFQMR scheme represented the worst choice. The expense of this scheme is about 31.8% up to 48.8% higher than in the BiCGSTAB method. Concerning viscous as well as inviscid fluid flows, the results indicate that, in combination with an ILU-preconditioner, the BiCGSTAB should be preferred and the TFQMR should be avoided. Next we consider the results obtained by using scaling as a preconditioner. Compared to the ILU-preconditioner, we get a completely different order of the Krylov subspace methods. Instead of the BiCGSTAB algorithm which represents the optimal choice in combination

FIG. 21. Percentage comparison of the CPU-time.

FIG. 22. Convergence histories for the ILU-preconditioned schemes. Test case 5.2.

FIG. 23. Percentage comparison of the CPU-time.

FIG. 24. Convergence histories for the ILU-preconditioned schemes. Test case 5.3.

336

DIFFERENT KRYLOV SUBSPACE METHODS

337

FIG. 25. Percentage comparison of the CPU-time.

with the ILU-preconditioning for all investigated test cases, the TFQMR always yields the fastest numerical scheme when scaling is used. The CGS, BiCGSTAB, and QMRCGSTAB methods required between 25.7% up to 99.7% more CPU-time (see Figs. 18, 26, 34). Furthermore, the results indicate a dramatic dependency of the GMRES scheme on the chosen preconditioner. While the method showed an acceptable convergence behaviour in connection with an ILU-preconditioner the performance of the scheme extremely decreased in most cases when scaling was used. Finally, we consider the Figs. 20, 28, and 36 which represent the convergence histories of the iterative solvers without preconditioning. In general, the schemes required ten times more CPU-time for the simulation of the flow field. As we have seen in the case of scaling the GMRES method needs distinctly more CPU-time.

FIG. 26. Convergence histories for the schemes with scaling as a preconditioner. Test case 5.3.

FIG. 27. Percentage comparison of the CPU-time.

FIG. 28. Convergence histories for the schemes without preconditioning. Test case 5.3.

FIG. 29. Percentage comparison of the CPU-time.

338

FIG. 30. Convergence histories for the ILU-preconditioned schemes. Test case 5.4.

FIG. 31. Percentage comparison of the CPU-time.

FIG. 32. Convergence histories for the ILU-preconditioned schemes. Test case 5.5.

339

FIG. 33. Percentage comparison of the CPU-time.

FIG. 34. Convergence histories for the schemes with scaling as a preconditioner. Test case 5.5.

FIG. 35. Percentage comparison of the CPU-time.

340

DIFFERENT KRYLOV SUBSPACE METHODS

341

FIG. 36. Convergence histories for the schemes without preconditioning. Test case 5.5.

6. CONCLUSION

An efficient implicit finite volume method for the simulation of inviscid and viscous flow fields has been developed. The use of unstructured grids enables a specific discretization of arbitrarily complex geometries and the employment of time-accurate adaptation techniques. The described combination of isotropic triangulations with non-isotropic subgrids leads to a distinct improvement of the boundary layer resolution in the case of high Reynolds number flows. For every time step, the implicit formulation leads to a linear system of equations with a large, sparse, badly conditioned, and nonsymmetric matrix. The different preconditioned Krylov subspace methods prove to be efficient for the solution of the arising linear system of equations. In order to get a fast and stable scheme it is necessary to choose a suitable combination of a preconditioner and a Krylov subspace method in the context

FIG. 37. Percentage comparison of the CPU-time.

342

ANDREAS MEISTER

of the available resources. Carried over to the investigated flow problems, this means that in the case of sufficient storage resources, the ILU-preconditioned BiCGSTAB method has to be preferred. In the presented test cases the additional amount of storage of the ILUpreconditioner varies between 8 and 32 MB, but on finer grids or in three-dimensional calculations it can easily grow up to more than 500 MB. In order to minimize the amount of storage the numerical investigation indicated the superiority of a scaling in combination with the TFQMR algorithm. Hence, the proposed scheme represents a fast solver for inviscid and viscous flow problems and an excellent basis for the extension to turbulent [52] and chemical reacting flow fields [34] as well as moving grids [26, 50, 51]. A third-order essentially nonoscillatory (ENO) scheme based on this approach was described in [38]. ACKNOWLEDGMENTS The author thanks Oliver Friedrich for providing the grid generator [25] and Daniel Hempel for providing the adaptation algorithms [36, 35]. The author furthermore wishes to express his thanks to Thomas Sonar for suggesting several improvements during many helpful discussions.

REFERENCES 1. K. Ajmani, W.-F. Ng, and M. Liou, Preconditioned conjugate gradient methods for the Navier–Stokes equations, J. Comput. Phys. 110, 68 (1994). 2. A. Arnone, M.-S. Liou, and L. A. Povinelli, Integration of Navier–Stokes equations using dual time stepping and a multigrid, AIAA J. 33(6), 985 (1995). 3. R. E. Bank and T. F. Chan, An analysis of the composite step biconjugate gradient method, Numer. Math. 66, 295 (1993). 4. R. Barrett, M. Berry, T. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. van der Vorst, Templates for the Solution of Linear Systems (SIAM, Philadelphia, 1994). 5. T. J. Barth and D. C. Jesperson, The design and application of upwind schemes on unstructured meshes, AIAA Paper 89-0366, 1989. 6. J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring I, Math. Comp. 47, 103 (1986). 7. J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring II, Math. Comp. 49, 1 (1987). 8. J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring III, Math. Comp. 51, 415 (1988). 9. J. H. Bramble, J. E. Pasciak, and A. H. Schatz, The construction of preconditioners for elliptic problems by substructuring IV, Math. Comp. 53, 1 (1989). 10. A. Brandt, Guide to multigrid development, in Multigrid Methods, edited by W. Hackbusch and U. Trottenberg (Springer-Verlag, Berlin/Heidelberg/New York, 1981), p. 220. [Lecture Notes in Math., Vol. 960] 11. P. N. Brown, A theoretical comparison of the Arnoldi and GMRES algorithms, SIAM J. Sci. Stat. Comp. 12(7), 58 (1991). 12. X. Cai, D. E. Keyes, and V. Venkatakrishnan, Newton-Krylov-Schwarz: An Implicit Solver for CFD, ICASE Report 95-87, Hampton, Virginia, 1995. 13. T. F. Chan, E. Gallopoulos, V. Simoncini, T. Szeto, and C. H. Tong, A quasi-minimal residual variant of the bi-cgstab algorithm for nonsymmetric systems, SIAM J. Sci. Comput. 15(2), 338 (1994). 14. R. Courant, K. O. Friedrichs, and H. Lewy, Uber die partiellen Differentialgleichungen der mathematischen Physik, Math. Ann. 100, 32 (1928). 15. A. Dervieux, B. van Leer, J. Periaux, and A. Rizzi (Eds.), Numerical Simulation of Compressible Euler Flows (Vieweg, Braunschweig/Wiesbaden, 1989). [Notes on Numerical Fluid Mechanics, Vol. 26]

DIFFERENT KRYLOV SUBSPACE METHODS

343

16. J.-A. D´esid´eri, R. Glowinski, and J. P´eriaux (Eds.), Hypersonic Flows for Reentry Problems, number II (Springer, Antibes, 1990). 17. B. Edney, Anomalous Heat Transfer and Pressure Distributions on Blunt Bodies at Hypersonic Speeds in the Presence of an Impinging Shock, FFA Report 115, Stockholm, 1968. 18. W. S. Edwards, L. S. Tuckerman, R. A. Friesner, and D. C. Sorensen, Krylov methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 110, 82 (1994). 19. H. C. Elman, Iterative methods for linear systems, in Large-Scale Matrix Problems and the Numerical Solution of Partial Differential Equations, edited by J. Gilbert and D. Kershaw (Clarendon Press, Oxford, 1994), p. 69. 20. L. Fezoui, S. Lanteri, B. Larrouturou, and C. Olivier, Resolution Numerique des Equations de Navier–Stokes pour un Fluide Compressible en Maillage Triangulaire, INRIA Rapports de Recherche No. 1033, Sophia Antipolis, 1989. 21. R. W. Fletcher, Conjugate gradients methods for indefinite systems, in Dundee Biennial Conference on Numerical Analysis, edited by G. A. Watson (Springer, New York, 1975), p. 73. 22. R. W. Freund, Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices, SIAM J. Sci. Stat. Comp. 13(1), 425 (1992). 23. R. W. Freund, A transpose-free quasi-minimal residual algorithm for non-Hermitian linear systems, SIAM J. Sci. Comp. 14(2), 470 (1993). 24. R. W. Freund and N. M. Nachtigal, QMR: A quasi-minimal residual method for non-Hermitian linear systems, Numer. Math. 60, 315 (1991). 25. O. Friedrich, A new method for generating inner points of triangulations in two dimensions, Comp. Meth. Appl. Mech. and Eng. 104, 77 (1993). 26. O. Friedrich, D. Hempel, A. Meister, and Th. Sonar, Adaptive computation of unsteady flow fields with the DLR-T -code, in 77th AGARD Fluid Dynamics Panel Meeting and Symposium on Progress and Challenges in CFD Methods and Algorithms, Sevilla, 1995. [AGARD-CP-578] 27. A. L. Gaitonde, A dual-time method for the solution of the unsteady Euler equations, The Aeronautical Journal of the Royal Aeronautical Society, 98, 283 (1994). 28. S. K. Godunov, A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics, Mat. Sb. 47, 357 (1959). 29. G. Golub, A. Greenbaum, and M. Luskin (Eds.), Transpose-Free Quasi-Minimal Residual Methods for NonHermitian Linear Systems (Springer-Verlag, New York, 1994). [The IMA Volumes in Mathematics and its applications, Recent Advances in Iterative Methods, Vol. 60] 30. M. H. Gutknecht, Variants of BiCGSTAB for matrices with complex spextrum, SIAM J. Sci. Comp. 14(5), 1020 (1993). 31. W. Hackbusch, Multi-Grid Methods and Applications (Springer-Verlag, Berlin/Heidelberg/New York/Tokyo, 1985). [Computational Mathematics, Vol. 4] 32. D. H¨anel and R. Schwane, An implicit flux-vector splitting scheme for the computation of viscous hypersonic flow, AIAA Paper 89-0274, 1989. 33. D. H¨anel, R. Schwane, and G. Seider, On the accuracy of upwind schemes for the solution of the Navier–Stokes equations, AIAA Paper 87–1105, 1987. 34. V. Hannemann, Numerische Simulation von Stoss-Stoss-Wechselwirkungen unter Ber¨ucksichtigung von chemischen und thermischen Nichtgleichgewichtseffekten, Ph.D. thesis, Universit¨at G¨ottingen, 1996. 35. V. Hannemann, D. Hempel, and Th. Sonar, Dynamic adaptivity and residual control in unsteady compressible flow computation, Mathematical and Computer Modelling 20, 201 (1994). 36. D. Hempel, Isotropic refinement and recoarsening in two dimensions, Numerical Algorithms 13, 33 (1996). 37. M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, NBS J. Res. 49, 409 (1952). 38. D. Hietel, A. Meister, and Th. Sonar, On the comparison of four different implementations of an third-order ENO scheme of box type for the computation of compressible flow, Numerical Algorithms 13, 77 (1996). 39. A. Hilgenstock, The ViB-Code to Simulate 3-D Stator/Rotor Flow in Axial Turbines, DLR-Forschungsbericht 92-19, G¨ottingen, 1992.

344

ANDREAS MEISTER

40. C. Hirsch, Numerical Computation of Internal and External Flows, Vol. 2 (Wiley, Chicester/New York, 1988). 41. M. S. Holden, A. R. Wieting, J. R. Moselle, and C. Glass, Studies of aerothermal loads generated in regions of shock/shock interaction in hypersonic flow, AIAA 26th Aerospace Sciences Meeting, January 1988. [AIAA88-0477] 42. E. Issman and G. Degrez, Implicit iterative methods for a multidimensional upwind Euler/Navier–Stokes solver on unstructured grids. Von Karman Institute for Fluid Dynamics, Lecture Series (1995-14), RhodeSaint-Gen`ese, March 1995. 43. A. Jameson, Computational transonics, Comm. Pure Appl. Math. 41, 507 (1988). 44. W. Joubert, Lanczos methods for the solution of nonsymmetric systems of linear equations, SIAM J. Matrix Anal. Appl. 13(3), 926 (1992). 45. E. M. Kasenally, Gmback: A generalised minimum backward error algorithm for nonsymmetric linear systems, SIAM J. Sci. Comp. 16(3), 698 (1995). 46. R. J. LeVeque, Numerical Methods for Conservation Laws (Birkh¨auser, Basel, 1990). 47. M.-S. Liou and C. J. Steffen, A New Flux Splitting Scheme, NASA Technical Memorandum 104404, Cleveland, Ohio, May 1991. [J. Comput. Phys. 107, 23 (1993)] 48. A. Meister, Zur Verwendung der Mehrgittermethode bei der zeitgenauen L¨osung der instation¨aren Eulergleichungen, DLR-IB 221-93 A 20, G¨ottingen, 1993. 49. A. Meister, Ein Beitrag zum DLR-T -Code: Ein explizites und implizites Finite-Volumen-Verfahren zur Berechnung instation¨arer Str¨omungen auf unstrukturierten Gittern, DLR-IB 221-94 A 36, G¨ottingen, 1994. 50. A. Meister, Development of an implicit finite volume scheme for the computation of unsteady flow fields on unstructured moving grids, in Numerical Methods for Fluid Dynamics V, edited by K. W. Morton and M. J. Baines (Clarendon Press, Oxford, 1995), p. 481. 51. A. Meister, Zur zeitgenauen numerischen Simulation reibungsbehafteter, kompressibler, turbulenter Str¨omungsfelder mit einer impliziten Finite-Volumen-Methode vom Box-Typ, Ph.D. thesis, Technische Hochschule Darmstadt, 1996. 52. A. Meister and M. Oevermann, Computation of laminar and turbulent compressible flow fields on unstructured grids with an implicit finite volume scheme, in Proceedings of the 2nd Seminar on Euler and Navier–Stokes Equations, Prague, 1996. 53. N. M. Nachtigal, S. C. Reddy, and L. N. Trefethen, How fast are nonsymmetric matrix iterations? SIAM J. Matrix Anal. Appl. 13(3), 778 (1992). 54. M. Nagaoka and N. Horinouchi, Implicit compressible flow solvers on unstructured meshes, Comput. Fluid Dynamics J. 2(2), 169 (1993). 55. S. Osher and F. Solomon, Upwind schemes for hyperbolic systems of conservation laws, Math. Comput. 38, 339 (1982). 56. C. Pommerall and W. Fichtner, memory aspects and performance of iterative solvers, SIAM J. Sci. Comput. 15 (2), 460 (1994). 57. J. J. Quirk, A contribution to the great riemann solver debate, ICASE Report 92, Hampton, Virginia, 1992. 58. G. Radicati, Y. Robert, and S. Succi, Iterative algorithms for the solution of nonsymmetric systems in the modelling of weak plasma turbulence, J. Comput. Phys. 80, 489 (1989). 59. P. L. Roe, Approximate riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43, 357 (1981). 60. Y. Saad, Preconditioning techniques for nonsymmetric and indefinite linear systems, Journal of Comp. and Appl. Math. 24, 89 (1988). 61. Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM J. Sci. Comp. 14(2), 461 (1993). 62. Y. Saad, Krylov subspace techniques, conjugate gradients, preconditioning and sparse matrix solvers, von Karman Institute of Fluid Dynamics, Lecture series (1994-05), 1994. 63. Y. Saad, Iterative Methods for Sparse Linear Systems (PWS, Boston, 1996). 64. Y. Saad and M. H. Schultz, Conjugate gradient-like algorithms for solving nonsymmetric linear systems, Math. Comput. 44, 417 (1985).

DIFFERENT KRYLOV SUBSPACE METHODS

345

65. Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comp. (7), 856 (1986). 66. M. A. Saunders, H. D. Simon, and E. L. Yip, Two conjugate-gradient-type methods for unsymmetric linear equations, SIAM J. Numer. Anal. 25(4), 927 (1988). 67. D. C. Slack, D. L. Whitaker, and R. W. Walters, Time integration algorithms for the two-dimensional Euler equations on unstructured meshes, AIAA J. 32(6), 1158 (1994). 68. G. L. G. Sleijpen and D. R. Fokkema, Bicgstab(`) for linear systems involving unsymmetric matrices with complex spectrum, Electronic Transactions on Numerical Analysis 1, 11 (1993). 69. Th. Sonar, On the design of an upwind scheme for compressible flow on general triangulations, Numerical Algorithms 4, 135 (1993). 70. P. Sonneveld, CGS: A fast Lanczos-type solver for nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 10(1), 36 (1989). 71. J. L. Steger and R. F. Warming, Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods, J. Comput. Phys. 40, 263 (1981). 72. H. A. van der Vorst, The convergence behaviour of preconditioned CG and CG-S in the presence of rounding errors, in Lecture Notes in Mathematics, Vol. 1457 (Springer-Verlag, New York/Berlin, 1990), p. 126. 73. H. A. van der Vorst, Bi-cgstab: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13, 631 (1992). 74. H. A. van der Vorst and K. Dekker, Conjugate gradient type methods and preconditioning, J. Comp. Appl. Math. 24, 73 (1988). 75. B. van Leer, Flux-vector splitting for the Euler equations, in Eighth International Conference on Numerical Methods in Fluid Dynamics, edited by E. Krause (Springer-Verlag, Berlin, 1982), p. 507. [Lecture Notes in Physics, Vol. 170] 76. V. Venkatakrishnan, Implicit Schemes and Parrallel Computing in Unstructured Grid CFD, ICASE Report 95-28, Hampton, Virginia, 1995. 77. V. Venkatakrishnan and D. J. Mavriplis, Implicit solver for unstructured meshes, J. Comput. Phys. 105, 83 (1993). 78. Y. Wada and M.-S. Liou, A flux splitting scheme with high-resolution and robustness for discontinuities, AIAA Paper 94-0083 (1994). 79. W. Waidmann, F. Alff, M. B¨ohm, U. Brummund, W. Clauss, and M. Oschwald, Supersonic combustion of hydrogen/air in a SCRAMJET combustion chamber, 45th Conference of the International Astronautical Federation, 1994. [IAF-94-S.4.429] 80. H. F. Walker, Implementation of the GMRES method using householder transformations, SIAM J. Sci. Comp. 9(1), 152 (1988). 81. A. R. Wieting and M. S. Holden, Experimental study of shock wave interference heating on a cylindrical leading edge at mach 6 and 8, AIAA 22nd Thermophysics Conference, June 1987. [AIAA-87-1511] 82. J. Xu and X. Cai, A preconditioned GMRES method for nonsymmetric or indefinite problems, Math. Comp. 59(200), 311 (1992).