2. background#

2.1. Problem#

We consider the quadratic modal problem (QEP)

Find \((\lambda ,u)\) such as \(({\lambda }^{2}\mathrm{B}+\lambda \mathrm{C}+\mathrm{A})\mathrm{u}\mathrm{=}\mathrm{0,}\mathrm{u}\mathrm{\ne }0\) (2.1-1)

where \(\mathrm{A},\mathrm{B}\) and \(\mathrm{C}\) have square matrices of size \(n\), with real or complex coefficients, symmetric or not. In mechanics, this type of problem corresponds in particular to the study of the free vibrations of a damped and/or rotating structure. For this structure, we look for the eigenvalues \({\lambda }_{i}\) (and their associated eigenvectors \({\mathrm{u}}_{i}\)) that are closest, in the complex plane, to a given reference value (the « shift » \(\sigma\)) to know if an exciting force can create a resonance. In this standard case,

  • The matrix \(\mathrm{A}\) is the stiffness matrix, noted \(\mathrm{K}\), real symmetric (possibly increased by the complex symmetric matrix, noted \({\mathrm{E}}_{\mathit{hyst}}\), if the structure has hysteretic damping): \(\mathrm{A}\mathrm{=}\mathrm{K}+{\mathrm{E}}_{\mathit{hyst}}\). So \(\mathrm{A}\) is symmetric real or complex.

  • The \(\mathrm{B}\) matrix is the mass or inertia matrix denoted \(\mathrm{M}\) (real symmetric).

The matrix :math:`mathrm{C}` includes possible gyroscopic effects and those of viscous damping*via the combination:

\(\mathrm{C}\mathrm{:}\mathrm{=}{\mathrm{E}}_{\mathit{visq}}+\xi \mathrm{G}\) (2.1-2)

with \({\mathrm{E}}_{\mathrm{visq}}\), a (real symmetric) damping matrix induced by dissipative forces, \(\mathrm{G}\) gyroscopy matrix (real anti-symmetric) and \(\xi\) a real parameter representative of the rotation speed. So \(\mathrm{C}\) is potentially not symmetric real.

This type of problem is activated by the keyword TYPE_RESU =” DYNAMIQUE “. In the general case of QEP (2.1-1) and using only their property of symmetry and matrix arithmetic, we can decompose the perimeter of the Code_Aster modal operators as follows.

_images/Shape1.gif _images/Shape2.gif _images/Shape3.gif

\(\mathrm{A}\) / \(\mathrm{C}\)

Symmetric real

Non-symmetric real

Complex

Symmetric real

Most common case SIMULT/INV with no restrictions on methods; Real modes or complexes \((\lambda ,\stackrel{ˉ}{\lambda })\).

SIMULTavec “SORENSEN” /”QZ “;

Real or complex modes \((\lambda ,\stackrel{ˉ}{\lambda })\).

Untreated case

Non-symmetric real

SIMULTavec “SORENSEN” /”QZ “; Real modes or complexes \((\lambda ,\stackrel{ˉ}{\lambda })\).

SIMULTavec “SORENSEN” /”QZ “

Real or complex modes \((\lambda ,\stackrel{ˉ}{\lambda })\).

Untreated case

Symmetric complex

SIMULTavec “SORENSEN” /”QZ “ Real, complex, any modes, or \((\lambda ,\stackrel{ˉ}{\lambda })\).

Untreated case

Untreated case

Other complexes (Hermitian, non-symmetric…)

Untreated case

Untreated case

Untreated case

Untreated case

Table 2.1-1. Scope of use of the Aster operator CALC_MODES according to the analysis methods, depending on the properties of the matries \(\mathrm{A}\) and \(\mathrm{C}\) of the QEP (2.1-1) . \(\mathrm{B}\) is any real.

In the absence of hysteretic damping, the QEP standard to be solved therefore consists of symmetric real matrixs. These eigenvalues are either real or complex conjugated in pairs \((\lambda ,\stackrel{ˉ}{\lambda })\). The eigenvectors then potentially have complex components.

In the presence of hysteretic damping, QEP becomes complex and loses the previous property [1] _ . These modes are then real, complex in pairs or mismatched.

Notes:

*Unlike GEP [Boi09] processed by Code_Aster, the keyword TYPE_RESU =” FLAMBEMENT “*is not legal within the framework of QEP. *

  • Computatively, the scope of option [“CENTRE”, “PLUS_PETITE”, “”, “TOUT”] has been extended to real non-symmetric matrices \(\mathrm{A},\mathrm{B}\text{et}\mathrm{C}\) of QEP. This is not yet the case for option [“PROCHE”, “AJUSTE”] which is limited, for the moment, to symmetric real matrices. The latter cannot therefore take into account gyroscopic effects.

