An accurate modal method for computing response to periodic excitation

PII: sow-7949(%)00367-7

CompuCrs & Sfrucrures Vol. 63, No. 3, pp. 625-631. 1997 CQ 1997 Published by Elsevier Science Ltd Printed in Great Britain. All rights reserved w4_5-7949/97 Sl7.00 + 0.00

AN ACCURATE MODAL METHOD FOR COMPUTING RESPONSE TO PERIODIC EXCITATION Cheng Huang, Zhong-Sheng Liu and Su-Huan Chen Department of Mechanics, Jilin University of Technology, Changchun 130022, People’s Republic of China (Received

31 March 1995)

Abstract-The modal superposition method is one of the methods often used to extract the response to periodic excitations. It cannot, however, give an exact solution, and may sometimes lead to meaningless results if the truncating of modes is serious. Fast and accurate computation of the contribution of the truncated higher frequency both of lower and higher frequency modes to the response, without using the unavailable truncated modes, is the objective of this paper. An explicit expression on the contribution from these unavailable modes to the solution is presented in this paper in order to improve the solution accuracy. Using this explicit expression, the response to periodic excitations can be computed accurately and conveniently with little extra computational cost. Numerical examples are given to illustrate the

effectiveness of the present method. 0 1997 Published by Elsevier Science Ltd. All rights reserved.

II. INTRODUCIION

