Model reduction methods for dynamic analyses of large structures

compurers & Stmcrurcs Vol. 47. No. 4/s. pp. m-149, Printed in Great Britam.

1993

00457949/93 s6.00 + 0.00 0 1993 Peqtmon Pres¶ Ltd

MODEL REDUCTION METHODS FOR DYNAMIC ANALYSES OF LARGE STRUCTURES B. H;icXXu_ADt$ and L. ERIKSWN~ tRoya

Institute of Technology, Department of Solid Mechanics, S-100 44 Stockholm, Sweden

SASEA Brown Boveri, Corporate Research, S-721 78 Vi&r&s, Sweden Airs&a&-User-friendly finite element pre-processors permit detailed moclelling of engineering structures with complicated topologies, often resulting in a much larger number of degrees of freedom than that motivated by the expected complexity of the structural response. Dynamic analyses of these large-scale models might require a prohibitive amount of CPU time. In particular, the uncoupling of the equations of motion for linear systems with non-proportional damping by means of a complex eigenvector basis is computationally very demanding. We present general model reduction procedures, applicable to whole structures or substructures, that produce a series of reduced models each spanning appropriate subspaces (general Krylov spaces) of the solution space. Crucial steps in such procedures are the selection of basis vectors for the (sub)structures and the control of reduction errors. The set of basis vectors is generated by using a simple physical approach for successively reducing the residual error in the governing equilibrium equations. A frequency window method (shifting) is used to capture the relevant frequency content. The residuals of the equilibrium equations, computed on fully expanded level, are used as reliable measures of the reduction errors. The efficiency of the techniques is demonstrated in some numerical experiments and on a large engineering problem, using a program based on SAP IV.

1. lNTRODUCl’ION

Important types of dynamic analyses are computationally very demanding and can only be effectively performed on smaller models. In fact, many successful algorithms for real eigenvalue analyses are based on efficient model reduction procedures (e.g. subspace iteration and the Lanczos method). Other typical engineering examples are large structures with non-proportional damping or unsymmetric system matrices. Recent developments of computer hardware have mainly focused on powerful engineering workstations, equipped with fast (parallel) processors and relatively large core memories. This fact tends to change the traditional picture of how an optimal finite element analysis should be organized. For instance, it is clear that one should avoid extensive use of slow external memories. Here we outline model reduction techniques that create a series of reduced models that each can be effectively analysed in the core memory. For very large structural models the model reduction should be performed in two steps, first using reduced models for the substructures and then a fmal reduced model for the assembled superstructure. Crucial activities in such a procedure (often called dynamic substructuring) are the creation of the subspaces for the description of structure/substructures and the control of the total error committed during the reduction process. In dynamic substructuring the structure is considered as composed of a number of substructures (components). The dynamic behaviour of each substructure is described by a set of modes that are typically non-local and very few in comparison to the

