Eigensensitivity analysis for symmetric nonviscously damped systems with repeated eigenvalues

1. Introduction

The dynamic behavior of a discrete structural system is characterized by its eigenpairs, i.e., eigenvalues and corresponding eigenvectors. Parameter changes in the system will result in variations in its dynamic characteristics. The calculation of eigenpair derivatives is essential in understanding and determining the impact of parameter changes upon the system's dynamic behavior without reanalysis. Using the derivatives also facilitates the tracking of eigenpairs in the presence of parametric variations. The eigensolution sensitivities of structural and mechanical systems with respect to structural design parameters provide a guidance in structural optimal design [1, 2], finite model updating [3], structural damage detection [4] and system identification [5]. Therefore, eigenpair sensitivity has played a central part in structural stability analysis and has emerged as an important area of research.

Many efforts in both theoretical and computational aspects for standard and generalized eigenvalue problems have been made in response to the demand of engineering applications. Basically theoretical results on existence of derivatives of eigenvalues and eigenvectors of a matrix or matrix pencil depending on one parameter have been considered by Lancaster [6], Kato [7], Rellich [8], Meyer and Stewart [9]. The analytic dependence of a semisimple eigenvalue of an analytic matrix-valued function depending on one parameter is analyzed in [10]. The sensitivity of simple eigenvalues and corresponding eigenvectors of a matrix depending on several parameters is presented by Sun [11]. The sensitivity analysis of multiple eigenvalues and associated invariant subspaces of a matrix or matrix pencil depending on several parameters have been investigated by Sun [12, 13], Chu [14], Xie and Dai [15]. The existence of derivatives of eigenvalues and corresponding eigenvectors of an analytic matrix-valued function depending on several parameters is considered in [16]. Moreover, a numerical method for computing eigenpair derivatives is presented, and the method admits the use of various normalization strategies for eigenvectors and may be used for general nonlinear eigenvalue problem. In practical numerical computation, Fox and Kapoor [17] have presented exact expressions for the derivatives of eigenvalues and eigenvectors of symmetric matrix pencils with respect to any design variables. Many authors have extended Fox and Kapoor's approach to compute eigenpair derivatives for nonsymmetric matrices or nonsymmetric matrix pencils, see, e.g., [18-22]. Nelson [23] presented an efficient method to compute the first-order derivatives of simple eigenvalue and corresponding eigenvector for standard eigenvalue problems, where only eigenpair under consideration is required. This was extended to symmetric generalized eigenvalue problems with repeated eigenvalues having well separated derivatives in [24]. The requirement that the repeated eigenvalues must have well separated derivatives was relaxed in [25]. Recently, many authors developed numerical methods for computing derivatives of eigenvalues and eigenvectors of matrices or matrix pencils, see, e.g., [26-28]. Interested readers are referred to survey papers [13, 29, 30, 31]. Some of the existing works on sensitivity of eigensolutions of undamped systems are reviewed by Adhikari in [2, 5].

To apply above-described results to systems with general nonproportional (viscous) damping, the equations of motion has to be linearized to standard or generalized eigenvalue problem, and the size of the problem is doubled. Adhikari [32], Lee et al. [33], Adhikari and Friwell [34], Choi et al. [35], and Guedria et al. [36] presented some numerical methods for computing derivatives of simple eigenvalues and corresponding eigenvectors of the quadratic eigenvalue problem without use of linearization. Friwell and Adhikari [37] extended Nelson’s method to the quadratic eigenvalue problem. The sensitivity of multiple eigenpairs of the quadratic eigenvalue problem are considered in [38-40], however, the methods given there require that the repeated eigenvalues have well separated derivatives. Multiple eigenvalues with repeated first derivatives and well separated second derivatives are considered in [41]. Qian et al. [42] extended the methods of [25] for generalized eigenvalue problems to quadratic eigenvalue problems, and proposed numerical methods for computing first and higher order derivatives of multiple eigenpairs of quadratic eigenvalue problems. Li et al. [43] proposed a method for sensitivity analysis of eigenvalues and eigenvectors bsaed on a combined normalization. Wang and Dai [44] developed a new algotithm for computing the repeated eigenvalue and eigenvector derivatives of viscous damped systems.

