An Unconditionally Stable Explicit Method for Structural Dynamics

A Available onlin ne at at www.sciencedirect.com www.sciencedirect.com Available online

Procedia Engiineering 00 (201 1) 000–000 Procedia Engineering 14 (2011) 2519–2526

Proceedia Engineeering www w.elsevier.com/loocate/procedia

The Twellfth East Assia-Pacific Conference C on Structurral Engineerring and Coonstruction

An Unnconditiionally Stable S E Explicit Method for Struuctural Dynaamics S. H. Yina* Deppartment of Civill Engineering, Naational Taipei University of Techno ology

Abstract A new unconnditionally stabble explicit timee-integration m method is propo osed herein, wh hich can inheritt the numericall characteristics of any exissting implicit Runge-Kutta aalgorithms. Th hus, the proposed method nnot only holdss computationaal efficiency buut also possessees good numeriical characterisstics. Moreoverr, the paper preesents an exactt derivation off the increment of mechanical energy of an uundamped system to investigaate the stabilityy condition and d algorithmic ddamping of the proposed algorrithm for linearr problems, wh hich is differentt from the analyysis of spectrall radii of integgration algorithhms in previou us studies. Thhe numerical characteristics c of the proposeed method aree computationaally examined from fr the viewpo oints of energy conservation off conservative systems. s

© 2011 Published by Elseevier Ltd Keywords: Uncconditionally stabble explicit metho od, time-integratioon method.

1.

Intordu uction

Numericaal integration of ordinary differential d equuations is the most importaant technique iin continuouss time dynam mics. For solvinng structural dynamic d probblems, many methods m have been proposeed such as thee Newmark-β method (New wmark 1959), Wilson-θ meethod (Wilson n et al. 1973), HHT-α methhod (Hilber ett WBZ-α methood (Wood et al. a 1981), HP--θ1 method (H Hoff and Pah hl 1988a, b), C CH-α method d al. 1977), W (Chung andd Hulbert 19993). These methods m are implicit, seccond-order acccurate and A A-stable (i.e.. unconditionally stable) inn linear dynaamics, and caapable of con ntrolling numerical dissipaation in high-modes when addequate param meters are usedd. frequency m

000 ; fax: +0-000 -000-0000 . * Corresponnding author. Tell.: +0-000-000-00 E-mail addr dress: [email protected] .

1877–7058 © 2011 Published by Elsevier Ltd. doi:10.1016/j.proeng.2011.07.317

2520

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Procedia Engineering 00 (2011) 000–000

2

In contrast to implicit methods, explicit methods do not involve any iterative procedure and thus need far less computational effort in each time step. However, explicit methods can only have conditional stability, which means that a very small time step is required to obtain a stable integration result. To improve this shortcoming, some explicit methods with unconditional stability have been proposed (Chang 2002, 2007). It will be very favorable if an explicit method holds computational efficiency and also possesses the same numerical characteristics as the fourth order L-stable implicit Runge–Kutta methods (Hairer et al. 1993, 1996) for linear dynamic problems. In this study, the parameters used in the proposed explicit method are determined by equating its map (difference equation) derived from a differential equation to the corresponding one to the existing implicit Runge–Kutta algorithm such that the proposed explicit method can inherit the numerical characteristics of the emulated algorithm. Moreover, the paper presents an exact derivation of the increment of mechanical energy of an undamped system from the nth to the (n+1)th time step for the Runge–Kutta family to investigate their stability condition and algorithmic damping for linear problems, which is different from the analysis of spectral radii of integration algorithms in previous studies. In addition, the numerical characteristics of the proposed method are computationally examined based on energy conservation of conservative systems. 2.

Proposed Explicit Method

After the mathematical discretization of a linear structural system, its equation of motion can be written as

&& + Cx& + Kx = Mx f,

(1)

where M, C, and K are the mass, viscous damping, and stiffness matrices, x , x& , and &x& are the displacement, velocity, and acceleration vectors, and f is the external force vector. The general formulations of the proposed explicit method can be expressed as

xi +1= xi + Δtβ1x& i + Δt 2β 2&& xi x& i +1= x& i + Δtγ1&& xi + Δtγ 2&& xi +1

(2) ,

where the subscript i denotes the ith time step,

Δt is the time increment, and β1 , β 2 , γ1 and γ 2 are the