The damping modeling [Lev96] in Code_Aster can be broken down into two classes: proportional viscous damping (called Rayleigh damping :math:`{mathrm{E}}_{mathit{visq}}mathrm{:}mathrm{=}alpha mathrm{K}+beta mathrm{M}`) or hysteretic damping (:math:`{mathrm{E}}_{mathit{hist}}mathrm{:}mathrm{=}(1+ieta )mathrm{K}`). Each is available at the global level of the structure or by means of localized amortization (groups of cells, ad hoc discrete elements). The latter make it possible to better represent the heterogeneity of the structure in relation to damping. However, there is still the delicate question of identifying the coefficients*:math:`(alpha ,beta eta )`*and their influences on the result. *

2.2. Properties of clean modes#

The preceding table 2.2-1 does not take into account particular cases (pure gyroscopy, over-damped structure, etc.). For the record, we will list some of them that are often found in literature. But it should be remembered that the context of Code_Aster models (almost systematic use of double Lagranges) and the sorting of eigenmodes programmed in modal operators do not take these specificities into account. Currently in QEP, we only retain coupled modes \((\lambda ,\stackrel{ˉ}{\lambda })\) and we keep the one with a positive imaginary part. Depending on the case, the presence of eigenvalues of other types (real or complex non-conjugated) is indicated by a ALARME or for information purposes.

We therefore have (using the elements and the notations from §3.3 of [Boi09] and the paper [TM01]), we have the distribution of the following properties. These are sometimes cumulative. Figures 2.2-1 and 2.2-2 illustrate, on canonical cases, some cases.

Table 2.2-1. Properties of the modes of a standard QEP (of size \(n\) ) according to those of its matrices. Possible links with Code_Aster models.

_images/1000000000000252000000EF30C4AF50A8D898E4.png

Figure 2.2-1. Examples of the spectrum of a QEP \((n=8)\) of a simplified nuclear enclosure model [TM01]. Figure on the left, we are in case No. 3 of table 2.2-1. The one on the right, we add hysteretic damping to \(\mathrm{K}\) so mismatched modes appear (general case no. 1) .

_images/10000000000001BD00000156B6E0569B23778247.png _images/10000000000001D7000001572BEECCCCE8EF3828.png

Figure 2.2-2. Examples of the spectrum of a QEP \((n=50)\) of a damped mass-spring system [TM01]. Figure on the left, we are in case No. 5 of table 2.2-1. The one on the right, we add the over-damping constraint and the spectrum becomes real negative in two packets (case no. 6).

2.3. Specificities of QEP and elements of theory#

A major algebraic specificity between QEP and the other lower-level modal problem classes (GEP and SEP) is the existence of twice as many eigenmodes ( \(\mathrm{2n}\) ) as the size of the discretized problem. Eigenvectors can no longer be linearly independent, in addition, eigenvalues can be finite or infinite and all « this little world » generally lives in the complex plane [3] _ . It is then necessary define robust heuristics to filter the eigenvalues desired by the user and distinguish the real ones, the complex ones, the mismatched complexes from the coupled ones…

Another complication inherent in QEP is algorithmic. There is no Schur decomposition (or generalized Schur) as for SEP (resp. GEP) on which the resolution algorithm will be able to rely. For example, for SEP \(\mathrm{A}\mathrm{u}\mathrm{=}\lambda \mathrm{u}\), this decomposition assures us of the existence of a unitary matrix (therefore well-conditioned and easily invertible) \(\mathrm{U}\), allowing the rewriting of the working matrix \(\mathrm{A}\) in a form that is easier to handle [4] _ : the upper triangular matrix \(\mathrm{T}\).

\(\mathrm{U}\mathrm{A}{\mathrm{U}}^{\text{*}}\mathrm{=}\mathrm{T}\) (2.3-1)

The QEP are therefore a very particular and important class [5] _ nonlinear problems. Solving them is less routine than for other classes of modal problems. In particular, few methods make it possible to solve the problem directly, it is often necessary to go through a GEP (via a *ad hoc linearization) then a SEP (spectral transformation) which induces a loss of spectral information and numerical instabilities (propagation of rounding errors).

In theoretical literature, rather than QEP (2.1-1), we are interested in the matrix polynomial (quadratic in \(\lambda\))