Determination of the response of a complex structure to periodic excitations, p(t), is frequently encountered in structural analysis. Usually, direct integration method and/or mode superposition method are used to obtain the response. The number of operations required in the direct integration method is directly proportional to the number of time steps used in the analysis. Therefore, in general, the use of direct integration method such as Wilson 0 method and Newmark method can be expected to be effective when the response for a relatively short duration (i.e. for a few time steps) is required. However, if the integration must be carried out for many time steps, it may be very expensive. In fact, it is this aspect that makes mode superposition analysis of some structures feasible with regard to cost, whereas direct integration would be prohibitively expensive [ 11. Using the mode superposition method, to avoid using direct integration method for a single-degreefreedom equation of motion, the periodic excitation, p(t), may be represented by a convergent series of harmonic functions whose frequencies are integral multiples of a certain fundamental frequency wo. The frequencies representing integral multiples of the fundamental frequency are called harmonics, with the fundamental frequency being the first harmonic. Such series of harmonic functions are known as Fourier series. Since calculation of the transient structural response for such large-scale systems is computational expensive, computational methods that can significantly reduce the problem size and computational cost and still retain the solution accuracy are highly desirable. This fact has been widely recognized

and some techniques for reducing the degrees of freedom have been developed. In a recent literature, Noor [2] presented an excellent review about recent advances and applications of the reducing degrees of freedom which are referred to as reduced-basis techniques. The reduced-basis method that uses the natural vibration modes of the structures, referred to as modal methods, is widely used. Two of the most widely used modal methods are the mode-displacement method [MDM] and modal acceleration method [MAM] [3]. Recently, a higher-order modal method, entitled the force derivative method [FDM] [4], has been developed to increase the convergence rate of the MDM and MAM. It was found [5] that the FDM obtains converged solutions using the smallest numer of basis vectors and the smallest amount of CPU time of the three methods MDM, MAM and FDM. Liu et al. [6] presented, very recently, a very effective method, entitled a hybrid method, to improve the solution accuracy of the MDM and MAM. This hybrid method expresses the solution into two parts: the solution from MDM plus the contribution of the higher-frequency truncated modes to the response. In this hybrid method, the contribution of the higher-frequency truncated modes to the solution is, exactly and explicitly, expressed in the form of a convergent power series whose items can be computed in succession by solving static responses of the structure to different static forces. In Refs [7, 81, Liu et al. used this hybrid method to compute the eigenvector derivatives in structural dynamics and the numerical examples therein indicate that this hybrid method converges to 625

Cheng Huang et al.

626

an accurate solution using fewer modes than MDM with little extra computational effort. However, this hybrid method developed by Liu et al. in Refs [6-S] does not contain the situation that both of the lower-frequency vibration modes and the higher-frequency vibration modes are truncated with the middle-frequency vibration modes available. Such situations are often encountered for example, in the dynamic analysis of coupled acoustic-structural systems [9]. This paper tries to further extend this hybrid method presented by Liu et al. to include the situation of truncating both of the lower-frequency and higher-frequency modes. In this extended hybrid method, the contribution from the unavailable lower-frequency modes and higher-frequency modes to the dynamic response is expressed exactly and explicitly as a fast convergent power series. Each term of this power series can be easily evaluated in succession through solving static problems with different force distributions. This paper uses the hybrid method presented by Liu et al. and the extended hybrid method to compute the response to periodic excitations, and gives a qualitative discussion about the advantage of these methods. A numerical example is given to illustrate the effectiveness of the presented method.

2. TECHNICAL

BACKGROUNDS

Consider a multi-degree-of-freedom structural system forced by a periodic excitation Fp(t) governed by MX(t) + KX(t) = Fp(t)

(1)

where KoRN x ’ and ME RN x N are real symmetric stiffness and mass matrices, respectively, FoRN is a force distribution column vector, X(t)E RN is response column vector, a dot denotes differentiation with respect to time. p(t) is a periodic function with period To, that is, p(t + To) =p(t). p(t) may be separated into its harmonic components by means of a Fourier series expansion, which may be defined as

p(t) = a0+ 2 a, cos(m~,t) + 5 b, sin(noot) (2) “=I I=, where w. = 2x/T, is the fundamental

frequency, and

a,, and b, are the coefficients of the nth harmonic. The coefficients a~, a,, and b,, are related to p(t) by the

equations I + To p(t) dt = average value of p(t)

(3)

p(t)cos(nwot)

s 1 +

b,=+

dt, n # 0

Cl r

(4)

To

p(t)sin(nw,,t) dt

(3

where z is an arbitrary time [IO]. Now we have to determine the response of the system to harmonic excitation. Consider the multidegree-of-freedom structural system force by a harmonic excitation F sin(nwt) governed by M%(t) + KX(t) = F sin(noot).

(6)

Setting w = nwo and X(t) = X sin wt, one can obtain from eqn (6), X = (K - o*M)-‘F

(7)

where XER represent the response distribution column vector. From eqn (7) one can obtain the solution X by solving a set of linear equations with coefficient matrix equal to (K - IM), i.e. X = (K - w2M)-‘F. When many o values are of concern, the linear equation solution may take much computational effort; however, modal summation method may be less computationally burdensome when many o values are of concern or the quantitative effects of w on the response need to be evaluated. Such cases where many w values are of concern are often encountered in dynamic analysis of structures. For example, periodic excitation forces are often expressed as its Fourier series representation. Therefore, the modal summation method has been widely used in many structural analysis programs such as NASSTRAN. When the modal summation method is used, eqn (7) needs to be expressed using its modes as

where A, and q~,,o’= 1, 2, . . , N) are the jth eigenvalue (the power ofjth natural frequency) and the jth eigenvector (jth vibration mode) associated with 1, of the following eigenproblem, (K - l,M)+, = 0

(9)

$;M4, = 1

(10)

where 11is in order of ascent, i.e. I, < A2< . . . < AN. For the dynamic analysis of the complex structures by finite element method, there may be thousands of degrees of freedoms. Therefore, in most cases, only a few of the lower modes are calculated while the other higher modes are truncated in order to save

Computing response to periodic excitation

computational effort required for obtaining eigen pairs. Assuming that the first L(L<
(11) Equation (11) is called modal displacement method (MDM). When L <
4P4.

(12)

The first term in the above equation is the pseudo static response, while the second term gives the method, namely, the mode acceleration method (MAM). The presence of 1j in the denominator of this term improves the convergence of the method as compared to the MDM.

3. BRIEF REVIIEW OF THE HYBRID METHOD 14

Separating the right side of eqn (8) in terms of the lower available modes and the higher unavailable modes, eqn (8) ca.n be rewritten as

x = x, +

XH

(13)

where XLo R denotes the contribution from the available lower modes to the solution X,

627

yj can be obtained procedure,

using the following

iterative

yl = K-IF yj = K-‘(My,_

1).

This hybrid method has been shown to converte the exact solution using very small number of vibration modes than MDM even if only two or three first items of the power series given by eqn (16) are retained with other terms truncated. However, this hybrid method presented in Ref. [6] does not contain the situation that both of the lower-frequency modes and the higher-frequency modes are truncated with the middle-frequency modes retained. This paper tries to extend this hybrid method to include this situation.

4. EXTENDED HYBRID METHOD

Consider the situation that the first LT lower-frequency modes and the last N - HT + 1 higher-frequency modes are truncated with the middle HT - LT - 1 modes retained. In terms of the unavailable first LT lower-frequency modes, the last (N - HT + 1) unavailable higher-frequency modes, and the available middle-frequency modes, eqn (8) can be rewritten as (19) where XLTrepresents the contribution from the first L unavailable truncated modes to X,

(14)

from and XH represents the contribution truncated higher modes to the solution X,

the

In terms of the hybrid method [6], XHcan be obtained without using the higher truncated modes,

XI,=

5 K-'(JJj-Zj)

(16)

where XR represents the contribution available modes to X,

from

where XHT represents the contribution unavailable higher truncated modes,

from the

N 44’

xHT=j_;Ts-‘02F.

the

(22)

j-l

where (17)

In order to compute (X,r + XHT) accurately without using those unavailable modes, let us try to express (XL, + XHr) into a fast convergent power

Cheng Huang et al.

628

series. For the purpose of succinctness, eqns (20) and (22) can be written in matrix form

expressed as &,I

XLr = &r(ALr - w~I)-‘@;~F

(23)

XHr = &r&T.

(24)


16 - W21< I&r+ I - QI - ~‘1)-‘@;rF

IQ - w21< l&r+ I - 61. where OLr is an unavailable mode matrix, Arr is an unavailable diagonal matrix with the jth eigenvalue as its jth diagonal, @“r is an unavailable mode matrix, AHP is an unavailable diagonal matrix with the (j + HT - 1)th eigenvalue as its jth diagonal,

(32)

Using the above power series expression and substituting eqns (29) and (30) into eqns (23) and (24), respectively, one has T, (33)

(25) (26) xHT= F (w’ - u)~U~F k=O

(27) &r = diag&,

LT+ I, . . ,&I.

(28)

where

In fact, these two diagonal matrices (ALT- 0*1)-’ and (A,, - 021)-’ can be expressed into the following power series:

(Au -

w*1)-’ =

kzo (w’ -

(34)

TA.= @L&P,, - ~l)-‘-‘@;~,

(k = 0, 1,2, . . .)

(35)

Uk = u%,(A”T - ur)-“- ‘@jrfT, u)~(A~~- uI)-~ - ’ (k=0,1,2

I... ).

(36)

(29) Therefore, XLr + XHTcan be written as, (A,, - d1)-’

= f

(w’ - u)yA.,

- al)-“-’

