AN ASSESSMENT OF TIME INTEGRATION SCHEMES FOR NON-LINEAR DYNAMIC EQUATIONS

Journal of Sound and Vibration (1996) 192(1), 321–331

AN ASSESSMENT OF TIME INTEGRATION SCHEMES FOR NON-LINEAR DYNAMIC EQUATIONS Y. M. X Department of Civil and Building Engineering, Victoria University of Technology, Melbourne, Victoria 3000, Australia (Received 20 March 1995, and in final form 24 July 1995) In the transient analysis of engineering structures by the finite element method there are several commonly used time integration schemes such as the method of central difference, the Newmark method, the Wilson-u method, the Houbolt method and the a-method of Hilber–Hughes–Taylor. In this paper these time integration schemes are applied to simple non-linear dynamic systems and the performances of these schemes are compared. It is found that there are significant differences in the accuracy of these schemes. The result of some of these schemes are so inaccurate that it becomes clear that they are not suited for time integration over a long time duration. The behaviour of these schemes is also compared with that of the fourth order Runge–Kutta method, which has been the preferred integration scheme for non-linear systems of a single or a few degrees of freedom. 7 1996 Academic Press Limited

1. INTRODUCTION

In the analysis of an engineering structure under impact, earthquake or similar dynamic loading the finite element method is usually used for the space discretization which results in a set of ordinary differential equations Mu¨ + Cu˙ + p(u) = f(t),

(1)

where M is the mass matrix, C is the damping matrix, f is the external force vector dependent on time t, and p is the internal force vector, which is generally a non-linear function of the displacement vector u. To obtain the time history of variables such as displacement and acceleration from equation (1), a step-by-step numerical time integration is often used. Numerous time integration schemes have been proposed. Among them the most commonly used time integration schemes in finite element applications are the method of central difference, the Newmark method [1], the Wilson-u method [2], the Houbolt method [3] and the a-method of Hilber–Hughes–Taylor [4]. Excellent reviews of these schemes can be found in the books by Hughes [5] and Argyris and Mlejnek [6]. For linear systems, the stability and accuracy of these time integration schemes have been thoroughly studied and well understood (see, for example, the book by Wood [7]). However, for non-linear systems the behaviour of these schemes needs to be carefully examined. It has been shown by Xie and Steven [8] that even for a simple non-linear system of an unforced periodic oscillator the application of the Newmark method may lead to chaotic predictions. 321 0022–460X/96/160321 + 11 $18.00/0

7 1996 Academic Press Limited

. . 

322

The main purpose of this paper is to compare the performances of the time integration schemes for non-linear dynamic systems as represented by equation (1). In particular, the following two equations are considered: u¨ + S1 u(1 + S2 u 2 ) = 0

(2)

which represents the well-known unforced and undamped Duffing oscillator, and u¨ + S tanh (u) = 0.

(3)

In these equations, S1 , S2 and S are constants. When S1 q 0, S2 q 0 and S q 0, equation (2) represents a hardening elastic spring and equation (3) a softening elastic spring. Both equations (2) and (3) are examples of conservative systems which maintain a constant total energy E in their exact solutions. For equation (2) E = 12 u˙ 2 + 12 S1 u 2 + 14 S1 S2 u 4

(4)

E = 12 u˙ 2 + S ln [cosh (u)].

(5)

and for equation (3)

Both equations (2) and (3) have solutions representing periodic oscillations. The exact period and time history of the vibration governed by equation (2) or equation (3) can be found analytically by using methods described in references [9, 10] or obtained numerically by using time integration with very small time steps. Although the Runge–Kutta method [11, 12] is generally not used in finite element codes for transient dynamic analysis, it has been the preferred time integration scheme for solving non-linear systems of a single or a few degrees of freedom. For the completeness of this study the fourth order Runge–Kutta method is also included for comparison. 2. A BRIEF DESCRIPTION OF VARIOUS TIME INTEGRATION SCHEMES

Details of the time integration schemes considered in this paper can be found in many texts. Therefore only a brief description of the Newmark method, the Wilson-u method, the Houbolt method and the a-method of Hilber–Hughes–Taylor and the fourth order Runge–Kutta method is given below. 2.1.    Let un , u˙ n and u¨ n be the known displacement, velocity and acceleration at time tn . The displacement and velocity at the next time station tn + 1 = tn + Dt are related to un , u˙ n and u¨ n by un + 1 = un + Dtu˙ n + Dt 2( 12 − b)u¨ n + Dt 2bu¨ n + 1 ,

(6a)