coefficient matrices determined as the following procedure. First, for convenience, we assume that the system is undamped and free of external force. Eq. (1) can be reduced to a system of first-order ordinary differential equations by introducing new variables which are usually made to be derivatives of the original variables as follows

= y& F= ( y , t ) Ay , where

I⎤ ⎡ 0 ⎡x ⎤ A = ⎢ −1 and y = ⎢ ⎥ . ⎥ ⎣-M K 0⎦ ⎣ x& ⎦

(3)

The general description of an s-stage Runge–Kutta method for solving Eq. (3) can be given by the formulas

2521

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Procedia Engineering 00 (2011) 000–000 s

y i +1= y i + Δt ∑ bmk m , m =1

s

s

k m = F(ti + cm Δt , y i + Δt ∑ aml= k l ) A(y i + Δt ∑ aml k l )

.

=l 1 =l 1

where the coefficients bm , cm , and aml in the Butcher tableau determine the particular method and its accuracy and stability properties. Then,

y i +1= y i + Δtb1k1 + Δtb2k 2 + L + Δtbs k s ⎡ R1 R 2 ⎤ = ⎢= ⎥ y i Ry i ⎣R3 R 4 ⎦

.

(4)

Eq. (4) is the so-called map (difference equation) of the Runge–Kutta method in which R is known as the amplification matrix. In addition, in order to construct the map of the proposed scheme, Eq. (2) is rearranged by substituting && x i = −M −1Kx i and && x i +1 = −M −1Kxi +1 into Eq. (2) as

⎤ ⎡ xi ⎤ Δtβ1 I − Δt 2β 2 M −1K ⎡ xi +1 ⎤ ⎡ = ⎥⎢ ⎥. ⎢ x& ⎥ ⎢ 2 2 −1 −1 −1 −1 ⎣ i +1 ⎦ ⎣ −Δtγ1M K − Δtγ 2 M K (I − Δt β 2 M K ) I − Δt γ 2 M Kβ1 ⎦ ⎣ x& i ⎦

(5)

To make the proposed explicit method possess the same numerical characteristics as the emulated Runge–Kutta methods, we equate the amplification matrix of the proposed method to the one of the Runge–Kutta algorithm. Then, the coefficient matrices β1 , β 2 , γ1 and γ 2 in Eq. (2) can be obtained as

β1 = 1

R Δt 2 β 2 1 2 (I − R1 )K −1M = Δt . −1 −1 − 1 γ1 = ( R + ( I − R 4 ) R 2 R1 ) K M Δt 3 γ 2 1 (I − R 4 )R 2 −1K −1M = Δt

(6)

Note that, unlike a general time integration method, it is interesting to find that the coefficient matrices

β1 , β 2 , γ1 and γ 2 depend on the structural properties and the size of the integration time step, rather than some constants. It should be mentioned that, once β1 , β 2 , γ1 and γ 2 are determined before time

marching, they are unchanged in a complete step-by-step integration procedure for solving linear dynamics. In this study, three Implicit Runge–Kutta methods are emulated, including the implicit midpoint method (Newmark average acceleration method), the modified extended backward differentiation formulae (ME-BDF) method, and the fourth order L-stable singly-diagonally-implicit Runge Kutta (SDIRK) method. To clearly demonstrate the implementation of the proposed method, a single degree-offreedom system is considered here. The matrix A in Eq. (3) can be expressed as

3

2522

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Procedia Engineering 00 (2011) 000–000

4

1⎤ ⎡ 0 = A ⎢= ⎥, ω 2 ⎣ −ω 0 ⎦

k . m

Thus, the amplification matrix R is reduced to be a 2 × 2 matrix in which R1, R2, R3, and R4 are reduced as scalars. The coefficient matrices β1 , β 2 , γ1 and γ 2 used in the proposed explicit method are also reduced as scalars and obtained using Eq. (6) as in Table 1. Note that the above procedure can be applied to any kind of Runge–Kutta algorithms and produces the equivalent explicit method. Table 1: The coefficient matrices The implicit midpoint method

β1

β2

β1 , β 2 , γ 1 and γ 2 used in the proposed explicit method The 4th order L-stable SDIRK method

The ME-BDF method

4 + ω 2 Δt 2

(2 + 7ω 2 Δt 2 + ω 4 Δt 4 ) 2(1 + ω 2 Δt 2 )3