XLT + XHT= f (a’ - u)‘(Tk + G)F. k=O

k=O

(37)

(30) where u E R is a shift value which is required to satisfy the following conditions in order to make the series

It can be proven (see Appendix) that To + U, = (K - uM)-’ - ‘@,&lR- aI)-‘@;

convergent

T, + u, = (K - uM)-‘(M(K

0’=1,2 ,...,

LT,HT,HT+l,...,

Tk + u, = (K - uM)-‘[M(K \

N).

- uM)-‘][M(K

- aM)-‘)

- @,Q(A,Q - al)-*@:

(31)

- uM)-‘1 Y t

(38)

. [M(K - uM)-‘I - @R(& - u~)-~-‘@R ,

(39)

(40)

from eqn (40) we can obtain Since the available eigenvalues, generally speaking, contain the value of o2 within their range, i.e. lZ~r< w2 < L, the above conditions can be further

X‘, + X/r = $0 (a2 - u)%k - rk)

(41)

where

yk = (K - uM)-‘[M(K - uM)-‘][M(K

- uM)-‘1 . . . [M(K - uM)]F

(42)

Computing response to

1

4

7

10 Fig.

(43) In practical computations, yk can be obtained using the following iterative procedure, y,, ==(K - aM)-‘F yk = (K - aM)-‘Myk _ ,

w

q

and only the first few terms in eqn (41) are to be calculated because of its fast convergent rate. It should be noted that the shift value, cr, is required to satisfy the condition, O#Aj,

(j=LT+l,LT+2

,...,

periodic

excitation

13

629

16

19

22

1.

For the purpose of comparison, taking the 27th DOF (node 11 along z direction) as example, the response of this truss to this excitation is computed using the three methods. Exact solution (MDM without modal truncation), MDM (with modal truncation) and the present method, respectively. Since MAM cannot be used directly to compute the response in the situation that the lower-frequency mode are also truncated, the solution from MAM is not included herein. In fact, if we let Q = 0.0 and only the first term of the power series be retained, this extended hybrid method reduces to MAM. As a measure of the error of the approximate methods, the error, denoted as e, is defined as

HT-1) e = (9 - 3)‘(ti - $)/4’+

in order to make matrix (K - aM) non-singular. 5. NJMERICAL

EXAMPLE

Consider a truss structure, shown in Fig. 1, with its given by: elasticity moduparameters lus = 6.9 x 1O’OPa, cross-section area of each rod: 4 x lO-4 mr, mass density: 7.8 x lo3 kg m-‘. The length between node 1 and node 22 is 14 m, between node 1 and node 3 is 3 m and the height of the structure is 2 m. The x axis is along the direction from node 1 to node 2.2, the y axis is along the direction from node 1 to node 3. It is fixed at node 1,3,22 and 23. The period excitation force, illustrated in Fig. 2, with its amplitude equal to 1.0 x 10SN acts at node 11 along z direction, node 5 along y direction and node 17 along - y direction.

where r#~is exact solution, $ is approximate solution. On the first situation, the first 3 modes (i.e. L = 3) are retained while the others are truncated. The first three modes retained are AL = diag(4724,8445, 15,151). For the excitation frequency values w = 125.6637, the response is calculated using the exact solution (MDM without modal truncation), MDM (with modal truncation) and the present method with the first 3 terms retained in eqn (16), respectively. Three results are shown in Fig. 3. The error of the MDM is 11.52759%, since the error of the present method is very small (e = 0.022586%), the response curves using exact solution method and the present method are not much different.

0.01

0.02

0.03

0.04

Period T = 0.05 (se@ Fig. 2.

Fig. 3.

0.0s

Cheng Huang er al.

630

4E-007

5 2E-007 2 g !2

0

k%-2E-007 -4E-007

o

0.005 0.010.0150.020.0250.030.035

Period T = 0.034 (set) Fig. 4.

On the second situation, the first 3 modes (1-3rd) and the last 36 modes (22257th) are truncated with the vibration modes from the 4-21st available. In using this extended hybrid method, a shift value u must be carefully chosen so that the condition set by eqn (31) is satisfied. It should be emphasized that the value of CJaffects significantly the convergent rate of the solution. Figure 4 shows the response results computed by three methods in case that ~7satisfies the above condition, 0 = 16,745.878. The curve in Fig. 4 confirms that the solution accuracy of the present method is much higher than that of MDM which is not only inaccurate sometimes but also misleading. The relative error of the solution from MDM is 5.501562%; however, the relative error of the solution from this extend hybrid method, with one to six terms of the power series retained, is 1.576862%.

unavailable modes to response. The contribution is given by a convergent series and the terms of series can be obtained in simple iterative form. Thus it can give accurate solution with little extra computational cost. In fact, the idea that takes advantage of the power series can be used not only in the problem this paper deals with. One of such applications addressed in Refs [ll-141 is to use this idea to improve the accuracy of substructure synthesis. The authors of this paper are trying to utilize this idea to reduce the size of large structures for the purpose of vibration control and the results are going to be published in forthcoming paper. Acknoi&dgement-This project is supported by National Natural Science Foundation of China. The authors are grateful to Dr Da-tong Song and Dr Wan-zhi Han for their helpful discussions during the preparation of this paper.

REFERENCES 1. K. J. Bathe and E. L. Wilson, Numerical Methods in Finite Elemenf Analysis, pp. 308-337. Prentice-Hall,

Englewood Cliffs, NJ (1976). 2. A. K. Noor, Recent advances and applications of reduction methods. ASME Appl. Mech. Rev., 47 (5), 125-146 (1994).

3. D. Wiliams, Dynamic loads in aeroplanes under given impulsive loads with particular reference to loading and gust loads on a large flying boat. Great Britain RAE Reports SME 3309 and 3316 (1945). 4. C. J. Camarda, R. T. Haftka and Riley, An evaluation of higher-order modal methods for calculation transient structural response. Comput. Strut?. 27 (l), 89-l (1987). 5. D. M. McGowanm and S. W. Bostic, Comparison of advanced reduced-basis methods for transient structural analys. AIAA J. 31 (9) 1712-1719 (1993). 6. Zhong-sheng Liu, Su-huan Chen, Wen-tong Liu and You-qun Zhao, An accurate modal methods for computing response to harmonic excitations Inr. J. Anal. Exp. Modal Anal. 9 (I), 1-13 (1994).

6.CONCLUDING REMARKS

An accurate modal method for computing response to periodic excitation is proposed in this paper: First, for the situations that the first lower modes truncated in order to save computational cost (e.g. L = 50, N = 3000). The advantage of the present method is that it can enhance drastically the accuracy with little extra computation effort. After the periodic excitation being expanded into the Fourier series expansion, each harmonic component can use the same zk, yk obtained by using eqns (17) and (18). If the [L][U] factorization of the stiffness matrix K has been obtained during solving the eigenproblem, the extra effort required for solving eqn (18) becomes much less. Second, on the situations that the lower-frequency and higher-frequency modes are both truncated and only the middle-frequency modes are kept, an effective method to deal with such modal truncation problem is presented in this paper. This method can compute the contribution from these

I. Zhong-sheng, Liu, Su-huan Chen, Min Yu and You-qun Zhao, Contribution of the truncated modes to eigenvector derivatives. AIAA J. 32 (7), ISSI-1553 (1994).

8. Zhong-sheng Liu and Su-huan Chen, An accurate method for computing Eigenvector derivatives for free-free structures. Comput. Struct. 52 (6), 1 (1994). 9. Zhong-sheng Liu, Su-huan Chen and You-qun Zhao, Computing eigenvector derivatives in structural dynamics. Acta Mechanica Solida Sinica 16 (3), 291-300 (1993).

Roy R. Craig, Jr, Structural Dynamics. An Iniroduction to Computer Methods. Wiley, New York pp. 163-167 (1981). Y. K. Cheung and A. Y. T. Leung, Finite Element Methods in Dynamics, p. 224. Science Press, Beijing (1991). Hu Haichang, Vibration Theory of Multi-degree-of-freedom Structures. Science Press, Beijing, pp. 103-110 (1987). 13. Hagiwara Ichiro and Z. D. Ma, Development of new mode-superposition technique for truncating lowerand/or higher-frequency modes. JSME Int. Series C 37 (I), 14-20 (1994).

14. Song Da-tong, Han Wan-Zhi, Chen Su-huan and Qiu Zhi-ping, Simplified calculation of eigenvector derivatives with repeated eigenvalues. AIAA J. (in press).

Computing response to periodic excitation

631

APPENDIX Using the equations of the eigenproblem, K8 = M@A

(AlI

eTM@ = I one can obtain, (K - aM)Q = M@(A - uf).

(A3)

Inverting the two sides of eqn (A3) yields, (K - aM)-’ = @(,I - al)-W’M-‘.

(A4)

Through inverting the two sides of eqn (A3), one can get, M-’ = BcPT. Substituting eqn (A5) into the right side of eqn (4) yields (K - aM)-’ = @(/4 - UT)-WT.

(‘46)

Using eqn (A6), it is easier prove that, (K - uM)-‘[M(K - uM)-‘][M(K - uM)-‘I, . . . , [M(K - uM)] = @(A - ul)-*-“FT.

L

J

7 k

(A7)

Partitioning the right sides of eqn (A7) in terms of available modes and unavailable modes, i.e. (K - uM)-‘[M(K - uM)-‘][M(K - uM)-‘I, . . , [M(K - uM)] = OLT(ALT- or)-*-‘@:, c I * x + @,Q(& - ur)-“-‘@: + e,,(n,,

- uL)-‘-‘cl$r

(A8)

and using the nomenclature of eqns (35) and (36), TX= @,,(A,,

- uI)-~ - ‘@:r,

(k = 0, 1,2,

. . .) (A9)

w = @HT(A”T- ul)-* - ‘a’,,,

(k = 0, 1,2, . .) (AlO)

one can get, ok + uk = (K - aM)-‘[M(K - uM)-‘][M(K - uM)-‘I,

. , [M(K - CM)-‘]

”

c

- @,(n, - al)-‘-‘@:.

(All)