\(\mathrm{Q}(\lambda )\mathrm{:}\mathrm{=}({\lambda }^{2}\mathrm{B}+\lambda \mathrm{C}+\mathrm{A})\) (2.3-2)

We also call it \(\lambda \text{-matrice}\) or dynamic matrix and we are looking for its spectrum denoted \(\mathit{spectre}(\mathrm{Q})\mathrm{:}\mathrm{=}\mathrm{\{}\lambda \mathrm{\in }C\mathrm{/}\mathit{det}\mathrm{Q}(\lambda )\mathrm{=}0\mathrm{\}}\). When this determinant is zero regardless of \(\lambda\), \(\lambda \text{-matrice}\) is said to be singular. Otherwise, it is said to be regular (we generally place ourselves in this framework!).

The characteristic polynomial of the problem is therefore written in the form \(\mathit{det}\mathrm{Q}(\lambda )\mathrm{=}\mathit{det}(\mathrm{B}){\lambda }^{\mathrm{2n}}+\) lower-order terms in \(\lambda\). Thus, as soon as \(\mathrm{B}\mathrm{=}\mathrm{M}\) is singular, this polynomial admits \(r<\mathrm{2n}\) finite roots to which must be added \(\mathrm{2n}-r\) infinite. On the other hand, these distinct eigenvalues can share the same eigenvectors. This is the « dual situation » of multiple eigenvalues in SEP where the same eigenvalue is attached to several eigenvectors. So for the canonical QEP

\(\mathrm{A}\mathrm{=}{\mathrm{I}}_{\mathrm{3,}}\mathrm{B}\mathrm{=}\left[\begin{array}{ccc}0& 6& 0\\ 0& 6& 0\\ 0& 0& 1\end{array}\right]\text{et}C\mathrm{=}\left[\begin{array}{ccc}1& \mathrm{-}6& 0\\ 2& \mathrm{-}7& 0\\ 0& 0& 0\end{array}\right]\) (2.3-3)

We have the following spectrum

\(\begin{array}{ccccccc}i& 1& 2& 3& 4& 5& 6\\ {\lambda }_{i}& 1\mathrm{/}3& 1\mathrm{/}2& 1& i& \mathrm{-}i& \mathrm{\infty }\\ {\mathrm{u}}_{i}& \left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]& \left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]& \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]& \left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]\end{array}\) (2.3-4)

From the uniqueness of \(\mathrm{B}\) comes the infinite eigenvalue. We also note the sharing of eigenvectors by the first two eigenvalues and the 4th/5th.

To more easily understand the infinite eigenvalues of \(\mathrm{Q}(\lambda )\), we associate them with the zero eigenvalues of the inverse polynomial

\(\mathit{rev}\mathrm{Q}(\lambda )\mathrm{:}\mathrm{=}{\lambda }^{2}\mathrm{Q}(\frac{1}{\lambda })\mathrm{=}(\mathrm{B}+\lambda \mathrm{C}+{\lambda }^{2}\mathrm{A})\) (2.3-5)

These finite eigenvalues of \(\mathit{rev}\mathrm{Q}(\lambda )\) are easier to calculate and share the same spectral characteristics (geometric and algebraic multiplicities, eigenspace…) as the initial infinite values.

For more information on these theoretical aspects one can consult, for example, the works of P.Lancaster [GLR82] [LT85] and the thesis of D.S.Mackey [Mac06] (and its numerous bibliographical references).

2.4. Linearization strategy#

2.4.1. Introduction#

The standard approach to solving a modal problem is to transform it to reveal canonical matrices revealing the modes sought. For example, the Schur decomposition for SEP or the generalized one for GEP. Unfortunately, this practice cannot be generalized to higher-order modal problems. Two solutions are then available:

  • Solve the non-linear problem directly (polynomial root, factorization of \(\lambda \text{-matrice}\), Newtons-type methods…).

  • Linearize it into a GEP, then process it via using the appropriate methods [Boi09].

The first strategy is available in the options [“PROCHE”, “AJUSTE”] (Müller-Traub method) to capture a first estimate of eigenvalues. We then mix it up with the following strategy, because to refine the estimates and accelerate convergence, we use an inverse power method (Jennings variant) on the linearized QEP. Some authors offer comprehensive direct approaches. They are often based on Newton’s variants. But their convergence is slow (one mode at a time) and not always assured (cf. Kublanovskaya 1970, G.Peters & J.H.Wilkinson 1979 and A.Ruhe 1973).