original number of dynamic degrees of freedom in the substructure. Hence, taking full advantage of modem computer architecture, in this way it is possible to perform fast and accurate modal analyses of ‘superlarge’ structure models. Currently, this means models with greater than, let us say, 100,000 original degrees of freedom. As is well known, in addition to the computational aspects there are several practical advantages for engineering organizations of using substructuring. It makes it easier to do data preparation, verification of parts, reanalyses after modifications, etc. Much has been written in the engineering literature on dynamic substructuring (appearing under such names as Guyan reduction, substructuring, component mode synthesis etc), see, for example [14]. Different ad hoc ways of coupling substructures and selecting component modes especially have been investigated [5,6]. Relatively recently more. systematic approaches to model reduction of dynamic systems have appeared. Wilson et al. [7-lo] have used an orthogonal set of specially selected Ritz (Krylov) vectors for the solution of dynamic problems by modal superposition. In a series of papers NoorOmid et al. [ 1l-1 51and Chen and Taylor [ 16-181 have proven the usability of especially orthogonal&d Krylov vectors in various dynamic analyses, resulting in different kinds of Lanczos or Amoldi methods. In general, excellent performance is shown for small problems. However, the numerical implementation of these algorithms should be made carefully to avoid loss of orthogonality due to round-off. Better efficiency is obtained by using block Lanczos [ 12,191. In the context of substructuring Craig et al. [20-231

736

B. H~GBLAD and L. ERIKSSON

used block Krylov techniques and they demonstrated good performance on small problems. In this paper we present some systematic and accurate model reduction techniques (based on the research in [26]) that are efficient for large structures of high complexity, typical for engineering problems in industry. A formalism is introduced where the set of basis vectors is generated by using a simple physical approach based on successive reduction of the residual errors in the governing equilibrium equations. The recurrence equations for the generation of the vector space (a general Krylov space) and the reduced equations of motion are analogous to the ones used in the Lanczos or Arnoldi algorithms. Efficiency is gained by a choice of starting vectors that is adapted to the mass distribution or coupling conditions on the boundary of the region, that might be the whole structure or a part of it (a substructure). Since the basis vectors are non-local the conditioning is in general improved and the irrelevant, crude approximations to the higher modes in a FE mesh are kept out. The organization of the paper is as follows. First we present the general idea behind the model reduction procedure and construct recurrence relations for classically and non-classically damped structures. The use of frequency windows is demonstrated. By using a shift in the frequency domain we can construct ‘high frequency’ substructures suitable for high frequency response calculations [25]. Then we review the algorithmic recurrence relations that must be used for constructing the reduced equations of motion and discuss some inherent computational difficulties. For a reliable model reduction procedure it is important to have good tools for error control, since in general much larger approximation errors are likely to occur than what is expected in the use of conventional finite elements. We concentrate on reliable a posteriori error bounds and present some results from simple numerical tests. Finally a large engineering example is analysed. 2. MODEL REDUCTION

The quality of the reduced model strongly depends on the properties of the subspace to which the approximate solution is restricted by the reduction procedure. Consider the dynamic behaviour of a linear structure or a non-linear structure in the neighbourhood of a chosen reference state. By a suitable FE discretization it is assumed that the (linearized) governing equations can be written on the partitioned form shown in eqn (1), where the mass matrix is denoted by M, the damping matrix by C and the stiffness matrix by K

Here index i refers to (interior) degrees of freedom where the applied loads are prescribed. Index b refers to (boundary) degrees of freedom with prescribed (boundary) motions. f, and f, are the corresponding forces. We want to describe the interior dynamic behavior of the (sub)structure by sets of modes that are few in comparison to the original number of dynamic degrees of freedom in the structure. For efficiency, the choice of these sets should be made adaptively, governed by inertia, damping and external loadings in the interior and prescribed motions on the boundary of the structure. If the system (1) is the result of a linearization of a non-linear system around a certain configuration the matrices in eqn (1) might well be unsymmetric. This fact will not change the general principle behind the reduction procedures. However, the feasible size of the final reduced system will be much smaller. Problems of this kind arise in rotor or vehicle (multibody) dynamics and aircraft flutter, for instance. In studying the behaviour of ui [interior to the (sub)structure] we consider u*(t) (i.e. even ti6 and iib) to be a known, but otherwise arbitrary quantity. For any combination of ub, ri6 and ilb the forces f,, are determined by the second row in eqn (1). From the first row of eqn (1) we obtain K,, u, = - K,bub- M,biib - M,ii, - C&ii6 - Ciii, + f,

or for non-singular II, = -K,;‘K,u,

(2)

K,, - K,;‘M,ii,

- K,;‘M,,ii,

- K,y’CibQb- K,‘C,,I,

+ K,‘f,.

(3)

This relation tells us that the solution can be written as a linear combination of vectors belonging to the union of the spaces that are spanned by the vector K;‘f, and the column vectors of the matrices K,‘K,b, K, ‘Ma, K,; ‘C,, K,; ‘C, and K,; ‘M,, . But for nonsingular M,i and/or C,, the column vectors of K,;‘M, and/or K,‘C,i trivally span the whole solution space. In a reduced model this last set is replaced by a much smaller set that still gives good accuracy. To avoid truncation in the description of the boundary motions and also guarantee compatibility at the interfaces of substructures, all boundary (physical) degrees of freedom will be retained in our reduced models. The general loading situation in eqn (2) can be considered as a superposition of two cases: (I) homogeneous boundary conditions with f, # 0, ub = 0 (usually applicable to ‘main’ (complete) structures) and (II) inhomogeneous boundary conditions with f, = 0, ub # 0 (applicable to substructures, for example). In the following we will motivate our model reduction procedure by using intuitive (physical) arguments,

Model reduction methods for analyses of large structures

starting with undamped systems. This will lead us to a natural choice of approximating subspaces: Krylov subspaces for properly chosen system matrices. 3.UNDAMPED

SYSTEMS

Homogeneous boundary conditions First , we consider the undamped case (C = 0), for homogeneous boundary conditions (case I). For simplicity, we suppose that the load f, can be written on separable form

f,=F,,Fz,h,...,F,l x k,(t), &O)rg3(t)9

* **v

g,w

= Fg,

(4)

where F is a time independent matrix with column vectorsF,,F,,F, ,,.., F, describing the spatial distribution of the force f, and g is a vector containing the corresponding time functions. Normally r (e n,, where ni is the number of interior degrees of freedom. Equation (3) then reduces to

737

where K,,,= [Go]m(Y’o) (m = 0, . . . ,j)

(10)

This is a series of block Krylov subspaces (231 for the matrix Go and normally j can be chosen such that r xjeq. In theory, this motivates the following expansion of the interior solution for the reduced model U,= KoUb+ K,qj')+ K2qj2'+. ..+

K,q,

(J’,

(11)

where the degrees of freedom are the boundary displacements II* and the interior generalized coordinate vectors 4” (r x 1). K,, (ni x nb) describes the static response and a,,, (n, x nb) (m = 1, . , . , j) contain the corrections due to dynamics. However, in practice with increasing j the Krylov blocks tend to become linearly dependent. The numerical realization must rely on a robust procedure for construction of a basis of linearly independent vectors belonging to the span(Ko,K1,K2,W3,...,Kj).

Inhomogeneous boundary conditions ui = -K,;‘M,,g,

+ K,;‘Fg(t).

(5)

We can formally decompose the solution into components a, = “lo’ + u!” I + &’ I + . . . + $’ + . . . ,

For the inhomogeneous eqn (3) 4”’ = -K,‘K,u,

(6)

where each component is successively computed by balancing the inertia forces originating from the previous component in the series. This means that the residual errors in the governing equilibrium equations are successively reduced. Thus, starting with the static part we obtain by inserting eqn (6) into (5)

case II we obtain from

n$” = -K-‘M,g, II

- KT’M_~!~’ I, II I

= -[l(ilMib II{‘)= -I(,

+ K,-‘M,(K,‘Klb)]Ub

‘M,,gj”

etc. or recursively (n > 1)

u{“’= K,y‘Fg(t )

n$“‘= -K;‘M,#-“.

u$” = -K,lM,,ti$o’ = -K,‘M,,K,‘Fg(t)

(12)

Introducing

a? = -K,‘M,,iij”

F b. = -K;‘K&

(n, x nb)

etc. and or recursively (n > 1) pb’ = -P,‘M, nl;“‘= -K,‘M,,g)“-“.

+

K?M,t(K;‘Kdl

(7)

(4 x 4)

(13)

Introducing Go = -K,‘M,,

and

r. = K;‘F

(8)

we see that the successive improvements of the solution are obtained by enlargements of the solution space with blocks of vectors (n, x r) %KI,IC2,K3,...,K/

CA.5 47/4/s--p

(9)

we see now that the blocks in the sequence (9) are built up according to Ko=ybo a,,, = [Go]‘“-‘Fbl (m = 1,. . . , j).

(14)

3. H~GBLAD and L. ERIKSSON

738

Again we see that this is a series of block Krylov subspaces for the matrix G,,. Normally j can be chosen such that nb x j 6 n,. Remark 1. For the eigenvalue problem (case I) we have f, = 0 and we can use the same scheme as given by eqns (9) and (10) but the starting block of vectors must be chosen differently from eqn (8). Instead of eqn (5) we obtain (assuming time dependence e’“‘) Ki,u, = w~M,,u,.

(15)

Here we use a block of np starting vectors P,, (n, x n,), chosen either as a set of random vectors or as a set of unit vectors with unit entries in components that correspond to the least dominant diagonal entries of the matrix M; ‘K,,, see [27]. Remurk 2. For the inhomogeneo~ case (case II) we see that an efficient alternative to the scheme (13,14) might be to enlarge the solution space with blocks of nb + n, column vectors belonging to the matrices

K, = [G~-‘Y’b,~G~!f’O]

(m = 1,. . . j)

(16)

where !POis constructed as in Remark 1 and rc,,,is an n, x (nc + n,) matrix. Normally j can be chosen such that (nb+n,) xj en,.

::)(:;)=w2(::

~~)(~~)+(t”,)~ (17)

where a2 = the square of the global eigenfrequency and the first row can be written K,ui = -Koub + w2Mibub + 02M,ui.

(18)

Here we can enrich the Krylov spaces with components belonging to eigenspaces associated with eigenvalues in a certain interval by choosing a shift p2 in that interval in the generation of the Krylov blocks. Thus, instead of eqns (12) we obtain (Ii, - ~~ji)u~o~= - K fb”b (Ii,, - /~~M,,)uj’)= 02M,~ub f (w* - p2)M II14” 1 (K, - p2M,)uj2’ = -(--a2

-(-o*-t/~~)M~,u~-‘).

(19)

Hence, introdu~ng G, = -o(;, - P%)-‘M,,

(n, x n,)

!Pg = -(K,, - prMi,)-‘K,,

(ni x n,)

tpftt,= (K,t - $MJ-‘M~

+ sG,,P#

(nn,x n,,) (20)

we obtain the following block Krylov subspaces that replace the Ones in sequence (g) rc$), i@, K$@,KY),. . . ) rc?’ with KP) = r&J

+p2)M,uj”,

(m = 1, . . . ,i).

(21)

The relations in eqn (16) are transformed

into

KP) = yi# rck)= [G~-‘Y”lf;)lG~Y’,,]

SHIITINC

We exemplify frequency shifting on an undamped system and consider a vibrational problem, an eigenvalue problem or a forced response problem,

i’e* (:Ij

(K,,-$M,,)ul”)=

~~)=[G~~-l~~)

rco=‘y’,

4. FREQUENCY

or recursively

(m = 1,. . . ,j).

(22)

In eqn (20) the factor s is given by f = (02 - /f2)/w2. For zero shift s = 1 and if we are looking for eigenvalues close to the shift p2 one should keep Is14 1. Assuming a diagonal mass matrix (M, = 0) we find that for frequency response calculations with the exciting frequency equal to the shift frequency F, the reduced substructure model is exact when only the first Krylov block r# is used and rc$) = 0 (m > 0). 5.

DAMPED

SYSTRMS

Propor~~onQ~damping

For systems with proportional damping the same Krylov blocks are applicable as for undamped systems according to the previous section. A common practice in the analysis of systems with smail nonproportional damping is to project the damping matrix on a subspace of the undamped eigensystem and then ignore the offdiagonal terms (‘the decoupling approximation’). However, for structures with close eigenfrequencies this can cause substantial errors in the response if the system is excited near these

139

Model reduction methods for analyses of large structures frequencies [24]. In this and other more typical cases, e.g. when strong localized dampers are used, we have to resort to efficient procedures for computation of complex eigensystems. Non -proportional damping, homogeneous boundary conditions

For structure with non-proportional damping and homogeneous boundary conditions (case I, applicable to whole structures, for example) we use the substitution v, = i, and rewrite the system M,,ii, + C,,L, + K,,u, = fi(f)

of the symmetric and sparsity of the matrices M,, C,, and K,,. Non -proportional dantping, inhomogeneous bounakry conditions

To be able to sort out different contributions to the solution space we suppose that the time dependence is given by the factor $’ so that the solution can be written as a formal power series in L. For non-homogeneous boundary conditions (applicable to substructures, for example) and non-proportional damping we obtain

(23) Kliui = -K,bu, - lCibub - 12Mbu,

on state form - LC,u, - l*Myu,. (;,

T)(;;)+(;

_\,](;;)=(fiz))

