- Email: [email protected]
A highly stable and accurate computational method for eigensolutions in structural dynamics
Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059 www.elsevier.com/locate/cma
A highly stable and accurate computational method for eigensolutions in structural dynamics Zhaohui Qi a, D. Kennedy a
b,*
, F.W. Williams
c
Department of Engineering Mechanics, Dalian University of Technology, Dalian 116023, PR China b Cardiff School of Engineering, Cardiff University, Cardiff CF24 3AA, Wales, United Kingdom c Department of Building and Construction, City University of Hong Kong, Kowloon, Hong Kong Received 18 October 2002; accepted 26 July 2005
Abstract A new computational method for the linear eigensolution of structural dynamics is proposed. The eigenvalue problem is theoretically transformed into a specific initial value problem of an ordinary differential equation. Based on the physical meaning of the sign count of the dynamic stiffness matrix, a stability control device is designed and combined with the fourth-order Runge–Kutta method. The resulting method finds the eigenvalues and eigenvectors at the same time, with high accuracy and high stability. Numerical examples show that the proposed method still gives high accuracy solutions when there is a great difference in magnitude among the eigenvalues, and also when some eigenvalues are very close to each other. 2005 Elsevier B.V. All rights reserved. Keywords: Structural dynamics; Eigenvalues; Eigenvectors; Computational methods
1. Introduction Since the natural frequencies and modes of vibration take an essential role in structural dynamics, it is a very important aim to have reliable and accurate methods for solving the eigenvalue problem of a structure. However, the computation of structural eigenvalues and eigenvectors remains a difficult task even now. Difficulties arise from both the discretization of the structure (e.g., by the finite element method) and the procedures of numerical analysis. When the higher natural frequencies of a structure are desired, smaller elements should be used [1–4], resulting not only in an increased order of the matrices representing the structure but also in an increased propensity of these matrices for ill-conditioning. Therefore, highly stable and high accuracy computational methods for eigensolutions of structures are necessary to obtain the eigenpairs with good precision. Among existing numerical methods for the eigensolution, the subspace iteration method has been widely used [5–11]. In spite of its many capabilities, sometimes this method is sensitive to the initial trial values of eigenvalues and eigenvectors and some eigenpairs of interest may be missed as a result of numerical instability. As an alternative to finding the eigenvalues directly, a few methods [12,13], such as the Wittrick–Williams algorithm [13– 17], try to find bounds on the eigenvalues by counting the number of eigenvalues exceeded at each of a sequence of trial frequencies. The procedure to locate the bounds employs only the signs of the diagonal elements of the upper triangular matrix resulting from Gauss elimination of the dynamic stiffness matrix. As a result, the method is so stable that it has been regarded as infallible [13]. However, if the sequence of trial frequencies is chosen by the bisection method, only linear *
Corresponding author. Tel.: +44 29 2087 5340; fax: +44 29 2087 4939. E-mail address: [email protected] (D. Kennedy).
0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.08.007
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4051
convergence on the eigenvalues can be achieved, and the solution of eigenvectors does not benefit from the method. Although efforts to improve the convergence and the computation of eigenvectors have been made by some authors [16–18], further work is needed in this area. Such root-counting methods have generally been applied to the transcendental eigenvalue problems arising from analytical solutions of the governing differential equations of the structural members. The object of this paper is to exploit the high stability of root-counting methods by applying them to the conventional linear eigenvalue problem, in order to obtain eigenvalues and eigenvectors with high accuracy. 2. Sign count of dynamic stiffness matrix The linear eigenvalue problem in structural dynamics can be stated mathematically as attempting to find a positive real parameter k in the equation KD ðkÞ ðK kMÞx ¼ 0;
ð1Þ
so as to make the n-dimensional vector x non-trivial, where the n · n dimensional real, symmetric, non-negative definite matrices K and M are the static stiffness and mass matrices, respectively. It is well known that there exists a non-singular matrix X whose ith column is the eigenvector xi associated with the eigenvalue ki, which satisfies ¼ Diagðk1 k; . . . ; ki k; . . . ; kn kÞ. ð2Þ Xt KD ðkÞX By virtue of SylvesterÕs law of inertia [9], if a series of non-singular matrices P1, . . . , Pm is found such that kÞPi ¼ Di ði ¼ 1; . . . ; mÞ; Pt KD ð i
ð3Þ
where each of the Di (i = 1, . . . , m) is a diagonal matrix, then each Di has the same number of negative elements, called the sign count of the dynamic stiffness matrix KD ð kÞ and denoted sfKD ðkÞg. Moreover, Eq. (2) shows that sfKD ðkÞg equals the number of eigenvalues exceeded by k. Knowing sfKD ðkÞg enables bounds on the eigenvalues to be obtained. For example, a trial frequency kk1 is a lower bound on the kth eigenvalue kk if sfKD ð kk1 Þg < k, or an upper bound if sfKD ðkk1 Þg P k. Once lower and upper bounds have been found, convergence on kk is achieved by successively evaluating the sign count at trial frequencies lying between the two bounds. Since only the signs, and not the values, of the elements of the corresponding diagonal matrix are needed when evaluating sfKD ð kÞg, the sign count is quite insensitive to numerical errors. It is this feature of the sign count that makes the algorithms based on it almost infallible. 3. Reduction of dynamic stiffness matrix to diagonal A generalized inner product of two vectors pi and pj of order n is introduced as kÞp . ðp ; p Þ ¼ pt KD ð i
j
ð4Þ
j
i
ð1Þ ð1Þ p1 ; p2 ; . . . ; pð1Þ n
Any set of independent vectors of order n can be transformed to an orthogonal set p1, p2, . . . , pn, in the spirit of the method of modified Gram–Schmidt orthogonalization, as follows. The first desired vector is simply chosen ð1Þ as p1 ¼ p1 . Then, for k = 2, . . . , n, after performing the transformations ðk1Þ pk1 ; pi ðkÞ ðk1Þ p ði ¼ k; . . . ; nÞ; ð5Þ pi ¼ pi ðpk1 ; pk1 Þ k1 ðkÞ
ðkÞ
ðkÞ
ðkÞ
ðkÞ
ðk1Þ
ðk1Þ
is exchanged with pk the kth vector is chosen as pk ¼ pk . If ðpk ; pk Þ ¼ 0, but ðpk0 ; pk0 Þ 6¼ 0 for some k 0 > k, then pk0 to make the choice valid. ð1Þ ð1Þ It can be seen from Eq. (5) that pk is a linear combination of pi (i = 1, . . . , k 1). If each pi is selected as the ith column of the unitary matrix of order n, then the components of pk have the feature that pkk = 1 and pkj = 0 for j > k. The matrix P whose kth column is pk is an upper triangular matrix whose diagonal elements are all unity, and satisfies Pt KD ð kÞP ¼ D;
ð6Þ
where D is a diagonal matrix. Being real and symmetric, the dynamic stiffness matrix can be uniquely decomposed into the form KD ð kÞ ¼ Ld DLt ; ð7Þ d
where Ld is a lower triangular matrix. Comparison between Eqs. (6) and (7) gives
4052
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
f12
f13
2
f23
f14
1 3
2
f24
1
4
3
f11
4
f22
3
f34
1
4 1
2
4
f33
2
3
f44
Fig. 1. Modes of displacements and forces for reduction of dynamic stiffness matrix to diagonal form.
P ¼ Lt d ;
ð8Þ
so that the matrix of congruence transformation can be obtained by the usual Gauss elimination process. If pi is regarded as a mode of general displacements of the structure whose stiffness matrix is KD ðkÞ, then kÞp f i ¼ KD ð i
ð9Þ
can be viewed as a required mode of forces to produce the displacement pi. Defined as the matrix whose kth column is fk, the force matrix F ¼ Pt D
ð10Þ
is a lower triangular matrix. Therefore, if j < i, the jth component of the ith mode of forces fij = 0, while the corresponding component of the mode of displacements pij 5 0. From a physical point of view, pi is generated by clamping the (i + 1)th, . . . , nth degrees of freedom of the structure and applying a force fii at the ith degree of freedom to make pii = 1; the forces fij (j = i + 1, . . . , n) are simply the corresponding constraint forces, as illustrated by Fig. 1. If fii = 0 at a specific trial frequency k, the next mode of displacements is evaluated by replacing k with a new trial fre quency. When k is an eigenvalue of the structure, it is not necessary to evaluate the remaining modes of displacements because pi is the eigenvector being sought. In the resulting diagonal matrix D, the kth diagonal element dk, being defined as the kth energy norm, is the input energy necessary to generate the mode of displacements pi, and is related to the forces through d k ¼ fkk .
ð11Þ
4. Qualitative behaviour of the energy norms Being expressed as d i ¼ pt KD ð kÞp ði ¼ 1; . . . ; nÞ; i
i
ð12Þ
the energy norms are definite functions of the trial frequency k, if the order of freedoms in the structure is fixed. The derivatives of these functions can be expressed in terms of the modes of displacements by the equations dd i dpti dKD ðkÞ ð13Þ pi ði ¼ 1; . . . ; nÞ. ¼ 2 f i þ pti dk dk dk The first term on the right-hand side of Eq. (13) is zero, because the first i 1 elements of fi and the last n i elements of pi are zero, and pii is always equal to unity. Thus, according to the definition of the dynamic stiffness matrix, Eq. (13) can be re-written as dd i ¼ pti Mpi dk
ði ¼ 1; . . . ; nÞ.
ð14Þ
The right-hand side of Eq. (14) must be negative because M is a positive definite matrix. Therefore, di (i = 1, . . . , n) are monotonic decreasing functions of the trial frequency k. Let KkD , Pk and Dk be the principal submatrices obtained by deleting the (k + 1)th, . . . , nth rows and columns from the matrices KD, P and D, respectively. These matrices are related to each other by Ptk KkD ð kÞPk ¼ Dk .
ð15Þ
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4053
Fig. 2. Energy norms as a function of trial frequency.
The kth energy norm can be expressed as D k ð kÞ ; dk ¼ Dk1 ð kÞ
ð16Þ
where Dk ð kÞ, the determinant of KkD ð kÞ, is equal to the determinant of Dk. From the analysis of Section 2, the number of negative elements in Dk is equal to the number of the roots of the equation D k ð kÞ ¼ 0 that are exceeded by k, because the matrix KkD can be regarded as the dynamic stiffness matrix of the structure with the (k + 1)th to nth degrees of freedom fixed. If Dk1 ð kÞ ¼ 0 while Dk ð kÞ 6¼ 0, there must be a small positive real number e such that the sign count of Dk1(k) increases by one while the sign count of Dk(k) is invariant as k increases from k e to k þ e, so that dk must jump from 1 to +1 at k. In general, the curve of the kth energy norm dk has a number of branches separated at the trial frequencies k where Dk1 ð kÞ ¼ 0 while Dk ð kÞ 6¼ 0, and in each of the branches dk decreases monotonically from +1 to 1, as illustrated in Fig. 2. 5. Criterion for eigenvalues in terms of energy norms From the analysis of Section 3, the last energy norm dn is the force required at the nth degree of freedom of the structure in order to produce the mode of displacements pn, the last column of the transformation matrix P. As no forces act at the other degrees of freedom, dn is the only possible non-zero force necessary to produce pn. Therefore, if kÞ ¼ 0; ð17Þ d n ð k must be an eigenvalue of the structure and pn must be the corresponding eigenvector. This criterion for eigenvalues applies to most practical structures except in the rare case where there exists an eigenvalue such that the last component of the corresponding eigenvector happens to be zero. This difficulty can be overcome by changing the order of the degrees of freedom to make Eq. (17) a valid criterion. Suppose a proper lower bound kli and
4054
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
a proper upper bound kui on the ith eigenvalue have been found. Then, by definition, the corresponding sign counts are sfKD ðkli Þg ¼ i 1 and sfKD ðkui Þg ¼ i, respectively, and the criterion expressed by Eq. (17) is valid if the conditions d n ðkli Þ > 0;
d n ðkui Þ < 0
ð18Þ
are satisfied simultaneously. The physical interpretation of the inverse of the dynamic stiffness matrix is well known; i.e., its jth column is the mode of displacements generated by imposing a unit force at the jth degree of freedom while no forces are imposed at the other degrees of freedom. Therefore, when the former jth degree of freedom is re-numbered as the last degree of freedom, the corresponding last energy norm d ðjÞ n can be obtained as d ðjÞ n ¼
1 ; ðK 1 D Þjj
ð19Þ
where ðK 1 D Þjj is the jth diagonal element of the inverse of matrix KD ðkÞ. If a matrix P has been found such that t P KD ðkÞP ¼ D, then the inverse of KD ðkÞ can be obtained easily from 1 t K1 D ðkÞ ¼ PD P .
ð20Þ
Combination of Eqs. (19) and (20) gives 1 d ðjÞ ; ð21Þ n ¼ X 2 n pjk k¼j d k where pij is the element of P at the jth column and the ith row. The components of the corresponding mode of displacements pðjÞ n are expressed as Xn pik pjk k¼minðj;iÞ dk ðjÞ pin ¼ ði ¼ 1; . . . ; nÞ. ð22Þ Xn p2jk k¼j d k Based on the above analysis, when Eq. (17) is employed as the criterion for eigenvalues, one can determine which degree of freedom should be the last one as follows: (1) When a proper lower or upper bound of the required eigenvalue is known, calculate all the d ðjÞ n (j = 1, . . . , n). (2) Check if Eq. (18) is satisfied. (3) If none of the d ðjÞ n (j = 1, . . . , n) satisfy the requirement, then set the trial frequency to the average of the lower and upper bounds and go back to step (1). (4) Any one of the freedoms j that satisfies Eq. (18) can be numbered last.
6. Computational method for the eigensolution For simplification of the notation, and without loss of generality, it is assumed that Eq. (17) is the criterion for all the eigenvalues. Then each eigenvalue corresponds to the zero point of a specific branch of dn. Within the interval ðkli ; kui Þ, dn is a monotonic decreasing function of k, which means that the map from k to dn is one to one. Therefore, k can be viewed mathematically as a function of dn, and the derivative of this function with respect to k is obtained by inverting the righthand side of Eq. (14), i.e., dk 1 ¼ t . ð23Þ dd i pi Mpi In order to find the ith eigenvalue ki, a trial value k0 is chosen within the interval ðkli ; kui Þ and the corresponding last energy 0 norm d n is evaluated. The point with coordinates ðk0 ; d 0n Þ serves as an initial point on the curve generated by Eq. (23). The eigenvalue ki can be found by locating the point (ki, 0). In other words, finding an eigenvalue amounts to solving an initial value problem of the ordinary differential Eq. (23). This transformation is of great value because computational methods for the eigensolution can be derived from a variety of highly developed numerical methods for the solution of ordinary differential equations. Denoting the right-hand side of Eq. (23) as f ð kÞ, the one-step Euler method gives t t p ðK k0 MÞpn p Kp k¼ k0 þ hf ð k0 Þ ¼ k0 þ n t ¼ tn n ; ð24Þ pn Mpn pn Mpn where the step of integration
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
h ¼ d 0n .
4055
ð25Þ
Eq. (24) shows that the method of RayleighÕs quotient can be regarded as the result of the application of the one-step Euler method in solving Eq. (23). Because of its high accuracy and simplicity, the fourth-order Runge–Kutta method is preferred, so that the eigenvalue is estimated as follows: k¼ k0 þ 16f ð k0 Þ þ 13f ð k1 Þ þ 13f ð k2 Þ þ 16f ð k3 Þ h; ð26Þ where k1 ¼ k0 þ 12f ð k0 Þh;
h ¼ d 0n ;
k2 ¼ k0 þ 12f ðk1 Þh;
k3 ¼ k0 þ f ðk2 Þh.
ð27Þ
If the corresponding last energy norm dn(k) is larger than a given tolerance, a new estimate of the eigenvalue is obtained from Eq. (26) with k0 ¼ k and h = dn(k). This process is executed repeatedly until dn(k) is small enough. In the meantime, it is also possible to improve the estimates of the bounds on the eigenvalues, since sfKD ðkÞg can be easily obtained during the evaluation of f ð kÞ. If sfKD ð kÞg ¼ k then k must be an upper bound on the kth eigenvalue kk as well as a lower bound on the (k + 1)th eigenvalue kk+1. If k > klkþ1 then k is a better lower bound on kk+1, and if k < kuk then k is a better upper bound on kk. It is obvious that the derivatives employed in Eq. (26) should be obtained from the same curve. This requires that, when the ith eigenvalue is being sought, sfKD ð kk Þg should be equal to either i or i 1 for each kk (k = 0, 1, 2, 3). Also it is necessary that d n ðkk Þ < 0 when sfKD ðkk Þg ¼ i and d n ðkk Þ > 0 when sfKD ðkk Þg ¼ i 1. In practice, these requirements are usually met as long as the starting trial frequency k0 lies between the proper bounds. If these requirements are not met, a new trial frequency l
u
k þ ki k¼ i 2
ð28Þ
is chosen. The replacement of Eq. (26) by Eq. (28) can only reduce the efficiency of the solution, but can never cause it to fail. Two factors are known to affect the efficiency of the proposed method for eigensolution, as shown in Fig. 3. The first factor, denoted by b, is the local bandwidth of the curve being employed, which is defined as the distance between the two boundary points at which the values of energy norm become infinite. Another factor, denoted by e, is the eccentricity of the curve, which is defined as the least distance between the zero point and the boundary points of the curve, expressed as a fraction of the bandwidth. The larger b and e are, the less effort is needed to determine which degree of freedom is suitable for selection as the last degree of freedom, and the faster the convergence should be. While finding the suitable branch of the curve, as described in Section 5, there may be more than one degree of freedom that satisfies Eq. (18). In this case, the one that minimizes d n ðkli Þ d n ðkui Þ is chosen. Although such a choice cannot guarantee that the branch of the selected curve has the largest b and e, it can prevent the most badly behaved curves from being selected.
Fig. 3. Definition of the factors b and e.
4056
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
7. Simplification of the structure The simplest structure for analysis is a chain structure. Fortunately, any complex structure can be transformed into an equivalent chain structure by applying the Householder transformation [7]. For example, a structure having n nodes with nd degrees of freedom per node can be simplified into a chain structure having n · nd nodes with one degree of freedom per node, as illustrated in Fig. 4. Being symmetric and positive definite, the mass matrix M of the structure can be decomposed into the form M ¼ Lm Ltm ;
ð29Þ
where Lm is a lower triangular matrix. Therefore, the eigenvalue problem of Eq. (1) can be re-written as t t Lm ðL1 m KLm kIÞLm u ¼ 0;
ð30Þ
where I is a unit matrix with the same order as that of the mass matrix. Since any real symmetric matrix can be transformed into a tri-diagonal matrix by the Householder transformation, Eq. (30) can be further simplified to ðK0 kIÞu0 ¼ 0;
ð31Þ
where 2
u0 ¼ Pth Ltm u;
a1
6b 6 1 6 6 6 6 0 t 1 t K ¼ Ph ðLm KLm ÞPh ¼ 6 6 6 6 6 6 4
3
b1 a2
b2
b2
a3
..
.
..
.
..
.
an1 bn1
7 7 7 7 7 7 7. 7 7 7 7 7 bn1 5
ð32Þ
an
t The orthogonal matrix Ph is the matrix of the Householder transformation for L1 m KLm . A great reduction in the computational work can be achieved by dealing with the simplified structure instead of the original structure. Because the stiffness matrix of the simplified structure is a tri-diagonal matrix and its mass matrix is a unit matrix, the corresponding energy norms are obtained as follows:
d 1 ¼ a1 ;
d iþ1 ¼ aiþ1
b2i di
ði ¼ 1; 2; . . . ; n 1Þ;
ð33Þ
and the corresponding modes of displacements are obtained by ð1Þ
p1 ¼ p 1 ;
ð1Þ
piþ1 ¼ piþ1
bi p di i
ði ¼ 1; 2; . . . ; n 1Þ;
ð34Þ
ð1Þ
where pi is the ith column of the unit matrix. Moreover, for a chain structure, when the jth degree of freedom is chosen as the last degree of freedom, in order to obtain d ðjÞ n it is not necessary to calculate the energy norms and the corresponding modes of displacements in the original order. Instead, a simpler algorithm is as follows:
Structure A having 3 nodes with 2 degrees of freedom per node
Chain B having 6 nodes with 1 degree of freedom per node
Fig. 4. Simplification of a structure A into a chain B.
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4057
b2i1 ði ¼ n; n 1; . . . ; j þ 1Þ; di b2 ¼ aiþ1 i ði ¼ 1; 2; . . . ; jÞ di
d n ¼ an ;
d i1 ¼ ai1
ð35Þ
d 1 ¼ a1 ;
d iþ1
ð36Þ
and pn ¼ pnð1Þ ; ð1Þ
p1 ¼ p 1 ;
bi1 p ði ¼ n; n 1; . . . ; j þ 1Þ; di i b ð1Þ ¼ piþ1 i pi ði ¼ 1; 2; . . . ; jÞ. di ð1Þ
pi1 ¼ pi1
ð37Þ
piþ1
ð38Þ
By comparing the computational work necessary for a general structure with that for a chain structure, it is seen that such simplification of the structure is highly recommended. 8. Numerical examples The following examples were solved using a computer program incorporating the proposed method. The simple Examples 1–4 demonstrate the numerical stability of the method in extreme cases that are prone to ill-conditioning. The larger Example 5 demonstrates the potential for applying the method to problems of industrial scale. Example 1. Considering a structure with dynamic stiffness matrix 1 0 2 1 KD ðkÞ ¼ k ; 1 2 0 1
ð39Þ
it can be seen that the two diagonal elements in KD(2) are both zero. However, taking a new trial frequency k0 ¼ 2 1018 , and reducing KD ð k0 Þ to a diagonal matrix gives d1 = 1018 and d2 = 1018. Since the eigenvalues of this structure are 1 and 3, the sign count of KD ð k0 Þ should be equal to 1, which is implied by the values of d1 and d2. Example 2. In the case of KD ð kÞ ¼ Diagð5; 3; 2; 6; 7Þ k Diagð1; 1; 1; 1; 1Þ, the eigenvalues do not always make the last diagonal element of the dynamic stiffness matrix equal to zero. This problem was solved to ensure that the proposed method works properly in this case. From the results shown in Table 1, it is confirmed that the proposed method can automatically recognise which freedom should be the last one and gives the correct answer. Example 3. Let the mass matrix be a unit matrix and the stiffness matrix K = RtVR, where V = Diag(0.001, 20, 500, 3, 7.1 · 108) and R is an orthogonal matrix. Being equal to the diagonal elements of the matrix V, the eigenvalues for this structure have a great variation in magnitude, which indicates that the dynamic stiffness matrix has a propensity for illconditioning. The results obtained by the proposed method are shown in Table 2, where e1 is the error of the last energy norm, e2 is the maximal component in the force vector KD(ki)vi, and ki and vi are approximated eigenvalues and eigenvectors, respectively. The results are very accurate even for such an ill-conditioned problem. Moreover, the eigenvalues and eigenvectors were obtained at the same time, with high accuracy. The physical meanings of e1 and e2 imply that they should be of roughly the same magnitude, which is indeed revealed by the results. Example 4. Consider a structure whose mass matrix M = RtDiag(1, 2, 3, 4, 5)R and stiffness matrix K = Rt Diag(50, 2 · 49.9, 3 · 50.1, 4 · 49.99, 5 · 50.01)R, where R is an orthogonal matrix. The eigenvalues of this structure are all very close to 50. The numerical solution of this problem is shown in Table 3, where the very small magnitudes of e1 and e2, defined as in Example 3, show that both the eigenvalues and eigenvectors were correctly obtained.
Table 1 Results for Example 2, showing the eigenvalues ki and eigenvectors vi i
ki
vi
1 2 3 4 5
2.00000000 3.00000000 5.00000000 6.00000000 7.00000000
0.00000000 0.00000000 1.00000000 0.00000000 0.00000000
0.00000000 1.00000000 0.00000000 0.00000000 0.00000000
1.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.00000000 1.00000000 0.00000000
0.00000000 0.00000000 0.00000000 0.00000000 1.00000000
4058
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
Table 2 Results for Example 3, showing the eigenvalues ki, eigenvectors vi and measures of mode accuracy e1 and e2 i 1 2 3 4 5
ki
e1
vi
0.001 3 20 500 7.1 · 108
0.6208 0.6095 0.3490 0.1210 0.3267
0.2525 0.7172 0.3796 0.1454 0.5066
0.4128 0.0293 0.0498 0.7578 0.5020
0.5490 0.3366 0.1780 0.5957 0.4457
e2 8
1.125 · 10 2.850 · 109 4.766 · 109 7.449 · 108 5.960 · 108
0.2812 0.0014 0.8366 0.1873 0.4312
5.215 · 108 3.210 · 108 1.490 · 108 5.588 · 108 6.706 · 108
Table 3 Results for Example 4, showing the eigenvalues ki, eigenvectors vi and measures of mode accuracy e1 and e2 i
ki
vi
1 2 3 4 5
49.90 49.99 50.00 50.01 50.10
0.1934 0.8505 0.0930 0.4716 0.0910
0.1450 0.2677 0.5804 0.4162 0.6303
0.1574 0.2141 0.3001 0.5255 0.7505
1
0.1443 0.3973 0.7513 0.4753 0.1761
3
0.9466 0.0366 0.0052 0.3199 0.0172
e1
e2
8.700 · 1015 4.051 · 1015 1.323 · 1014 3.530 · 1014 4.783 · 1014
8.743 · 1015 1.904 · 1015 1.141 · 1014 1.562 · 1014 3.204 · 1014
7
5
3m
4
2
8
6 4m
4m
4m
Fig. 5. A frame composed of beams.
Table 4 Comparison between the proposed method and the inverse power method for Example 5, showing the eigenfrequencies ki (in Hz), the number of iterations n and measure of mode accuracy e2 i
Proposed method
Inverse power method
ki
e2
n
ki
e2
n
1 2 3 4 5 6
15.2828042 38.7297544 44.4735753 48.3991123 50.5691556 57.7034310
0.002413 0.123841 0.001349 0.082575 5.332224 0.001154
4 4 4 4 4 4
15.2837170 38.7990895 44.5553287 48.4878425 50.7001757 57.8424109
0.745206 3.663227 4.991209 5.618068 4.694577 9.630830
7 31 72 89 45 62
Example 5. The lowest six eigenfrequencies and corresponding modes of vibration of the frame shown in Fig. 5 were found. The frame is fully built-in at node 1 and horizontal movement is prevented at node 2. All the beams in the frame have identical cross-sections and are made of the same material, with flexural rigidity EI = 5 MN m2, extensional rigidity EA = 900 MN and mass per unit length l = 35 kg m1. Each beam is approximated by Bernoulli elements of length 0.2 m, so that the frame has 250 nodes and provides a demonstration of the method on a realistically sized problem. The results are shown in Table 4, where e2 is a measure of mode accuracy as defined in Example 3. It is seen that the proposed method is more accurate and more efficient than the commonly used inverse power method.
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4059
9. Conclusions A new computational method for the eigensolution of structural dynamics has been proposed. In contrast to most existing methods, the proposed method theoretically transforms the eigenvalue problem into a specific initial value problem of an ordinary differential equation. This transformation enables many advanced numerical methods for ordinary differential equations to be used in the computation of the eigensolution. Based on the physical meaning of the sign count of the dynamic stiffness matrix, a stability control device has been designed and combined with the fourth-order Runge–Kutta method. The resulting method finds the eigenvalues and eigenvectors at the same time, with high accuracy and high stability. Numerical examples have shown that even in very unfavourable conditions, e.g., when there is great variation in magnitude among the eigenvalues or some eigenvalues are very close to each other, the proposed method still gives high accuracy solutions. Acknowledgements This work is supported by the China Association for International Exchange of Personnel. The first author is especially grateful for the financial support of Cardiff University and the National Science Foundation of China. The third author has returned to his chair at Cardiff University upon completion of his appointment at City University of Hong Kong. References [1] P. Ladeveze, J.P. Pelle, Accuracy in finite element computation for eigenfrequencies, Int. J. Numer. Methods Engrg. 28 (1989) 1929–1949. [2] N.-E. Wiberg, R. Bausys, P. Hager, Adaptive h-version eigenfrequency analysis, Comput. Struct. 71 (1999) 565–584. [3] P. Hager, N.-E. Wiberg, Adaptive eigenfrequency analysis by superconvergent patch recovery, Comput. Methods Appl. Mech. Engrg. 176 (1999) 441–462. [4] M. Duran, M. Miguez, J.-C. Nedelec, Numerical stability in calculation of eigenfrequencies using integral equations, J. Comput. Appl. Math. 130 (2001) 323–336. [5] K.-J. Bathe, E.L. Wilson, Large eigenvalue problems in dynamic analysis, J. Engrg. Mech. Div. ASCE 98 (1972) 1471–1485. [6] K.-J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, 1996. [7] L. Meirovitch, Computational Methods in Structural Dynamics, Sijthoff and Noordhoff, Netherlands, 1980. [8] A.R. Collar, A. Simpson, Matrices and Engineering Dynamics, Ellis Horwood, Chichester, 1987. [9] L. Mirsky, An Introduction to Linear Algebra, Constable, London, 1982. [10] K.-J. Bathe, S. Ramaswamy, An accelerated subspace iteration method, Comput. Methods Appl. Mech. Engrg. 23 (1980) 313–331. [11] S. Rajendran, M.V. Narasimhan, An accelerated subspace iteration method, Int. J. Numer. Meth. Engrg. 37 (1994) 141–153. [12] K.K. Gupta, Eigenproblem solution by a combined Sturm sequence and inverse iteration technique, Int. J. Numer. Methods Engrg. 7 (1973) 17–42. [13] W.H. Wittrick, F.W. Williams, A general algorithm for computing natural frequencies of elastic structures, Q. J. Mech. Appl. Math. 24 (1971) 263– 284. [14] F.W. Williams, M.S. Anderson, Inclusion of elastically connected members in exact buckling and frequency calculations, Comput. Struct. 22 (1986) 395–397. [15] F.W. Williams, W.H. Wittrick, Exact buckling and frequency calculations surveyed, J. Struct. Engrg. ASCE 109 (1983) 169–187. [16] F.W. Williams, D. Kennedy, Reliable use of determinants to solve non-linear structural eigenvalue problems efficiently, Int. J. Numer. Meth. Engrg. 26 (1988) 1825–1841. [17] A. Simpson, On the solution of S(x)x = 0 by a Newtonian procedure, J. Sound. Vib. 97 (1984) 153–164. [18] C.T. Hopper, F.W. Williams, Mode finding in nonlinear structural eigenvalue calculations, J. Struct. Mech. 5 (1977) 255–278.
A highly stable and accurate computational method for eigensolutions in structural dynamics Zhaohui Qi a, D. Kennedy a
b,*
, F.W. Williams
c
Department of Engineering Mechanics, Dalian University of Technology, Dalian 116023, PR China b Cardiff School of Engineering, Cardiff University, Cardiff CF24 3AA, Wales, United Kingdom c Department of Building and Construction, City University of Hong Kong, Kowloon, Hong Kong Received 18 October 2002; accepted 26 July 2005
Abstract A new computational method for the linear eigensolution of structural dynamics is proposed. The eigenvalue problem is theoretically transformed into a specific initial value problem of an ordinary differential equation. Based on the physical meaning of the sign count of the dynamic stiffness matrix, a stability control device is designed and combined with the fourth-order Runge–Kutta method. The resulting method finds the eigenvalues and eigenvectors at the same time, with high accuracy and high stability. Numerical examples show that the proposed method still gives high accuracy solutions when there is a great difference in magnitude among the eigenvalues, and also when some eigenvalues are very close to each other. 2005 Elsevier B.V. All rights reserved. Keywords: Structural dynamics; Eigenvalues; Eigenvectors; Computational methods
1. Introduction Since the natural frequencies and modes of vibration take an essential role in structural dynamics, it is a very important aim to have reliable and accurate methods for solving the eigenvalue problem of a structure. However, the computation of structural eigenvalues and eigenvectors remains a difficult task even now. Difficulties arise from both the discretization of the structure (e.g., by the finite element method) and the procedures of numerical analysis. When the higher natural frequencies of a structure are desired, smaller elements should be used [1–4], resulting not only in an increased order of the matrices representing the structure but also in an increased propensity of these matrices for ill-conditioning. Therefore, highly stable and high accuracy computational methods for eigensolutions of structures are necessary to obtain the eigenpairs with good precision. Among existing numerical methods for the eigensolution, the subspace iteration method has been widely used [5–11]. In spite of its many capabilities, sometimes this method is sensitive to the initial trial values of eigenvalues and eigenvectors and some eigenpairs of interest may be missed as a result of numerical instability. As an alternative to finding the eigenvalues directly, a few methods [12,13], such as the Wittrick–Williams algorithm [13– 17], try to find bounds on the eigenvalues by counting the number of eigenvalues exceeded at each of a sequence of trial frequencies. The procedure to locate the bounds employs only the signs of the diagonal elements of the upper triangular matrix resulting from Gauss elimination of the dynamic stiffness matrix. As a result, the method is so stable that it has been regarded as infallible [13]. However, if the sequence of trial frequencies is chosen by the bisection method, only linear *
Corresponding author. Tel.: +44 29 2087 5340; fax: +44 29 2087 4939. E-mail address: [email protected] (D. Kennedy).
0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.08.007
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4051
convergence on the eigenvalues can be achieved, and the solution of eigenvectors does not benefit from the method. Although efforts to improve the convergence and the computation of eigenvectors have been made by some authors [16–18], further work is needed in this area. Such root-counting methods have generally been applied to the transcendental eigenvalue problems arising from analytical solutions of the governing differential equations of the structural members. The object of this paper is to exploit the high stability of root-counting methods by applying them to the conventional linear eigenvalue problem, in order to obtain eigenvalues and eigenvectors with high accuracy. 2. Sign count of dynamic stiffness matrix The linear eigenvalue problem in structural dynamics can be stated mathematically as attempting to find a positive real parameter k in the equation KD ðkÞ ðK kMÞx ¼ 0;
ð1Þ
so as to make the n-dimensional vector x non-trivial, where the n · n dimensional real, symmetric, non-negative definite matrices K and M are the static stiffness and mass matrices, respectively. It is well known that there exists a non-singular matrix X whose ith column is the eigenvector xi associated with the eigenvalue ki, which satisfies ¼ Diagðk1 k; . . . ; ki k; . . . ; kn kÞ. ð2Þ Xt KD ðkÞX By virtue of SylvesterÕs law of inertia [9], if a series of non-singular matrices P1, . . . , Pm is found such that kÞPi ¼ Di ði ¼ 1; . . . ; mÞ; Pt KD ð i
ð3Þ
where each of the Di (i = 1, . . . , m) is a diagonal matrix, then each Di has the same number of negative elements, called the sign count of the dynamic stiffness matrix KD ð kÞ and denoted sfKD ðkÞg. Moreover, Eq. (2) shows that sfKD ðkÞg equals the number of eigenvalues exceeded by k. Knowing sfKD ðkÞg enables bounds on the eigenvalues to be obtained. For example, a trial frequency kk1 is a lower bound on the kth eigenvalue kk if sfKD ð kk1 Þg < k, or an upper bound if sfKD ðkk1 Þg P k. Once lower and upper bounds have been found, convergence on kk is achieved by successively evaluating the sign count at trial frequencies lying between the two bounds. Since only the signs, and not the values, of the elements of the corresponding diagonal matrix are needed when evaluating sfKD ð kÞg, the sign count is quite insensitive to numerical errors. It is this feature of the sign count that makes the algorithms based on it almost infallible. 3. Reduction of dynamic stiffness matrix to diagonal A generalized inner product of two vectors pi and pj of order n is introduced as kÞp . ðp ; p Þ ¼ pt KD ð i
j
ð4Þ
j
i
ð1Þ ð1Þ p1 ; p2 ; . . . ; pð1Þ n
Any set of independent vectors of order n can be transformed to an orthogonal set p1, p2, . . . , pn, in the spirit of the method of modified Gram–Schmidt orthogonalization, as follows. The first desired vector is simply chosen ð1Þ as p1 ¼ p1 . Then, for k = 2, . . . , n, after performing the transformations ðk1Þ pk1 ; pi ðkÞ ðk1Þ p ði ¼ k; . . . ; nÞ; ð5Þ pi ¼ pi ðpk1 ; pk1 Þ k1 ðkÞ
ðkÞ
ðkÞ
ðkÞ
ðkÞ
ðk1Þ
ðk1Þ
is exchanged with pk the kth vector is chosen as pk ¼ pk . If ðpk ; pk Þ ¼ 0, but ðpk0 ; pk0 Þ 6¼ 0 for some k 0 > k, then pk0 to make the choice valid. ð1Þ ð1Þ It can be seen from Eq. (5) that pk is a linear combination of pi (i = 1, . . . , k 1). If each pi is selected as the ith column of the unitary matrix of order n, then the components of pk have the feature that pkk = 1 and pkj = 0 for j > k. The matrix P whose kth column is pk is an upper triangular matrix whose diagonal elements are all unity, and satisfies Pt KD ð kÞP ¼ D;
ð6Þ
where D is a diagonal matrix. Being real and symmetric, the dynamic stiffness matrix can be uniquely decomposed into the form KD ð kÞ ¼ Ld DLt ; ð7Þ d
where Ld is a lower triangular matrix. Comparison between Eqs. (6) and (7) gives
4052
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
f12
f13
2
f23
f14
1 3
2
f24
1
4
3
f11
4
f22
3
f34
1
4 1
2
4
f33
2
3
f44
Fig. 1. Modes of displacements and forces for reduction of dynamic stiffness matrix to diagonal form.
P ¼ Lt d ;
ð8Þ
so that the matrix of congruence transformation can be obtained by the usual Gauss elimination process. If pi is regarded as a mode of general displacements of the structure whose stiffness matrix is KD ðkÞ, then kÞp f i ¼ KD ð i
ð9Þ
can be viewed as a required mode of forces to produce the displacement pi. Defined as the matrix whose kth column is fk, the force matrix F ¼ Pt D
ð10Þ
is a lower triangular matrix. Therefore, if j < i, the jth component of the ith mode of forces fij = 0, while the corresponding component of the mode of displacements pij 5 0. From a physical point of view, pi is generated by clamping the (i + 1)th, . . . , nth degrees of freedom of the structure and applying a force fii at the ith degree of freedom to make pii = 1; the forces fij (j = i + 1, . . . , n) are simply the corresponding constraint forces, as illustrated by Fig. 1. If fii = 0 at a specific trial frequency k, the next mode of displacements is evaluated by replacing k with a new trial fre quency. When k is an eigenvalue of the structure, it is not necessary to evaluate the remaining modes of displacements because pi is the eigenvector being sought. In the resulting diagonal matrix D, the kth diagonal element dk, being defined as the kth energy norm, is the input energy necessary to generate the mode of displacements pi, and is related to the forces through d k ¼ fkk .
ð11Þ
4. Qualitative behaviour of the energy norms Being expressed as d i ¼ pt KD ð kÞp ði ¼ 1; . . . ; nÞ; i
i
ð12Þ
the energy norms are definite functions of the trial frequency k, if the order of freedoms in the structure is fixed. The derivatives of these functions can be expressed in terms of the modes of displacements by the equations dd i dpti dKD ðkÞ ð13Þ pi ði ¼ 1; . . . ; nÞ. ¼ 2 f i þ pti dk dk dk The first term on the right-hand side of Eq. (13) is zero, because the first i 1 elements of fi and the last n i elements of pi are zero, and pii is always equal to unity. Thus, according to the definition of the dynamic stiffness matrix, Eq. (13) can be re-written as dd i ¼ pti Mpi dk
ði ¼ 1; . . . ; nÞ.
ð14Þ
The right-hand side of Eq. (14) must be negative because M is a positive definite matrix. Therefore, di (i = 1, . . . , n) are monotonic decreasing functions of the trial frequency k. Let KkD , Pk and Dk be the principal submatrices obtained by deleting the (k + 1)th, . . . , nth rows and columns from the matrices KD, P and D, respectively. These matrices are related to each other by Ptk KkD ð kÞPk ¼ Dk .
ð15Þ
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4053
Fig. 2. Energy norms as a function of trial frequency.
The kth energy norm can be expressed as D k ð kÞ ; dk ¼ Dk1 ð kÞ
ð16Þ
where Dk ð kÞ, the determinant of KkD ð kÞ, is equal to the determinant of Dk. From the analysis of Section 2, the number of negative elements in Dk is equal to the number of the roots of the equation D k ð kÞ ¼ 0 that are exceeded by k, because the matrix KkD can be regarded as the dynamic stiffness matrix of the structure with the (k + 1)th to nth degrees of freedom fixed. If Dk1 ð kÞ ¼ 0 while Dk ð kÞ 6¼ 0, there must be a small positive real number e such that the sign count of Dk1(k) increases by one while the sign count of Dk(k) is invariant as k increases from k e to k þ e, so that dk must jump from 1 to +1 at k. In general, the curve of the kth energy norm dk has a number of branches separated at the trial frequencies k where Dk1 ð kÞ ¼ 0 while Dk ð kÞ 6¼ 0, and in each of the branches dk decreases monotonically from +1 to 1, as illustrated in Fig. 2. 5. Criterion for eigenvalues in terms of energy norms From the analysis of Section 3, the last energy norm dn is the force required at the nth degree of freedom of the structure in order to produce the mode of displacements pn, the last column of the transformation matrix P. As no forces act at the other degrees of freedom, dn is the only possible non-zero force necessary to produce pn. Therefore, if kÞ ¼ 0; ð17Þ d n ð k must be an eigenvalue of the structure and pn must be the corresponding eigenvector. This criterion for eigenvalues applies to most practical structures except in the rare case where there exists an eigenvalue such that the last component of the corresponding eigenvector happens to be zero. This difficulty can be overcome by changing the order of the degrees of freedom to make Eq. (17) a valid criterion. Suppose a proper lower bound kli and
4054
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
a proper upper bound kui on the ith eigenvalue have been found. Then, by definition, the corresponding sign counts are sfKD ðkli Þg ¼ i 1 and sfKD ðkui Þg ¼ i, respectively, and the criterion expressed by Eq. (17) is valid if the conditions d n ðkli Þ > 0;
d n ðkui Þ < 0
ð18Þ
are satisfied simultaneously. The physical interpretation of the inverse of the dynamic stiffness matrix is well known; i.e., its jth column is the mode of displacements generated by imposing a unit force at the jth degree of freedom while no forces are imposed at the other degrees of freedom. Therefore, when the former jth degree of freedom is re-numbered as the last degree of freedom, the corresponding last energy norm d ðjÞ n can be obtained as d ðjÞ n ¼
1 ; ðK 1 D Þjj
ð19Þ
where ðK 1 D Þjj is the jth diagonal element of the inverse of matrix KD ðkÞ. If a matrix P has been found such that t P KD ðkÞP ¼ D, then the inverse of KD ðkÞ can be obtained easily from 1 t K1 D ðkÞ ¼ PD P .
ð20Þ
Combination of Eqs. (19) and (20) gives 1 d ðjÞ ; ð21Þ n ¼ X 2 n pjk k¼j d k where pij is the element of P at the jth column and the ith row. The components of the corresponding mode of displacements pðjÞ n are expressed as Xn pik pjk k¼minðj;iÞ dk ðjÞ pin ¼ ði ¼ 1; . . . ; nÞ. ð22Þ Xn p2jk k¼j d k Based on the above analysis, when Eq. (17) is employed as the criterion for eigenvalues, one can determine which degree of freedom should be the last one as follows: (1) When a proper lower or upper bound of the required eigenvalue is known, calculate all the d ðjÞ n (j = 1, . . . , n). (2) Check if Eq. (18) is satisfied. (3) If none of the d ðjÞ n (j = 1, . . . , n) satisfy the requirement, then set the trial frequency to the average of the lower and upper bounds and go back to step (1). (4) Any one of the freedoms j that satisfies Eq. (18) can be numbered last.
6. Computational method for the eigensolution For simplification of the notation, and without loss of generality, it is assumed that Eq. (17) is the criterion for all the eigenvalues. Then each eigenvalue corresponds to the zero point of a specific branch of dn. Within the interval ðkli ; kui Þ, dn is a monotonic decreasing function of k, which means that the map from k to dn is one to one. Therefore, k can be viewed mathematically as a function of dn, and the derivative of this function with respect to k is obtained by inverting the righthand side of Eq. (14), i.e., dk 1 ¼ t . ð23Þ dd i pi Mpi In order to find the ith eigenvalue ki, a trial value k0 is chosen within the interval ðkli ; kui Þ and the corresponding last energy 0 norm d n is evaluated. The point with coordinates ðk0 ; d 0n Þ serves as an initial point on the curve generated by Eq. (23). The eigenvalue ki can be found by locating the point (ki, 0). In other words, finding an eigenvalue amounts to solving an initial value problem of the ordinary differential Eq. (23). This transformation is of great value because computational methods for the eigensolution can be derived from a variety of highly developed numerical methods for the solution of ordinary differential equations. Denoting the right-hand side of Eq. (23) as f ð kÞ, the one-step Euler method gives t t p ðK k0 MÞpn p Kp k¼ k0 þ hf ð k0 Þ ¼ k0 þ n t ¼ tn n ; ð24Þ pn Mpn pn Mpn where the step of integration
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
h ¼ d 0n .
4055
ð25Þ
Eq. (24) shows that the method of RayleighÕs quotient can be regarded as the result of the application of the one-step Euler method in solving Eq. (23). Because of its high accuracy and simplicity, the fourth-order Runge–Kutta method is preferred, so that the eigenvalue is estimated as follows: k¼ k0 þ 16f ð k0 Þ þ 13f ð k1 Þ þ 13f ð k2 Þ þ 16f ð k3 Þ h; ð26Þ where k1 ¼ k0 þ 12f ð k0 Þh;
h ¼ d 0n ;
k2 ¼ k0 þ 12f ðk1 Þh;
k3 ¼ k0 þ f ðk2 Þh.
ð27Þ
If the corresponding last energy norm dn(k) is larger than a given tolerance, a new estimate of the eigenvalue is obtained from Eq. (26) with k0 ¼ k and h = dn(k). This process is executed repeatedly until dn(k) is small enough. In the meantime, it is also possible to improve the estimates of the bounds on the eigenvalues, since sfKD ðkÞg can be easily obtained during the evaluation of f ð kÞ. If sfKD ð kÞg ¼ k then k must be an upper bound on the kth eigenvalue kk as well as a lower bound on the (k + 1)th eigenvalue kk+1. If k > klkþ1 then k is a better lower bound on kk+1, and if k < kuk then k is a better upper bound on kk. It is obvious that the derivatives employed in Eq. (26) should be obtained from the same curve. This requires that, when the ith eigenvalue is being sought, sfKD ð kk Þg should be equal to either i or i 1 for each kk (k = 0, 1, 2, 3). Also it is necessary that d n ðkk Þ < 0 when sfKD ðkk Þg ¼ i and d n ðkk Þ > 0 when sfKD ðkk Þg ¼ i 1. In practice, these requirements are usually met as long as the starting trial frequency k0 lies between the proper bounds. If these requirements are not met, a new trial frequency l
u
k þ ki k¼ i 2
ð28Þ
is chosen. The replacement of Eq. (26) by Eq. (28) can only reduce the efficiency of the solution, but can never cause it to fail. Two factors are known to affect the efficiency of the proposed method for eigensolution, as shown in Fig. 3. The first factor, denoted by b, is the local bandwidth of the curve being employed, which is defined as the distance between the two boundary points at which the values of energy norm become infinite. Another factor, denoted by e, is the eccentricity of the curve, which is defined as the least distance between the zero point and the boundary points of the curve, expressed as a fraction of the bandwidth. The larger b and e are, the less effort is needed to determine which degree of freedom is suitable for selection as the last degree of freedom, and the faster the convergence should be. While finding the suitable branch of the curve, as described in Section 5, there may be more than one degree of freedom that satisfies Eq. (18). In this case, the one that minimizes d n ðkli Þ d n ðkui Þ is chosen. Although such a choice cannot guarantee that the branch of the selected curve has the largest b and e, it can prevent the most badly behaved curves from being selected.
Fig. 3. Definition of the factors b and e.
4056
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
7. Simplification of the structure The simplest structure for analysis is a chain structure. Fortunately, any complex structure can be transformed into an equivalent chain structure by applying the Householder transformation [7]. For example, a structure having n nodes with nd degrees of freedom per node can be simplified into a chain structure having n · nd nodes with one degree of freedom per node, as illustrated in Fig. 4. Being symmetric and positive definite, the mass matrix M of the structure can be decomposed into the form M ¼ Lm Ltm ;
ð29Þ
where Lm is a lower triangular matrix. Therefore, the eigenvalue problem of Eq. (1) can be re-written as t t Lm ðL1 m KLm kIÞLm u ¼ 0;
ð30Þ
where I is a unit matrix with the same order as that of the mass matrix. Since any real symmetric matrix can be transformed into a tri-diagonal matrix by the Householder transformation, Eq. (30) can be further simplified to ðK0 kIÞu0 ¼ 0;
ð31Þ
where 2
u0 ¼ Pth Ltm u;
a1
6b 6 1 6 6 6 6 0 t 1 t K ¼ Ph ðLm KLm ÞPh ¼ 6 6 6 6 6 6 4
3
b1 a2
b2
b2
a3
..
.
..
.
..
.
an1 bn1
7 7 7 7 7 7 7. 7 7 7 7 7 bn1 5
ð32Þ
an
t The orthogonal matrix Ph is the matrix of the Householder transformation for L1 m KLm . A great reduction in the computational work can be achieved by dealing with the simplified structure instead of the original structure. Because the stiffness matrix of the simplified structure is a tri-diagonal matrix and its mass matrix is a unit matrix, the corresponding energy norms are obtained as follows:
d 1 ¼ a1 ;
d iþ1 ¼ aiþ1
b2i di
ði ¼ 1; 2; . . . ; n 1Þ;
ð33Þ
and the corresponding modes of displacements are obtained by ð1Þ
p1 ¼ p 1 ;
ð1Þ
piþ1 ¼ piþ1
bi p di i
ði ¼ 1; 2; . . . ; n 1Þ;
ð34Þ
ð1Þ
where pi is the ith column of the unit matrix. Moreover, for a chain structure, when the jth degree of freedom is chosen as the last degree of freedom, in order to obtain d ðjÞ n it is not necessary to calculate the energy norms and the corresponding modes of displacements in the original order. Instead, a simpler algorithm is as follows:
Structure A having 3 nodes with 2 degrees of freedom per node
Chain B having 6 nodes with 1 degree of freedom per node
Fig. 4. Simplification of a structure A into a chain B.
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4057
b2i1 ði ¼ n; n 1; . . . ; j þ 1Þ; di b2 ¼ aiþ1 i ði ¼ 1; 2; . . . ; jÞ di
d n ¼ an ;
d i1 ¼ ai1
ð35Þ
d 1 ¼ a1 ;
d iþ1
ð36Þ
and pn ¼ pnð1Þ ; ð1Þ
p1 ¼ p 1 ;
bi1 p ði ¼ n; n 1; . . . ; j þ 1Þ; di i b ð1Þ ¼ piþ1 i pi ði ¼ 1; 2; . . . ; jÞ. di ð1Þ
pi1 ¼ pi1
ð37Þ
piþ1
ð38Þ
By comparing the computational work necessary for a general structure with that for a chain structure, it is seen that such simplification of the structure is highly recommended. 8. Numerical examples The following examples were solved using a computer program incorporating the proposed method. The simple Examples 1–4 demonstrate the numerical stability of the method in extreme cases that are prone to ill-conditioning. The larger Example 5 demonstrates the potential for applying the method to problems of industrial scale. Example 1. Considering a structure with dynamic stiffness matrix 1 0 2 1 KD ðkÞ ¼ k ; 1 2 0 1
ð39Þ
it can be seen that the two diagonal elements in KD(2) are both zero. However, taking a new trial frequency k0 ¼ 2 1018 , and reducing KD ð k0 Þ to a diagonal matrix gives d1 = 1018 and d2 = 1018. Since the eigenvalues of this structure are 1 and 3, the sign count of KD ð k0 Þ should be equal to 1, which is implied by the values of d1 and d2. Example 2. In the case of KD ð kÞ ¼ Diagð5; 3; 2; 6; 7Þ k Diagð1; 1; 1; 1; 1Þ, the eigenvalues do not always make the last diagonal element of the dynamic stiffness matrix equal to zero. This problem was solved to ensure that the proposed method works properly in this case. From the results shown in Table 1, it is confirmed that the proposed method can automatically recognise which freedom should be the last one and gives the correct answer. Example 3. Let the mass matrix be a unit matrix and the stiffness matrix K = RtVR, where V = Diag(0.001, 20, 500, 3, 7.1 · 108) and R is an orthogonal matrix. Being equal to the diagonal elements of the matrix V, the eigenvalues for this structure have a great variation in magnitude, which indicates that the dynamic stiffness matrix has a propensity for illconditioning. The results obtained by the proposed method are shown in Table 2, where e1 is the error of the last energy norm, e2 is the maximal component in the force vector KD(ki)vi, and ki and vi are approximated eigenvalues and eigenvectors, respectively. The results are very accurate even for such an ill-conditioned problem. Moreover, the eigenvalues and eigenvectors were obtained at the same time, with high accuracy. The physical meanings of e1 and e2 imply that they should be of roughly the same magnitude, which is indeed revealed by the results. Example 4. Consider a structure whose mass matrix M = RtDiag(1, 2, 3, 4, 5)R and stiffness matrix K = Rt Diag(50, 2 · 49.9, 3 · 50.1, 4 · 49.99, 5 · 50.01)R, where R is an orthogonal matrix. The eigenvalues of this structure are all very close to 50. The numerical solution of this problem is shown in Table 3, where the very small magnitudes of e1 and e2, defined as in Example 3, show that both the eigenvalues and eigenvectors were correctly obtained.
Table 1 Results for Example 2, showing the eigenvalues ki and eigenvectors vi i
ki
vi
1 2 3 4 5
2.00000000 3.00000000 5.00000000 6.00000000 7.00000000
0.00000000 0.00000000 1.00000000 0.00000000 0.00000000
0.00000000 1.00000000 0.00000000 0.00000000 0.00000000
1.00000000 0.00000000 0.00000000 0.00000000 0.00000000
0.00000000 0.00000000 0.00000000 1.00000000 0.00000000
0.00000000 0.00000000 0.00000000 0.00000000 1.00000000
4058
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
Table 2 Results for Example 3, showing the eigenvalues ki, eigenvectors vi and measures of mode accuracy e1 and e2 i 1 2 3 4 5
ki
e1
vi
0.001 3 20 500 7.1 · 108
0.6208 0.6095 0.3490 0.1210 0.3267
0.2525 0.7172 0.3796 0.1454 0.5066
0.4128 0.0293 0.0498 0.7578 0.5020
0.5490 0.3366 0.1780 0.5957 0.4457
e2 8
1.125 · 10 2.850 · 109 4.766 · 109 7.449 · 108 5.960 · 108
0.2812 0.0014 0.8366 0.1873 0.4312
5.215 · 108 3.210 · 108 1.490 · 108 5.588 · 108 6.706 · 108
Table 3 Results for Example 4, showing the eigenvalues ki, eigenvectors vi and measures of mode accuracy e1 and e2 i
ki
vi
1 2 3 4 5
49.90 49.99 50.00 50.01 50.10
0.1934 0.8505 0.0930 0.4716 0.0910
0.1450 0.2677 0.5804 0.4162 0.6303
0.1574 0.2141 0.3001 0.5255 0.7505
1
0.1443 0.3973 0.7513 0.4753 0.1761
3
0.9466 0.0366 0.0052 0.3199 0.0172
e1
e2
8.700 · 1015 4.051 · 1015 1.323 · 1014 3.530 · 1014 4.783 · 1014
8.743 · 1015 1.904 · 1015 1.141 · 1014 1.562 · 1014 3.204 · 1014
7
5
3m
4
2
8
6 4m
4m
4m
Fig. 5. A frame composed of beams.
Table 4 Comparison between the proposed method and the inverse power method for Example 5, showing the eigenfrequencies ki (in Hz), the number of iterations n and measure of mode accuracy e2 i
Proposed method
Inverse power method
ki
e2
n
ki
e2
n
1 2 3 4 5 6
15.2828042 38.7297544 44.4735753 48.3991123 50.5691556 57.7034310
0.002413 0.123841 0.001349 0.082575 5.332224 0.001154
4 4 4 4 4 4
15.2837170 38.7990895 44.5553287 48.4878425 50.7001757 57.8424109
0.745206 3.663227 4.991209 5.618068 4.694577 9.630830
7 31 72 89 45 62
Example 5. The lowest six eigenfrequencies and corresponding modes of vibration of the frame shown in Fig. 5 were found. The frame is fully built-in at node 1 and horizontal movement is prevented at node 2. All the beams in the frame have identical cross-sections and are made of the same material, with flexural rigidity EI = 5 MN m2, extensional rigidity EA = 900 MN and mass per unit length l = 35 kg m1. Each beam is approximated by Bernoulli elements of length 0.2 m, so that the frame has 250 nodes and provides a demonstration of the method on a realistically sized problem. The results are shown in Table 4, where e2 is a measure of mode accuracy as defined in Example 3. It is seen that the proposed method is more accurate and more efficient than the commonly used inverse power method.
Z. Qi et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4050–4059
4059
9. Conclusions A new computational method for the eigensolution of structural dynamics has been proposed. In contrast to most existing methods, the proposed method theoretically transforms the eigenvalue problem into a specific initial value problem of an ordinary differential equation. This transformation enables many advanced numerical methods for ordinary differential equations to be used in the computation of the eigensolution. Based on the physical meaning of the sign count of the dynamic stiffness matrix, a stability control device has been designed and combined with the fourth-order Runge–Kutta method. The resulting method finds the eigenvalues and eigenvectors at the same time, with high accuracy and high stability. Numerical examples have shown that even in very unfavourable conditions, e.g., when there is great variation in magnitude among the eigenvalues or some eigenvalues are very close to each other, the proposed method still gives high accuracy solutions. Acknowledgements This work is supported by the China Association for International Exchange of Personnel. The first author is especially grateful for the financial support of Cardiff University and the National Science Foundation of China. The third author has returned to his chair at Cardiff University upon completion of his appointment at City University of Hong Kong. References [1] P. Ladeveze, J.P. Pelle, Accuracy in finite element computation for eigenfrequencies, Int. J. Numer. Methods Engrg. 28 (1989) 1929–1949. [2] N.-E. Wiberg, R. Bausys, P. Hager, Adaptive h-version eigenfrequency analysis, Comput. Struct. 71 (1999) 565–584. [3] P. Hager, N.-E. Wiberg, Adaptive eigenfrequency analysis by superconvergent patch recovery, Comput. Methods Appl. Mech. Engrg. 176 (1999) 441–462. [4] M. Duran, M. Miguez, J.-C. Nedelec, Numerical stability in calculation of eigenfrequencies using integral equations, J. Comput. Appl. Math. 130 (2001) 323–336. [5] K.-J. Bathe, E.L. Wilson, Large eigenvalue problems in dynamic analysis, J. Engrg. Mech. Div. ASCE 98 (1972) 1471–1485. [6] K.-J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, 1996. [7] L. Meirovitch, Computational Methods in Structural Dynamics, Sijthoff and Noordhoff, Netherlands, 1980. [8] A.R. Collar, A. Simpson, Matrices and Engineering Dynamics, Ellis Horwood, Chichester, 1987. [9] L. Mirsky, An Introduction to Linear Algebra, Constable, London, 1982. [10] K.-J. Bathe, S. Ramaswamy, An accelerated subspace iteration method, Comput. Methods Appl. Mech. Engrg. 23 (1980) 313–331. [11] S. Rajendran, M.V. Narasimhan, An accelerated subspace iteration method, Int. J. Numer. Meth. Engrg. 37 (1994) 141–153. [12] K.K. Gupta, Eigenproblem solution by a combined Sturm sequence and inverse iteration technique, Int. J. Numer. Methods Engrg. 7 (1973) 17–42. [13] W.H. Wittrick, F.W. Williams, A general algorithm for computing natural frequencies of elastic structures, Q. J. Mech. Appl. Math. 24 (1971) 263– 284. [14] F.W. Williams, M.S. Anderson, Inclusion of elastically connected members in exact buckling and frequency calculations, Comput. Struct. 22 (1986) 395–397. [15] F.W. Williams, W.H. Wittrick, Exact buckling and frequency calculations surveyed, J. Struct. Engrg. ASCE 109 (1983) 169–187. [16] F.W. Williams, D. Kennedy, Reliable use of determinants to solve non-linear structural eigenvalue problems efficiently, Int. J. Numer. Meth. Engrg. 26 (1988) 1825–1841. [17] A. Simpson, On the solution of S(x)x = 0 by a Newtonian procedure, J. Sound. Vib. 97 (1984) 153–164. [18] C.T. Hopper, F.W. Williams, Mode finding in nonlinear structural eigenvalue calculations, J. Struct. Mech. 5 (1977) 255–278.