Multiparameter spectral analysis for aeroelastic instability problems
MMultiparameter spectral analysis for aeroelastic instability problems
Arion Pons and Stefanie Gutschmidt Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, UK. Department of Mechanical Engineering, University of Canterbury, Christchurch
Key words multiparameter eigenvalue problems; instability; aeroelasticity.
AMS classification
This paper presents a novel application of multiparameter spectral theory to the study of structural stability, with particular emphasis on aeroelastic flutter. Methods of multiparameter analysis allow the development of new solution algorithms for aeroelastic flutter problems; most significantly, a direct solver for polynomial problems of arbitrary order and size, something which has not before been achieved. Two major variants of this direct solver are presented, and their computational characteristics are compared. Both are effective for smaller problems arising in reduced-order modelling and preliminary design optimization. Extensions and improvements to this new conceptual framework and solution method are then discussed. Introduction
Predicting and controlling aeroelastic instability forms a major part of the discipline of aeroelasticity. Many different physical systems show flutter instability, and many models exist to describe them. In a linear system, or the linearization of a nonlinear system, the onset of flutter or divergence can be formulated in the well-known stability criterion:
Im(๐) > 0 for stability (1) where ๐ are the time-eigenvalues of the system according to the Fourier transform ๐(๐ก) =๐ฬ๐ ๐๐๐ก for the system coordinate ๐ [1]. These eigenvalues may be nondimensionalised. Flutter occurs when the system parameters are such that the system is on the stability boundary, Im(๐ f ) = 0 . The flutter point may then be described as a tuple which includes the set of system parameters (particularly, an airspeed parameter) and the modal frequency of instability, ๐ f . Typically only one or two flutter or divergence points are of industrial relevance. Eq. 1 is not however the only criterion that can be used to characterize instability. This leads us into the study of aeroelastic methods. Apart from the four established methods โ the p-method , classical flutter analysis , the k-method (or U-g method , or
V-g method ) and the p-k method ; all of which are detailed and discussed in a number of reference works [1โ3] โ recent years have seen a proliferation of new aeroelastic methods. Several authors have refined the existing established methods for particular scenarios or applications [4โ6]. The application of concepts from robust control theory have yielded a series of methods, including the ๐ -method by Lind and Brenner [7,8], the ๐ -k method by Borglund [9โ12], and others [13โ15]. The prime advantage of these ๐ -type methods is that they facilitate the propagation of uncertainty * Corresponding author, [email protected], tel. +44 755 366 3296 ultiparameter analysis for instability problems Page of distributions through the system, allowing a worst-case flutter speed estimate to be made in a system with high uncertainty. Other developments have come from other fields: Afolabi [16,17] characterized coupled-mode flutter as a loss of eigenvector orthogonality, using methods from catastrophe theory. Irani and Sazesh [18] used stochastic methods, while Gu et al. [19] devised a genetic algorithm, and a number of authors [20โ23] have applied neural networks to the detection of flutter points. It is in the context of these developments that we propose our method of analyzing flutter problems. The central methodological contribution of this paper is the concept that the solution of an aeroelastic system for its flutter points is nothing other than a multiparameter eigenvalue problem. We will show the simple link between the aeroelastic stability problem and multiparameter spectral theory, and how this allows for direct solution of a variety of flutter problems. The purpose of this paper is to detail this link and to explore the mechanisms by which this direct solution may be accomplished. For this purpose we will apply our analysis to simple example problems โ however, the method does extend to problems that are significantly more complex; these are considered in Pons and Gutschmidt [24, 25]. Aeroelastic flutter as a multiparameter eigenvalue problem (MEP)
Consider a linear finite-dimensional system with eigenvector ๐ฑ โ โ ๐ and arbitrary continuous dependence on both an eigenvalue parameter ๐ โ โ , and another structural or environmental parameter ๐ โ โ : A(๐, ๐)๐ฑ = ๐ (2) where
A โ โ ๐ร๐ . Any complex structural parameter can of course be split into two real parameters. The stability problem for this system (with respect to parameter ๐ ) is to find ๐ such that an eigenvalue of the problem ( ๐ ) has zero imaginary part. This point is the โstability boundaryโ: for a system with multiple structural parameters, the stability boundary may be a line or other higher-dimensional surface. We then note that the condition Im(๐) = 0 is equivalent to modifying the original definition of the problem such that ๐ โ โ and not ๐ โ โ . Such a manoeuvre does not seem to be immediately useful: under ๐ โ โ , a solution to Eq. 2 only exists on the stability boundary, and nowhere else. In order to develop, for example, iterative methods for flutter point calculation, we need to be able to define some form of solution in the subcritical and supercritical areas (above and below the stability boundary, respectively). There is an easy way of doing this. Following [26], we take the complex conjugate of Eq. 2 as another equation:
A(๐, ๐)๐ฑ = ๐,
Aฬ (๐, ๐)๐ฑฬ = ๐. (3) As ๐ โ โ and ๐ โ โ are unaffected by the conjugation, this operation enforces these conditions. This procedure has been utilized before in the analysis of delay differential equations [26], and (in a limited form) in the context of Hopf bifurcation prediction [27], but has never been applied to aeroelastic or other structural stability problems. Equation 3 is nothing other than a multiparameter eigenvalue problem (MEP): an eigenvalue problem in which the eigenvalue point is not simply defined by a scalar and an eigenvector, but by an ๐ -tuple and an eigenvector. A number of methods of analysis have been developed for such problems, and in this paper we will explore some of these. However, as the methods that are available depend strongly on the structure of matrix A , we will first define a system to work with. ultiparameter analysis for instability problems Page of An Example Section model 3.1. Formulation
Consider first the simple section model shown in Figure 1. This model has two degrees of freedom: plunge โ and twist ๐ . The governing equations for this model are easy to derive; they are: ๐โฬ + ๐ โ โฬ + ๐ โ โ โ ๐๐ฅ ๐ ๐ฬ = โ๐ฟ(๐ก), ๐ผ ๐ ๐ฬ + ๐ ๐ ๐ฬ + ๐ ๐ ๐ โ ๐๐ฅ ๐ โฬ = ๐(๐ก), (4) where ๐ and ๐ผ ๐ are the section mass and polar moment of inertia, ๐ โ and ๐ ๐ are the section plunge and twist stiffnesses, ๐ โ and ๐ ๐ are the section plunge and twist damping coefficients, and ๐ฅ ๐ is the sectionโs static imbalance โ defined as the distance along the ๐ฅ -axis from the pivot point to the centre of mass. Taking the Fourier transform, [โ(๐ก), ๐(๐ก)] = [โฬ, ๐ฬ]๐ ๐๐๐ก , of this model, we obtain: (โ๐๐ + ๐๐ โ ๐ + ๐ โ )โฬ + ๐๐ฅ ๐ ๐ ๐ฬ = ๐ฟ(๐, โฬ, ๐ฬ), ๐๐ฅ ๐ ๐ โฬ + (โ๐ผ ๐ ๐ + ๐๐ ๐ ๐ + ๐ ๐ )๐ฬ = ๐(๐, โฬ, ๐ฬ). (5) To model the aerodynamic loads in the frequency domain we use Theodorsenโs unsteady aerodynamic theory [3]: ๐ฟ = โ๐ (๐ฟ โ โฬ + ๐ฟ ๐ ๐ฬ) , ๐ = ๐ (๐ โ โฬ + ๐ ๐ ๐ฬ) . (6) The aerodynamic coefficients {๐ฟ โ , ๐ฟ ๐ , ๐ โ , ๐ ๐ } are complex functions of ๐ โ the reduced frequency; an aerodynamic parameter related to the airspeed ( ๐ ) by ๐ = ๐๐ ๐โ . Other structural and environmental parameters involved are the air density ( ๐ ), the semichord ( ๐ ) and the distance along the ๐ฅ -axis from the midchord to the pivot point, as a fraction of the semichord ( ๐ ). See Pons [24] or Hodges and Pierce [3] for details. We assume that the flow over the airfoil is quasisteady: that is, that Theodorsenโs function takes a value of universally [1,3]. We will deal with general Theodorsen aerodynamics in a later paper, as this requires iterative or approximate multiparameter solution methods which we will not cover here. We also assume without loss of generality a lift-angle of attack coefficient of ๐ถ ๐ฟ๐ผ = 2๐ , as per thin-airfoil theory [28]. It is then customary to nondimensionalise Eq. 5. Further details of this are given in Pons [24]. The final result is a flutter problem of the form ((M + G + G
1๐ + G ) ๐ โ D ๐ โ K ) ๐ฑ = ๐. (7) with dimensionless parameters defined as in Table 1 and the nomenclature, and the matrix coefficients G = 1๐ [1 ๐๐ ( + ๐ )], (8) G = 1๐ [ โ2๐ 2๐(1 โ ๐)โ๐(1 + 2๐) ๐๐(1 โ 2๐)], ultiparameter analysis for instability problems Page of Table 1 : Dimensionless parameter values for the section model
Parameter Value mass ratio โ ๐ radius of gyration โ ๐ bending nat. freq. โ ๐ โ torsional nat. freq. โ ๐ ๐ bending damping โ ๐ โ torsional damping โ ๐ ๐ static imbalance โ ๐ ๐ โ0.1 pivot point location โ ๐ โ0.2 Figure 1:
Diagram of section model G = 1๐ [0 20 1 + 2๐] , M = [ 1 โ๐ ๐ โ๐ ๐ ๐ ], (8) D = [2๐๐ โ ๐ โ
00 2๐๐ ๐ ๐ ๐ ๐ ] , K = [๐ โ2
00 ๐ ๐ ๐2 ], Note that this nondimensionalisation is a convenience and not necessary prerequisite for using multiparameter solution methods. However something that will be of great use is to rearrange Eq. 7 into two polynomial forms by defining new eigenvalue parameters:
ฮฅ = ๐ ๐โ , ๐ = 1 ๐ โ , and ๐ = 1 ๐โ . These two forms are then ((M + G )๐ + G ฮฅ๐ + G ฮฅ โ D ๐ โ K )๐ฑ = ๐ , (9) which we term the ฮฅ - ๐ form, and ((M + G ) + G ๐ + G ๐ โ D ๐ โ K ๐ )๐ฑ = ๐ , (10) which we term the ๐ - ๐ form. Both are preferable to the ๐ - ๐ form, Eq. 7. We should note at this point that, although we have here derived these two forms for a small 2-DOF section model, these forms arise in many other models; often with significantly larger matrices. The solution methods we develop for this form will thus be applicable to a wide variety of problems. In many aeroelastic systems structural damping is negligible, in which case D = 0 . While this may seem to be a trivial case of the preceeding systems, the omission of the structural damping term does allow us to change the structure of the system, leading to a faster computation time (as we will show later). In the case of Eq. 10, we have ((M + G ) + G ๐ + G ๐ โ K ฮ)๐ฑ = ๐, (11) where
ฮ = ๐ . This system is now linear in ฮ , whereas it was previously quadratic in ๐ . Note that the ๐ - ๐ form is the only polynomial form presented which allows us to make one parameter linear in this situation. ultiparameter analysis for instability problems Page of The three polynomial systems that we have introduced in this section (Eq. 9, 10 and 11) are not yet fully constrained MEPs, as they only consist of one equation. To constrain the system we employ the method introduced in Section 2 โ the addition of an equation representing the conjugate of the initial equation. We obtain: ((M + G ) + G ๐ + G ๐ โ K ฮ)๐ฑ = ๐, ((Mฬ + Gฬ ) + Gฬ ๐ + Gฬ ๐ โ Kฬ ฮ)๐ฑฬ = ๐, (12) ((M + G ) + G ๐ + G ๐ โ D ๐ โ K ๐ )๐ฑ = ๐, ((Mฬ + Gฬ ) + Gฬ ๐ + Gฬ ๐ โ Dฬ ๐ โ Kฬ ๐ )๐ฑฬ = ๐, (13) and ((M + G )๐ + G ฮฅ๐ + G ฮฅ โ D ๐ โ K )๐ฑ = ๐, ((Mฬ + Gฬ )๐ + Gฬ ฮฅ๐ + Gฬ ฮฅ โ Dฬ ๐ โ Kฬ )๐ฑฬ = ๐. (14) These systems are all polynomial multiparameter eigenvalue problems [29]. A number of solution methods are known for such systems. Direct solution via linearization 4.1. Linearization
Any polynomial MEP can be made into a linear one (consisting of a zeroth-order term and first-order terms in each eigenvalue) via the process of linearization [30]. This process bears resemblance to the linearization of single-parameter polynomial eigenvalue problems; a process which is well-known [31]. Multiparameter linearization is relevant because a number of direct solution methods exist for linear MEPs. Consider first Eq. 12. This equation contains only one quadratic variable ( ๐ ), as ฮ is already linear. To linearize this system, we define a new eigenvector which contains a factor of ๐ : ๐ช = [๐ฑ; ๐๐ฑ] . By expanding the systemโs coefficient matrices and using this factor of ๐ in the eigenvector to reduce the order of the quadratic term, we obtain a linear problem of double the size: ([M + G
00 โ๐ผ ๐ ] + [โK
00 0] ฮ + [G G ๐ผ ๐ ([Mฬ + Gฬ
00 โ๐ผ ๐ ] + [โKฬ
00 0] ฮ + [Gฬ Gฬ ๐ผ ๐ (15) The upper rows of Eq. 15 represents Eq. 12 directly, and the lower row represents the identity ๐ผ ๐ (๐๐ฑ) = ๐๐ผ ๐ ๐ฑ . Note that other linearizations of similar form are possible. However, irrespective of the exact linearization used, the resulting system will be of the form: (A + Bฮ + C๐)๐ช = ๐, (Aฬ + Bฬ ฮ + Cฬ ๐)๐ชฬ = ๐. (16) Note that it is not actually necessary to linearize the second (i.e. conjugate) equation of the aeroelastic system, because the linearized conjugate equation will be the conjugate of the linearized first equation. This is a property of this method of linearization. This same linearization process can be applied to any quadratic MEP. Any problem of the form (A + B๐ + C๐ + D๐๐ + E๐ + F๐ )๐ฑ = ๐, (17) ultiparameter analysis for instability problems Page of can be linearized as ([A B C0 โ๐ผ ๐
00 0 โ๐ผ ๐ ] + [ 0 D E๐ผ ๐ + [ 0 0 F0 0 0๐ผ ๐ (18) We have A = M + G B = G C = โD D = โK E = 0
F = G (19) for Eq. 13, and changing ๐ โ ๐ and ๐ โ ฮฅ , A = โK B = โD C = 0
D = M + G E = G F = G (20) for Eq. 14. These linearised systems have matrix coefficients triple the size of the original problem. Consider Eq. 16. Post-multiplying the first equation in this system by
Cฬ ๐ฒ and premultiplying the second by C๐ฑ , we have (A + Bฮ + C๐)๐ฑ โ (Cฬ ๐ฒ) = 0, (C๐ฑ) โ (Aฬ + Bฬ ฮ + Cฬ ๐)๐ฒ = 0. (21) These two equations are both equal to zero so we may equate them. After cancelling the terms in ๐ , the expression can be manipulated into: ฮ ๐ณ = ฮฮ ๐ณ (22) and an enlarged eigenvector ๐ณ = ๐ฑ โ ๐ฒ and the operator determinants ฮ = B โ Cฬ โ C โ Bฬ ฮ = C โ Aฬ โ A โ Cฬ ฮ = A โ Bฬ โ B โ Aฬ . (23) Eq. 22 is a generalized eigenvalue problem (GEP), in the single parameter ๐ . Solvers for the generalized eigenvalue problem are very widely available. If the linear system has square coefficients of size ๐ then the operator determinants are of size ๐ . The operator determinants can also be used to define a second GEP in ๐ . By multiplying the first and second equations of Eq. 16 by Bฬ ๐ฒ and B๐ฑ respectively, we can also show that: ฮ ๐ณ = ๐ฮ ๐ณ. (24) However, it is only necessary to solve one of Eq. 22 or Eq. 24: once one has been solved, then its solutions can be substituted back into the original polynomial system, which yields another smaller GEP. Alternatively, if the eigenvalue is simple then it may be computed more cheaper via Rayleigh quotients in ๐ณ or (decomposing ๐ณ = ๐ฑ โ ๐ฒ ), ๐ฑ and ๐ฒ . The system is thus completely solved for its flutter points. This direct solution method is known in mathematical literature as the operator determinant method . Its computational complexity is ๐ช(๐ ), irrespective of whether ๐ is the size of the linear coefficients ( A , B , etc.) or the polynomial coefficients ( M , G , etc.) [30,32,33]. This large complexity arises from solving the generalized eigenvalue problem, an ๐ช(๐ ) process by the QZ algorithm [34], with operator ultiparameter analysis for instability problems Page of determinants of size ๐ = ๐ช(๐ ) . The operator determinant method has not previously been used in aeroelasticity, and has only rarely seen engineering application in the study of dynamic model updating [35,36]. One important caveat of the operator determinant approach is that the matrix ฮ must not be singular. If it is, then the eigenvalues of the original polynomial system (e.g. Eq. 12-14) will not generally coincide with those of the GEPs of the linearized problem (Eq. 22 and 24) [30,37,38]. A linear MEP with singular ฮ is said to be singular MEP. A large proportion of the linear flutter problems that arise in the study of aircraft aeroelasticity are singular. This is largely because the linearization of polynomial problems tends to generate singular linear problems, even if the original polynomial problem has all its matrix coefficients at full rank โ compare Eq. 18. Some smaller linearizations are not necessarily singular: for example Eq. 15. However, here we find that (for our aerodynamic model) the coefficient matrix G is not full rank and so the linearized problem is singular anyway. In many circumstances we will be dealing with a singular problem, and for a long time the lack of a working solver for such singular problems has been a major obstacle to the application of multiparameter methods to real-world problems. However, recently a solution has been proposed. Muhiฤ and Plestenjak [29] proved that the eigenvalues of a polynomial system are equivalent to the finite regular eigenvalues of the pair of singular operator determinant GEPs constructed via linearization. The finite regular eigenvalues of a linear two-parameter system are the pairs ( ๐ , ๐ ) such that ๐ โ {1,2} [37]: rank(A ๐ + B ๐ ๐ + C ๐ ๐) < max (๐ ,๐ก)โโ rank(A ๐ + B ๐ ๐ + C ๐ ๐ก), (25) that is, the finite regular eigenvalues are the set of points that cause the singular problem to have its maximum rank, even though this is not full rank. On the basis of this proof, Muhiฤ and Plestenjak [29] devised a set of algorithms which would extract the common regular part of the singular matrix pencils (i.e. matrix-valued functions polynomial in a variable ๐ ) ฮ โ ๐ฮ and ฮ โ ๐ฮ . This common regular part is represented by two smaller nonsingular matrix pencils ( ฮ โ ๐ฮ and ฮ โ ๐ฮ ), the eigenvalues of which are the finite regular eigenvalues of the singular problem and thus the eigenvalues of the polynomial problem. In practical terms, this can be seen as a compression of the singular operator determinant matrices into smaller full-rank matrices. The operator determinant method, as presented in Section 4.2, can be applied to these compressed pencils. The algorithms involved in the extraction of the common regular part are presented in [29] and published also in MATLAB code [39]. We will not detail them here as they are complex. Direct solution via Quasi-linearization 5.1. Quasi-linearization
Hochstenbach et al. [30] recently presented another method of linearization for polynomial MEPs. Instead of increasing the size of the coefficient matrices, this method of linearization increases the number of parameters and equations in the system. To differentiate it from strict linearization, Hochstenbach et al. [30] term this new method quasi-linearization . The process is as follows. Considering again Eq. 12, we define a new eigenvalue parameter ๐ผ = ๐ . The equations then become ultiparameter analysis for instability problems Page of ((M + G ) + G ๐ + G ๐ผ โ K ฮ)๐ฑ = ๐, ((Mฬ + Gฬ ) + Gฬ ๐ + Gฬ ๐ผ โ Kฬ ฮ)๐ฑฬ = ๐. (26) This is now a linear three-parameter eigenvalue problem, but with only two equations. We need a third equation to constrain the system. We have the relation ๐ผ โ ๐ = 0 , but this is nonlinear. However, noticing that we can write this relation as det ([๐ผ ๐๐ 1]) = 0 (27) we can recast it as a new MEP to add to Eq. 26 (using an arbitrary eigenvector ๐ฒ ): ((M + G ) + G ๐ + G ๐ผ โ K ฮ)๐ฑ = ๐, ((Mฬ + Gฬ ) + Gฬ ๐ + Gฬ ๐ผ โ Kฬ ฮ)๐ฑฬ = ๐, ([0 00 1] + [0 11 0] ๐ + [1 00 0] ๐ผ) ๐ฒ = ๐. (28) In a similar way, we can linearize Eq. 13 with the definitions ๐ผ = ๐ and ๐ฝ = ๐ : ((M + G ) + G ๐ + G ๐ผ โ D ๐ โ K ๐ฝ)๐ฑ = ๐, ((Mฬ + Gฬ ) + Gฬ ๐ + Gฬ ๐ผ โ Dฬ ๐ โ Kฬ ๐ฝ)๐ฑฬ = ๐, ([0 00 1] + [0 11 0] ๐ + [1 00 0] ๐ผ) ๐ฒ = ๐, ([0 00 1] + [0 11 0] ๐ + [1 00 0] ๐ฝ) ๐ฒ = ๐, (29) and Eq. 14 with the definitions ๐ผ = ๐ , ๐ฝ = ฮฅ , ๐พ = ฮฅ๐ : ((M + G )๐ผ + G ๐ฝ + G ๐พ โ D ๐ โ K )๐ฑ = ๐, ((Mฬ + Gฬ )๐ผ + Gฬ ๐ฝ + Gฬ ๐พ โ Dฬ ๐ โ Kฬ )๐ฑ = ๐, ([1 00 0] ๐ผ + [0 11 0] ๐ + [0 00 1]) ๐ฒ = ๐, ([0 10 0] ๐ผ + [0 01 0] ๐ฝ + [1 00 1] ๐พ) ๐ฒ = ๐, (30) These systems can be solved via the operator determinant method, in a more general form than that presented in Section 4.2. The general form of a nonhomogeneous MEP is ๐ ๐ (๐ผ) = A ๐0 ๐ฑ ๐ + โ ๐ ๐ A ๐๐ ๐ฑ ๐๐๐=1 , ๐ = 1, โฆ , ๐ (31) where ๐ผ is the vector of eigenvalues ( ๐ ๐ being the individual eigenvalues), A ๐๐ are the coefficient matrices, which can be complex and of different sizes for each equation, and ๐ฑ ๐ are the eigenvectors. Eq. 31 can be visualized as a non-square array of matrices: [A A A โฏ A A A A โฏ A โฎ โฎ โฎ โฑ โฎA ๐0 A ๐1 A ๐2 โฏ A ๐๐ ]. (32) ultiparameter analysis for instability problems Page of In a process analogous to that presented in Section 4.4, we can construct a series of operator determinants for this system. There are
๐ + 1 such operator determinants ( ฮ through to ฮ ๐ ), of size ๐ ๐ for constant coefficient size ๐ . They correspond to taking determinants of Eq. 32, with certain columns removed or inserted, and with the normal scalar multiplication operation replaced by a Kronecker product between matrices. Definitions and derivations of these determinants may be found in [40โ44], and one software implementation in [39]. In the two-parameter case this analysis collapses into that of Section 4.2. These general operator determinants allow us to compute the eigenvalues of the Eq. 31 via generalized eigenvalue problems. Providing ฮ is nonsingular, it holds that ฮ ๐ ๐ฑ = ๐ ๐ ฮ ๐ฑ. (33) In this way we are able to solve the quasi-linearized systems given in Section 5.1 (we will discuss singularity in Section 5.4). However, note that already by using this quasi-linearization process we have made the operator determinants smaller than they were with standard linearization. Using Eq. 12, the operator determinants are square and of size , and using Eq. 13 or 14 they are of size . With standard linearization they are of size and , respectively. In the two-parameter case, the computation of the operator determinants is trivial. However in the case of larger systems we must find some other approach. In the general case, the operator determinants may be computed by modifying the Leibniz formula for the determinant [40]:
ฮ(M) = โ sgn(๐ฌ) ๐ฌ โ ๐ ๐ โ ๐๐=1 M ๐,๐ฌ ๐ , (34) where ๐ ๐ is the set of permutations of the set {1, โฆ , ๐} , sgn(๐ฌ) is the sign of the permutation vector ๐ฌ โ ๐ ๐ . M is the square matrix array (a modification of Eq. 32) corresponding to the operator determinant desired. The notation โ ๐๐ denotes the repeated application of the Kronecker product. Note that the tensor determinant definition in [40] is slightly erroneous as the factor (โ1) sgn(๐) in their tensor determinant expressions should be either sgn(๐) or (โ1) ๐(๐) , where
๐(๐) is the number of inversions in ๐ . Alternatively, Muhiฤ and Plestenjak [39] devised an operator determinant Laplace expansion [45] that is able to compute the operator determinants of a system of arbitrary size, via a process of recursion: ฮ(M) = โ(โ1) ๐+1 M ๐1 โ ๐ ๐1๐๐=1 , (35) where ๐ ๐1 denotes the minor entry corresponding to M ๐1 . The summation in Eq. 35 follows the first column of M , though, of course, many other summation paths could be used. It should be noted that the use of the Laplace or Leibniz methods for computing the determinant of an ordinary matrix have very high non-polynomial computational complexity costs โ ๐ช(๐!) and
๐ช(๐! ๐) , respectively [46,47].
One major advantage of the direct solver based on quasi-linearization is that it generates linearized problems that are not always singular โ when the coefficient matrices are of full rank, Eq. 28-30 are in general nonsingular. However in our case the coefficient G is singular, and so it happens that all these three equations are singular anyway. The theory of the ultiparameter analysis for instability problems Page of compression process noted in Section 4.3 has never been extended to the general multiparameter case, and so we have no rigorous method of solving singular MEPs with ๐ > 2 . However, we do have a practical method of doing so. Numerical experiments (some of which are detailed in Section 6) indicate that, although the compression algorithm takes only ฮ , ฮ and ฮ as an input, it will successfully compress these three operator determinants for the general multiparameter problem (ignoring all the others). In the case of Eq. 28-30, with these three operator determinants we can solve for all the eigenvalue parameters of the system, irrespective of the order in which we arrange the eigenvalues. We should note that there is no justification for this procedure other than the experimental evidence that the algorithm works โ evidence confirmed also by the nonrigourous ๐ > 2 compression process utilized in [39]. The common regular part relationship proved by [29] applies only to two-parameter eigenvalue problems arising from linearization, and indeed it is not clear how the concept of the common regular part would extend to an
๐ > 2 problem with more than two pencils. That the two-parameter algorithm does work for the first three operator determinants of an
๐ > 2 problem implies that the algorithm will generalize โ it is a question of working out what this generalization is. This is an interesting area for further research. For the purposes of this paper, however, we will use the compression algorithm for the quasi-linearized ฮ , ฮ and ฮ without proof. Numerical experiments 6.1. Undamped model
We are now in a position to compute the flutter points of the models we introduced in Section 3. We have devised two direct solvers to do so: a direct solver with linearization and a direct solver with quasi-linearization. Consider first the undamped section model (Eq. 12), with parameter values as per Table 1. To validate our direct solution methods, we first produce a modal damping plot (Figure 2). Three points of neutral stability may be observed ๐ = ยฑ1.00 and ๐ = 0 . The point ๐ ๐น = 1.00 is the only physical flutter point. We can link the modal damping curve for this flutter point with the lower modal frequency path in the Re(๐) plot (this requires data not observed directly on Figure 2). We thus can estimate that ๐ ๐น =1.32 rad/s ( ฮ ๐น = 0.57 s /rad ). The point at ๐ = โ1.00 along the same mode corresponds to a nonphysical flutter event occurring at negative airspeed. Finally, at ๐ = 0 both modes have neutral stability (at ๐ = 0.548 rad/s and ๐ = 1.426 rad/s or ฮ = 3.33 s /rad and ฮ =0.492 s /rad ). This, however, merely represents the fact that the structure is undamped at zero airspeed due to the lack of any structural damping in the model. The divergence point of this system does not appear on this plot, as it occurs at infinite ๐ and ฮ . However, the computation of divergence points does not require multiparameter methods anyway, as if the frequency ๐ is assumed to be zero in the governing equation (e.g. Eq. 9) then the problem becomes a single parameter eigenvalue problem in the airspeed parameter (e.g. ฮฅ ). This problem can then be solved with single-parameter solvers. We then compute the flutter points of the system via our two direct solution methods, yielding flutter points identical to those detailed above for the presented accuracy. With the solver based on linearization, the uncompressed operator determinants are of size and the compressed ones of size . With the solver based on quasi-linearization the size of the uncompressed determinants can be reduced significantly โ down to โ while producing compressed determinants of the same size. The quasi-linearization method thus reduces the time required for the compression process, as the finite-regular part extraction algorithm ultiparameter analysis for instability problems Page of works by successively increasing the rank of the system (not all at once). However, for this problem the computation time is too small for any meaningful assessment: we will investigate computation time more fully in Section 6.3. Figure 3 shows the system flutter points superimposed on a contour plot [24]. As can be seen, there is an exact agreement between the direct solutions and the contour plot solutions (the intersections of Re(๐) = 0 and
Im(๐) = 0 ). Figure 2:
Modal damping plot for the undamped section model.
Figure 3:
Contour plot for the undamped section model, showing identical solutions from the two direct solvers. ultiparameter analysis for instability problems Page of We now simulate the damped system, in both forms (Eq. 13 and 14). Figure 4 shows a contour plot for the ๐ - ๐ form of this system (Eq. 16), with direct solutions from both solution methods. There is a single physical flutter point at ๐ ๐น = 1.650 , ๐ ๐น = 0.8336 s/rad , and a nonphysical flutter point at a small negative ๐ . The divergence point occurs at infinite ๐ and ๐ . As expected, the solutions from both direct solution methods are identical to each other and the solutions from the contour plot. Figure 5 shows a contour plot of the ฮฅ - ๐ form (Eq. 14) with direct solutions from both solution methods. This physical flutter point can be located at ฮฅ ๐น = 1.98 Hz , ๐ ๐น = 1.20 rad/s and the divergence point at ฮฅ ๐ท = 3.99 Hz . Again, the flutter points computed with direct solvers agree exactly with those seen on the contour plot, and the results from the system as a whole agree with those of the ๐ - ๐ form ( ๐ ๐น =1 ๐ ๐น โ = 0.83 s/rad and ๐ ๐น = ฮฅ ๐น ๐ ๐น โ = 1.65 ). Given the high computational complexity of these direct solution methods โ
๐ช(๐ ) โ we are interested in the maximum system size for which a direct solution is practical. We have already found that for ๐ = 2 the computational effort required is tiny, and so at the very least these direct solvers are useful for small reduced-order models as might be used in a preliminary design analysis. This alone is of use, as the directness of these solvers makes them ideal for use in optimization routines or other applications in which a large number of computations must be performed with limited user guidance. To gain a better understanding of the computational complexity characteristics of our methods, we simulate a series of systems of increasing matrix coefficient size. We generate random complex-valued matrices for the polynomial coefficients ( M , G , etc.), of size ๐ = 2 ๐ with ๐ โ {1, 2, โฆ , 5} . For robustness, we average the results for ๐ = 1 over random matrices, for ๐ = 2 over random matrices, for ๐ = 3 over random matrices; and for systems larger than this we generate only one matrix. Figure 6 shows the wall-clock solution time against system size for the direct solver with linearization, for a 64-bit Intel i7-4770 with 3.4 GHz processor and 16 GB RAM, running MATLAB R2014b. Note that the ๐ solution time denotes the time required to compute the ๐ components of the solution via a series of one-parameter eigenvalue problems. As can be seen, the compression process is the most expensive component of the algorithm, making up a constant fraction of about 65% of the total computation time over the entire range of ๐ . The GEP solution time is initially completely negligible but becomes more significant as system size increases: by ๐ = 32 it makes up 34% of the total computation time. The ๐ -solution and setup process are never of much significance. We expect from computational complexity theory that the gradient of the GEP-solution time curve in log-log units will be approximately 6.0 (Section 4.2). Fitting a linear curve through this data yields a gradient of 5.46, with an ๐ of . However, the GEP solution-time curve is slightly concave, with a maximum gradient of 6.47. The algorithm can then be said to have complexity ๐ช(๐ ) overall. Figure 7 shows the wall-clock solution time against system size for the direct solver with quasi-linearization. This simulation was run on the same Intel i7-4770 platform with MATLAB R2014b. The random coefficient matrices used are identical to those in Figure 6.And although now the quasi-linearized system in no longer singular, we still apply the compression algorithm to capture some of its overhead costs.. As can be seen, these costs are not large but are still more significant than the GEP and ๐ solution times for some system sizes. ultiparameter analysis for instability problems Page of Figure 4: Contour plot for the damped section model ( ๐ - ๐ form), showing identical solutions from the two direct solvers. Figure 5: Contour plot for the damped section model ( ฮฅ - ๐ form), showing identical solutions from the two direct solvers. ultiparameter analysis for instability problems Page of Figure 6: Wall-clock solution time against system size for the direct solver with linearization.
Figure 7: Wall-clock solution time against system size for the direct solver with quasi-linearization. ultiparameter analysis for instability problems Page of These GEP and ๐ solution times are effectively identical to those of the direct solver with linearization, as the algorithm is solving a GEP of exactly the same size, yielding exactly the same eigenvalues. However, the most striking aspect of Figure 7 is the fact that the setup time occupies the majority of the required computation time right up to ๐ = 16 (after which it is surpassed by the GEP solution time). The setup process involved using defining the system array (Eq. 32) in a cell array, and computing the necessary operator determinants of this array ( ฮ and ฮ ) using the modified Leibniz formula (Eq. 34). On investigating these two procedures we find that it is the computation of the Kronecker products within the operator determinant which occupies the greatest time. However, despite this, the actual computational complexity of the algorithm as a whole is lower than that of the direct solver with linearization. GEP solution time is still approximately ๐ช(๐ ) , though the concavity that was noted in Figure 6 is still present. Finally, Figure 8 shows a comparison of the total computation times for the two algorithms. As can be seen, below a system coefficient size of approximately 10, linearisation is faster than quasi-linearisation (at ๐ = 2 , over an order of magnitude faster). Above this coefficient size, quasi-linearisation is faster (at ๐ = 32 , twice as fast). To keep computation times below 10s, the system size must be below 11; to keep it below 1s, 8, and below 0.1s, 5. As it stands, the operator determinant method is not suitable for use with finite-element models or any other systems with a large number of degrees of freedom. It is, however, useful for the solution of simple reduced order models, e.g. in a preliminary design optimization. Figure 8: Total wall-clock solution time against system size for the two direct solver variants. ultiparameter analysis for instability problems Page of Discussion 7.1. Alternative solution methods
We have so far discussed the operator determinant method โ a solution method which is direct but computationally expensive. However, two other classes of method are also available. The Sylvester-Arnoldi type methods are a set of closely-related algorithms that are only valid for two-parameter problems, but are fast and can handle large systems [48]. They still use operator determinants to reformulate the system into a set of generalized eigenvalue problems, but instead of solving this set of GEPs directly, the solution procedure is optimized based on the given knowledge that the operator determinants consist of a difference of two Kronecker products. For example, each step of a Krylov subspace procedure solution procedure for the GEP reduces to the solution of a Sylvester equation. Several related algorithms may be devised this way, and in some cases the computational complexity of the solution process can be reduced all the way down to
๐ช(๐ ) [48]. However, this technique shows little potential for generalization to ๐ -parameter systems, as it relies on being able to reformulate the operator determinant GEP into a simple and well-known matrix equation. As the operator determinant definition becomes more complex, the resulting matrix equation changes and efficient solvers may not be available. There is also no easy extension for singular problems. Subspace methods for one-parameter eigenvalue problems are based around generating a series of linear spaces that eventually approximate one of the systemโs eigenspaces (the linear space of eigenvectors corresponding to a given eigenvalue). The Jacobi-Davidson and Rayleigh-Ritz methods are well-known one-parameter subspace methods, which can be generalized to apply to two-parameter systems [26,32,49โ51], though this is not without difficulties [49]. These methods do not invoke the operator determinants, and show potential for generalization both to ๐ -parameter and to polynomial systems. The Jacobi-Davidson is applicable to singular systems and has previously been tested on singular linearized aeroelastic stability problems [25]; but its performance was observed to be poorer than the operator determinant method. However the solution of singular MEPs has only been achieved very recently (2010 [29]), and so it is likely that the next few years will bring significant developments in this area. The core methodology presented in this paper can be extended to a wide variety of problems. Not only more complex two-parameter flutter models โ considered in Pons [24] and Pons and Gutschmidt [25] โ but also stability problems from a wide variety of fields, including systems with entirely different eigenvalue definitions. Indeed, any combination of model parameters in a stability model can indeed be treated as multiparameter eigenvalues; and scalar constraints on these parameters can cast as eigenvalue problems by introducing a scalar eigenvector [24]. For example, in aeroelastic model, given the location of a flutter point, we could solve for the sets of parameter values that could generate such a flutter point โ allowing us to perform model identification based on flutter point information. This could pave the way for a least-squares approach for overconstrained multiparameter systems. Alternatively, we could introduce a flight altitude / air density / Mach number parameter and compute points on the aeroelastic flutter envelope of the aircraft. The multiparameter formulation provides a versatile way of analysing stability problems in any combination of parameters. ultiparameter analysis for instability problems Page of Conclusion
In this paper we have demonstrated and discussed the use of multiparameter solution techniques for the solution of aeroelastic and related stability problems. We have introduced the link between multiparameter spectral theory and stability analysis, and we showed how this link can be used to reformulate stability problems with a complex-valued stability metric and a pertinent environmental parameter into a two-parameter eigenvalue problem. We demonstrated that this allows the direct solution of stability problems that are linear or polynomial in these parameters, and we discussed aspects of the solution process, including the linearization and quasi-linearization of polynomial problems, general N-parameter problems, computational costs and approaches to problem singularity. We discuss extensions to these methods, including the generalization to more complex problems; and further applications, including parameter identification and flight envelope computation. The application of multiparameter methods to stability problems โ in aeroelasticity and in other disciplines โ has the potential to provide a wide variety of new methods for stability analysis.
References [1] R. L. B ISPLINGHOFF , H. A SHLEY , R. L. H ALFMAN . Aeroelasticity . Addison-Wesley: Reading, MA, 1957. [2] H. J. H ASSIG . An approximate true damping solution of the flutter equation by determinant iteration.
Journal of Aircraft (11): 885โ889. DOI: 10.2514/3.44311. [3] D. H. H ODGES , G. A. P IERCE . Introduction to Structural Dynamics and Aeroelasticity . 2nd ed. Cambridge University Press: New York, 2011. [4] H. H
ADDADPOUR , R. D. F IROUZ -A BADI . True damping and frequency prediction for aeroelastic systems: The PP method.
Journal of Fluids and Structures (7): 1177โ1188. DOI: 10.1016/j.jfluidstructs.2009.06.006. [5] A. N AMINI , P. A LBRECHT , H. B OSCH . Finite element-based flutter analysis of cable-suspended bridges.
Journal of Structural Engineering (6): 1509โ1526. DOI: 10.1061/(ASCE)0733-9445(1992)118:6(1509). [6] P. C. C HEN . Damping perturbation method for flutter solution: the g-method.
AIAA Journal (9): 1519โ1524. DOI: 10.2514/2.1171. [7] R. L IND , M. B RENNER . Robust aeroservoelastic stability analysis . Springer-Verlag: London, 1999. [8] R. L IND . Match-Point Solutions for Robust Flutter Analysis.
Journal of Aircraft (1): 91โ99. DOI: 10.2514/2.2900. [9] D. B ORGLUND . Robust aeroelastic stability analysis considering frequency-domain aerodynamic uncertainty.
Journal of Aircraft (1): 189โ193. DOI: 10.2514/2.3074. [10] D. B ORGLUND . The ฮผ -k method for robust flutter solutions. Journal of Aircraft (5): 1209โ1216. DOI: 10.2514/1.3062. [11] D. B ORGLUND . Upper Bound Flutter Speed Estimation Using the ฮผ -k Method. Journal of Aircraft (2): 555โ557. [12] D. B ORGLUND , U. R INGERTZ . Efficient computation of robust flutter boundaries using the ฮผ -k method. Journal of Aircraft (6): 1763โ1769. [13] Y. G U , Z. Y ANG , W. W ANG , W. X IA . Dynamic Pressure Perturbation Method for Flutter Solution: the Mu-Omega Method. In Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural ultiparameter analysis for instability problems Page of Dynamics, and Materials Conference, Structures, Structural Dynamics, and Materials Conference , Palm Springs, CA, 2009. DOI: 10.2514/6.2009-2312. [14] Y. G U , Z. Y ANG . Generalized Mu-Omega Method with Complex Perturbation to Dynamic Pressure. In
Proceedings of the 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference , Orlando, FL, 2010. DOI: 10.2514/6.2010-2799. [15] D. B ORGLUND . Robust Aeroelastic Analysis in the Laplace Domain: The ยต-p Method. In
Proceedings of the 2007 International Forum on Aeroelasticity and Structural Dynamics (IFASD) , Stockholm, Sweden, 2007. [16] D. A FOLABI . Flutter analysis using transversality theory.
Acta mechanica (1โ4): 1โ15. DOI: 10.1007/BF01180214. [17] D. A FOLABI , R. M. V. P IDAPARTI , H. T. Y. Y ANG . Flutter Prediction Using an Eigenvector Orientation Approach.
AIAA Journal (1): 69โ74. DOI: 10.2514/2.353. [18] S. I RANI , S. S AZESH . A new flutter speed analysis method using stochastic approach.
Journal of Fluids and Structures : 105โ114. DOI: 10.1016/j.jfluidstructs.2013.03.018. [19] Y G U , X. Z HANG , Z. Y ANG . Robust flutter analysis based on genetic algorithm.
Science China Technological Sciences (9): 2474โ2481. [20] P. W. B LYTHE , I. H. H ERSZBERG . The Solution of Flutter Equations Using Neural Networks. In , Institution of Engineers, Australia, 1993; 415. [21] C. H. C HEN . Determination of flutter derivatives via a neural network approach.
Journal of Sound and Vibration (4): 797โ813. DOI: 10.1016/S0022-460X(02)01279-8. [22] C. H. C HEN , J. C. W U , J. H. C HEN . Prediction of flutter derivatives by artificial neural networks.
Journal of Wind Engineering and Industrial Aerodynamics (10): 1925โ1937. DOI: 10.1016/j.jweia.2008.02.044. [23] A. N ATARAJAN , โAeroelasticity of morphing wings using neural networks,โ Doctoral Dissertation, Virginia Polytechnic Institute and State University, Virginia, USA, 2002. [24] A P ONS , โAeroelastic flutter as a multiparameter eigenvalue problem,โ Masterโs Thesis, University of Canterbury, New Zealand, 2015. [25] A P ONS AND S G UTSCHMIDT . Multiparameter Solution Methods for Semi-Structured Aeroelastic Flutter Problems.
AIAA Journal , article in advance, 2017. DOI: 10.2514/1.J055447. [26] K. M EERBERGEN , C. S CHRรDER , H. V OSS . A Jacobi-Davidson method for two-real-parameter nonlinear eigenvalue problems arising from delay-differential equations.
Numerical Linear Algebra with Applications (5): 852โ868. DOI: 10.1002/nla.1848. [27] K. M EERBERGEN , A. S PENCE . Inverse Iteration for Purely Imaginary Eigenvalues with Application to the Detection of Hopf Bifurcations in Large-Scale Problems.
SIAM Journal on Matrix Analysis and Applications (4): 1982โ1999. DOI: 10.1137/080742890. [28] F. M. W HITE . Fluid mechanics . 6th ed. McGraw-Hill: New York, 2009. [29] A. M UHIฤ , B. P LESTENJAK . On the quadratic two-parameter eigenvalue problem and its linearization.
Linear Algebra and its Applications (10): 2529โ2542. DOI: 10.1016/j.laa.2009.12.022. [30] M. E. H OCHSTENBACH , A. M UHIฤ , B. P LESTENJAK . On linearizations of the quadratic two-parameter eigenvalue problem.
Linear Algebra and its Applications (8): 2725โ2743. DOI: 10.1016/j.laa.2011.07.026. ultiparameter analysis for instability problems Page of [31] F. T ISSEUR , K. M EERBERGEN . The Quadratic Eigenvalue Problem.
SIAM Review (2): 235โ286. DOI: 10.1137/S0036144500381988. [32] M. E. H OCHSTENBACH , B. P LESTENJAK . A Jacobi-Davidson Type Method for a Right Definite Two-Parameter Eigenvalue Problem.
SIAM Journal on Matrix Analysis and Applications (2): 392โ410. DOI: 10.1137/S0895479801395264. [33] B. P LESTENJAK . A continuation method for a weakly elliptic two-parameter eigenvalue problem.
IMA Journal of Numerical Analysis (1): 199โ216. DOI: 10.1093/imanum/21.1.199. [34] A. Q UARTERONI , R. S ACCO , F. S ALERI . Numerical mathematics . 2nd ed. Springer-Verlag Berlin: Berlin, Germany, 2007. [35] N. C OTTIN . Dynamic model updating โ a multiparameter eigenvalue problem.
Mechanical Systems and Signal Processing (4): 649โ665. DOI: 10.1006/mssp.2000.1362. [36] N. C OTTIN , J. R EETZ . Accuracy of multiparameter eigenvalues used for dynamic model updating with measured natural frequencies only.
Mechanical Systems and Signal Processing (1): 65โ77. DOI: 10.1016/j.ymssp.2004.10.005. [37] A. M UHIฤ , B. P LESTENJAK . On the singular two-parameter eigenvalue problem.
Electronic Journal of Linear Algebra : 420โ437. [38] M. E. H OCHSTENBACH , B. P LESTENJAK . Backward error, condition numbers, and pseudospectra for the multiparameter eigenvalue problem.
Linear Algebra and its Applications : 63โ81. DOI: 10.1016/S0024-3795(03)00613-X. [39] B. P LESTENJAK , A. M UHIฤ . MultiParEig 1.0 . MATLAB File Exchange: Ljubljana, Slovenia, 2015. [40] T. K Oล IR . Finite-dimensional multiparameter spectral theory: the nonderogatory case.
Linear Algebra and its Applications โ : 45โ70. DOI: 10.1016/0024-3795(94)90396-4. [41] B. B INDING , P. J. B ROWNE . Two parameter eigenvalue problems for matrices.
Linear Algebra and its Applications : 139โ157. DOI: 10.1016/0024-3795(89)90293-0. [42] F. V. A TKINSON . Multiparameter spectral theory.
Bulletin of the American Mathematical Society (1): 1โ28. DOI: 10.1090/S0002-9904-1968-11866-X. [43] F. V. A TKINSON . Multiparameter Eigenvalue Problems: Matrices and Compact Operators . Academic Press: London, 1972. [44] B. P. R YNNE . Multiparameter spectral theory and Taylorโs joint spectrum in Hilbert space.
Proceedings of the Edinburgh Mathematical Society (01): 127. DOI: 10.1017/S0013091500006635. [45] H. A NTON . Elementary Linear Algebra . 2nd ed. John Wiley & Sons: New York, New York State, USA, 1977. [46] H. A. E ISELT , C. L. S ANDBLOM . Linear Programming and its Applications . Springer-Verlag: Berlin, Germany, 2007. [47] L. N. T REFETHEN , D. B AU . Numerical linear algebra . Society for Industrial and Applied Mathematics: Philadelphia, PA, 1997. [48] K. M EERBERGEN , B. P LESTENJAK . A Sylvester-Arnoldi type method for the generalized eigenvalue problem with two-by-two operator determinants . Department of Computer Science, Katholieke Universiteit Leuven: Leuven, Belgium, 2014. [49] M. E. H OCHSTENBACH , T. K OSIR , B. P LESTENJAK . A Jacobi-Davidson Type Method for the Two-Parameter Eigenvalue Problem.
SIAM Journal on Matrix Analysis and Applications (2): 477โ497. DOI: 10.1137/S0895479802418318. ultiparameter analysis for instability problems Page of [50] M. E. H OCHSTENBACH , B. P LESTENJAK . Harmonic RayleighโRitz extraction for the multiparameter eigenvalue problem.
Electronic Transactions on Numerical Analysis : 81โ96. [51] M. E. H OCHSTENBACH , A. M UHIฤ , B. P LESTENJAK . JacobiโDavidson methods for polynomial two-parameter eigenvalue problems.
Journal of Computational and Applied Mathematics : 251โ263. DOI: 10.1016/j.cam.2015.04.019.: 251โ263. DOI: 10.1016/j.cam.2015.04.019.