2 4 + ω 2 Δt 2

(1 + 7ω Δt + 2ω Δt ) 2(1 + ω 2 Δt 2 )3

4

2

2

4

4

4(786432+114688ω 2 Δt 2 − 4352ω 4 Δt 4 − 1184ω 6 Δt 6 + 7ω 8 Δt 8 ) 3(16+ω 2 Δt 2 )5 1572864 + 360448ω 2 Δt 2 + 22272ω 4 Δt 4 − 352ω 6 Δt 6 + 3ω 8 Δt 8 3(16+ω 2 Δt 2 )5

8 × (147456 + 33792ω 2 Δt 2

γ1

1 2

2 + 3ω 2 Δt 2 2(2 + 7ω 2 Δt 2 + ω 4 Δt 4 )

+ 2072ω 4 Δt 4 − 31ω 6 Δt 6 ) 3(786432+114688ω 2 Δt 2 − 4352ω 4 Δt 4 − 1184 + 7ω 8 Δt 8 )

(1572864 + 360448ω 2 Δt 2 + 22272ω 4 Δt 4

γ2

1 2

1 + 7ω 2 Δt 2 + 2ω 4 Δt 4 2 + 7ω 2 Δt 2 + ω 4 Δt 4

− 352ω 6 Δt 6 + 3ω 8 Δt 8 ) 4(786432+114688ω 2 Δt 2 − 4352ω 4 Δt 4 − 1184ω 6 Δt 6 + 7ω 8 Δt 8 )

3.

Numerical stability from viewpoints of Energy conservation

For linear or nonlinear conservative systems, energy conservation should be satisfied. Theoretically, the dynamics of the map resulting from different step-by-step integration schemes for solving such systems should closely obey this law. In other words, the law can be employed to investigate stability condition and algorithmic damping of these schemes. To achieve this aim, first, the equation of motion of a single degree-of-freedom undamped system is considered as

mx&& + kx = 0, and the energy conservation relation can be obtained by multiplying Eq. (7) with the velocity

d 1 2 1 2 ( mx& + kx ) = 0. dt 2 2

(7)

x& as

2523

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Procedia Engineering 00 (2011) 000–000

The increment of mechanical energy of an conservative system from the ith to the (i+1)th time step can be expressed as

1 1 1 1 1 1 ( mx& 2 + kx 2 ) ii +1 = ( mx&i +12 + kxi +12 ) − ( mx&i 2 + kxi 2 ) 2 2 2 2 2 2 . 1 1 = m( x&i +1 + x&i )( x&i +1 − x&i ) + k ( xi +1 + xi )( xi +1 − xi ) 2 2

(8)

Recall the map of Runge–Kutta methods for a single degree-of-freedom system as

⎡ xi +1 ⎤ ⎡ R1 R2 ⎤ ⎡ xi ⎤ ⎡ R1 = ⎢ x& ⎥ ⎢= ⎥⎢ ⎥ ⎢ 2 ⎣ i +1 ⎦ ⎣ R3 R4 ⎦ ⎣ x&i ⎦ ⎣ −ω R2

R2 ⎤ ⎡ xi ⎤ . R1 ⎥⎦ ⎢⎣ x&i ⎥⎦

(9)

By substituting Eq. (9) into Eq.(8), Eq.(8) is written as

1 1 2 i +1 1 ( mx& 2 + = kx ) i m[ω 4 R2 2 xi 2 − 2 R1 R2ω 2 xi x&i + ( R12 − 1) x&i 2 ] 2 2 2 1 + k[( R12 − 1) xi 2 + 2 R1 R2 xi x&i + R2 2 x&i 2 ] 2 1 1 ( R12 + ω 2 R2 2 − 1)( kxi 2 + mx&i 2 ) = 2 2

(10)

Finally, we can obtain the incremental rate of mechanical energy of the system as

( Ei +1 − Ei ) / Ei = ( R12 + ω 2 R2 2 − 1) ,

(11)