(24)

or Hdv + Aw = d(r),

0 “1

.

V,

Setting d = FIOl[slOIT =

DklW,

(26)

where D is a time-independent matrix and introducing a shift P, eqn (25) might be written in the form w = -(A - ,uH)-‘H& + (A - pH)-‘D[gIO]‘.

With a ‘shift’ p for enriching the Krylov space with components belonging to eigenspaces near fl, eqn (30) is changed

(25)

where H and A are the coefficients matrices, d(t) is the right-hand side in eqn (24) and w=

(30)

(27)

-I%&ub

-(A* - /J*)M~u~. (31)

- (A -p&u,

(It should be noticed that p is not a proper shift for 1.) But for a complex (imaginary) P the left-hand side of eqn (31) will be complex. Since we want to work with real quantities on substructure level we prefer a shifted stiffness matrix that is real for a purely imaginary frequency shift: (IL, + p*M,). Hence instead of eqn (31) we use: (K,i + /J*M..)u!” ” I = -K,q (K..I, + p*M ”+I!‘) I = - rlC,u, - lC,,uj”’ = J [ - G + GQ,,

Introducing r,,=(A-pH)-‘H !P,,=(A-pH)-‘D

+P*W-‘&h

(2nix2n,) (2nix2r)

(28)

(K, + ~*M,)uj*’ = -1*M,q

- lC,iuj”

-(a* - /~*)M,,ujOj

we obtain the Krylov blocks

= -r12M&

I@“,KY', K$“),KY’, . . . , Kj”’ with

- J*C,,(K, + p*M,)-‘[-Ca

Kp= !P#

+ C,i& + ~%)-‘KtiPb + (A*- P*)M,,&

(m = l,...

,j)(normally

2r xj e 2n,).

(29)

For efficiency the shift p should be chosen complex or at least imaginary @ = iwO). Hence complex arithmetic is needed. However, it is not necessary to form the matrices A and H explicitly. On partitioned form it is possible to take full advantage

+ P*M,,)-‘Ktiub or recursively (n > 2) (K,+~*M,,)uj”‘=

-LC,,u)“-” -(A* - fi*)M,,u~-*‘.

(32)

B. H~GBLAD and

740

From eqn (32) it is clear that the contributing component uQ) is associated with the power I’ (at least for p = 0). Thus formally the solution should be sought in the span(@), KY),KY),KY),. . . , K,“‘), where the generalized block-Krylov subspaces are given by the column vectors of the matrices @‘I:

-(K,, + $M,)-‘Kib

~‘l”‘= -(K,, + p*M,,)-‘(Cib + CfiKg)] K$“)= -(K,, + p*M,)-‘[C,rc~

+ sM~,Ic~) - M,]

L. ERIKSSON

For the substructures we demand that the reduced equations of motion shall conform to the conventional way of describing finite elements, i.e. it must be possible to assemble the matrices by simple addition of contributions belonging to common degrees of freedom. The conventional Ritz projection fulfils that condition and we use it to formulate the reduced equations of motion for substructures. However, for main structures consisting of a whole structure without substructures or an assembled superstructure, we use Lanczos method for symmetric/nonsymmetric and Amoldi method for non-symmetric systems. Substructures For a substructure we employ Ritz method and project the substructure equations on a subspace spanned by the set of orthogonal&d Krylov vectors. We consider the model reduction as a transformation