The mentioned previous studies are only restricted to the case of undamped or viscously damped systems. However, a physically realistic model of damping mechanism within the scope of linear dynamic analysis (e.g., the damping in rubbers or composite materials) may not be a viscous damping model. Possibly the most general damping model is to take into account the hereditariness, hysteretic or frequency-dependent behavior, which often assume that the damping forces depend on the past history of motion via convolution integrals over some decaying kernel functions. Increasing the use of mechanical structures with nonviscous damping, such as composite structural materials, active control and damage tolerant systems in rocket, spacecraft, satellite, ships and automobiles has led to renewed demand for dynamic analysis of nonviscously damped systems. In this study, we consider a class of nonviscous damping models in which the damping forces depend on the past history of motion via convolution integrals over some kernel functions. The equations of motion describing free vibration of a $N$-degree-of freedom linear system with such damping can be expressed by:

where $\mathbf{M}$ and $\mathbf{K}\in {\mathbf{R}}^{N\times N}$ are symmetric mass and stiffness matrices, respectively, $\mathcal{G}\in {\mathbf{R}}^{\mathrm{N}\times \mathrm{N}}$ is the symmetric matrix of kernel functions. In the special case when $\mathcal{G}\left(t-\tau \right)=C\delta (t-\tau )$, Eq. (1) reduces to the case of viscously damped systems. Therefore the damping model of this kind is a further generalization of the familiar viscous damping.

During the past one decade, some studies have been presented for this nonviscous eigenproblem. Adhikari [45] extended the method of Fox and Kapoor to systems with nonviscous damping in the form of Eq. (1). In this method the eigenvector derivative was obtained as a linear combination of all of the eigenvectors. Adhikari and Friswell [46] extended Nelson’s method to symmetric and asymmetric nonviscously damped systems. However, Nelson's method is lengthy and clumsy for programming. Recently, Mirzaeifar et al. [47] proposed a new method based on a combination of algebraic and modal methods for generally asymmetric viscously damped systems. More recently, Li et al. [48] proposed a method of design sensitivity analysis of asymmetric viscously damped systems with distinct and repeated eigenvalues, which can compute the left and right eigenvector derivatives separately and independently. Li et al. [49-50] discussed the eigensensitivity of asymmetric viscously damped systems by two different normalization to eigevectors. Unfortunately, all these methods for eigensensitivity with nonviscously damped systems are only restricted to the case of distinct eigenvalues. In reference [51], Li et al. presented a method to calculate the eigensensitivity of general nonlinear eigenproblems with repeated eigenvalues by dividing the derivatives of eigenvectors into a particular solution and a homogeneous solution, but the process of finding particular solution is complicated.

In this paper, we present an improved algorithm for calculation of eigenpair derivatives for symmetric nonviscously damped systems with repeated eigenvalues based on the method in [51]. In the proposed method, the derivatives of eigenvectors are computed in $N$-space and the mode shape derivatives are divided into a particular solution and a homogeneous solution. A simplified method is given to calculate the particular solution by solving a linear equation with non-singularity coefficients, the homogeneous solution are computed by the second order derivative of eigenequation. The proposed method requires only the eigenvalues and eigenvectors under consideration, gives exact solution and guarantees numerical stability.

Throughout this paper we use the following notation. ${\mathrm{}\mathbf{C}}^{m\times n}$ denotes the set of complex matrices, ${\mathrm{}\mathbf{C}}^n={\mathbf{C}}^{n\times 1}$, ${\mathrm{}\mathbf{C}=\mathbf{C}}^1$. ${\mathbf{I}}_n$ is the identity matrix of order $n$, $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_1,\dots ,a_n)$ stands for the diagonal matrix with diagonal elements $a_1$, ... , $a_n$. ${\mathbf{A}}^T$ denotes the transpose of a matrix $\mathbf{A}$.