Theoretically, the incremental rate should be zero because energy conservation should be satisfied. However, due to the temporal discretization, this rate may be positive, equal to zero, or negative, depending on R1 and R2 for different Runge–Kutta algorithms. Figure 1 shows the incremental rate of mechanical energy versus ωΔt in which ω is the natural frequency of the system and Δt is the integration time step. Based on the observations on these two plots, there are some discussions. For the 4th order explicit Runge–Kutta method, the condition of stability of 2 2 2 the method can be obtained by equating ( R1 + ω R2 − 1) to zero. The critical value of ωΔt is 2 2 4 4 2 3 calculated as 2 2 in the case where R1 =1 − ω Δt / 2 + ω Δt / 24 and R2 =Δt − ω Δt / 6 . When ωΔt is more than 2 2 , ΔEi /Ei is more than zero. Thus, the energy grows with time, and the method becomes unstable. In contrast, for the implicit midpoint method (Newmark average acceleration method), ΔEi /Ei is equal to zero and the energy can be conserved for all ωΔt. It is a well-known unconditional stable (A-stable) algorithm with no algorithmic damping. The ME-BDF method and the 4th order SDIRK method have two features in common. ΔEi /Ei is less than zero for all ωΔt more than zero and is equal to -1 when ωΔt goes to infinity. The former means that the two methods are not only unconditional stable but also provide algorithmic damping, and the latter implies that the energy Ei+1 is completely damped out (i.e. the two methods are so-called L-stable). In particular, the properties of algorithmic damping of the two methods are compared. The ME-BDF method is seen to damp the low modes too strongly since ΔEi /Ei drops drastically when ωΔt is more than 0.2. As a result, the ME-BDF method may not only dissipate spurious high-frequency response due to the spatial discretization but also suppress real and dominant low-frequency response. In contrast, the 4th order SDIRK method performs

5

2524

6

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Proceedia Engineering g 00 (2011) 000–0 000

