A relaxation method for large eigenvalue problems, with an application to flow stability analysis
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 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 formwhere 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)
- et al.
Global two-dimensional stability measures of the flat plate boundary-layer flow
Eur. J. Mech. B-Fluid
(2008) - et al.
Multifrontal parallel distributed symmetric and unsymmetric solvers
Comput. Method Appl. Mech. Eng.
(2000) - et al.
High-order low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems
J. Comput. Phys.
(2007) - et al.
Krylov methods for the incompressible Navier–Stokes equations
J. Comput. Phys.
(1994) - et al.
Jacobian-free Newton–Krylov methods: a survey of approaches and applications
J. Comput. Phys.
(2004) Compact finite difference schemes with spectral-like resolution
J. Comput. Phys.
(1992)- et al.
A preconditioned Krylov technique for global hydrodynamic stability analysis of large-scale compressible flows
J. Comput. Phys.
(2010) - et al.
Boundary conditions for direct simulations of compressible viscous flows
J. Comput. Phys.
(1992) - et al.
Fast eigenvalue calculations in a massively parallel plasma turbulence code
Parallel Comput.
(2010) - et al.
Steady solutions of the Navier–Stokes equations by selective frequency damping
Phys. Fluids
(2006)
Closed-loop control of an open cavity flow using reduced-order models
J. Fluid Mech.
Three-dimensional non-reflective boundary conditions for acoustic simulations: far field formulation and validation test cases
Acta Acoust.
Cited by (16)
The verification of multiplicity support of a defective eigenvalue of a real matrix
2022, Applied Mathematics and ComputationCitation 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 PhysicsCitation 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.
Effectivity and efficiency of selective frequency damping for the computation of unstable steady-state solutions
2018, Journal of Computational PhysicsComputation of frequency responses for linear time-invariant PDEs on a compact interval
2013, Journal of Computational PhysicsCitation 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].
Transition and turbulence in horizontal convection: Linear stability analysis
2017, Journal of Fluid MechanicsA novel optimization algorithm for the selective frequency damping parameters
2022, Physics of Fluids