u˙ n + 1 = u˙ n + Dt(1 − g)u¨ n + Dtgu¨ n + 1 ,

(6b)

where b and g are parameters. The Newmark method requires that the equilibrium equation (1) be satisfied at tn + 1 : i.e., Mu¨ n + 1 + Cu˙ n + 1 + p(un + 1 ) = f(tn + 1 ).

(7)

By substituting the relationships (6a) and (6b) into equation (7) u¨ n + 1 becomes the only remaining unknown in the non-linear equation, which may be solved by using the Newton–Raphson method. The Newmark method with b = 1/4 and g = 1/2 is known as the average acceleration method, which has been most frequently used. The Newmark method with b = 0 and g = 1/2 is equivalent to the method of central difference.

   

323

2.2.  -u  The Wilson-u method is a variation of the Newmark method. Here the equilibrium equation (1) is required to be satisfied at tn + u = tn + uDt (with u e 1); i.e., Mu¨ n + u + Cu˙ n + u + p(un + u ) = f(tn + u ).

(8)

The displacement and velocity at tn + u are related to un , u˙ n and u¨ n by equations similar to equation (6a) and (6b): un + u = un + uDtu˙ n + (uDt)2( 12 − b)u¨ n + (uDt)2bu¨ n + u ,

(9a)

u˙ n + u = u˙ n + uDt(1 − g)u¨ n + uDtgu¨ n + u .

(9b)

By substituting the relationships (9a) and (9b) into equation (8), u¨ n + u can be found by solving the non-linear equation. The acceleration at tn + 1 is then deduced from u¨ n and u¨ n + u by linear interpolation, u¨ n + 1 = (1 − 1/u)u¨ n + (1/u)u¨ n + u ,

(10)

from which the displacement and velocity at tn + 1 can be obtained by using the standard Newmark formulae (6a) and (6b). In the Wilson-u method one assumes b = 1/6 and g = 1/2. The third parameter u is often chosen to be 1·4. 2.3.    The Houbolt method is one of the first time integration schemes developed for the computer analysis of aircraft structures. As for the Newmark method, the equilibrium equation is satisfied at tn + 1 = tn + Dt as given in (7). However, the velocity and acceleration at tn + 1 are related to the previous displacements by u˙ n + 1 = (11un + 1 − 18un + 9un − 1 − 2un − 2 )/(6Dt),

(11a)

u¨ n + 1 = (2un + 1 − 5un + 4un − 1 − un − 2 )/Dt 2.

(11b)

By substituting these relationships into equation (7) one obtains a non-linear equation for un + 1 in terms of the known displacements un , un − 1 and un − 2 . The Houbolt method is a three-step algorithm, which is not convenient to implement. Apart from the extra memory needed for storing the displacements at several steps, the Houbolt method requires an initialization to specify u−1 and u−2 at the starting point at t0 . The initialization scheme as described in reference [6] is used in this study. 2.4.  a-  –– Compared with other methods discussed in this section, the a-method is a relatively new integration scheme which appears to have improved accuracy characteristics. In recent years, the a-method has been implemented in general finite element analysis codes such as LUSAS (London University Stress Analysis System) [13]. In the a-method again one uses the standard Newmark formulae (6a) and (6b). The only variation from the Newmark method is that it operates with a modified equilibrium relationship, Mu¨ n + 1 + Cu˙ n + 1 + (1 + a)p(un + 1 ) − ap(un ) = f(tn + 1 ), where a is an additional parameter.

(12)

. . 

324

2.5.  –  There are many different versions of the Runge–Kutta method. The one used in this study is a fourth order scheme as given in reference [12]. First of all the second order differential equation system (1) needs to be transformed to the following first order system u˙ = v,

v˙ = M−1 [f(t) − Cv − p(u)],

(13a, b)

which is of the general form y˙ = F(t, y),

(14)

where y is a vector including both the displacement u and the velocity v, and the function F is known from the right side of equation (13). With known un and vn , and therefore known yn , the displacement and velocity at the next time station is given as yn + 1 = yn + (k1 /6) + (k2 /3) + (k3 /3) + (k4 /6),

(15)

in which k1 = DtF(tn , yn ),

k2 = DtF(tn + Dt/2, yn + k1 /2),

k3 = DtF(tn + Dt/2, yn + k2 /2),

k4 = DtF(tn + Dt, yn + k3 ).

(16a, b) (16c, d)

3. COMPARISONS OF SEVEN TIME INTEGRATION SCHEMES