better as a ““numerical low w-pass filter””, suppressingg high-frequen ncy response and a leaving loow-frequency y response inttact. Therefore, from this viewpoint, v it iis valuable to produce the equivalent exxplicit method d that possessses the same numerical n chaaracteristics ass the 4th orderr SDIRK meth hod. Note thaat the stability y ( −1 ≤ ΔEi / Ei ≤ 0 ) derived in thhis study is eequivalent too criteria of ttime integratiion scheme (i.e. 0 ≤ ρ (ωΔt ) ≤ 1 in which ρ is the sppectral radiuss based on eigenvalue analysis of the ccorrespondingg ( 12 + ω 2 R2 2 ) . amplificatioon matrix and is equal to (R

Figure 1: The inncremental rate of o mechanical eneergy of a linear SD DOF system for different time integration algorithm ms including the 4th order expliccit Runge–Kutta method, m the impliicit midpoint methhod, the ME-BDF F method, and the 4th order SDIRK K method.

4.

Applicaation to nonliinear systemss

For lineaar dynamics, a rigorous stability anallysis of time-stepping alg gorithms baseed on energy y conservationn is feasible. However, for fo nonlinear dynamics, itt is difficult to analyticallly obtain thee expression oof ΔEi /Ei of these algorith hms. Consequuently, the nu umerical stability of these algorithms in n dynamics off an undampeed nonlinear system s needs to be compu utationally inv vestigated by oobserving thee time variatioon of the mechanical enerrgy of the sysstem. Here, we w consider a governing eqquation of an n && + kx + ax 3 = undamped nnonlinear systtem as mx mark explicitt 0 , and solvee for it by using the Newm method and the proposedd explicit meth hods which reespectively in nherit the New wmark averagee acceleration n DIRK method. For convenience, these thhree proposed d method, the ME-BDF meethod, and thee 4th order SD med the NAA explicit methhod, the ME-BDF explicit method, andd the 4th orderr explicit metthods are nam SDIRK expllicit method. In this exxample, m, k, and a a are chosen as 1 kg, 1 Nt/m, and 0.1Nt/m3, the in ntegration tim me step is used d as 1.0 sec, and the iniitial condition ns are considdered as x (0 0) = 0.5 andd x& (0) = 0 . Ideally, thee 2 2 4 a 0.1266 for tthis nonlinearr mechanical energy ( mx& / 2 + kx / 2 + ax / 4 ) sshould be kept a constant as conservativee system at any time. However, duue to the teemporal disccretization annd numericall characteristiics of differennt algorithms, the mechaniccal energy maay oscillate orr decay with titime. Figure 2 shows the tiime series of mechanical m en nergy of the ssystem solved by different time t integratioon algorithmss including thhe Newmark explicit e method, the NAA exxplicit method d, the ME-BD DF explicit meethod, and thee 4th order SD DIRK explicit method. Onee may find thaat the NAA ex xplicit method and the 4th order SDIRK K explicit metthod perform better than the t other meth thods since th he correspond ding calculatedd mechanicall energy oscilllates around 0.1266 0 with very small ampplitude. It is no ot reasonable that the mechhanical energy y computed byy the ME-BDF F explicit metthod decays v ery fast due to o its too strong g algorithmic damping.

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Proceedia Engineering g 00 (2011) 000–0 000

o an nonlinear syystem with cubic stiffness (a=0.1) for different timee integration Figure 2: The tiime series of mecchanical energy of algorithms inclluding the Newmark explicit method, the NAA expplicit method, the ME-BDF expliciit method, and thhe 4th order SDIRK explicitt method.

5.

Conclu usions

In this paper, we proppose a new unconditionall u ly stable expllicit time-integration methood which can n o any existinng implicit Ru unge-Kutta allgorithms. In addition, thee inherit the nnumerical chaaracteristics of paper presennts an exact derivation off the incremennt of mechan nical energy of o an undampped system to o investigate tthe stability coondition and algorithmic a daamping of the proposed algorithm for linnear problems,, which is diff fferent from thhe analysis of spectral radii of integration n algorithms in n previous stuudies. Finally, the numericcal characteristics of the prooposed methodd for nonlineaar systems are computationaally examinedd from the vieewpoints of ennergy conservaation of conseervative system ms. Acknowledgments The author wishes to acknowledge a the National Science Coun ncil, Taiwan for f the generoous support off this work.

REFERE ENCES [1] [2] [3] [4] [5] [6]

Butcherr JC. Numerical methods m for ordin nary differential eequations. John Wiley W & Sons. ISB BN 0471967580,, 2003 Chang SY. Explicit Pseeudodynamic Algorithm with Uncconditional Stabillity. J. of Engine eering Mechanicss, ASCE. 128(9),, 2002, pp. 935-947. Chang SY. Enhanced, Unconditionally Stable Explicit P Pseudodynamic Algorithm. A J. off Engineering Meechanics, ASCE.. 133(5), 2007, pp. 541-5554. Chung, J and Hulbert GM. G A time integrration algorithm ffor structural dyynamics with improved numericall dissipation: Thee generallized-α method. J. J Appl. Mech. 60 0(6), 1993, pp. 3771–375. Hairer E, Nørsett SP, and a Wanner G. Solving ordinaryy differential equ uations I: Nonstiff problems. Beerlin, New York:: Springeer-Verlag. ISBN 978-3-540-56670 9 0-0, 1993. Hairer E and Wanner G. G Solving ordina ary differential eqquations II: Stifff and differential-algebraic probleems. Berlin, New w York: S Springer-Verlag, ISBN I 978-3-540--60452-5; 1996.

2525

7

2526

8 [7] [8] [9] [10] [11] [12] [13]

S.H. Yina / Procedia Engineering 14 (2011) 2519–2526 Author name / Procedia Engineering 00 (2011) 000–000 Hilber HM, Hughes TJR, and Taylor RL. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Eng. Struct. Dyn. 5, 1977, pp. 283–292. Hoff C and Pahl PJ. Development of an implicit method with numerical dissipation for time integration algorithms in structural dynamics. Comput. Meth. Appl. Mech. Eng. 67, 1988a, pp. 367–385 Hoff C and Pahl PJ. Practical performance of h1-method and comparison with other dissipative algorithms in structural dynamics. Comput. Meth. Appl. Eng. 67, 1988b, pp. 87–110 Newmark NM. A method of computation for structural dynamics. J. Engrg. Mech. Div. 85, 1959, pp. 67–94. Park KC. An improved stiffly stable method for direct integration of nonlinear structural dynamic equations. J. Appl. Mech. 42, 1975, pp. 464–470. Wilson EL, Farhoomand I, and Bathe KJ. Nonlinear dynamic analysis of complex structures. Earthquake Eng. Struct. Dyn. 1, 1973, pp. 241–252. Wood WL, Bossak M, and Zienkiewicz OC. An alpha modification of Newmark’s method. Int. J. Numer. Methods Eng. 15, 1981, pp. 1562–1566.