The linearization strategy is therefore often preferred. It doubles the size of the modal problem to be treated, changes its nature (at the risk of losing spectral information) [6] _ ) and introduces numerical instabilities (sensitivity to rounding errors [7] _ and to the numerical parameters [8] _ ). But it has the immense advantage of reusing all the « digital artillery » already deployed for the GEP [9] _ .

It is thus possible to simultaneously capture the entire spectrum (cf. the QZ method usable with the options [“CENTRE”, “PLUS_PETITE”, “TOUT”]) or only the part centered on an area of interest (cf. IRAM or Lanczos) using the appropriate spectral transformation.

2.4.2. Principle#

Now that the context has been outlined, to conform to the terminology used in the literature, we will now simplify the notations by considering the QEP associated with the \(\lambda \text{-matrice}\) \(\mathrm{Q}(\lambda )\mathrm{:}\mathrm{=}({\lambda }^{2}\mathrm{M}+\lambda \mathrm{C}+\mathrm{K})\). This matrix is also called a « dynamic matrix » in structural dynamics jargon (and in Code_Aster computer documents/sources). It is the counterpart to the one described in the note dedicated to GEP: \(\mathrm{Q}(\lambda )\mathrm{:}\mathrm{=}(\mathrm{K}\mathrm{-}\lambda \mathrm{M})\).

Definition 1

From a theoretical point of view, linearizing \(\mathrm{Q}(\lambda )\) is the same as finding another \(\lambda \text{-matrice}\) (square complex of size \(\mathrm{2n}\)) of the type \(\mathrm{A}\mathrm{-}\lambda \mathrm{B}\) such that there are two other \(\lambda \text{-matrices}\) (\(\mathrm{E}(\lambda )\) and \(\mathrm{F}(\lambda )\)) with a non-zero constant determinant, verifying the relationship

\(\left[\begin{array}{cc}\mathrm{Q}(\lambda )& {0}_{n}\\ {0}_{n}& {\mathrm{I}}_{n}\end{array}\right]\mathrm{=}\mathrm{E}(\lambda )(\mathrm{A}\mathrm{-}\lambda \mathrm{B})\mathrm{F}(\lambda )\) (2.4-1)

The matrices \(\mathrm{A}\) and \(\mathrm{B}\) are called the « companion » matrices for the linearization of QEP.

Clearly, the spectra of the initial QEP \((\mathrm{Q}(\lambda ))\) and that of its linearization \((\mathrm{A}\mathrm{-}\lambda \mathrm{B})\) coincide. However this decomposition is not unique and it is necessary to choose, if possible, the one that:

  • Preserves the properties of the initial matrices (arithmetic, symmetry…),

Has the**least sensitivity to rounding errors* (« stable backward »: balanced, well-conditioned, invertible matrices…),

Manipulates**canonical matrix* (identity, triangular…),

Adapts to a**large perimeter of QEP * (viscous damping, gyroscopic, over-damping…).

A more intuitive way to understand linearization techniques is to introduce a variable change of the type \(\mathrm{v}\mathrm{=}\lambda \mathrm{u}\) [10] _ in (2.1-1) and rearranging it in matrix form:

\((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{0}_{n}& {\mathrm{I}}_{n}\\ \mathrm{-}\mathrm{K}& \mathrm{-}\mathrm{C}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{\mathrm{I}}_{n}& {0}_{n}\\ {0}_{n}& \mathrm{M}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \mathrm{v}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-2)

This « very classical » linearization falls within the framework of definition 1 if we use \(\lambda \text{-matrices}\):

\(\mathrm{E}(\lambda )\mathrm{:}\mathrm{=}\left[\begin{array}{cc}\mathrm{-}(\mathrm{C}+\lambda \mathrm{M})& \mathrm{-}{\mathrm{I}}_{n}\\ {\mathrm{I}}_{n}& {0}_{n}\end{array}\right],\mathrm{F}(\lambda )\mathrm{:}\mathrm{=}\left[\begin{array}{cc}{\mathrm{I}}_{n}& {0}_{n}\\ \lambda {\mathrm{I}}_{n}& {\mathrm{I}}_{n}\end{array}\right]\) (2.4-3)

In fact, the linearization of a QEP follows the general framework of that of matrix polynomials. By setting \(\mathrm{Q}(\lambda )\mathrm{:}\mathrm{=}\mathrm{\sum }_{i\mathrm{=}1}^{k}{\lambda }_{i}{\mathrm{X}}_{i}\) (here \(k\mathrm{=}\mathrm{2,}{\mathrm{X}}_{i}\mathrm{=}\mathrm{M}\) etc), we are always sure to find two linearizations, called « companions », of the form \({\mathrm{A}}_{1}\mathrm{-}\lambda {\mathrm{B}}_{1}\) and \({A}_{2}-\lambda {B}_{2}\) with

\(\begin{array}{cc}{\mathrm{B}}_{1}\mathrm{=}{\mathrm{B}}_{2}\mathrm{:}\mathrm{=}\left[\begin{array}{cc}\mathrm{-}{\mathrm{X}}_{k}& {0}_{(k\mathrm{-}1)n}\\ {0}_{(k\mathrm{-}1)n}& \mathrm{-}{\mathrm{I}}_{(k\mathrm{-}1)n}\end{array}\right],& \\ {\mathrm{A}}_{1}\mathrm{:}\mathrm{=}\left[\begin{array}{cccc}{\mathrm{X}}_{k\mathrm{-}1}& {\mathrm{X}}_{k\mathrm{-}2}& \text{...}& {\mathrm{X}}_{0}\\ {\mathrm{I}}_{n}& {0}_{n}& \text{...}& {0}_{n}\\ \text{...}& \text{...}& \text{...}& \text{...}\\ {0}_{n}& \text{...}& \mathrm{-}{\mathrm{I}}_{n}& {0}_{n}\end{array}\right],& {\mathrm{A}}_{2}\mathrm{:}\mathrm{=}\left[\begin{array}{cccc}{\mathrm{X}}_{k\mathrm{-}1}& \mathrm{-}{\mathrm{I}}_{n}& \text{...}& {0}_{n}\\ {\mathrm{X}}_{k\mathrm{-}2}& {0}_{n}& \text{...}& \text{...}\\ \text{...}& \text{...}& \text{...}& \mathrm{-}{\mathrm{I}}_{n}\\ {\mathrm{X}}_{0}& {0}_{n}& \text{...}& {0}_{n}\end{array}\right]\end{array}\) (2.4-4)

In practice, linearizations that are more adapted to the case being treated are preferred (cf. next paragraph).

2.4.3. A few examples#

In the literature, the following linearizations** are often mentioned (with \(N\) a complex regular square matrix of size \(n\)):

(\({L}_{1}\)) \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{0}_{n}& \mathrm{N}\\ \mathrm{-}\mathrm{K}& \mathrm{-}\mathrm{C}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{N}& {0}_{n}\\ {0}_{n}& \mathrm{M}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-5)

(\({L}_{2}\)) \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{K}& {0}_{n}\\ {0}_{n}& \mathrm{N}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{C}& \mathrm{M}\\ \mathrm{N}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-6)

The choice depends on the potential singularity of matrices \(\mathrm{K}\) and \(\mathrm{M}\). We want to be able to invert the diagonal terms, so if \(\mathrm{K}\) is singular (resp. \(\mathrm{M}\)), we use (\({L}_{1}\)) (resp. (\({L}_{2}\))). To balance the terms of the matrices and facilitate the manipulations, we often take \(\mathrm{N}\mathrm{=}\alpha {\mathrm{I}}_{n}\) with \(\alpha\) real balancing factor (cf. [Boi09] annex 1) constructed from the matrix norms of the other terms \(\alpha \mathrm{=}\alpha (\mathrm{\parallel }\mathrm{K}\mathrm{\parallel },\mathrm{\parallel }\mathrm{M}\mathrm{\parallel },\mathrm{\parallel }\mathrm{C}\mathrm{\parallel })\). To work with a symmetric GEP (or even positive definite or semi-definite), you can also initialize the auxiliary matrix \(\mathrm{N}\text{à}\alpha \mathrm{K}\text{ou}\beta \mathrm{M}(\alpha ,\beta \mathrm{\in }\mathrm{\Re })\).

To solve QEP with gyroscopic effects more effectively, the literature recommends the « Hamiltonian » counterparts [11] _ «  previous linearizations. Let the « companion » matrices be

(\({L}_{3}\)) \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{K}& {0}_{n}\\ \mathrm{C}& \mathrm{N}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{0}_{n}& \mathrm{K}\\ \mathrm{-}\mathrm{M}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-7)

(\({L}_{4}\)) \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{0}_{n}& \mathrm{-}\mathrm{K}\\ \mathrm{M}& \mathrm{N}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{M}& \mathrm{C}\\ {0}_{n}& \mathrm{M}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-8)