The accuracy and the relative computational costs of different time integration schemes need to be compared so that the user may have some idea as to which scheme might be used and which scheme should be avoided. In this study seven schemes are considered. These are the average acceleration method (the Newmark method with b = 1/4 and g = 1/2); the method of central difference (the Newmark method with b = 0 and g = 1/2); the Newmark method with b = 0·3025 and g = 0·6; the Wilson-u method with b = 1/6, g = 1/2 and u = 1·4; the Houbolt method; the a-method with a = −0·1, b = 0·3025 and g = 0·6; and the fourth order Runge–Kutta method. All the parameters chosen here are those which have been frequently used in the literature. For the conservative systems represented by equations (2) and (3), the exact solution maintains a constant energy at any time in the response history as given in equations (4) and (5). In these cases the inaccuracy of the numerical solution can be conveniently measured by the deviation of the total energy from its initial value. The error of the time integration may be defined as e = =(E − E0 )/E0 = × 100%,

(17)

in which E and E0 are calculated from equation (4) or equation (5) by using the current displacement and velocity, and the initial displacement and velocity, respectively. In what follows the seven integration schemes as described above are applied to equations (2) and (3). 3.1. u¨ + S1 u(1 + S2 u 2 ) = 0 (   S1 q 0  S2 q 0) Here S1 = 100 and S2 = 10, and the initial conditions are u0 = 1·5 and v0 = 0. The period of this oscillator T is equal to 0·15. The phase portrait of the exact solution is given in Figure 1(a). Also shown in Figure 1 are the phase portraits of the numerical results from different integration schemes with the same time step size Dt = T/20 over a time duration

   

325

Figure 1. Phase portraits of the hardening spring u¨ + 100u(1 + 10u 2 ) = 0, u0 = 1·5, v0 = 0, Dt = T/20, time duration = 100T. (a) Exact solution; (b) Newmark b = 0·25, g = 0·5; (c) Newmark b = 0, g = 0·5; (d) Newmark b = 0·3025, g = 0·6; (e) Wilson-u, b = 16 , g = 12 u = 1·4; (f) Houbolt; (d) a-method, a = −0·1, b = 0·3025, g = 0·6; (h) Runge–Kutta.

. . 

326

T 1 Maximum percentage errors of the total energy over the time duration of 100T for the hardening spring u¨ + 100u(1 + 10u 2 ) = 0, u0 = 1·5, v0 = 0

Dt T/15 T/20 T/25 T/50 T/100 T/200 T/1000

Average acceleration method

Method of central difference

Newmark with b = 0·3025 and g = 0·6

Wilson

Houbolt

a-Method

Runge– Kutta

7·1 4·1 2·6 0·7 0·2 0·0 0·0 (1·48 s)

7·7 4·3 2·7 0·7 0·2 0·0 0·0 (0·99 s)

99·8 99·6 99·3 97·7 93·1 82·9 39·5 (1·51 s)

93·9 87·6 79·6 36·1 7·1 1·0 0·0 (1·68 s)

98·8 97·3 95·3 75·3 30·0 5·4 0·1 (1.92 s)

67·8 49·2 34·3 6·7 0·9 0·1 0·0 (1·59 s)

44·5 18·4 7·3 0·3 0·0 0·0 0·0 (2·00 s)

of 100T. It is clear that the results vary greatly with different schemes. Some of the results have deviated far from the exact solution. The percentage error of the total energy as defined in equation (17) is increasing with time, with small oscillations in some cases. The maximum values of the percentage errors for the time duration of 100T are summarized in Table 1 for the seven schemes with various time step sizes. The numbers in the last row of Table 1 are the CPU time in seconds after 100 000 steps of time integration with Dt = T/1000. All the computations are carried out on a Pentium 90 MHz Personal Computer. It is seen that the results of the method of central difference are almost as accurate as those of the average acceleration method. Nonetheless, the computational cost for an explicit method such as that of central difference is, as expected, much lower than that for the implicit average acceleration method. The worst scheme among the seven is obviously the Newmark method with b = 0·3025 and g = 0·6. Even with the step size as small as T/1000, it gives an error as high as 39·5%. For the a-method, although the parameters b and g are exactly the same as those in the previous Newmark method, the introduction of the extra parameter a has improved the accuracy significantly. For the fourth order Runge–Kutta method when the time step is large it is not as accurate as the average acceleration method or the method of central difference, but when the time step is being reduced the convergence rate of the Runge–Kutta method is higher than those of any other schemes considered here. The grave errors in the results of the Newmark method with b = 0·3025 and g = 0·6, the Houbolt method and the Wilson-u method are caused by the strong numerical damping possessed by these schemes. In particular, the numerical damping in the Newmark method with b = 0·3025 and g = 0·6 is unacceptably high. Even for a moderate time duration of 100T, what is supposed to be a periodic vibration is rapidly damped towards the static equilibrium point u = 0 and v = 0. Although numerical damping also exists in the a-method, it is not as severe as that of the Houbolt method or the Wilson method. 3.2. u¨ + S tanh (u) = 0 (   S q 0) Here S = 100 and the initial conditions are u0 = 4·0 and v0 = 0. The period of this oscillator is equal to 1·14. The phase portrait of the exact solution is given in Figure 2(a). Also shown in Figure 2 are the phase portraits of the numerical results from different integration schemes with the same time step size of T/20 over the time duration of 100T. Although a softening spring is considered now, the behaviour of the integration schemes is similar to that for the previous hardening spring. Strong numerical damping is again

   