The remainder of this paper is arranged as follows. In Section 2, we presents theoretical background of the eigenvalue derivatives of nonviscously damped systems with repeated eigenvalues. In Section 3, we focus on dealing with the derivatives of eigenvectors, and derive the particular solutions to the governing equation of the eigenvector derivatives by solving a linear equation with non-singularity coefficients. Numerical stability of the proposed method was discussed in Section 4. In Section 5, one numerical example is performed by our proposed method. Finally, we make some concluding remarks in Section 6.

2. Theoretical background

The determination of eigenvalues and eigenvectors of nonviscously damped systems has been discussed by Adhikari [43]. Here we briefly outline the topics required for further developments. Taking the Laplace transform of Eq. (1) gives:

here $\stackrel{-}{\mathbf{u}}\left(s\right)=\mathbf{L}\left[\mathbf{u}\right(t\left)\right]\in C^N$, $\mathbf{G}\left(s\right)=\mathbf{L}\left[\mathcal{G}\right(t\left)\right]\in {\mathbf{C}}^{N\times N}$ and $\mathbf{L}[\cdot ]$ denotes the Laplace transform. The eigenvalues ${\lambda }_i$ are the roots of the characteristic equation:

whose order $m$ of the characteristic polynomial is generally higher than $\text{2}N$, that is, $m=2N+p$, $p\ne 0$. Thus, although the system has $N$ degrees of freedom, the number of eigenvalues is more than $\text{2}N$. This is a major difference between the nonviscously damped systems and the viscously damped systems where the number of eigenvalues is exactly $\text{2}N$, including any multiplicities. The eigenvector ${\mathbf{u}}_i$ corresponding to the eigenvalues ${\lambda }_i$ can be defined:

The normalized condition of ${\mathbf{u}}_i$ is:

where:

Suppose the system matrices in Eq. (1) are functions of some design parameter $p$, and at $p=p_0$, the eigenvalues problem has a semisimple eigenvalue ${\lambda }_0$ with multiplicity $n>1$. In this paper, we assume that the eigenvalue is usually multiple only at the isolated point $p=p_0$ and $\lambda \left(p\right)$ and $\mathbf{u}\left(p\right)$ are sufficiently differentiable at $p=p_0\text{.}$ Under these assumptions, the eigenvectors of Eq. (1) are uniquely determined for all $p\ne p_0$ in the neighborhood of $p_0$. We will develop a new method for computing the derivatives of $\lambda \left(p\right)$ and $\mathbf{u}\left(p\right)$. Denote the corresponding eigenvectors of ${\lambda }_0$ by $\mathbf{X}=(x_1,...,x_n)$, then:

where $\mathrm{\Lambda }={\lambda }_0{\mathbf{I}}_n$. To simplify notation, here and henceforth “$p_0$” is omitted for variables evaluated at $p=p_0$. Since the multiplicity is $n$, the orthonormal condition for the matrix $\mathbf{X}$ is as follows:

As $p$ varies, correspondingly, one will have the continuously varying diagonal matrix $\mathrm{\Lambda }\left(p\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\lambda }_1\right(p),...,{\lambda }_n(p\left)\right)$ and eigenvector subsets $\mathbf{Z}\left(p\right)=\left(z_1\right(p),...,z_n(p\left)\right)$ which satisfy:

The concern here is with the sensitivity of the above mentioned eigenpair, i.e., the derivatives of $\mathrm{\Lambda }\left(p\right)$ and $\mathbf{Z}\left(p\right)$ at $p=p_0$.