Here the objective is not to respect a property of symmetry but rather those « Hamiltonian » [12] _ «  of the problem. If \(\mathrm{K}\) is singular (resp. \(\mathrm{M}\)), we use (\({L}_{4}\)) (resp. (\({L}_{3}\))).

In the case of a QEP called « T-palindromic » [Mac06] \((\mathrm{K}\mathrm{=}{\mathrm{M}}^{T}\text{et}\mathrm{C}\mathrm{=}{\mathrm{C}}^{T})\) we prefer linearization

(\({L}_{5}\)) \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{K}& \mathrm{K}\\ \mathrm{C}\mathrm{-}\mathrm{M}& \mathrm{K}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{M}& \mathrm{K}\mathrm{-}\mathrm{C}\\ \mathrm{-}\mathrm{M}& \mathrm{-}\mathrm{M}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-9)

2.4.4. Strategies selected in Code_Aster#

We therefore have a very wide variety of choices to adapt the linearization of QEP to each situation. For the moment, Code_Aster does not offer the user a parameter that allows the user to modulate only the linearization. This is fixed for a given modal solver:

  • With OPTIONparmi [“PROCHE”, “”, “AJUSTE”] (2ndphase) (2ndphase) the transition QEP/GEP is carried out*via*the companion matrices of (L1), by asking :math:`mathrm{N}mathrm{=}{mathrm{I}}_{n}mathrm{.}` The linearized system is not symmetric but it is not a prerequisite for**applying an inverse power method.

(L1) “\((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{0}_{n}& {\mathrm{I}}_{n}\\ \mathrm{-}\mathrm{K}& \mathrm{-}\mathrm{C}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}{\mathrm{I}}_{n}& {0}_{n}\\ {0}_{n}& \mathrm{M}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-10)

  • With OPTIONparmi [“CENTRE”, “”, “PLUS_PETITE”] and METHODE =” TRI_DIAG “or “SORENSEN “, the passage QEP/GEP is carried out via the companion matrices of (L2) by setting \(\mathrm{N}\mathrm{=}\mathrm{M}\). If the stiffness, mass and damping matrices are real symmetric, we fall back into the perimeter of the standard GEP (real symmetric). The methods of Lanczos and Sorensen are then accessible. If \(\mathrm{K}\) becomes complex or if one of them is not symmetric, only IRAM remains in the race (with QZ see next item)!

(L2) “\((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{K}& {0}_{n}\\ {0}_{n}& \mathrm{M}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{C}& \mathrm{M}\\ \mathrm{M}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-11)

However, when \(\mathrm{M}\) is not invertable (due to blocking Lagranges for example), the inversion of the diagonal term of \(\mathrm{A}\) is no longer assured. To overcome this problem, a « regularization » of the matrix \(\mathrm{M}\), noted \({M}_{R}\), is introduced into the diagonal term. This has the property of cancelling, by multiplication, all the components of the core of \(\mathrm{M}\mathrm{:}{\mathrm{M}}_{R}\mathrm{u}\mathrm{=}0\text{si}\mathrm{u}\mathrm{\in }\mathit{Ker}(\mathrm{M})\). Thus, we can define \({\mathrm{M}}_{R}^{\mathrm{-}1}\) which formally verifies the \({\mathrm{M}}_{R}^{\mathrm{-}1}\mathrm{M}\mathrm{\approx }\mathrm{M}{\mathrm{M}}_{R}^{\mathrm{-}1}\mathrm{\approx }{\mathrm{I}}_{n}\) property on the « non-core components of \(\mathrm{M}\) » of the manipulated vectors. Because of the particular structure of the dualized matrix \(\mathrm{M}\) (cf. [Boi09] §3.2), this manipulation amounts to proscribing components related to Lagranges and blocked degrees of freedom (the \({\mathrm{v}}_{i}^{\mathrm{2,3}}\) of the demonstration of property no. 4 [Boi09]). This does not hinder the process since we are only interested in physical modes and not in these modeling artifacts [13] _ .

(L2) » \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{K}& {0}_{n}\\ {0}_{n}& {\mathrm{M}}_{R}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{C}& \mathrm{M}\\ \mathrm{M}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-12)

  • With OPTIONparmi [“CENTRE”, “”, “PLUS_PETITE”, “TOUT”] and METHODE =”QZ”,* we preferred to avoid these complications. The objective is to **emphasize numerical robustness to deal with small cases (\(<{10}^{3}\) degrees of freedom). The scenario (L2) was empirically selected by posing \(\mathrm{N}\mathrm{=}\alpha {\mathrm{I}}_{n}\), with, \(\alpha \mathrm{:}\mathrm{=}\frac{1}{\mathrm{3n}}(\mathrm{\parallel }\mathrm{K}\mathrm{\parallel }+\mathrm{\parallel }\mathrm{M}\mathrm{\parallel }+\mathrm{\parallel }\mathrm{C}\mathrm{\parallel })\)

