# A combination of the lanczos algorithm with the substructure technique

## A combination of the lanczos algorithm with the substructure technique

Journal of Sound and Vibration (1995) 186(4), 607–616 A COMBINATION OF THE LANCZOS ALGORITHM WITH THE SUBSTRUCTURE TECHNIQUE I.-W. L Department of ...

Journal of Sound and Vibration (1995) 186(4), 607–616

A COMBINATION OF THE LANCZOS ALGORITHM WITH THE SUBSTRUCTURE TECHNIQUE I.-W. L Department of Civil Engineering, Korea Advanced Institute of Science and Technology, Science Town, Taejon 305-701, Korea  G.-H. J Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Science Town, Taejon 305-701, Korea (Received 5 July 1993, and in final form 13 October 1994) A modified Lanczos method combined with a substructure technique was used for calculating natural frequencies and mode shapes of large structural systems. The method does not require generation and storage of stiffness and mass matrices of the original structure. It uses only the stiffness and mass matrices of each substructure. No approximating assumptions are needed other than the usual assumption of a linear elastic system modelled by finite elements. Thus, natural frequencies and the mode shapes for the finite element model employed are the same as those with or without the substructuring algorithm. The method is demonstrated by calculating the first ten natural frequencies and the corresponding mode shapes for an open truss helicopter tail-boom structure.

1. INTRODUCTION

When the dynamic analysis of a structure is performed by using the mode superposition method, free vibration modes and the corresponding natural frequencies for the finite element model of the structure must be calculated first. Because most of the analysis time is spent on finding natural frequencies and the mode shapes, it is very important to calculate them as efficiently as possible. Complex structures are usually analyzed by a mathematical model such as the finite element method. This model has typically a very large number of degrees of freedom, often exceeding 1000 to 10 000. Thus, the substructure technique of dividing the structure into a number of substructures and calculating natural frequencies and the corresponding mode shapes has been developed. An augmented fixed interface mode set was first proposed by Hurty  on the 1960s in connection with the well-known component mode synthesis (CMS) method. Hurty’s work opened up a new area of research on structural dynamics as a number of new CMS methods have since appeared in the literature [2–4]. A CMS method is a Rayleigh–Ritz based approximation method born out of the need to analyze linear structural dynamics problems of unusually high order. In particular, two new mode sets proposed in the early 1970s were shown to have excellent convergence properties in the sence of CMS. They are the MacNeal–Rubin mode set [5, 6] and the Benfield–Hruda mode set . They have the merit that the dimensions of the operating matrix are of low order. However, they give an approximate solution according to the assumed boundary condition. Also a number of 607 0022–460X/95/390607+10 \$12.00/0

.-.   .-. 

608

researchers have developed and modified the CMS method [9–12]. In the 1980s, Arora  proposed an analytical method using the substructure technique and the subspace iteration method. This paper presents an analytical method which combines the Lanczos algorithm and the substructure technique. The method does not require generation and storage of the stiffness and mass matrices of the original structure. It requires only generation of the stiffness and mass matrices of the substructure. The method gives an exact solution of the finite element model because there are no approximating assumptions except the usual assumption of a linear elastic system. Because matrix operations are of the same type, the analysis time for calculating natural frequencies and the mode shapes can be reduced by using the parallel processing capabilities of computers. Thus the method proposed, of combining the substructure technique with the Lanczos algorithm, is an effective solution method for the eigenvalue problem. 2. LANCZOS METHOD

Consider the following generalized eigenvalue problem such as (1)

KF=lMF,

where K and M are the stiffness matrix and mass matrix, l and F are the eigenvalue and eigenvector, respectively, (a list of nomenclature is given in the Appendix). The Lanczos method was originally proposed to tridiagonalize the matrices [14–16]. Once the coefficient matrices of a generalized eigenvalue problem have been tridiagonalized, eigenvalues and eigenvectors can be effectively calculated by using a transformation method such as Jacobi’s method. The method is as follows [17–19]. Let X0 be an arbitrary starting vector for finding the eigenvalues and eigenvectors. Normalize this vector with respect to the matrix M to obtain X1 , with g being the normalization factor: X1=X0/g,

g=(X0TMX0 )1/2 .

(2)

Consider next a general procedure for calculating the vectors, X2 , . . . Xq . Let b1=0, and then calculate the vectors X2 , . . . Xq by using the following iteration scheme for i=2, . . . q: KX i=MXi−1 ,