and recursively (n > 2) K$‘)= -(K, + ~zMii)-‘[C,,rc~: I + sM&‘!~]. (33)

(35) In eqn (33) the factor s is given by s = (12 - /l*)/n*.

(34)

For zero shift s = 1, and if we are looking for the eigenvalues close to the shift p one should put s = (A*- p*)/n* x 2. The precise value is not critical. Remark 3. For a well-defined substructure, physics will secure linear independence of the vectors in the blocks K$‘),KY), w$‘),etc. However, for the rest of the blocks the situation is different. The damping matrix for localized dampers contains very few elements different from zero. Then it is likely that the block rep) adds few new linearly independent vectors to the Krylov block ~8). However, it is easily seen that the rank of C, and C, is a simple function of the number of damping elements. For few damping elements the existing independent columns of the matrix (C, + C,,r#)) are best created on damping element level. The blocks KY), K’@, a$“), etc. are treated analogously. In this way it is possible to contruct an efficient procedure for creating a vector basis belonging to the span (rep), a?), rc$‘),K$“),. . . , rcf’)). 6. REDUCED EQUATIONS OF MOTION

In the previous sections we have outlined recursive techniques for generating a Krylov space with presumably good approximation properties. In practical numerical work we need an algorithm for generation of linearly independent basis vectors that span the Krylov space. This is usually accomplished by incorporation of some sort of (mass) orthogonalization in the generation procedures. This orthgonalization and the subsequent formulation of the reduced system equations can be done in a number of alternative ways. In our approach we have worked with some modified variants of the well-known Ritz and Lanczos (Amoldi) procedures.

where the span (Q) = the span&, , x2, x3, . . . , IC,)and lpbo= rcb). We obtain

Model nduction methods for aimlyses of large structures

Equation (35) defines the substructure modes by enlarging the vectors in the Krylov blocks and supplying them with values in the boundary degrees of freedom. Each vector in the block K&“)is enlarged with a unit value in an associated boundary degree of freedom. When the shift frequency p is equal to zero the columns of

741

Introducing right (Q,) and left (P,) Lanczos vectors (w = Q,x) we obtain PTA-‘HQ,% + Pfa,x = PTA-‘g(t) or with obvious notations

aI Bz 1 T,%+U,x=b,,

(43)

where T, is tri-diagonal and U, in general complex diagonal when a complex shift is used:

FKl

() I

can be interpreted as static interface modes (shape function), since they give the correct static solution (in discrete sense) due to interface displacements. All vectors in the rest of the Krylov blocks are set to zero on the boundary. For the vectors in the blocks r$ (m = 1,. . .) we used modified Gram-Schmidt orthogonal&&ion so that each vector is mass orthogonal&d against each other vector within the different blocks. If the mass matrix is singular, in general the procedure will still work (for semidefinite inner product spaces, see [32]). For further details, see [26].

72

T, =

recurrence efficients are The

a2 Y3

:

83

.

a3

* * .

. . .

. . .

. . .

. . .

. . .

relations

For proportional damping the supersystem has real eigenmodes. Assuming the general transformation

for generating

MK-‘MB + MK-‘Cii + Mu = MK-‘f(t).

these co-

(45)

S,=(A-‘~‘p,-pj~-p,-,~

(46) J

vjt

I

Vi+

I

sTr

pJ+l’sj-?

whereq,,e,e, ***, are mass orthogonalized Krylov vectors we rewrite the equations of motion for Lanczos method

WI

0’ xi)

r,=(A-‘H)q,-q:-q._,: I

Main structure

j

j

jJ+l~J+l, vi+1

(47) where aj = p;A-lHa~,

(38)

yJ= p,?-, A-‘Hq,

Hence after projection QWK-‘MQSI

(42)

+ QWK-‘CQ$

and

+ Q?MQy = Q%lK-‘f(t)

(39)

bj = e’_ , A- ‘HP,.

(48)

Right and left Lanczos vectors satisfy a biorthogonality condition

or T,g + (CrTj+ flI,)$ + Ijy = Q?<-‘Mf(t)

(4)