(L2) * \((\underset{\mathrm{A}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{-}\mathrm{K}& {0}_{n}\\ {0}_{n}& \alpha {\mathrm{I}}_{n}\end{array}\right]}}\mathrm{-}\lambda \underset{\mathrm{B}}{\underset{\underbrace{}}{\left[\begin{array}{cc}\mathrm{C}& \mathrm{M}\\ \alpha {\mathrm{I}}_{n}& {0}_{n}\end{array}\right]}})\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]\mathrm{=}{0}_{\mathrm{2n}}\) (2.4-13)

2.5. Implementation in Code_Aster#

By design, the process of a quadratic modal calculation is very similar to that of a generalized problem. We will therefore refer to the functional diagram for solving GEP described in [Boi09] (cf. §3.8). Its stages are broken down as follows:

2.5.1. Pre-treatments (linearization, shift calculation)#

Determination of the linearized form and the associated spectral transformation according to the method selected (power, Lanczos, IRAM and QZ) and the search option “AJUSTE”/”PROCHE” or “PLUS_PETITE”/”CENTRE”.

Note that since the natural modes are in the complex plane, the Sturm test currently implemented is inoperative. So the heuristic associated with OPTION “PROCHE” and “AJUSTE” is not based on modal positions, and the “BANDE” option is forbidden. Likewise for post-verification based on the Sturm test [Boi12].

For the subspace method (OPTION among [“CENTRE”, “”, “PLUS_PETITE”] and METHODE =” TRI_DIAG “or” SORENSEN “), once the linearized problem \(({L}_{2})’’\) is updated, a « shift-and-invert » spectral transformation is applied (cf. [Boi09] §3.6/§5). This will allow us to complete the transformation into a SEP [14] _ \(({S}_{1})\), to direct the search for the spectrum in an area of interest to the user and to better separate the eigenvalues.

\(({S}_{1})\) \(\underset{{\mathrm{A}}_{\sigma }}{\underset{\underbrace{}}{{(\mathrm{A}\mathrm{-}\sigma \mathrm{B})}^{\mathrm{-}1}}}\mathrm{B}\underset{\mathrm{W}}{\underset{\underbrace{}}{\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]}}\mathrm{=}\underset{\mu }{\underset{\underbrace{}}{\frac{1}{\lambda \mathrm{-}\sigma }}}\underset{\mathrm{W}}{\underset{\underbrace{}}{\left[\begin{array}{c}\mathrm{u}\\ \lambda \mathrm{u}\end{array}\right]}}\) (2.5-1)

Does the inverse power method apply a specific spectral transformation to the linearized problem \(({L}_{1})’\).

As for the QZ method (OPTION among [“CENTRE”, “”, “PLUS_PETITE”, “TOUT”] and METHODE =”QZ”), this spectral transformation is superfluous. The algorithm directly treats the linearized problem \({({L}_{2})}^{\text{*}}\) and calculates all of its modes. Depending on the user’s wishes (OPTION, CALC_FREQ =_F (NMAX_FREQ…)), it then filters the results.

Notes:

  • As with linearization techniques, no user parameter yet allows changing spectral transformations. In quadratic, it would be interesting to use transformations that are more adapted to the complex plane than the traditional « shift-and-invert ». For example, a Cayley transformation to make conjugate eigenvalue pairs converge \((\lambda ,\stackrel{ˉ}{\lambda })\) simultaneously

\(\underset{{\mathrm{A}}_{\sigma }}{\underset{\underbrace{}}{(\mathrm{A}\mathrm{-}\sigma \mathrm{B}){\mathrm{B}}^{\mathrm{-}1}(\mathrm{A}\mathrm{-}\stackrel{ˉ}{\sigma }\mathrm{B}){\mathrm{B}}^{\mathrm{-}1}}}\mathrm{w}\mathrm{=}\underset{\mu }{\underset{\underbrace{}}{\frac{1}{(\lambda \mathrm{-}\sigma )(\stackrel{ˉ}{\lambda }\mathrm{-}\sigma )}}}\mathrm{W}\) (2.5-2)

  • This Cayley strategy is being generalized to treat gyroscopy (see V.Mehrman and D.Watkins 2001) and simultaneously converge quadruplets \((\lambda ,\stackrel{ˉ}{\lambda },-\lambda ,-\stackrel{ˉ}{\lambda })\)

  • Some authors also propose to invert the linearization phase and the spectral transformation phase [Mee07].