ai−1=X Ti MXi−1 , Xi=X i /bi ,

X i=X i−ai−1 Xi−1−bi−1 Xi ,

bi=(X MX i ) . T i

1/2

(3–5) (6, 7)

Vectors Xi , i=1, . . . q, generated by using these relations are M-orthogonal: that is, XTi MXj=dij ,

(8)

where dij is the Kronecker delta. The matrix X=[X1 X2 X3 · · · · Xq ] satisfies the relation Tq=XT(MK−1M)X.

(9)

where Tq is the tridiagonal matrix of order q called the Hessenberg-matrix: a1 b2 K G b a b G 2 2 3 b3 a3 b4 G Tq=G · · G · G G k

· · bq−1

· aq−1 bq

L G G G G. G G bq G aq l

(10)

   

609

By using equation (9), one can relate the eigenvalues and eigenvectors of Tq , when q=n, to those of the eigenvalue problem of equation (1), which may be written as (1/l)MF=MK−1MF.

(11)

F=XF ,

(12)

Tn F =(1/l)F .

(13)

With the relation

Hence, the eigenvalues of Tn of equation (13) are the reciprocals of those of the eigenvalue problem of equation (1), and the eigenvectors of the two problems are related, as given in equation (12). If, qQn, equation (13) can be written as Tq F =(1/l)F ,

(14)

which will give approximate solutions of equation (1). When the number of degrees of freedom of the original of the structure is very large, the dynamic analysis cannot be performed by using a small computer, as described in section 1. Thus, the substructure technique is needed to calculate the natural frequencies and mode shapes of the original structure. According to the above iteration scheme, the inversion of the stiffness matrix must be calculated to find X i in equation (3). To calculate X i , in the proposed method one uses the substructure technique described in the following section. The natural frequencies and mode shapes of the original structure thus are calculated by using only the stiffness and mass matrices of the substructure.

3. SUBSTRUCTURE ANALYSIS BY THE DISPLACEMENT METHOD

When the displacements are calculated by using the displacement method, one first analyses the structure with fixed boundaries of all substructures. Next the external loads acting on the interior degrees of freedom are eliminated. The displacements are then obtained by considering the reaction forces needed to fix the substructure boundaries and the external forces acting on the boundary degrees of freedom with relaxation of the boundary conditions for the original structure. Finally, actual displacements of the structure can be obtained by superposition of the above two solutions for the displacements . The complete set of equilibrium equations for the structure, regarded as a free body, may be written in matrix form as (15)

KU=P,

where K is the stiffness matrix of the original structure, and U represents a column matrix of displacements corresponding to external forces P. The stiffness matrix of the rth substructure, regarded as a free body, can conveniently be partitioned into K(r)=

\$

K(r) bb K(r) ib

%

K(r) bi , K(r) ii

(16)

.-.   .-. 

610

where the superscript r denotes the rth substructure and subscripts b and i refer to the boundary and interior displacements, respectively. Naturally, because of the symmetry of (r) the stiffiness matrix, K(r) bi is a transpose of Kib . By using the above stiffness matrix, the (r) substructure displacements U can be related to the external forces P(r) by K(r)U(r)=P(r) ,

(17)

or

\$

K(r) bb K(r) ib

%\$ % \$ %

K(r) bi K(r) ii

U(r) P(r) b b . (r) = Ui P(r) i

(18)

3.1.   When the substructure boundaries on the complete structure are fixed, the boundary fixing must be sufficient to restrain rigid body degrees of freedom on each substructure considered separately. The substructure interior displacements and boundary reactions due to P(r) i when U(r) b =0 can be determined from the equations boundaries −1 =[K(r) Pi [U(r) i ] ii ]

(r) (r) −1 (r) R(r) Pi , b =Kbi [Kii ]

and

(19, 20)

fixed

where the matrix inversion of K(r) ii is permissible because the boundary fixing restrains all rigid body degrees of freedom. Before considering ‘‘matching’’ of displacements on common boundaries, it is necessary to evaluate the substructure stiffness matrix associated with the displacements U(r) b : (r) (r) (r) −1 (r) K(r) b =Kbb −Kbi [Kii ] Kib .

(21)