with the damping matrix proportional to the mass or stiffness matrix, C = aM + flK, resulting in a symmetric, tridiagonal matrix T,. The recurrence relations for the coefficients in T, might be found in several references, e.g. [1 11,[19], [26] and [31], or they could be derived as special cases of the more general relations valid for the case with non-proportional damping. For non-proportional damping the eigenvectors might be non-real and we work with the Lanczos form [cf. eqn (2511 A-‘H*

+ w = A-‘g(t).

(41)

P,?QJ= Uj.

(49)

In case of symmetry we obtain s, = Hr, and we choose j?J=lrJWrJ11/2and

v,+,=g.

(50) 1

7. COMPUTATIONAL

J

ASPECTS

For large-scale structures the model reduction should be performed in two steps, first on substructure level and then on final assembled superstructure level. In this way it is possible to keep the need for computational power in the final dynamic analyses within reach for commonly available computer

742

B. H&J~BLA~I and L. ERIKSSON

where the Krylov vectors for K,;‘M, appear naturally, which supports our choice of these vectors as a basis that is adapted to the coupling conditions of the substructure. In Sec. 9 it is shown that the efficiency of the Krylov vectors is superior when high accuracy is needed. For the reduction of a whole system or an assembled superstructure we use Lanczos or Amoldi methods. Different Lanczos methods differ in using complete, selective or no orthogonalization. All the methods generate a sequence of tridiagonal matrices with eigenvalues that are progressively better estimates of the extremal eigenvalues of the unreduced system. The extremal parts of the spectrum are well approximated by Lanczos procedures since the Lanczos vectors are generated in the ‘direction’ of the gradient of the Rayleigh quotient (for the real symmetric case) [3 I]. In theory mass orthogonalization need only be done against the two latest Lanczos vectors. In practice severe loss of orthogonality will occur and full or partial orthogonalization against all earlier vectors is needed for every new vector. By using suitable shifts other subsets of the spectrum ui = -K,;‘K,ub + u,,,, (51) can be captured. In order to cope with multiple eigenvalues block Lanczos must be used. This version of Lanczos is capable of capturing multiplicities up to where the interior dynamic part u,~ can be written the block size. Unfortunately these algorithms are difficult to implement due to the sensitivity to round(52) utd = Aid&et off errors. For example, in solving large, real eigenvalue problems, (K - w*M)u = 0, we use block and the with the load fine, = (&M,)(-K;‘K,u*) Lanczos according to Matthies [19] ending up with a operator AU = (K,, - w*M,)-I. If we choose a set of substructure modes w,,, block tridiagonal matrix T,. To get rid of the roundlinearly independent off effects double Gram-Schmidt orthogonalization is @=I,..., n) to construct an approximation to the operator Aid we can write the error in the solution as needed (see [26] for details). For the case of nonproportional damping the matrix H is indefinite and (53) the orthogonalization is made in an indefinite inner u,d - (khappr = Wm 7 n )f,nert product space [32]. Even in this case severe loss of with the ‘error’ operator E(w,, n) = [A, - (A,,&Pr,,X]. orthogonality will occur and full re-orthogonalization against all earlier vectors is needed for every new Using a suitable energy norm we obtain vector. A source of trouble in engineering models is the frequent occurrence of singular mass matrices. A d IIE(wm, n)ll II&t II. (54) singular mass matrix makes the M-orthogonality difficult to obtain with the Gram-Schmidt process on substructure level. It also causes H to become singuIf we, for a given number of modes n, try to minimize lar and makes the H-orthogonality difficult to obtain the corresponding operator norm IIE(w,,,,n)ll we find on superstructure level. Therefore we avoid that that it will be minimized by the set of n eigenmodes situation by replacing all zero masses with small that span the least dominant eigenspace of the masses. Although this will result in a high spectral fixed-fixed interior substructure eigenvalue problem condition number for the M- or H-matrix, the precise of the (shifted) operator Ard[28]. Hence, in this implications of this on the accuracy is not clear. But general sense the set of natural modes is optimal. if the substructure is physically well defined, no rank However one should notice that the special character of the load is not considered. Better efficiency is deficiency will occur in the Krylov blocks and any obtained if we minimize the norm ]]E(w,, n)l,, I], i.e. difficulties are not likely to develop. For structures with a complex eigensystem there we adapt our choice of modes to the load. In fact a are additional sources of trouble. For undamped formal expansion of the solution of eqn (52) gives structures there are efficient tools for determining the number of eigenvalues in an interval. However, for q, = -K;‘K,*u~ + (w2Kn’Mi,)( -K,iK,& non-proportionally damped structures there are no + (02K,7 ‘M,,)‘( - K, IK,bub) + . . . , (55) simple Sturm sequencing technique [31] which allows

resources. For complex eigensystems the assembled system is expanded twice the original size and for ill-conditioned cases the subsequent computations might have to be performed in four-fold precision, resulting in a demand for computer power that is eight times that needed for the corresponding real case. A key feature in our approach is the systematic model reduction of the substructures that makes it possible to obtain an accurate substructure description within a limited frequency interval. Hence, by repeated ‘spectral slicing’ it is possible to obtain overall accurate modal properties with limited computer resources. For the substructure description there is at least one serious, alternative candidate to our choice of Krylov vectors: the set of natural modes for the fixed-fixed interior substructure eigenvalue problem. The motive for choosing Krylov vectors is mainly computational efficiency, measured as the computational effort needed for achieving a certain level of accuracy. Indeed we can rewrite eqn (18) (assuming for simplicity Mlb = 0)

Model reduction methods for analyses of large structures the user to determine if there are any eigenvalues in any region in the complex plane. Special attention should be paid to overcritically damped modes since they can cause trouble in the complex modal superposition if they are not considered properly (i.e. in pairs!). These eigenvalues are located on the negative real axis in the complex plane and occur in pairs since they can be considered as originating from ‘degenerated’ complex conjugate pairs. One root in the pair is located near origin and the other further away in such a way that every largest root in any pair is greater than every smallest root in any pair. A Sturm check test with a real negative shift, using the second order system matrix in eqn (59), gives a number of sign changes that is equal to the number of real eigenvalue pairs in which the shift is located between the roots in the pair [29]. That means that if the shift varies from zero and to very large negative values the number of sign changes will start from zero, pass a maximum and then decrease to zero. This is easily seen from the structure of the diagonal&d system matrix [left-hand side of eqn (31)] where each term in diagonal is the product of two factors containing the roots of a pair. Remark 4. For the time history response of the ‘Lanczos-reduced’ supersystem we can use two alternatives, for example 1. Complete eigen decomposition [cf. eqn (43)] of the system ~~“*TjU~i’z

- l/d&x,

= 0

(56)

(with U;‘RT,U,-‘~z symmetric and complex) followed by superposition of the uncoupled modal equations. This is a safe variant since the unphysical ‘ghost’ eigenvalues can be sorted out. There exist efficient eigenvalue solvers for the tridiagonal systems in eqn (56), see [33]. It should be noted that it is possible to define Ti such that U,~il~,U,~i~z is complex and symmetric even for systems with unsymmetric coefficient matrices [33]. 2. Direct time-stepping scheme applied on the tridiagonal system using the trapezoidal method (Crank-Nicholson), for instance Cr, + Ar,‘2U,)x,+, = (T, - AtDJ,h, +(hn-t-hn+1)/2.

743

not excite the unphysical components in the shifted Lanczos system and the computed result can safely be multipled by the factor e”’ (a is the shift applied on the real part) to give the correct answer. However, in general direct numerical integration of the reduced Lanczos system is not a robust procedure and is not recommended. Remark 5. For unsymmetric systems we can use the unsymmetric Lanczos method with a complex final tridiagonal matrix[33] that can be made symmetric. However, for strongly unsymmet~c systems Amoldi’s method is recommended fl5]. Here T will be on upper Hessenberger form. Amoldi’s method does not seem to be much tested for indefinite systems. 8. ERROR MEASURES

The total error is a sum of the theoretical reduction errors (characteristic for the chosen subspace) and round-off committed in the process of generation of the reduced system (resulting in loss of orthogonality, etc.). For the Lanczos procedures there are many theoretical results avaifable, mainly for error estimators of a priori type. However, these are valid only in exact arithmetic and are not so useful in practice. For a specific problem an a posteriori error measure should be used that takes full advantage of the special features of the problem. For this purpose residuals are employed that arise when the approximate solution is substituted into the unreduced system equations. For substructures we have a further level of truncation errors originating from the substructure description. They are also captured by our residual error measures. It is important to relate the relative residual error to the accuracy of the solution. Here we give some formulas for the accuracy of the solution in terms of the residual errors. The relations are derived in the Appendix. In the following we consider the complete model reduction as a transformation U=QY=(q,,q*r...,~)Y3 where Q contains the orthonormal vectors.

(57)

In theory this is very attractive, since complex arithmetic is avoided at all stages and the time history evolution of the residual (error) vector (due to the Lanczos vector reduction) may be computed to all steps. However, Lanczos systems of different orders contain unphysical eigenvalue components (e.g. due to round-off) with positive real parts. These components will inevitably be triggered by a numerical integrator of the type (57). A remedy that makes it possible to formally perform the numerical integration in (57) is to shift all eigenvalues into the negative complex plane so that all will have negative real parts. Normally, the applied external loads do

