A relaxation method for large eigenvalue problems, with an application to flow stability analysis

https://doi.org/10.1016/j.jcp.2012.01.038Get rights and content

Abstract

Linear stability analysis of fluid flows usually involves the numerical solution of large eigenvalue problems. We present a spectral transformation allowing the computation of the least stable eigenmodes in a prescribed frequency range, based on the filtering of the linearized equations of motion. This “shift-relax” method has the advantage of low memory requirements and is therefore suitable for large two- or three-dimensional problems. For demonstration purposes, this new method is applied to compute eigenmodes of a compressible jet.

Introduction

The stability of fluid flow is a fundamental question in fluid dynamics, which has significant implications on the design, operation and control of flow devices. Consequently, hydrodynamic stability theory has taken a central role in fluid dynamics research, and remarkable progress has been made over the past decades. Early investigations of generic flow configurations, such as channel flows or boundary layers, have recently given way to more complex two- and three-dimensional geometries and more complex flow physics. The resulting eigenvalue problems from simple configurations were sufficiently small to allow their solution by direct techniques (such as the QR-algorithm). For more complex stability problems, however, direct techniques no longer provide a feasible solution; iterative eigenvalue algorithms have to be employed to compute a subset of the full spectrum that effectively describes the essential dynamics of small disturbances superposed on a steady base flow. Many of these algorithms for large-scale eigenvalue problems have been developed within the linear-algebra community and are readily available through several public domain libraries [12]. Among them, the two most commonly used in fluid-flow problems are the Arnoldi [16] and the Krylov–Schur [23] algorithms. These methods can, in principle, extract any portion of the full spectrum through the repeated application of a restarting step. In practice, however, only the dominant modes (i.e. those associated with the eigenvalues of largest magnitude) can be computed in many applications; yet, these modes may not provide relevant information about the physics of the problem and, in some cases, may even be spurious. As far as the asymptotic stability behavior is concerned, the least stable modes (i.e., the modes with the largest exponential growth rate) are far more important. They can be computed by coupling an iterative eigenvalue algorithm to a time-stepping routine for the linearized equations of motion (abbreviated by the linear operator L) over a given time interval Δt [10]. This way, the iterative eigenvalue solver will efficiently extract the modes that are most amplified over a time interval Δt, that is to say, it extracts the least stable modes of L. The propagation time Δt, in general, affects the speed of convergence of the iterative solver, but does not influence which modes the solver will converge to.

This technique is generally sufficient for bounded flows that are governed by a limited number of (or even a single) dominant instability mechanism, since the associated spectrum consists of eigenvalues that are well separated. In this case, the principal eigenvalues are easy to isolate by the iterative algorithm. When multiple and competing mechanisms are at play, the spectrum is far more complicated, and physically relevant modes are more difficult to extract. In particular, eigenvalue clusters near the neutral axis, stemming from continuous branches or even numerical artifacts, pose a great challenge to the convergence of the iterative algorithm. Even though unstable modes may still be extracted, the stable part of the spectrum quickly moves beyond the reach of the iterative algorithm; in this case, a different strategy is called for.

The region of convergence may be manipulated and adjusted by a rational transformation of the complex eigenvalue plane. The “shift–invert” method [16] allows the computation of the modes whose eigenvalues are closest to a complex shift parameter σ. But at each iteration of the eigenvalue solver, the method requires the solution of a linear system of the form (L  σI)x = b. Most studies accomplish the latter solution using a direct LU decomposition which has to be performed once at the start and is used for all successive iterations, until the shift parameter σ is changed to access different parts of the spectrum. A variant of the shift–invert method, known as the Cayley transformation, yields better convergence, if an iterative solution of the linear system is chosen [18].

The LU decomposition is based on a matrix representation of the linear stability operator. Some global stability investigations used spectral discretization methods which resulted in a dense matrix of moderate size [2]. Later studies took advantage of a sparse representation, in particular, when the operator arises from a finite-element or finite-difference discretization [5]. In this case, the number of non-zero elements is proportional to the number of degrees of freedom N, making a sparse matrix representation convenient to handle computationally. Highly efficient multi-frontal LU solvers for large-scale sparse matrices are readily available (see e.g. [9], [3], [4]), but the sparsity of the output matrices is not always guaranteed. Even though the bandwidth of the factorization is not greater than that of the original matrix [11, p. 152], all elements between the upper and lower band may be non-zero. Computing and storing the decomposition thus results in substantial memory requirements. For example, for the discretization of a two-dimensional problem with N degrees of freedom on a structured mesh, the bandwidth scales with N1/2, in which case the LU decomposition would contain up to O(N3/2) non-zero elements. In three dimensions, as the bandwidth increases to N2/3, the memory requirements go up as N5/3. For a discretization of the compressible Navier–Stokes equations on a two-dimensional domain with 256 × 512 points using a finite-differences scheme with a six-points stencil, one can estimate that storing the LU decomposition requires about 80 GB of memory (work space requirements during the decomposition tend to be even larger), thus illustrating the limitation of this method for larger-scale problems. In some cases, appropriate reordering of the matrix entries can somewhat alleviate the problem by improving the sparsity of the factorization; this approach, however, does not provide a viable and extendable solution for large-scale problems.