Equation (21) will be used subsequently to assemble the boundary stiffness matrix Kb for the complete structure . 3.2.   When the boundaries are relaxed, boundary reactions and external forces applied on the boundaries will not be in balance and, therefore, relaxation will induce boundary displacements of such magnitude as to satisfy equilibrium at each joint on the boundary. To calculate these boundary displacements, the complete structure can be regarded as an assembly of substructures subjected to external loading. A typical resultant boundary load Q(r,r+1) between the rth substructure and the (r+1)th substructure is given by b (r+1) =−R(r) +P(r,r+1) . Q(r,r+1) b b −Rb b

(22)

The equation of equilibrium in terms of boundary displacements for the complete structure can now be written as Kb Ub=Qb ,

(23)

where Kb and Qb are (r) Kb=s (a(r) K(r) b a ) T

and

Qb=s (b(r) Q(r,r+1) ), b

r

T

(24, 25)

r

where a(r) and b(r) are Boolean transformation matrices for Kb and Qb . The boundary displacement can be obtained from equation (23) and written as Ub=K−1 b Qb .

(26)

    (r) i

611 (r) i

From equation (18), the substructure interior displacements U due to the forces P and boundary displacements are given by (r) −1 (r) (r) −1 (r) (r) U(r) i +[Kii ] Pi −[Kii ] Kib Ub .

(27)

When the boundary displacements are set equal to zero, the substructures are completely isolated from one another so that the application of an interior force causes displacements only in one substructure. It is evident, therefore, that the interior displacements U(r) i with the boundary fixed can be calculated for each substructure separately, by using equation (19). Although the determination of the boundary displacements U(r) b involves the complete structure, considerable computation advantage is derived from the fact that Kb is of much lower order than the complete stiffness matrix K . (r) Thus Ub and U(r) i can be obtained from equation (26) and equation (27), and Ub is obtained from Ub . Therefore, the displacements of the structure, U, can be obtained. 4. PROPOSED METHOD

Basically, the method is that the natural frequencies and the corresponding mode shapes are calculated by using the Lanczos algorithm and the substructure technique. When using the Lanczos algorithm, one must multiply the mass matrix and a vector as shown in equations (3), (4) and (7), and calculate K−1 to find X  i by using equation (3). The method uses the substructure technique instead of the direct calculation of K−1 to find X  i . That is to say, the inverse operation of the stiffness matrix of the original structure to find the displacement is replaced by the substructure technique discussed in section 3. These operations can be expressed as and

p=Mr

(28, 29)

Kq=r,

where p, q and r are arbitrary vectors and M and K are n×n matrices of the original structure. The substructure technique can be used to find q in equation (29). If equations (28) and (29) can be calculated by using the substructure technique, the substructure technique can easily be applied to the Lanczos method. When the structure is divided into a number of substructures, one does not have the stiffness and mass matrices of the original structure but those of each substructure. Thus, the calculation of equation (28) can be executed as l

l

p=Mr= s (A(r) M(r)A(r) )r= s {(A(r) M(r)A(r) )r}, T

T

r=1

(30)

r=1

where M(r) is the mass matrix of the rth substructure, A(r) is the Boolean transformation matrix for the rth mass matrix, and l is the number of substructures. Now, in respect of the calculation of equation (29), one also does not know the stiffness matrix of the original structure, and therefore one calculates r by using the substructure technique discussed in section 3. Namely, the boundary stiffness matrix Kb and the load (r) (r) vector Qb of the original structure are obtained by synthesizing K(r) b and Pb (=Qb ), and thus it follows immediately that l

(r) Kb= s b(r) K(r) b b , T

(31)

r=1

l

l

(r) −1 (r) (r) Qb= s e(r) [P(r) P(r) b −Kbi (Kii ) Pi ]= s e b , T

r=1

T

r=1

(32)

.-.   .-. 

612

where b and e are the Boolean transformation matrices. Equation (33) is then obtained from equations (31) and (32) as (r)

(r)

Kb Ub=Qb .

(33)

Then, the displacements of the original structure can be obtained by calculating Ub in equation (33) and U(r) i in equation (27), and by synthesizing Ub and Ui . Once equations (28) and (29) have been calculated, a general eigenvalue problem such as equation (1) can be transformed into a standard eigenvalue problem of a tridiagonal matrix. Dynamic analysis can then be performed without knowledge of the stiffness and mass matrices of the original structure when the structure is divided into several substructures.

5. NUMERICAL EXAMPLES