327

Figure 2. Phase portraits of the softening spring u¨ + 100 tanh (u) = 0, u0 = 4·0, v0 = 0, Dt = T/20, time duration = 100T; (a)–(h) as Figure 1.

observed in the Newmark method with b = 0·3025 and g = 0·6, the Houbolt method and the Wilson-u method. The average acceleration method and the method of central difference give fairly accurate results when the time step Dt is as large as T/20. The errors

. . 

328

T 2 Maximum percentage errors of the total energy over the time duration of 100T for the softening spring u¨ + 100 tanh (u) = 0, u0 = 4·0, v0 = 0

Dt T/15 T/20 T/25 T/50 T/100 T/200 T/1000

Average acceleration method

Method of central difference

Newmark with b = 0·3025 and g = 0·6

Wilson

Houbolt

a-Method

Runge– Kutta

8·8 4·7 3·0 0·7 0·2 0·1 0·0 (3·72 s)

14·0 7·3 3·9 0·9 0·2 0·1 0·0 (2·22 s)

100·0 100·0 100·0 100·0 100·0 98·1 32·4 (3·74 s)

100·0 100·0 100·0 36·0 5·5 0·8 0·0 (3·79 s)

100·0 100·0 100·0 99·9 27·7 4·1 0·1 (4·36 s)

97·3 59·8 33·7 5·6 0·9 0·2 0·0 (5·05 s)

60·8 1·0 8·3 0·3 0·0 0·0 0·0 (4·51 s)

of these integration schemes with various time step sizes are summarized in Table 2. The Newmark method with b = 0·3025 and g = 0·6 proves once again to be extremely inaccurate—even with a time step as small as T/1000 it gives an error as high as 32·4%. The numbers in the last row of Table 2 are the CPU time in seconds after 100 000 steps of time integration with Dt = T/1000.Compared to the previous example of hardening spring, the present example of softening spring requires more than twice the computational time. 4. DISCUSSION

From Figure 2 it appears that the Runge–Kutta method is the most accurate integration scheme of the seven schemes. In Table 2 it is shown that the maximum error in the total energy drops to 1·0% when Dt = T/20 is used for the Runge–Kutta method. However, when a smaller time step Dt = T/25 is used, the maximum error in the total energy is ‘‘increased’’ to 8·3%. Such behaviour of having a less accurate solution with a smaller time step is not uncommon with non-linear systems. Xie and Steven [8] have already shown that with the average acceleration method a smaller time step may even change a stable solution to an unstable one.

Figure 3. Instability of the a-method. u¨ + 100 tanh (u) = 0, u0 = 10·0, v0 = 0, Dt = T/8; a = −0·1, b = 0·3025, g = 0·6.

   

329

Among the seven schemes considered in this paper, five are unconditionally stable for linear systems. These are the average acceleration method, the Newmark method with b = 0·3025 and g = 0·6, the Houbolt method, the Wilson-u method and the a-method. When they are applied to nonlinear systems, such unconditional stability is no longer valid. Even the a-method, which has moderate numerical damping, may give an unstable solution. An example is shown in Figure 3. Here the softening spring as discussed in section 3.2 is considered, but with different initial conditions, u0 = 10·0 and v0 = 0. The period of this oscillator is equal to 1·8. The unstable solution shown in Figure 3 was obtained when Dt = T/8 was used for the a-method. When a large time step Dt is used, in addition to possible instability, other problems may arise from the numerical time integration. Consider equation (2) and take S1 = −0·5 and S2 = −1·0 and the initial conditions u0 = 0·5 and v0 = 0.The period of this oscillator is equal to 7·3. The phase portrait of the exact solution is shown in Figure 4(a). If