An alternative that avoids the computation of the LU decomposition of the operator L consists of iterative algorithms [24] to solve the linear system arising from the shift–invert or Cayley transformation. Together with ILU-type preconditioners, this approach has been applied to incompressible [15] and compressible flows [18]; in [22], the authors use un-preconditioned iterative solvers for the computation of unstable modes in plasma flows. These methods yield a reduction of computational costs associated with the solution of the linear system, but they do not provide the same level of versatility as direct methods do. Indeed, if one chooses the shift parameter σ close to the spectrum of L, then (L  σId) becomes ill-conditioned. In this case, the cost of preconditioning as well as the number of iterations for the linear solver to converge have to be assessed critically. In contrast, if σ is selected farther from the spectrum, the iterative linear solver is more likely to converge with a “cheap” preconditioner, but at the same time the focusing effect of the shift–invert transformation is rather weak; consequently, it may not be possible to extract the desired modes.

The present paper presents a method for selectively extracting modal information from a linear operator L without relying on the (iterative or direct) solution of a linear system. This approach has been inspired both by the “shift–invert” technique for the solution of eigenvalue problems [16] and by the selective frequency damping method of Akervik et al. [1] for the computation of unstable steady flow. Similar to the latter method, we propose to use a relaxation procedure to selectively stabilize parts of the spectrum away from a chosen frequency shift, after which a standard Krylov method is employed to obtain the least stable modes of the relaxed system. Although the spectral transformation involved in the present “shift-relax” technique is somewhat less flexible than the “shift–invert” technique, its low memory requirement and ease of implementation make it suitable and attractive for large-scale stability computations of two- or three-dimensional flows.

Section snippets

Definition of the problem

Let the dynamics of the problem under consideration be governed by a set of non-linear equations of the formq̇=F(q),where q is the state vector and F denotes a discrete integro-differential operator with appropriate boundary conditions. For simplicity, only finite-dimensional operators (which arise after spatial discretization) will be considered in this paper. We assume that this operator has a fixed point q0, such that F(q0) = 0. If this base state is stable to finite-amplitude perturbations,

Governing equations for local and global computations

We consider a compressible jet of radius R, with characteristic velocity U0, density ρ0 and temperature T0 measured on the centerline, discharging into a fluid at rest with density ρ and temperature T. These same quantities are used to make the problem non-dimensional. In a cylindrical coordinate system (x, r, θ), the nonlinear governing equations are expressed in terms of the conservative flow variables q = (ρ, ρux, ρur, ρuθ, ρE)T, where ρ is the density, u = uxex + urer + uθeθ is the flow velocity, and E

Non-parallel base flow

The SR method is now applied to compute global modes (18) of a non-parallel jet. The geometry of the computational domain is represented in Fig. 6. The jet exits from an idealized nozzle, modeled as an infinitely thin adiabatic wall at r = 1 and x  0. Only the upper half-plane r > 0 is resolved in the calculations, and appropriate symmetry conditions for axisymmetric flow are imposed on the jet axis r = 0. In order to control the jet profile at x = 0 (diffusive effects inside the pipe can be

Conclusion

A new numerical procedure for the solution of large eigenvalue problems has been presented. A relaxation technique using a first-order temporal bandpass filter is coupled to to the linearized equations of motion, such that the least stable eigenmodes of the filtered system lie in a prescribed frequency band of interest centered around a shift frequency. These modes are then recovered through propagation over a finite time interval, using standard eigenvalue extraction techniques. This

Acknowledgments

This work was supported by DGA Grant No. 2009.60.034.00.470.75.01. The authors are thankful to Miguel Fosas and Patrick Huerre for their comments.

References (26)

  • S. Balay, K. Buschelman, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, L. Curfman McInnes, B.F. Smith, H. Zhang,...
  • A. Barbagallo et al.

    Closed-loop control of an open cavity flow using reduced-order models

    J. Fluid Mech.

    (2009)
  • C. Bogey et al.

    Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases

    Acta Acoust.

    (2002)
  • Cited by (16)

    • The verification of multiplicity support of a defective eigenvalue of a real matrix

      2022, Applied Mathematics and Computation
      Citation Excerpt :

      The computation of eigenvalues of matrices is one of the fundamental problems in the area of numerical computation with a wide range of applications, including vibrational analysis of structures [33], stability analysis of fluid flows [8,41], optimal control problems [6], quantum mechanics [35], and delay differential equations [22,23].

    • Stability analysis and improvement of the solution reconstruction for cell-centered finite volume methods on unstructured meshes

      2019, Journal of Computational Physics
      Citation Excerpt :

      Hence, the direct calculation of the eigenvalues and their corresponding eigenvectors (e.g., the LQ and the QR factorization algorithms [35] which are embedded in off-the-shelf linear algebra packages such as LAPACK [36]) are time-consuming and in some cases fail due to the computational limitations. Instead, we use iterative methods such as Arnoldi's and Krylov's subspaces (see [37–39]) that are much more efficient in calculating the parts of the spectrum that are of interest (the rightmost eigenvalues in our case). Even so, better and faster convergence is achieved for well-separated eigenvalues.

    • Computation of frequency responses for linear time-invariant PDEs on a compact interval

      2013, Journal of Computational Physics
      Citation Excerpt :

      In particular, wall-bounded shear flows of both Newtonian and viscoelastic fluids have non-normal dynamical generators of high spatial order and the ability to accurately compute frequency responses for these systems is of paramount importance; additional examples of systems with non-normal generators, for which the tools developed in this paper are particularly well-suited, can be found in the outstanding book by Trefethen and Embree [8] and the references therein. The utility of non-modal analysis in understanding the dynamics of infinitesimal fluctuations around laminar flow conditions has been well-documented; see [1,9–15] for Newtonian fluids and [16–20] for viscoelastic fluids. In viscoelastic fluids with large polymer relaxation times, analysis is additionally complicated by the fact that pseudo-spectral methods exhibit spurious numerical instabilities [21,22].

    View all citing articles on Scopus
    View full text