As an example of calculating natural frequencies and the corresponding mode shapes with substructures, an open truss helicopter tail-boom structure  is considered. The geometry of the structure is shown in Figure 1. The structure is modelled by 108 truss elements and 28 joints as shown in Figures 1–3. The structure has 72 degrees of freedom, so the eigenvalue problem is of dimension 72×72. For example calculations, the cross-sectional area of all members is 1·0 in2, the material specific weight 0·10 lb/in3 and Young’s modulus is 10·5×106 psi. The tail-boom structure is divided into three substructures by partitioning it at nodes 9–12 and 17–20, as shown in Figure 3. Substructures 1 and 2 each have 24 boundary degrees of freedom, 12 interior degrees of freedom and 36 members. Substructure 3 has 36 members, and 12 boundary degrees of freedom and 24 interior degrees of freedom. The lowest ten natural frequencies and the corresponding mode shapes of the tail-boom structure were calculated with and without substructures. Of course, the natural frequencies and the corresponding mode shapes are exactly the same in these two cases (see Table 1). To calculate eigensolutions, the subroutine program needed to find the stiffness matrix and mass matrices, etc., was extracted from ADINA (A utomatic D ynamic I ncremental N on linear A nalysis) program. All the elements of the starting vector for both programs are 1. When one computes ten natural frequencies by using 20 Lanczos vectors, the computing time of the program based

Figure 1. Geometry of helicopter tail-boom.

   

613

Figure 2. Finite element model for the tail-boom structure (degrees of freedom for joints 5 and 28 are, 1, 2, 3, and 70, 71, 72 respectively).

on the algorithm with and without substructuring are 46·58 s and 42·42 s, respectively, for the SUN3-60 computer. These computing times do not include the times needed to calculate the stiffness and mass matrices. the back-up storage by the sky-line method without substructuring is 2286, with substructuring 2382 (substructures 1, 2 and 3 are 582, 900 and 900 respectively). These results are summarized in Tables 2 and 3.

Figure 3. Nodal numbering system for substructural formulation for the finite element model of the helicopter tail-boom structure. (a) An overall numbering for the boundary nodes; (b) numbering system for boundary and interior nodes for each substructure. Note that for clarity the diagonal members are not shown.

.-.   .-. 

614

T 1 Natural frequencies of helicopter tail-boom structure Natural frequency (no.) 1 2 3 4 5 6 7 8 9 10

Natural frequency Lanczos method Lanczos smethod (full structure) (substructure) 21·845 23·129 101·642 105·297 107·370 200·832 227·730 239·328 241·834 377·745

Exact solution

21·845 23·129 101·642 105·297 107·370 200·832 227·730 239·328 241·834 377·835

21·85 23·13 101·6 105·3 107·4 200·8 227·7 239·3 241·8 377·7

T 2 Information on the stiffness and mass matrices Substructure 1 Substructure 2 Substructure 3 Original

NEQ

NWK

NWM

BW (mean)

BW (max.)

24 36 36 72

291 450 450 1143

291 450 450 1143

13 13 13 16

21 21 21 21

NEQ, number of equations; NWK, number of stiffness matrix elements; NWM, number of mass matrix elements; BW, half bandwidth.

6. CONCLUSIONS

An efficient numerical procedure for calculating natural frequencies and the corresponding mode shape for large (or complex) finite element models of a structure by using substructures has been developed. With the proposed method, it seems disadvantageous that analysis time and backup storage are increased more than those without substructures. If, however, the number of degrees of freedom of the structure is very large (with relatively small boundary degrees of freedom), analysis time and backup storage do not cause problems because they increase very slightly. The method developed can also be practical for calculating natural frequencies and mode shapes for a large finite element model of a structure by using a small computer. This is the most important feature of this method. The method is based on the Lanczos algorithm and the substructure technique. The method does not use the component mode substitution concept nor approximating assumptions. Thus, for a given finite element model of the structure, natural frequencies and

T 3 Comparison of Lanczos method and the proposed method Method Lanczos method (full structure) Proposed method (substructure)

Space (NWK+NWM) 2286 582 900 (2382) 900

Time (CPU) 41·42 s 46·58 s

   

615

