Application of Ritz vectors for dynamic analysis of large structures

Compurm & Structurrs Printed in Great Britain.

Vol.21.

No.

5. pp. 901-907.

0045-7949/m

1985

$3.00

+ .Oa

PergamonPress Ltd.

APPLICATION OF RITZ VECTORS FOR DYNAMIC ANALYSIS OF LARGE STRUCTURES R. R. ARNOLDand R. L. CITERLEY Anamet Laboratories, Inc., San Carlos, CA 94070, U.S.A. and

NASA-Ames

M. CHARGIN and D. GALANT Research Center, Moffett Field, CA 94035, U.S.A. (Received 19 January

1984)

Abstract-The use of an orthogonal set of specially selected Ritz vectors is shown to be very effective in reducing the cost of dynamic analysis by modal superposition. Several mechanical structures are examined and the Ritz vector approach is compared to the classical eigenvector approach on the basis of cost, accuracy, and elapsed analysis (throughput) time. Mathematical proof of the completeness of orthogonal Ritz vectors is provided for the case of a positive definite mass matrix and a symmetric stiffness matrix.

INTRODUCTION Since the late 195Os, there have been many significant achievements in the field of dynamic analysis of large structures. Some of these achievements are well known and are recognized by such names as Guyan reduction [l], substructuring [2], and component mode synthesis [3]. Their timely implementation into large scale finite element programs has led to significant improvements in the computational efficiency of dynamic analysis. Without these improvements the most advanced large scale computers would have difficulty completing a dynamic analysis for a moderate to large size problem. The introduction of superelement technology into MSCI NASTRAN in 1980- 1981, along with the capability for handling experimental data, has provided the structural analyst an unrivaled capability for solving large complex problems. In spite of these accomplishments, the dynamic analysis of large structures (~10,000 DOF) requires relatively long computer runs, some of which necessitate labordebilitating overnight runs. A relatively recent paper by Wilson [4] discusses the use of Ritz vectors [5] for the solution of dynamic problems using modal superposition. Wilson derives a numerical algorithm that is used to analyze several small models (cl00 degrees of freedom) of typical steel frame building structures. His work demonstrated improved accuracy with fewer vectors as compared to the use of eigenvectors. Wilson’s algorithm has been implemented into the MSC/NASTRAN computer code for both the conventional and superelement dynamic analysis procedures. The results convincingly demonstrate savings of 90% for eigenvalue extraction and 50% for response calculations using relatively non-optimized MSC/NASTRAN DMAP algorithms. Additionally, it is noted that throughput was increased by a factor of three, with overnight runs virtually

eliminated. Such an achievement is consistent with Wilson’s original goal [6], which was the reduction of computational effort by a factor of 10 for the dynamic analysis of large structures. The research reported herein has been conducted to address several practical issues which arise from the use of a Rayleigh-Ritz procedure. On a practical level, the Wilson algorithm is used to demonstrate the applicability of Ritz vectors to larger multi-degree-of-freedom (DOF) problems. Specifically, accuracy, cost, and elapsed machine time (throughput) are discussed and real-life examples presented. On a theoretical level, the convergence and completeness of Ritz vectors are discussed and a formal proof of completeness is provided.