(58) Krylov (Lanczos)

Eigenvalue problems The residual of the fully expanded system corresponding to the computed eigenmode v, = Qy, and the eigenvalue xi is obtained from AWQY,

+

&CQy, + KQy, = r.

(5%

We define the a posteriori (physical) relative residual error as

II41

III
B. H~C~BLALIand L. ERIKSSON

744

Fig. I. Damped cantilever beam with four substructures.

For the undamped

estimator is the time varying residual

case we obtain

r(t)=Mii+Crk+Ku-f(t)

(61) where K(M) is the condition number of the mass matrix M. For the damped case (eigenvalue x and eigenmode wc) we obtain (using appropriate scalar products and R as the expanded version of r)

G&&%&x max(L&), (62) I

I

where K(H) is the condition number of the matrix H. For the mode shapes in the undamped case we obtain error estimates of the type (Iv, -

P,V,IIM <

computed on the fully expanded system. The formulas (61)-(64) indicate that the relative residual error in (60) could be used as an error estimator in all cases. 9.

loo

10-Z

&i@

“IIKQT, , mm I

(63)

d--(3;

+ I w,

NUMERICAL EXPERIMENTS

The model reduction procedures outlined in this paper have been implemented in a derivative program of SAP IV [35]. The simple cantilever beam in Fig. 1 is made up of four substructures (SUBl-SUB4) with 156DOF in each. We use dynamic substructuring with Krylov vectors and Lanczos reduction on supersystem level. Complex eigenvalues and eigenmodes are computed in two intervals for the frequency (imaginary) part of the eigenvalues, near 0 Hz and near 6850 Hz (which is used as an imaginary shift on supersystem level and as a real shift on substructure level). The diagrams in Figs 2 and 3 show the influence of the number of

JrM-‘r llrll

(64

I

(Pi is a projection operator onto the eigenspace of of). Equation (63) indicates that it is important that the approximate eigenvalues have good accuracy when they are closely spaced.

$

10-j

1

10-4

e

10-S

3 i!

10-6 lo-7

Time history analysis For a Lanczos system there is special a posteriori error estimators, but they do not in general consider the effect of round-off (see, e.g. [15]). When substructures are used a safe but computationally expensive

lo-9

1

2

3

4

5

6

number of blocks Fig. 2. Relative residual error near 0 Hz.

7

745

Model reduction methods for analyses of large structures

and the eigenvalue problem is also solved using an EISPACK subroutine [34]. The size of the problem sticiently small (6OODGF, expanded was 12OODGF) for solving it on a workstation. The computations were done in fourfold precision and the results are considered as the ‘correct’ ones. The relative error in the eigenvalues is defined as

number of blocks Fig. 3. Relative residual error near 6850 Hz.

Krylov blocks per substructure on the achievable accuracy, measured by the relative residual error in eqn (60). Each Krylov block contains 12 vectors (nb = 12, ni = 144). When increasing the number of blocks from one (corresponding to Guyan reduction [l]) to six in Fig. 2, the relative residual errors for the five modes closest to zero pass through a minimum when three blocks are used. The achievable accuracy is typically half the machine precision. In Fig. 3 the same study is made for the five modes ‘closest’ to the shift @ = i6850). For seven blocks the results were unusable with the numerical precision employed (double precision). The results in the diagrams can be interpreted in such a way that there is an optimum number of Krylov blocks above which the enhancement of the approximation properties is destroyed by the increased round-off The reliability of the residual error estimator according to eqn (60) is demonstrated in Fig. 4. The same damped beam structure is analysed as in Fig. 3 101

where 5 is computed by the EISPACK subroutine and 1, by the same technique as in the last problem. The diagram shows that the relative residual error, eqn (60), is somewhat conservative for the modes considered. However, one should remember that the residual error considers errors both in the eigenmodes and eigenvalues. The level of the errors are over-all low for the six modes around the shift frequency (6850 Hz). In Fig. 5 we compare the use of Krylov vectors versus eigenmodes for the description of the interior motion of a substructure in the beam in Fig. 1. The relative residual errors of the three ‘lowest’ eigenmodes are plotted. The superior performance of the Krylov vectors is clearly demonstrated for this structure. The optimum number of Krylov vectors is 36 (three blocks). The same number of eigenmodes gives quite unsatisfactory results (relative error of about 0.1). Only when all interior eigenmodes (144 modes) are considered the error comes down to the level obtained by using Krylov vectors. 10. ENGINEERING

PROBLEM

As a representative engineering problem we have a chosen a bogie for traction vehicles, consisting of a bogie frame, wheels with wheel shafts and brake disks, primary suspension and vertical dampers (see Fig. 6). The complete model contains 78,000 DOF distributed on three substructures, bogie frame

1 100,

i;

I

estimated error

10-7 -

1

2

3

4

5

6

10-s -

mode no Fig. 4. Comparison between relative residual error and real relative error for the modes near the shift frequency 6850 Hz.

number of modes Fig. 5. Krylov vectors versus eigenmodes.

B. H~G~LNI and L. E&KSSON

746

Fig. 6. Bogie for traction vehicles. Table 1. List of substructure data Unrcdttced DOF

Krylov blocks

RedUlXd DOF

Bogie frame Wheel pair 1 Wheel pair 2

29,000 24,000 24,000

3 4 4

270 144 144

Total

78.000

substructure

1400

(29,000 DOF) and two wheel pairs (24,OOODOF each). We compute 20 complex eigenfrequencies (damped) in a frequency window around 300 Hz using first complex Lanczos reduction on the whole system and then on a supersystem (1400 DOF) consisting of reduced and assembled substructures (see Table 1).

The problem was run on a workstation IBM RISC/6000, model 550, equipped with 256 Mbyte core memory and with an average performance of about 25 MFLOPS. The resulting accuracy for the two alternatives is shown in Table 2, where all digits for the listed eigenfrequencies are the same, only the relative residual errors are somewhat larger for the substructuring alternative. However, when using the full model the CPU time was around 35 hr. When substructuring was used the CPU time was reduced to about 60 min (reduction factor of 35!) and the disc space needed was reduced by a factor of five. The bogie frame substructure is reduced by using shifted Krylov blocks (shift frequency 300 Hz). The accuracy of the obtained substructure description is studied by computing the frequency response due to

Table 2. List of residual errors for modes around 300 Hz

Mode number

I 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

Angular frequency (rad/sec)

Frequency (Hz) 2.91738E + 2.91846E + 2.93374E + 2.937758 + 2.94224E + 294526E + 294733E + 2.95039E + 2.95515E + 2.95757E + 2.%379E + 2.96857E + 2.98136E + 2.98539E + 2.98867E + 299335E + 299416E + 299877E + 3.00374E + 3.01306E +

02 02 02 02 02 02 02 02 02 02 02 02 02 02 02 02 02 02 02 02

1.83304E + 1.83373E + 1.84332E + 1.84584E + 1.84866E + 1.85056E + 1.85186E + 1.85378E + 1.85677E + 1.85829E + 1.86220E + 1.8652lE + 1.87324E + I .87577E + 1.87784E + 1.88077E + 1.88128E + 1.88418E + 1.8873lE + 1.89316E +

03 03 03 03 03 03 03 03 03 03 03 03 03 03 03 03 03 03 03 03

Real part

Eq. model damping

- 1.85977E + 00 -2.54666E + 00 -7.5954lE - 02 -1.691288-01 - 1.07473E - 01 -3.7694lE - 01 - 164829E - 01 -8.852788 - 02 -2.37179E - 01 -2.81529E - 01 -2.81146E-01 -2.11162E-01 -5.04804E - 01 - 1.21367E + 00 -1.39146E+OO -1.381llE+OO - 1.12475E + 00 - 1.12777E + 00 -2.61503E-01 -9.67194E - 01

1.01458E - 03 1.38879E - 03 4.12050E - 05 9.1626lE-05 5.81355E - 05 2.0369OE - 04 890072E - 05 4.7755lE - 05 1.27737E - 04 1.51499E - 04 1.50975E -04 l.l32llE-04 2.6948lE - 04 6.47022E - 04 74099lE-04 7.343328-04 5.97866E - 04 5.98544E - 04 1.38559E-04 5.108898 -04

Full model residual error 2.21143E 148258E 2.10777E 1.83983E 2.14194E 2.52022E 2.18824E 2.65284E 3.78875E 5.00528E 2.32702E 4.0708lE 6.95922E 4.5149lE 8.22345E 1.20920E 144292E 597329E 3.77049351.723108 -

06 06 06 06 06 06 06 06 06 06 06 06 07 07 07 06 06 08 07 06

Substr. model residual error 4.86023E 2.62797E 4.0156OE 2.1394OE 1.09920E 3.24438E 2.22882E 1.00023E 5.57848E 1.562248 1.087638 1.42564E 1.59107E 9.24644E 9.39437E 1.0704OE360292E 2.09057E 9.85705E 4.68287E -

05 05 05 05 05 06 06 06 07 07 07 06 06 06 06 05 05 05 05 05

Model reduction methods for analyses of large structures 10-Z

workstations. Currently the procedures are not fully automated: the user has to choose the shiRs or optimal number of Krylov blocks for the substructures. In practice some engineering judgement is needed to choose, e.g. the (complex) shifts elWently. A good rule is to base the choice on a series of Sturm checks for the corresponding real CaseyFuture activities will include further studies on unsymmetric systems.

shit frequency:300 Hz

10-3 -

141

10-4 -

10-8 Acknowledgement-The authors would like to thank ASEA Brown Boveri for the permission to publish this work.

10-9 -

RMERENCES

210

280

290

300

310

320

330

frequency (Hz) Fig. 7. Accuracy of a shifted bogie substructure frequency window around 3OOHz.

in a

an applied load to the substructure (properly constrained) and then compare it with the response of the unreduced substructure model. The frequency response is computed by solving equation systems for each frequency in both cases. The relative response error of the reduced model is acceptable in a frequency window of about 30 Hz around the shiA frequency as illustrated in Fig. 7. In this case, exceptionally good accuracy is obtained very near the shifting frequency.

Il. CONCLUDING REMARKS In this paper we have presented model reduction methods for structural dynamics that have shown to be accurate and computationally very efficient for large scale problems. The methods can be used systematically to capture the essential parts of the subspace that contains the solution of the unreduced problem. This is obtained by adaption to external loads, mass distribution and coupling conditions (for substructures) when the reduced basis functions are selected. For broad banded excitations frequency slicing may be utilized both on substructure and supersystem level. A general (and physically motivated) procedure for selecting basis vector (Krylov vectors) for substructures with (non-) proportional damping is introduced. In general the a posteriori residual errors for the fully expanded system are shown to be reliable estimates for controlling the reduction and round-off errors. A very competitive technique is obtained for the generation of complex eigenspaces of large structures, if shifted Krylov blocks on substructure levels are combined with the use of the corresponding complex (imaginary) shift on superstructure level. So far the introduced methods have shown good performance on complex engineering problems and it __ __ seems likely that new classes of large problems are within reach for efficient analysis on powerful

1. R. J. Guyan, Reduction of stiff&

and mass matrices. AIAA Jnl3, 380 (1965). 2. W. C. Hurty, Dynamic analysis of structural systems using component modes. AIAA Jnl3, 678-685 (1965). 3. R. R. Craig and M. C. C. Bampton, Coupling of substructures for dynamic analysis. AIAA Jnl 6, 1313-1319 (1968). R. L. Goldman, Vibration analysis by dynamic partitioning. AIAA Jnl7, 1152-l 154 (1969). R. H. b&Neal, A hybrid method of component mode synthesis. Cornput. Struct. 1, 581-601 (1971). R. M. Hintz, Analytical methods in component modal synthesis. AIAA Jnl 13, 1007-1016 (1975). E. L. Wilson, M. Yan and J. M. Dickens, Dynamic analysis by direct superposition of Ritz vectors. Earthquuke Engng Struct. Dyn. 10, 813-821 (1982). 8. E. L. Wilson and E. P. Bayo, Use of special Ritz vectors in dynamic substructure analysis. ASCE, J. Struct. Engng 112, 1944-1954 (1986). 9. K-J. Joo and E. Wilson, An adaptive linite element technique for structure dynamic analysis. Comput. Struct. 30, 1319-1339 (1988). 10. Kuan-Jung Joo, E. L. Wilson and P. Leger, Ritz vectors and generation criteria for mode superposition analysis. Earthquake &gng Struct. Dyn. 18, 149-167 (1989). 11. B. Nour-Gmid and R. W. Clough, Dynamic analyses of structures using Lanczos w-ordinates. Eurthquake E@ng Struct. Dyn. 12, 565-577 (1984). 12. B. Nour-Gmid and R. W. Clough, Block Lanczos method for dynamic analyses of st&tures. Eurthquuke Enmuz Struct. Dvn. 13. 271-275 (1985). 13. B. “N&u-Gmid, -Lanc&s algoritdm fdr transient heat conduction analysis. Int. J. Numer. Meth. Eitgng 24, 251-262 (1987). 14. B. Nour-Gmid and M. E. Regelbrugge, Lanczos method for dynamic analysis of damped structural systems. Earth&&e Engng- Strut. Dyn: 1% 1091-l 104 (1989). 15. B. Nour-Gmid. W. S. Dunbar and A. D. Woodbury. Lanczos and krnoldi methods for the solution of convection-diffusion equations. Comp. Meth. Appl. Mech. Enana %s. 75-95 (1991). 16. H. Chen an; R.. L. Taylor, Solution of eigenproblems for damped structural systems by the Lanczos algorithm. Comput. Struct. 30, 151-161 (1988). 17. H. Chen and R. L. Taylor, Solution of viscously damped linear systems using a set of load-dependent vectors. Earthquake Enmra Struct. Dyn. 19, 653-665 (1990). - 18. A. Ibrahimbegovic, H. Chen, E. L. Wilson and R. L. Taylor, Ritz method for dynamic analysis of large discrete systems with non-proportional damping. Earthquake Engng Strut. Dyn. 19, 887-889 (1990). 19. H. G. Matthies, A subspace Lanczos method for the generalized symmetric eigenproblem. Comput. Struct.

21, 319-325 (1985).

B. H;iccBun

748

20. R. R. Craig, A study of modal coupling methods. In Proceedings of the International Conference on Spacecraft Structures, European Space Agency, Noordwijk. The Netherlands (1985). 21. A. L. Hale, The block-Krylov method in component synthesis: an alternative to component modes. AIAA Paper 86-0925 (1986). 22. R. R. Craig and Z. Ni, Component mode synthesis for model order reduction of nonclassically-damped systems. In AIAA Guidance, Navigation and Control Conference, Monerey, CA (1987). 23. R. R. Craig and A. L. Hale, Block-Krylov component synthesis method for structural model reduction. AIAA J. Guiaknce Control Dyn. 11, 562-570 (1988). 24. S. Park, I. Park and F. Ma, Decoupiing approximation of nonclassically damped structures. AIAA Jni 30, 2348-2351 (1992). 25. K.-W. Min, T. Igusa and J. D. Achenbach, Frequency window method for strongly coupled and multiply connected structural systems: multiple-mode windows. J. Appt. Me&. 59, 244-252 (1992). 26. L. Eriksson, Modei suction methods for structu~l engineering dynamics. Licentiate thesis, KTH, Stockholm (1993). 27. K. J. Bathe, Finite Element Procedures in Engineering Analysb. Prentice-Hall, Englewood Cliffs, NJ (1982). 28. J. Babuska, M. Prager and E. Vitasek, Numerical Processes in Dtyerential Equations. Wiley Interscience, Prague 09% 29. P. Lancaster, ~~~-~trices and Vibrating Systems. Pergamon Press (1966). 30. B. N. Parlett, The Symmetric Eigenvalue Problem. Prentice-Hall, Englewobd Cliffs, NJ-(1980). 31. G. H. Golub and C. F. van Loan. Matix Cowmutations. Johns Hopkins University Press,’ Baltimore fi989). 32. J. Bognar, &i&nite Inner Product Spaces. Springer, New York (1974). 33. J. K. Cullum and R. A. WilloulIhby, Lancros Algorithms for Large Symmetric Eigenvaiue Computations, Volume 1, Theory; Volume 2, Programs; Progress in Scientific Computing, Volumes 3 and 4. (Edited by S. Abarbanel, R. Glowinski, G. Golub, P. Henrici and H. 0. Kreiss). Birkhauser, Base1 (1985). 34. B. S. Garbow, J. M. Boyle, J. J. Dongarra and C. 8. Moler (Bds), EISPACK, matrix eigensystem routines, Lecture Notes in Computer Science, 6 and 51. Springer, New York (1976/s). 35. K. J. Bathe, E. L. Wilson and F. E. Peterson, SAP IV-a structral analysis program for static and dynamic response of linear systems. Report EERC 73-l 1, College of Engineering, University of California, Berkeley, June (1973), revised April (1974). 36. G. Still, Computable bounds for eigenvalues and eigenfunctions of elliptic differential operators. Numer. Math. 54, 201-223 (1988).

and L.

bIKssoN

may be written (g - tifm)s, = 0

(f,, ae,) = a,, (Kronecker delta). (A2)

gandRarcsym~t~cstiffnessandmassmatricrs(m and e,, 9, am eigenvectors belonging to V,,,. We obtain the Rayleigh quotient

where v, is the expanded appro~mation

xm)

to $, . Introduce the (A4)

Assume the expansion ”

=,r;se,.

VP

Q - @M)v, = c c,(c$ - cj;)M+, 1-I

(As)

(A6)

and

I =

k-l

min (e$ - @)‘cv,, Mv,).