the corresponding mode shapes are exactly the same with or without the substructuring method. When the structure is divided into several substructures, there are matrix operations of the same types. Also, it has been pointed out that the analysis time for natural frequencies and mode shapes can be reduced by using the computational parallel processing procedure for each substructure. REFERENCES 1. W. C. H 1965 American Institute of Aeronautics and Astronautics Journal 3(4), 678–685. Dynamical analysis of structural systems using component modes. 2. G. M. L. G 1964 Journal of Sound and Vibration 1, 41–59. Branch mode analysis of vibrating systems. 3. R. R. C, . and C.-J. C 1977 Proceedings of the AIAA/ASME 18th Structure, Structural Dynamics, and Materials Conference, B, 89–99. On the use of attachment modes in substructure coupling for dynamic analysis. 4. R. M. H 1975 American Institute of Aeronautics and Astronautics Journal 13(8), 1007–1016. Analytical methods in component modal synthesis. 5. R. H. MN 1975 Computers and Structures 1(8), 581–601. A hybrid method of component mode synthesis. 6. S. R 1975 American Institute of Aeronautics and Astronautics Journal 13(8), 995–1006. Improved component-mode representation for structural dynamic analysis. 7. W. A. B and R. F. H 1971 American Institute of Aeronautics and Astronautics Journal 9(7), 1255–1261. Vibration analysis of structures by component mode substitution. 8. G. A. M 1984 AIAA/AAS Astrodynamics Conference, Seattle, Washington. A modal reduction method for use with nonlinear simulations of flexible multibody spacecraft. 9. P. C. H 1980 Journal of Applied Mechanics 47, 177–184. Modal identities for elastic bodies with application to vehicle dynamics and control. 10. W. S and S. D 1981 Journal of Mechanical Design 103, 643–651. The application of finite element methods to the dynamic analysis of flexible spatial and co-planar linkage systems. 11. W. S. Y and E. J. H 1986 Journal of Structural Mechanics 14(1), 105–126. Dynamics of articulated structures, part I: theory. 12. W. S. Y and E. J. H 1986 Journal of Structural Mechanics 14(2), 177–189. Dynamics of articulated structures, part II: computer implementation and applications.. 13. J. S. A and D. T. N 1980 International Journal for Numerical Methods in Engineering 15, 333–341. Eigensolution for large structural system with substructures. 14. T. J. R. H 1987 The Finite Element Method. Englewood Cliffs, New Jersey, Prentice-Hall. 15. B. N-O and R. W. C 1984 Earthquake Engineering and Structural Dynamics 12, 567–577 Dynamic analysis of structures using Lanczos coordinates. 16. J.H. W 1965 The Algebraic Eigenvalue Problem. Oxford: Clarendon Press. 17. K. J. B 1982 Finite Element Procedures in Engineering Analysis. Englewood Clifis, New Jersey, Prentice-Hall. 18. B. N-O and B. N. P 1983 Interntional Journal for Numerical Methods in Engineering 19, 859–871. Lanczos versus subspace iteration for eigenvalue problems. 19. V. I. W, R. K. R and C. N. C 1983 Computers and Structures 16(4), 253–257 Lanczos eigenvalue algorithm for large structures on a minicomputer. 20. J. S. P 1968 Theorey of Matrix Structural Analysis. New York: McGraw-Hill . 21. R. J. G 1965 American Institute of Aeronautics and Astronautics Journal 3(2), 380. Reduction of stiffness and mass matrices. 22. R. U 1966 Journal of Sound and Vibration 4(2), 149–155. Reduction of the number of unknowns in the displacement method applied to kinetic problems.

APPENDIX: NOMENCLATURE X 1 , X 2 , . . . Xq X , X T q a 1 , bi

Lanczos vectors working vectors defined in equations (3)–(7) tridiagonal matrix components of tridiagonal matrix, Tq

616 l F F U Ub Ui U(r) U(r) b U(r) i P Pb Pi P(r) P(r) b P(r) i Qb Q(r,r+1) b R(r) b K K(r) Kb K(r) b Kbb , Kbi , Kib , Kii , (r) (r) (r) K(r) bb , Kbi , Kib , Kii M M(i)

.-.   .-.  eigenvalue defined in equation (1) eigenvector defined in equation (1) eigenvector defined in equation (12) [Ub Ui ], structure displacements boundary displacements interior displacements [U(r) U(r) b i ], substructure displacements boundary displacements for rth substructure interior displacements for rth substructure [Pb Pi ], externally applied forces boundary forces interior forces [P(r) P(r) b i ], externally applied forces on the rth substructure boundary forces on the rth substructure interior forces on the rth substructure resultant boundary forces for original structure resultant boundary forces between rth substructure and (r+1)th substructure boundary reactions due to P(r) i stiffness matrix for original structure stiffness matrix for rth substructure boundary stiffness matrix for original structure boundary stiffness matrix for rth substructure component submatrices of K component submatrices of K(r) mass matrix of original structure mass matrix for ith substructure