Figure 4. Phase portraits of the spring u¨ − 0·5u(1 − u 2 ) = 0, u0 = 0·5, v0 = 0, T = 7·3, Dt = 2·5, time duration 2500. (a) Exact solution; (b) Newmark, b = 0·25, g = 0·5; (c) Newmark, b = 0·3025, g = 0·6; (d) Wilson-u, b = 16 , g = 12 , u = 1·4; (e) Houbolt; (f) a-method, a = −0·1, b = 0·3025, g = 0·6.

330

. . 

Dt = 2·5 (1T/3) is used for the time integration over the time duration of 2500, the solution by the average acceleration method becomes chaotic, as shown in Figure 4(b). Instead of a periodic oscillation, the numerical results jump irregularly between the system’s two static equilibria u = 1 and u = −1. For the two explicit schemes, i.e., the method of central difference and the Runge–Kutta method, a converged solution of the non-linear equilibrium equation cannot be obtained with the normal Newton–Raphson method when Dt = 2·5 is used. For the Newmark method with b = 0·3025 and g = 0·6, the Houbolt method, the Wilson-u method and the a-method the results are shown in Figures 4(c)–(f). It is clear that when large time steps are used the numerical damping in these four schemes is so high that the results bear no resemblance to the exact solution shown in Figure 4(a). 5. CONCLUSIONS

Time integration schemes commonly used in transient finite element analysis have been examined together with the fourth order Runge–Kutta method. It is found that the accuracy of these schemes varies greatly depending on the scheme, the parameters and the time step sizes. The average acceleration method and the method of central difference give reasonable results when Dt Q T/20. The Newmark method with b = 0·3025 and g = 0·6, the Houbolt method and the Wilson-u method possess strong numerical damping, which makes them unsuitable for any practical applications involving integration over a long time duration. In particular, it becomes clear from the results of this study that the Newmark method with b = 0·3025 and g = 0·6 should never be used for practical purposes. The numerical damping of the a-method is not as strong as that of the Newmark method with b = 0·3025 and g = 0·6, the Houbolt method or the Wilson-u method. As far as the relative computational costs are concerned, the most efficient method is the central difference method. The Runge–Kutta method costs slightly more than the implicit Newmark methods. The Houbolt method, the Wilson-u method and the a-method are also more expensive than the implicit Newmark methods. Although the Runge–Kutta method is not usually used for transient finite element analysis, it is far more accurate than other schemes when the time step is small (Dt Q T/50). If the time step is too large (Dt q T/5) a converged solution of the non-linear equilibrium equation may not be found by the method of central difference and the Runge–Kutta method. For other schemes, numerical results obtained when using such a large time step are often meaningless and sometimes chaotic. REFERENCES 1. N. M. N 1959 Journal of the Engineering Mechanics Division, American Society of Civil Engineers, 85(EM3), 67–94. A method of computation for structural dynamics. 2. E. L. W, I. F and K. J. B 1973 Earthquake Engineering and Structural Dynamics 1, 241–252. Nonlinear dynamic analysis of complex structures. 3. J. C. H 1950 Journal of the Aeronautical Sciences 17, 540–550. A recurrence matrix solution for the dynamic analysis of aircraft. 4. H. M. H, T. J. R. H and R. L. T 1977 Earthquake Engineering and Structural Dynamics 5, 283–292. Improved numerical dissipation for time integration algorithms in structural dynamics. 5. T. J. R. H 1987 The Finite Element Method, Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs, New Jersey: Prentice-Hall. 6. J. A and H.-P. M 1991 Dynamics of Structures. Amsterdam: North-Holland. 7. W. L. W 1990 Practical Time-stepping Schemes. Oxford: Clarendon Press.

   

331

8. Y. M. X and G. P. S 1994 Communications in Numerical Methods in Engineering 10, 393–401. Instability, chaos, and growth and decay of energy of time-stepping schemes for non-linear dynamic equations. 9. W. L. W and M. E. O 1988 Communications in Applied Numerical Methods 4, 205–212. Stability properties of some algorithms for the solution of nonlinear dynamic vibration equations. 10. T. M. C and J. H. G 1989 Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics 56, 149–154. An alternative frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. 11. C. W. G 1971 Numerical Initial Value Problems in Ordinary Differential Equations. Englewood Cliffs, New Jersey: Prentice-Hall. 12. W. H. P, S. A. T, W. T. V and B. P. F 1992 Numerical Recipes. Cambridge: Cambridge University Press; second edition. 13. FEA Ltd 1993 LUSAS Finite Element System. Kingston upon Thames, Surrey, U.K.