When an eigenvalue is repeated, calculation of the eigenpair derivatives becomes more complicated and has been an active subject of research. In this situation, the associated eigenvectors are only defined up to a subspace with dimension equal to the geometric multiplicity of the eigenvalue. Adjacent eigenvectors should be first determined from the given eigenvectors $\mathbf{X}$. The adjacent eigenvectors $\mathbf{Z}$ can be expressed in terms of $\mathbf{X}$ by an orthogonal transformation such as:

where $\mathbf{\Gamma }=({\gamma }_1,...,{\gamma }_n)$ is an orthogonal transformation matrix of dimension $n$×$n$ which need to be determined. It is natural that the adjacent eigenvectors satisfy the eigenequation and the orthonormal condition too:

The next procedure is to find $\mathbf{\Gamma }$ and then to find $\mathit{Z}$ and $\partial \mathbf{\Lambda }/\partial p$. Differentiating Eq. (9) with respect to the design parameter $p$ and letting $p\to p_0$ gives:

where:

Premultiplying each side of Eq. (13) by ${\mathbf{X}}^T$ and substituting $\mathbf{Z}=\mathbf{X}\mathbf{\Gamma }$ into it gives a new subeigenvalue problem, such as:

where:

This is a new subeigenvalue problem, it implies that $\partial {\lambda }_i/\partial p$ ($i=$ 1,…, $n$) are $n$ eigenvalues of matrix $\mathrm{\Xi }$ and column vectors ${\mathbf{\gamma }}_i$ ($i=$ 1,…, $n$) of $\mathbf{\Gamma }$ are specified eigenvectors corresponding to $\partial {\lambda }_i/\partial p$ respectively. Assuming that $\partial {\lambda }_i/\partial p$ are distinct, an unique $n$×$n$ orthogonal matrix $\mathbf{\Gamma }$ can be generated that can be used to define unique eigenvectors $\mathbf{Z}$ for the repeated eigenvalues by $\mathbf{Z}=\mathbf{X}\mathbf{\Gamma }$.

3. The proposed method

Due to the singularity of $D\left({\lambda }_0\right)$, the eigenvector derivatives $\partial Z/\partial p$ cannot be computed directly using Eq. (13). Using knowledge of linear algebra, the general solution of Eq. (13) has the form:

where $\mathbf{Z}$ is the homogeneous solution, $\mathbf{c}$ is a $n$×$n$ constant vector, ${\mathbf{V}}_1$ is a particular solution of Eq. (13), i.e.:

In view of Eq. (15), particular solution ${\mathbf{V}}_1$ may be required to be $\partial \mathbf{D}\left({\lambda }_0\right)/\partial s$-orthogonal with respect to $Z$, i.e.:

Using $\mathit{Z}=\mathit{X}\mathbf{\Gamma }$ and Eq. (17), we have:

Combining Eq. (16) and (18), we can write the following single matrix equation:

where the order of coefficient matrix on the left side of Eq. (19) is $(N+n)\times (N+n)$, we will prove the non-singularity of the coefficient matrix in next section. Eq. (19) can be split into:

Premultiply Eq. (20) by ${\mathbf{X}}^T$ and using Eq. (8), we get:

From which Eq. (20) is simplified to Eq. (16). This means that the part ${\mathbf{V}}_1$ of the solution to Eq. (19) is the particular solution to Eq. (13). Using non-singularity of the coefficient matrix of Eq. (19), we may solve Eq. (19) to achieve the particular solution ${\mathbf{V}}_1$ of Eq. (13). Once the particular solution ${\mathbf{V}}_1$ is obtained, the next procedure is to compute the constant matrix $\mathbf{c}$.

In order to compute the $n$×$n$ matrix $\mathbf{c}$, differentiating Eq. (12) with respect to $p$ yields:

Substituting Eq. (15) and (17) into (23), one gets:

where:

So the diagonal elements of $\mathbf{c}$ can be determined by:

To extract the off-diagonal elements, employing the second-order differential information of Eq. (11), one obtains:

Premultiplying each side of Eq. (26) by ${\mathbf{Z}}^T$ and use the relation Eqs. (14) (15), yields:

where:

So the off-diagonal elements of $\mathbf{c}$ should be:

Then, the first-order eigenpair sensitivity is completed when the derivatives of eigenvalues are distinct. The procedure of the proposed method is summarized as follows.

Algorithm.

Step 1. Input: The $N\times N$ coefficient matrices $\mathbf{M}$, $\mathbf{C}$, $\mathbf{K}$ and their derivatives, of order 1 to 2, the eigenvalue ${\lambda }_0$ and matrices of corresponding right and left eigenvectors $X\in C^{N\times n}$ satisfying normalizing conditions.

Step 2. Calculate $\mathrm{\Xi }=-{\mathbf{X}}^T\left({\lambda }_0^2\frac{\partial \mathbf{M}}{\partial p}+{\lambda }_0\frac{\partial \mathbf{G}\left({\lambda }_0\right)}{\partial p}+\frac{\partial \mathbf{K}}{\partial p}\right)\mathbf{X}.$

Step 3. Solving the subeigenvalue problem $\mathrm{\Xi }\mathbf{\Gamma }=\mathbf{\Gamma }\frac{\partial \mathbf{\Lambda }}{\partial p}$ and normalize $\mathbf{\Gamma }$ so that ${\mathbf{\Gamma }}^T\mathbf{\Gamma }={\mathit{I}}_n$.

Step 4. Compute the adjacent eigenvectors $\mathit{Z}$ by $\mathit{Z}=\mathbf{X}\mathbf{\Gamma }$.

Step 5. Define ${\mathbf{A}}^{*}=\left[\begin{array}{cc}\mathbf{D}\left({\lambda }_0\right)& \frac{\partial \mathbf{D}\left({\lambda }_0\right)}{\partial s}\mathbf{X}\\ {\mathbf{X}}^T\frac{\partial \mathbf{D}\left({\lambda }_0\right)}{\partial s}& 0\end{array}\right]$.

Step 6. $\left[\begin{array}{l}\mathbf{F}\\ 0\end{array}\right]=\left[\begin{array}{c}-\frac{\partial \mathbf{D}\left({\lambda }_0\right)}{\partial s}\mathbf{Z}\frac{\partial \mathbf{\Lambda }}{\partial p}-\left({\lambda }_0^2\frac{\partial \mathbf{M}}{\partial p}+{\lambda }_0\frac{\partial \mathbf{G}\left({\lambda }_0\right)}{\partial p}+\frac{\partial \mathbf{K}}{\partial p}\right)\mathbf{Z}\\ 0\end{array}\right]$.

Step 7. Solving the equation ${\mathbf{A}}^{\mathbf{*}}\left[\begin{array}{c}{\mathbf{V}}_1\\ \mathbf{\mu }\end{array}\right]=\left[\begin{array}{l}\mathbf{F}\\ 0\end{array}\right]$ get the particular solution ${\mathbf{V}}_1$ by $\left[\begin{array}{c}{\mathbf{V}}_1\\ \mathbf{\mu }\end{array}\right]=[{\mathbf{A}}^{*}{]}^{-1}\left[\begin{array}{l}\mathbf{F}\\ 0\end{array}\right]$.

Step 8. Compute the matrix $\mathbf{c}$ by Eqs. (25) and (27).

Step 9. Obtain the derivatives of eigenvectors by $\partial \mathbf{Z}/\partial p={\mathbf{V}}_1+\mathbf{Z}\mathbf{c}$.

Step 10. Output: $\partial \mathbf{\Lambda }/\partial p$, $\partial \mathbf{Z}/\partial p$.

4. Numerical stability of the proposed method

Identifying the non-singularity of the coefficient matrix ${\mathbf{A}}^{*}$ in Eq. (19) may be used to prove the numerical stability of the proposed method in the case of multiple eigenvalues:

where $\mathbf{X}$ is an $N\times n$ eigenvectors matrix, and ${\lambda }_0$ is the repeated eigenvalue of multiplicity $n$. To show that the coefficient matrix ${\mathbf{A}}^{*}$ is always non-singular, consider the solution of following matrix equation:

If Eq. (30) has only zero solution, we can conclude that the coefficient matrix ${\mathbf{A}}^{*}$ is non-singular. Splitting Eq. (30) as following forms:

Premultiply Eq. (32) by ${\mathbf{X}}^T$, we get:

Substituting ${\mathbf{Z}}_2=0_n$ into Eq. (31), there is:

i.e., ${\mathbf{Z}}_1$ is the eigenvectors matrix corresponding to ${\lambda }_0$. Because $\mathbf{X}$ is also the eigenvectors matrix, then we have:

where $\alpha$ is constant matrix. Substituting the above equation into Eq. (32):

Utilizing normalization Eq. (8), we can obtain $\mathbf{\alpha }=0_n$, then ${\mathbf{Z}}_1=\mathbf{X}\mathbf{\alpha }=0_{N\times n}$. Therefore the Eq. (30) has the unique zero solution, it can be concluded that ${\mathbf{A}}^{*}$ is a full rank matrix, then the numerical stability of the proposed method can be guaranteed.

5. Numerical example

In this section, a four degrees-of-freedom mass-spring system with nonviscous damping [51], shown in Fig. 1, will be used to illustrate the procedures of the proposed method. As it can be observed, two distinct cases are considered for the viscoelastic damping. On the one hand, the dissipative element is connected between the two masses. In this case, the viscoelastic damping force depends on the relative motion of the two masses. On the other hand, the two masses have also dissipative elements connecting them to the ground. In this case, the nonviscous damping force only depends on the absolute motion of the individual masses. The mass, stiffness and $\mathbf{G}\left(s\right)$ matrices are given by

Fig. 1A four degrees-of-freedom mass-spring system with nonviscous damping

Setting $m=$ 1 kg, $k=$1000 N/m and $k_1=$1000 N/m. Clearly, this model is a Boit model. The stiffness $k$ is chosen as the design parameter $p$. The derivatives of system matrices with respect to the design parameter $p$ are:

We consider the 2-repeated eigenvalues –1.4282×10^-3$+\text{44.7277}i\text{,}$ the corresponding normalized eigenvectors are:

Solving Eq. (14), the eigenvalue derivatives $\partial \mathbf{\Lambda }/\partial p$ and $\mathbf{\Gamma }$ can be obtained:

The eigenvector derivatives can be found by Eqs. (15), (25) and (27):

In order to verify the accuracy of eigenpair derivatives, we perturb $p$ by $\mathrm{\Delta }p=0.1$, and compute the linear approximation of perturbed system by Taylor series:

Comparing the approximated eigenpairs and the actual eigenpairs, the result are listed in Table 1 and Table 2.

Tables 1 and 2 show that the proposed method is efficient for the calculation of derivatives of semisimple eigenpairs of symmetric nonviscous quadratic eigenvalue problem.

6. Conclusions

In this paper, a new method for computation eigenvector derivatives with repeated eigenvalues for symmetric nonviscously damped eigensystems has been proposed. The derivatives of eigenpairs are computed in $N$ space and the mode shape derivatives are divided into a particular solution and a homogeneous solution. The method requires only the information of those eigenvectors corresponding to the repeated eigenvalues to extend the governing equations of particular solutions so that the coefficient matrix is non-singular and symmetric. In contrast to the methods in the reference, this algorithm of eigenpair derivatives is accurate, and only involves those degenerate eigenvectors to be differentiated, this is the highlight of the algorithm. The algorithm is simple, and it is easy to be implemented on computers. Numerical examples have demonstrated the validity of the proposed method. Finally, we should point out that there are unfinished tasks in this study. Among these, the case when the first derivatives of eigenvalues has repeated roots is especially interesting topics that remain to be further investigated.