RITZ PROCEDURE An approximate method for calculating the fundamental frequency of a mechanical system was developed by Lord Rayleigh in 1894 [71. Rayleigh’s well-known assumed that the deflected shape of a vibrating body could be represented by a simple one term expression, such as X = aQ(x, co),

(1)

where a is an unknown amplitude, o is the lowest natural frequency, and @(x, o) is a shape function corresponding to the lowest natural frequency. By equating potential to kinetic energy, Rayleigh was able to estimate the fundamental frequency of the vibrating system. While being very useful for simple structures, the Rayleigh method remained relatively unexploited until Ritz, in 1911, provided a more generalized formulation [5]. Ritz assumed that the deflected shape could be represented by a series of functions a,,(~, 901

902

w) that satisfy the boundary

R. R. ARNOLDet al.

conditions.

subspace iteration with arbitrary starting vectors [lo], in which the lowest N eigenvectors are found by taking L. Ritz vectors with

Thus,

(2) Equating potential and kinetic energies leads to a frequency equation from which the frequencies of the various modes of vibration are calculated. Ritz’s method not only provided a better estimate of the fundamental frequency, but also allowed the approximate calculation of higher frequencies. Both the Rayleigh and Ritz approaches received limited use, primarily because of reasonably accurate estimate of mode shapes was possible only for rather simple systems. With the advent of finite elements and high-speed digital computers, the normal mode approach for characterizing dynamic behavior has been used almost exclusively in the time period 1960-1982. It should be noted that Meek [8] introduced a one-step Ritz procedure for the seismic analysis of frame structures, and Citerley [9] incorporated the Ritz procedure in a finite element program to predict the dynamic response of a structure composed of beams, plates, and shells. By developing an initial approximation to each mode shape, it was shown that accurate frequency and response predictions could be obtained by using Ritz vectors. However, the accuracy of any solution obtained was a direct function of how well the estimated and actual mode shapes compared. Of greater importance, these works demonstrated the usefulness of a Ritz approach for solving large DOF problems when combined with a digital computer and the finite element methodology. The transition from Ritz vectors to eigenvectors seemed quite natural, inasmuch as eigenvectors had a solid mathematical foundation and could be formulated on modern digital computers in a reasonably efficient manner. Additionally, eigenvectors could be calculated directly from the mass and stiffness properties of the structure, without the need to guess mode shapes. Unfortunately, and unrealized, was the extra cost of solving an inherently larger eigenvalue problem. Wilson, in 1982, developed an approach for the dynamic analysis of large structures using an orthogonal set of Ritz vectors which reduced the problem size dramatically. His approach is described in detail in [4] and is briefly reviewed here. Perhaps somewhat by accident, but certainly guided by tremendous insight, Wilson found an approach that would reduce the computational effort for eigenvalue extraction by an order of magnitude. Using special Ritz vectors calculated from the characteristics of the dynamic load, Wilson developed an algorithm that produced not only more accurate response predictions, but also provided an error estimate prior to eigenvalue extraction-the most expensive aspect of the modal method of dynamic analysis. The Wilson approach is very similar to

L = min (2N, N -t 8). Similarly, the MSC/NASTRAN procedure of generalized dynamic reduction (GDR) uses arbitrary Ritz vectors for starting the iteration process. The MSC~AST~N user’s manual recommends a minimum of 1.5N random starting vectors be used to achieve satisfactory prediction of eigenvalues. However, by using a rather special set of orthogonal Ritz vectors, it is shown by Wilson [4] that both fewer vectors and iterations are required to achieve convergence. In Wilson’s method, the first Ritz vector corresponds to the statically deflected shape, with the static load obtained from the time-independent part of the dynamic load. Using inverse iteration and Gram-Schmidt orthogonalization, additional orthogonal vectors are calculated. Once a set of vectors has been determined that reflects the influence of the dynamic load to be considered, a final set of L Ritz vectors is obtained. These vectors are orthogonalized with respect to a reduced stiffness and mass matrix. The entire procedure is outline in Table 1. By contrast, the normal mode approach calculates a set in which some of the eigenvectors are not orthogonal to the dynamic load and therefore do not contribute to the response. Further, it is observed that the use of these special Ritz vectors automatically includes the advantages of Guyan reTable I. Algorithm for generation of orthogonal Ritz vectors {4J 1. Given mass, stiffness matrices M, K and load vector f 2. Trianguiarize stiffness matrix K = LrDL 3. Solve for first vector, x1 Kxt = f x:Mx, = 1 4. Solve for additional vectors i = 2, . . . , L

Kxl’ = Mxi-, Cj =

x,TMx~

XT* = Xi* -

solve for XT normalized and find XI

solve for $ computeforj i-l 2

= I,. . . ,i

-1 Cj Xj

M-orthogonalized

j-l

xTMxi = 1

M-normalized

and find Q

5. O~hogona~ized Ritz vector with respect to stiffness matrix [K* - (!)TM]Zi = 0 X” = xz

Solveforzi(i= I,..., L), where z,r Mzi = 1 compute find Ritz vectors

Applicationof Ritz vectors for dynamic analysis of large structures duction (without the approximation of an assumed mass dis~bution to the residual structure’s boundaries), static condensation, and static correction due to higher mode truncation.

903

such that Mv = 0, is mapped into zero and can never be excited. Thus, eqn (5) can now be compared to the form suggested by Parlett fl I] for convergence of the Lanczos algorithm.

PROOF OF CONVERGENCE OF PROCEDURE

ESTIMATION OF ERRORS

From Step 5 in Table 1 we obtain a system of equations in which M is symmetric positive semidefinite and K is symmetric. In order to show convergence of the basic Ritz procedure, four basic steps must be taken:

The source of errors when using the Lanczos algorithm usually results from round-off errors of the computed vectors in the orthogonality check. Parlett has shown that the Lanczos step can be stated for eqn (5) as

(a) Change coordinates so that M is diagonal. (b) Relabel the coordinates and equations so that any zero eigenvalues occur along the lower corner of the diagonal. (c) Rescale the variables so that the eigenvalues of M are absorbed into the variables. (d) Discard the equations and variables corresponding to the corner of zero eigenvalues.

Zi - Q: Qi = 0,

(6a)

E Q< - Qi Ti s ri ii,

(6b)

Qrri = 0,

where for the ith step, the residual vector ri E ii 4: - qi C+ - qi-1 pi,

Step (a) is accomplished when the eigenvalues A and their corresponding eigenvectors $ are obtained such that M = JITA* +.

II k--l

qi

rj-r/Pi,

E

0-f E

qrI2

II,

qi,

Finally, Pz+1 = /Irj 11.The cancellation of orthogonality is measured by the ratio p?+J(l3? + of + p:+ ,). From Parlett’s discussion of Paige’s theorem, the residual error given in eqn (6~) is bounded by (7)

(4)

Substituting and pre- and post-multiplying by a permutation matrix results in *K+=JIX

Pi E

(3)

In Step 3 of Table 1 a force vector, f, is selected. The choice of this vector is arbitrary. The spatial distribution of the forcing function could be used as the appropriate vector. In addition, the three orthogonal loads (g) could be used for frame structures, and a Fourier harmonic distribution for shell structures. Thus, for static equilibrium Kx = Mg.

(6~)

where e is the round-off error

= h2JIg.

rij

=

Sjj

=

Pi+1

Saj>

eigenvectors Of Ti.

Let A* be such that With the above proof in mind, an error criteria for a specific dynamic load can be developed. If the modal mass is taken to be unity, then an error norm can be defined for the spatial distribution of load, f, as

Z 0 A* A = [ 0 0I and if htL = 0, A& = 0, then

e = fT ElfT f , E=g, where

(5)

(8)

where E =

f - 2 PjMOXj,

Pj =

“Xjf9

7z = A* JI K JITA, X = X+x, g=AJIg.

OX = XQ.

Note that any contribution to the solution that This error norm will vary from e = 1.Oif no vectors arises from the null space of M, i.e. a vector, V, are used, to e = 0.0 if all vectors are used.

904

R. R. ARNOLDet

As an estimate to the frequency content of the forcing function, the power spectra of the load needs to be examined. For a forcing function of the form

al.

Table 2. Computational accuracy of Ritz procedure for the FM2 Ritz vector Frequency no. (Hz)

Modal participation factor

Contribution of Ritz vector

215.5 - 6.37 21.975 -5.1125 13.334 - 14.795 5.4224 - 10.237 4.4256 - 11.565

0.719 0.266 0.026 - 0.015 0.002 0.028 - 0.00137 0.0013 -0.0144 0.00054

Error

F = f(s) F(t) 1

and using the Fourier integral, the power spectral density function of the wave form F(t) is defined by I rsl2 12 s,2 F(t) eeimr dt

(F(w))* = lim _Z

IJ

27Fs

.

For w < a, some cutoff frequency, the frequency content of the solution process can be established to some arbitrary and predetermined limit, y, such that

Thus, with the examination of two limits specified by eqns (8) and (9), a converged solution can be assured and the necessary number of Ritz vectors can be established.

APPLICATIONS

In support of the Air Force Weapons Laboratory, Advanced Radiation Technology Office (ARTO), it has been necessary to find new methods of performing dynamic analysis of super-large multi-degree-of-freedom models (>lOO,OOO DOFs). These

2 3 4 5 6 7 8 9 10

66 90 127 166 170 176 184 563 916 1943

0.280 0.014 -0.011 0.004 0.002 - 0.00096 - 0.00040 -0.00091 0.005 - 0.00001

Conclusion: Only two Ritz vectors are required to characterize response to within 1.4%.

large models occur as a result of the requirement to predict and control the errors that occur in pointing a laser beam. The Airborne Laser Laboratory (ALL), consisting of a KC-135 aircraft, a gas dynamic laser (GDL), a graphite-epoxy optical bench, an Airborne Dynamic Alignment System (ADAS), and an Airborne Pointing and Tracking System (APT), is used as an example for evaluating the application of the Ritz vector approach to dynamic analysis. Specifically, the ALL contains many systems and subsystems that incorporate optical mirrors and sensors. During actual use, the atmospheric environment creates static, dynamic, and thermal loads which cause pointing errors to occur as the laser and tracking beams travel between mirrors. Even relatively small structural deformations

l

- A-set

GRID

Points

Fig.1. Finite element model of the BCM with flexures and flexure mounts.

Application of Ritz vectors for dynamic analysis of large structures

to O-

Weldment

GRID Points

Hidden

Fig. 2. Undeformed

structural plot of the FM2 finite element model.

500.-

s

NASTRAN

SOL.

31

(EIGENVECTORS)

A

NASTRAN

SOL.

31

(RITZ

VECTORS)

400.-

2 g uw 300.. v) P

z 2 zoo.

I

c

100.

0. 0

I

I

2

4

l_._/-,

6 NUMBER

I

a

10

OF VECTORS

Fig. 3. Reduction in CPU time for mode extraction

using Ritz vectors.

1

12

R. R. ARNOLD etal.

906 lOOO.-

0

dASTRAN

SOL.

68

(EIGENVECTORS)

n

NASTRAN

SOL.

68

(RITZ

8

NASTRAN

SOL.

31

(EIGENVECTORS)

A

NASTRAN

SOL.

31

(RITZ

VECTORS)

NON-OPTIMIZED

VECTORS)

OPTIMIZED

800.-

0.1

0

4

2

6 NUMBER

10

a

12

OF VECTORS

Fig. 4. Reduction in total solution time using Ritz vectors.

can produce unacceptable high error in beam pointing. Consequently, extreme accuracy is required to predict the combined response of the models which represent the ALL. For this study, two different mirrors used in the ALL are examined: the fold mirror (FM2) and the beam combining mirror (BCM). The corresponding finite element models are shown in Figs. 1 and 2, and consist of approximately 1,000 DOFs as a result of combining beam, plate, and solid finite elements. The Wilson algorithm of [4] is used with minor modifications compatible with the MSC/NASTRAN computer program. The results concern: (1) (2) (3) (4)

CPU time for eigenvalue extraction. CPU time for total response solution. Accuracy of solution. Throughput (turn-around time).

Shown in Fig. 3 are the modal extraction times for both the Ritz and eigenvaiue procedures for the FM2 model. The Ritz procedure achieves a 10 to 1 advantage at the lower number of modes, and maintains an approximate 10 to 1 advantage for the

higher number of modes. Similarly, the total cost for a response analysis for a similar mirror is shown in Fig. 4, both with superelements and by conventional analysis. The total run cost is cut by approximately 50% for this example. The error associated with each Ritz mode is shown in Table 2. Note that the error used is based on the error norm presented in Wilson’s paper [41 and applies to the particular dynamic load chosen. Using the Ritz error norm, the analyst can determine, prior to the more expensive mode extraction and response calculation, the number of Ritz modes required in the analysis. Tables 3 and 4 summarize the results for both the BCM and FM2 and, as can be observed, the Ritz procedure clearly is superior to the use of eigenvectors. Finally, it is noted that in demonstration tests conducted by the authors on commercial and government computing facilities, throughput was increased by a factor of three. That is to say, three “Ritz jobs” could be submitted and returned in the same elapsed time as was required for one eigenvalue job. Overnight runs occurred for eigenvalue

Table 3. Comparison of eigenvectors

Extraction method Eigen Ritz Eigen Ritz

Number of modes 8 8 12 1.5

Total CPU (set) 625 288 776 331

Eigenvalue extraction module (se@ 386 2 535 16

and Ritz vectors for the BCM Ritz computation (see) 29 62

Elapsed runtime (min) 73 22 128 37

Response at t = .0013 for GP 4330 1.14 2.8 3.0 3.1

x x x x

10-6 IO-’ 10-5 10-s

907

Application of Ritz vectors for dynamic analysis of large structures Table 4. Comparison

Extraction method

Number of modes

Eigen Ritz Eigen Ritz

5 5 10 9

of eigenvectors

and Ritz vectors for the FM2

Total CPU (set)

Eigenvalue extraction module (set)

460 372 650 439

460 2 650 3

Ritz computation (set) 66 63

(1) Runs occurred during period of inconsistent sequently not comparable.

runs when 12 or more modes were requested, whereas the Ritz runs returned during the day even for

15 modes.

The

impact

to the

analyst’s

time

schedule is quite obvious. Finally, it is noted that the use of Ritz vectors is not restricted to structural problems. For any mathematical formulations that use modal superposition of eigenvectors, the use of Ritz vectors has potential application. REFERENCES

R. J. Guyan, Reduction of stiffness and mass matrices. A.Z.A.A. J. 3(2) (1965). NASTRAN Theoretical Manual, Section 4.3. (1972). R. H. MacNeal, Vibrations of Composite Systems.

OSR TN-55-120, Air Research and Development Command (1950). E. L. Wilson, M. Yuan and J. M. Dickens, Dynamic analysis by direct superposition of Ritz vectors.

Elapsed runtime (min) (1) (1) (1) (1)

Response at 60 cps for GP 5277 1.40 1.43 1.43 1.44

x x x x

10-3 10-3 10-3 lo-3

computer operation and are con-

Earthquake

Engineering

and Structural Dynamics

10,

813-821 (1982). 5. W. Ritz, Gesammelte Werke. Paris, 265 (1911). 6. Personal communication. 7. J. W. S. Rayleigh, Theory of Sound, 2nd ed., Vol. 1. Macmillan, London (1894) Sec. 88, reprinted by Dover Publications, New York (1945). 8. J. L. Meek and C. Felippa, DYjlDFR-Dynamic Analysis of 3-Dimensional Frame Structures. Div. of Structural Engineering and Structural Mechanics, University of California, Berkeley (1965). 9. R. L. Citerley, R. W. Clough a&S. Tl Yamahara, The Elastic Response of Structures of Three-Dimensional Frame and Panels in Three Dimensions to Dynamic Loads, Anamet Report No. 966.500 (1966). 10. K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Elemenf Analysis. Prentice-Hall, Englewood

Cliffs, New Jersey (1976). 11. B. N. Parlett, A new look at the Lanczos algorithm for solving symmetric systems of linear equations. Lin. Alg. Appl. 29, 323-346

(1980).