2.5.2. Factorization of dynamic matrices#

As already pointed out, these dynamic matrix factorizations like \(\mathrm{Q}(\lambda )\) are not used to implement a Sturm test. The latter being illegal for QEP. However, they remain a very present and very expensive ingredient [15] _ for:

  • Estimate the cost function of the heuristic part (Müller-Traub method) with OPTION =” AJUSTE “.

  • When you need to manipulate a work matrix with inverses. To be more effective, we then look for a formulation that only involves \(\mathrm{Q}{(\lambda )}^{\mathrm{-}1}\) (OPTION among [“PROCHE”, “AJUSTE”] + SOLVEUR_MODAL =_F (OPTION_INV +” DIRECT “), or else OPTION among [” CENTRE “,” “,” « ”,” PLUS_PETITE “] + =_F (SOLVEUR_MODAL METHODE =” TRI_DIAG “/” SORENSEN “)).

The global approach (OPTION among [“CENTRE”, “PLUS_PETITE”, “”, “TOUT”] + SOLVEUR_MODAL =_F (METHODE =”QZ”)) is not concerned with this preliminary factorization of the dynamic matrix.

2.5.4. Audit post-treatments#

This last part includes the only post-processing that verifies the smooth running of the calculation (in the absence of a Sturm test adapted to the complex case). We calculate the relative error on the norm [17] _ from the initial QEP residual

\(\begin{array}{c}\mathrm{u}\mathrm{\Leftarrow }\frac{\mathrm{u}}{{〚\mathrm{u}〛}_{\mathrm{\infty }}}\\ \text{Si}\mathrm{\mid }\lambda \mathrm{\mid }>\text{SEUIL\_FREQ}\text{alors}\\ \frac{{∥{\lambda }^{2}\mathrm{Mu}+\lambda \mathrm{Cu}+\mathrm{Ku}∥}_{2}}{{∥\mathrm{Ku}∥}_{2}}\text{?}<\text{SEUIL,}\\ \text{Sinon}\\ {∥{\lambda }^{2}\mathrm{Mu}+\lambda \mathrm{Cu}+\mathrm{Ku}∥}_{2}\text{?}<\text{SEUIL,}\\ \text{Fin si.}\end{array}\)

Algorithm 1. Residue norm test.

This sequence is parameterized by the keywords SEUIL and SEUIL_FREQ, belonging to the keyword factors VERI_MODE and CALC_FREQ respectively. On the other hand, this post-processing is activated by initializing STOP_ERREUR to “OUI” (by default value) in the keyword factor VERI_MODE. When this rule is activated and not respected, the calculation stops, otherwise the error is just signalled by an alarm. Of course, we can only highly recommend that you do not deactivate this pass-it-all setting!

2.5.5. Display in the message file#

In the message file, information relating to the modes selected is mentioned. For example, in the case of a QEP, we specify the modal solver used and the list of \({\lambda }_{i}\) eigenvalues retained in ascending order of imaginary part \(\frac{\text{Im}({\lambda }_{i})}{2\pi }\) (column FREQUENCE). Then we specify its damping \({\xi }_{i}=\frac{-\text{Re}({\lambda }_{i})}{{\lambda }_{i}}\) and its error norm (cf. algorithm No. 1).


CALCUL MODAL: METHODE GLOBALE FROM TYPE QR

ALGORITHME QZ_EQUI

NUMERO FREQUENCE (HZ) **** AMORTISSEMENT **** NORME D’ERREUR

1 1.23915E+02 1.55604E-08 9.93686E-02

2 1.24546E+02 -3.76772E-10 1.81854E-01

3 4.97033E+02 -5.88927E-10 1.41043E-03

4 4.99575E+02 1.41315E-12 2.14401E-03

5 1.11531E+03 4.75075E-12 1.17669E-04

Example 2. Impressions in the message file of the eigenvalues retained during QEP. Excerpt from the sdll123a test case.

Now that the context of QEP in Code_Aster has been outlined, we will focus on power-type methods and their implementation in OPTION “PROCHE” and “AJUSTE”.