(A7)

min Iwj - r# I 6

(A8)

f

Thus

I

where r=(K-#M)v,,

(AS)

Further using the expansion (AlO) g&S

(Al 1) For the approximate eigenvector we obtain

APPENDIX:

DERIVATION

Q - &;‘M)v, = (K - ~~M)(v, - P,v& + (K - r$M)P,v,

OF ERROR BOUNDS

(A121

First we derive the error bounds for real eigenvaiues PO]. In the full expanded space (V, , n = total number of degrees of freedom) we have (If - o:M)$&, = 0

d = (r, M-‘r)

(IL,,M$,) = a,, (Kronecker delta), (Al)

where (, .) denotes the usual Euclidian scalar product of two vectors and K and M are symmetric stiffness and mass matrices (n x n) and $,, $, are eigenvectors belonging to V,,. In the truncated subspace V, (m = retained number of degrees of freedom, in general m < n) the projected problem

2 min (o$ - rj:)2iivz - P,v& J*l or

Ilv,-P,v,n,c

f

(Al3)

(A14)

Model reduction methods for analyses of farge structures (P, is an projection operator onto the eigenspace of wf). An even better estimate can be derived [36]

llrl,- V,Ilu<6,(1 +cp, where 6, is the right-hand side of eqn (A14) For the complex eigensystem (symmetric matrices) analogous formulas are obtained if ‘complex’ scalar products are used. The final inequalities on the right-hand

749

sides of the formulas (61x63) are easily derived by considering the definition of a matrix condition number and annlvina the followinn relations (Al3 @v, Kv) llMll d (Kv, M-‘Kv) Q (Kv, Kv)llM-‘Il.

(A16)