IN'IERNATION,AI. JOUKNAL FOR NUMERICAI. METHODS IN ENGINEERING, VOL. 39, 4199-4214 (1996) TWO SIMPLE FAST INTEGRATION METHODS FOR LARGE-SCALE DYNAMIC PROBLEMS IN ENGINEERING SU M M A RY A new simple explicit two-step method and a new family of predictor-corrector integration algorithms are developed for use in the solution of numerical responses of dynamic problems. The proposed integration methods avoid solving simultaneous linear algebraic equations in each time step, which is valid for arbitrary damping matrix and diagonal mass matrix frequently encountered in practical engineeringdynamic systems. Accordingly, computational speeds of the new methods applied to large system analysis can be far higher than those of other popular methods. Accuracy. stability and numerical dissipation are investigated. Linear and nonlinear examples for verification and applications of the new methods to large-scale dynamic problems i n railway engineering are given. The proposed methods can be used as fast and economical calculatIon tools for solving large-scale nonlinear dynamic problems in engineering. time integration: explicit method; predictor-correctoralgorithm; Newmark method; structural dynamics; large-scale dynamic problem K E Y WORDS: 1. INTRODUCTION Step-by-step time integration algorithms are widely used to solve the equations of multidegree of freedom systems in structural dynamics, especially in non-linear structural dynamics. Discussions on these algorithms can be found in the book of Belytschko and Hughes.' However, for large-scale problems, as is frequently the case in modern dynamic analysis in practical engineering problems, the calculation speed becomes very low and the cost is high. Therefore, it is necessary to develop more efficient integration algorithms for large-scale complex system analysis. Basically. there are two general classes of algorithms for dynamic problems: implicit and explicit. When implicit algorithms, such as the Newmark-P method' and the Wilson-0 r n e t h ~ d , ~ are applied to large-scale structural dynamic problems, the cost and the computational time increase dramatically with the degrees of freedom of systems because large-scale simultaneous algebraic equations must be solved in each time step although large time steps are permitted. On the contrary, explicit schemes tend to be inexpensive. Hoff and Taylor4 have pointed out that if lumped mass and damping matrices are used, an explicit scheme probably consists of pure vector operations. This is very convenient for computers with vector processors and the disadvantage of conditional stability can be partially alleviated through a vectorized implementation. Therefore, these explicit algorithms become more competitive on large-scale problems compared to the more stable implicit algorithms. * Professor of Mechanical Engineering CCC 0029--5981/96/244199- 16 ,(? 1996 by John Wiley & Sons, Ltd. Receiced 9 June 1994 Recised 13 October 1995 4200 W.-M. ZHAI During the last ten years or so, advances have been made to improve explicit methods. A state-of-the-art in high-order explicit schemes for use in structural dynamic applications has been summarized in Reference 4 by Hoff and Taylor. As indicated by these authors, it seems that the second-order accurate central difference method still remains the most popular explicit algorithm although several new developed explicit methods such as Katona-Zienkiewicz5 method have their own advantages. These new developed methods usually have lower stabilities than the central difference method, and some of them need more than one function evaluation (calculation of the internal forces), that is the most expensive part in nonlinear problems which are discretized by finite elements, finite differences or boundary elements. When applied to nonlinear structural dynamics, the central difference method is still supposed to solve a set of linear algebraic equations in each step unless the problem is an ideal one which satisfied the following conditions: (a) the mass matrix is diagonal, and (b) the damping matrix can be neglected or is proportional to the mass matrix. However, condition (b) is difficult to satisfy for practical engineering problems, while condition (a) can usually be observed. The present work is first aimed to develop a new simple explicit method, which has at least the same stability limit as the central difference method, and needs only simple vector operations in each step no matter what form the damping matrix of the system is. In Section 2, the new explicit scheme and its stability, accuracy, and numerical dissipation are presented, respectively. Secondly, the paper will construct a new simple predictor-corrector algorithm taking the new explicit method as a predictor and the Newmark implicit method as a corrector, which will be given in Section 3. Because the advantage and disadvantages of both the implicit and the explicit methods are complementary, the combination of these two classes of algorithms can possibly provide access to satisfactory synthetic characteristics on accuracy, stability and cost. Hughes and Liu6 have presented a family of implicit-explicit algorithms for transient finite-element analysis, in which a predictor-corrector family, constructed from the Newmark family, showed the potential prospect of the combination algorithms. In order to examine the integration properties of the new methods, several numerical examples are analyzed in Section 4. As an example for demonstrating the effectiveness of the new methods in application to large-scale problems in practical engineering, the coupling dynamics of railway vehicle and track system is numerically simulated in Section 5. 2. NEW EXPLICIT INTEGRATION METHOD The matrix equation of linear structural dynamics is MX + GX + KX = F or MA + CV + KX = F (1) where M, C and K are the mass, damping and stiffness matrices, respectively, F is the vector of applied loads (a given function of time, F = F(t)),X, V and A (i.e. X,% and X) are the vectors of displacements, velocities and accelerations, respectively. The initial value problem consists of finding a function X = X(t) which satisfies equation (I), and the initial conditions X(0) = x, V(0) = vo Here X, and Vo are given vectors of initial displacements and velocities, respectively. 420 1 TWO SIMPLE FAST INTEGRATION METHODS 2.1. Integration scheme The new explicit scheme for approximate solutions of equations (1) and (2) is proposed as follows + V,At + (f + $)AnAt2 - $A,-1At2 V,+1 = V, + ( 1 + q)A,At - qA,-iAt Xn+l= X, (3) where X,, V, and A, are the approximations to X(t = nAt), V(t = nAt) and A(t = nAt), respectively, Af is the time step, and $ and cp are free parameters that control the stability and numerical dissipation of the algorithm. Obviously, the above scheme is similar in form to the well-known Newmark method. In fact, the construction of the new scheme is enlightened by the latter. Substituting equation (3) into equation (1) at time step t = (n 1)At + MAn+l + C V n + l + K X n + l = F n + l (4) and rearranging the terms yield A,+1 =M-lF',+, (5) where F,+l = F,+1 - KX, - (C + KAt)V, - [(l + q ) C + (i+ $)KAt]A,At + (vC + $KAt)An-lAt (6) in which F,, I = F[t = (n + l)At]. To start the integration procedure, one can easily let cp = I) the initial conditions (2) as well as A0 = M-'(Fo =0 at the first time step and use - CVo - KXo) (7) Therefore, the scheme is self-starting. If the mass matrix is diagonal, as is frequently the case in structural dynamics and will be assumed here, the new integration algorithm is explicit and need not solve any equations. 2.2. Stability To investigate the stability of the new explicit method, we need only consider the linear homogeneous form of equation (4) without damping for the single-degree-of-freedom case as a , + , +w2x,+1 =o (8) where w=Jiclm (9) The difference form of equation (8) is xn+2 + [(f + $)a2- 2]x,+1 + [(l+ cp - 2 ~ ) R +2 llx, + ($ - cp)R2x,- 1 =0 (10) where R = oAt (11) The eigenvalue equation of equation (10) takes the following form: EL3 + [(f + $)R2 - 2112 + [(+ + cp - 2$)sr2 + 1]E. + ($ - c p ) P = 0 (12) 4202 W.-M. ZHAI Table 1. Conditions of stabiiity of the new explicit method cp - __ * At .. .. - - - The requirement for stability is liLl < 1. By use of the transformation 2 l + z 1 -z = __ equation (12) becomes [4 + 2(cp - 2+)SzZ]Z3 + [4 + (4l+b- 4cp - 1 ) W ] Z 2+ 2cpSzZz + Q2 =0 (14) and then the requirement for stability is simply R,(Z) ,< 0. According to the Routh-Hurwitz criterion, stable steps can be derived, which have been given in Table I in detail. It can be seen from Table I that the range of stable steps is very wide. When cp = i,b = i, the stability limit is Ar < 2/w, which is the same as that of the central difference method. And if cp = j!j and I) = $, the stable step range is ( 2 / ~ ) \ , 4d At < (2/to) in which the maximum stable step is larger than 2/ro. & 2.3. Accuriicy Applying Taylor formula to A,- we have A,,-, = A,l - A,At + iAnAt2 - ... (1 5 ) where a superposed dot denotes a time derivative. Substituting equation (15) into equation (3) yields X,,, V,, + V,At + $A,At2 + $AnAt3 - i$A,At4 + O(At5) = V,, + A,At + cpAl,At2 - fipA,At3 + O(At4) = X, (16) TWO SIMPLE FAST INTEGRATION METHODS 4203 Local truncation errors can then be written as + (A + 411/)AnAt4 + O ( A t 5 ) E ( V ) = (+ - cp)A,At2 + (6 + f q ) A n A t 3+ O(At4) E(X) = (b - ll/)A,At3 (17) If rC, = & the order of accuracy of E(X) is O ( A t 4 ) ,and if = f the order of accuracy of E(V) will be O ( A t 3 ) .Otherwise the orders of accuracy decrease to O(At3)and O ( A t 2 ) ,respectively. Obviously, the new explicit integration method has the same order of accuracy as that of the Newmark's implicit method. 2.4. Numerical dissipation and dispersion In order to show numerical dissipation and dispersion of the new algorithm, the analytical procedure used by Hilber et aL7 is adopted here. Rewriting equation (10) as Xn+2 - 2P1Xn+1+ PZX, - P3xn-1 =0 (18) where The general solution of equation ( 1 8) is where ii(i = 1,2,3) is the root of eigenvalue equation (12) and cl is the constant defined by the initial data. If cp = II/, there is a zero root in equation (12) which is the so-called spurious root, ,13. And the other roots When Rc- 4 2ql+ 1 i., and i., become two complex conjugate eigenvalues of equation (12), called principal roots, which satisfy lil.21 d 1. Let us write as =P Qi = exp[Q( in which P =P, - %+i)] (23) 4204 W.-M. ZHAI Newmark (p=1/4) -16 0.00 .04 .08 (b) 0.00.04 .12 -16 .20 . 2 4 At/T .08 .12 .16 . 2 0 .24 At/T Figure I . (a) Algorithmic damping characteristics of the new explicit method: amplitude decay function for various cp( = I+);(b) Relative period errors of new explicit scheme (9= $), Houbolt and Newmark methods In this case, the solution of equation (18) could be written in the form x, = exp( - ztot,) [cl cos(ot,) + c2 sin (or,)] (25) where (I, = QjAt t, = nAt As measures of the numerical dissipation and dispersion, we consider the amplitude decay function AD and the relative period error TD, respectively, AD = 1 - exp( - 274) T D = ('T - T ) / T = Q/SZ - 1 where T = 2nJw and T = 2 4 8 . In Figure l(a) algorithmic amplitude decay values for various cp( = $) are shown. The continuous control of numerical dissipation in the explicit scheme is evident. When both cp and $ are larger than i or one of them is larger than f and the other equals to 4, there will be algorithmic damping to decrease amplitudes, and when cp = $ = i, there is no numerical dissipation. In Figure l(b) the relative period errors are plotted versus At/T for various cases. The relative period errors of the Houbolt and Newmark (B = 1/4) methods are depicted in Figure l(b) for comparison. The period obtained with the present explicit scheme is less than the actual period, whereas the calculated periods with the Houbolt, Newmark and Park methods are longer than the actual period. It is clear from Figure l(b) that the relative period error of the new explicit method is much less than those of the other methods. For example, when cp = $ = 4 the period error is about half of the error of the trapezoidal rule. 4205 TWO SIMPLE FAST INTEGRATION METHODS 3. NEW PREDICTOR-CORRECTOR INTEGRATION METHOD 3. I , Basic procedure Basically, the new predictor-corrector scheme for approximately solving equations (1) and (2) consists of the following four steps: I. Predict X and V by use of the above new explicit method: + V,At + (4+ $)A.At2 = V, + (1 + cp)A,At - v A , - Xp."+1 = X, Vp.,+1 where subscript p refers to the predicted value. 11. Evaluate A from equation (1) at the time step t = (n - $An1 lAt2 (28) At + 1)At: Ap,n+l = M-I(Fn+l - K X p . , + l - C V p , , + l ) 111. Correct X and V with the Newmark implicit method:2 X,+ 1 = X, + V,At + (i- P)AnAt2 + P A , , + 1 At2 V n + l = V f , + ( l -y)A,At +7Ap,,+1At (30) where fi and 7 are free parameters of the Newmark method. IV. Evaluate A as the final values from equation (4): An+* =M-'(Ff,+i- K X + i - C V , + , ) (31) Similar to the new explicit method, the above integration procedure can be automatically started by letting cp = I) = 0 at the first time step and using the initial conditions (2) as well as (7). 3.2. ModiJed procedure In order to enhance accuracy of the algorithm, some modification to the basic predictorcorrector algorithm can be made by the consideration of local truncation errors. By the use of Taylor formula, local truncation errors of equations (28) and (30) can be obtained as + (& + +rl/)A,At4 + O(Ar5) q)A,At2 + (b + fq)A,At3 + O(At4) E p ( X j = (b - 49A,,Ar3 Ep(V)= (f - and + (A - fP)A,At4 + O(At5) Ec(V ) = (4- y)A,,At2 + (f - fl/)A,At3 + O(At4) E,(X) = (f - fi)A,At3 (33) where a superposed dot denotes a time derivative, and subscripts p and c' represent the predictor and the corrector, respectively. Modification to predictor equation (28) can then be made by the error theory as + Exp(Xp,n Vm.n+l = Vp.n+l + 8 r . p V p . n - Vc,n) Xm,n+ I = Xp.n+ 1 &,n) (34) 4206 W.-M. ZHAl Table 11. Modified coefficients Parameter where subscript m means the modifier. Modification to corrector equation (30) is similar to equation (34), i.e. X,t+l =Xc,n+, + C x c ( X p . n + l -Xc.n+,) Vnt~= V c . n + l + ~ i . c ( V p . n +- ~V o , n + ~ ) (35) cxp, cVp,c , ~and crc in equations (34) and (35) indicate the modified coefficients relating to truncation errors, as shown in Table 11. Therefore the modified predictor-corrector algorithm will consist of six steps: I. Predict: X p , n + l ,V p , , + , =equation (28) 11. Modify: X,.,, V,,,, = equation (34) 111. Evaluate: Ap,n+l= M - ' ( F , + , - KX,.,+, - CV,,,,,+,) ,, IV. Correct: X c , n + l ,V c . n + l =equation (30) V. Modify: X,+ 1 , V,, = equation (35) VI. Evaluate: A,,,. = equation (31) In this case, cp = $ = cXp= c,, = 0 should be employed at the first time step to start the integration procedure automatically. Obviously, the basic and the modified predictor-corrector integration algorithms are explicit and need not solve any system of simultaneous equations as long as the mass matrix M is diagonal. 4. VERIFICATION EXAMPLES In this section, three examples including a linear, a piecewise linear and a non-linear problem are employed to examine the accuracy and stability of these new integration methods. Exumple I In the linear case, we consider the following initial-value problem i+x=O X(0) = 0 x(0) = 1 (37) 4207 TWO S I M P L E t AST I K T E G R A T I O YMETHODS ----Explicit( p = 0 . 5 0 ) ---Houbolt -.-Explicit ( VI = 0.55 ) Displacement(m) 0 9 6 3 Time (9) Figure 2. Some integrated results with time step At =0 5 s for linear case 4 2 0 -2 -4 0 5 10 15 20 25 Figure 3. Examination of stability of thc new explicit method with a critical step Atcr = 2 s when cp = =f When the time step is sufficiently small, e.g. Ar = 005 s, the numerical results with the new integration methods are almost the same as the exact solution. Comparison of the new explicit method with the trapezoidal rule and with the Houbolt method is shown in Figure 2, with a time step of 0.5 s. To examine the stability, we consider the explicit scheme with the same value, 0.5, both for cp and According to Table I, the critical stable step is At,, = 2 i o = 2s for this example, which is well proved in Figure 3. A little increment of At from 2.0 to 2.01 s causes a divergent result. +. Example I1 Consider the piecewise linear system such as .t + k.u = 0 i ( 0 ) = 10 x(0) = 0 4208 W.-M. ZHAI Figure 4. Piecewise linear spring characteristic used in Example I1 40 - 20 30 10 Ath) - c Figure 5. Calculated errors 6 of explicit and predictor-xorrector ( P C )algorithms plotted versus time step At in which the variation of parameter k is given in Figure 4. Numerical results obtained with different values of and b are compared under the same condition of cp = y = $. It is shown that the best accuracy of the predictor-corrector schemes is attained if the optimum parameters of the explicit algorithm (cp = $ = +) are matched with y = f and fl = which refer to a high- accuracy implicit member of the Newmark family.' In Figure 5, calculated errors of the predictor-corrector schemes are compared with those of the explicit algorithm for various time steps. It is clear that the accuracy of the pure explicit integration has been largely improved by use of the present predictor-corrector integrations. Table 111 lists the comparison of the best accuracy predictor-corrector method (cp = $ = y = 5, j? = with the explicit algorithm (cp = y = f) as well as with the exact solution for three kinds of time steps. + A, A) Example 111 Consider a non-linear system with cubically hardening spring quoted in Reference 9: + i ioox + 1 0 0 0 ~ 3= o X(0) = 60 x(0) = 0 (39) 4209 TWO SIMPLE FAST INTEGRATION METHODS Table 111. Comparison of present integrated results with the exact solution t A t , = 0-005s At 1 Exact solution 10 20 30 40 50 60 0.48965 0.91873 1.23416 1.40091 1-40829 1.25591 0.48975 0.91891 1.23439 1.40115 1.40852 1.25612 - 0.0204% - Explicit At, = 0.05 s 0,48965 0.91873 1.23416 1.40092 1.40830 1.25592 . Error Explicit P-C* At, = 0.1 s P-C* Explicit 0.50000 0.93750 1.25781 1.42523 1-43139 1,27598 0.48958 0.91825 1.23448 1.40111 1.40857 1.25607 - - 1~oooO 1.50000 1.34000 0.91667 1.38306 1.24702 2.1 Yo 0.0522% .__ ~ 0% P-C* 8.846% 1.274% *P-C means predictor-corrector algorithm Converged Explicit ....-.-- ---- Jm pl ici t -E 2 0. 0 0.1 0. 2 0.3 Timeh) 0. 4 0.5 Figure 6 . Comparison of responses obtained with new explicit, Newmark implicit, and new prediction+orrector algorithms The same time step, At = 0.015s, employed by Park’ is adopted here. Numerical response processes with new explicit, Newmark implicit and present predictor-corrector algorithms are plotted in Figure 6. It can be seen that the response by the best accuracy predictor-corrector scheme is located between those by the new explicit (cp = $ = i)and by the Newmark implicit ( f l = $) methods. Both responses by the explicit and by the predictor-corrector methods trace the actual behaviour more closely than that by the Newmark method, and therefore much more closely than those by the Park method and by the Houbolt method, see Reference 9. When time step increases to 0.027 s, the best accuracy predictor-corrector scheme will be unstable in this example. The same situation happens to the pure explicit-predictor method at its critical step, 0.025s. Numerical experiments have shown that the critical stable step for this non-linear example could be increased to 0-033s if cp = y = and II/ = /? = &, which is close to the unstable step, about 0.035 s, of the implicit trapezoidal rule used in this nonlinear problem. The above parameter values usually imply good stability of the predictor-corrector algorithms but lead to somewhat reduction of amplitudes. 4210 W.-M. Z H A l ---__._ Figure 7. Dynamic model of thc interactive vehicldtrack system 5. APPLICATION EXAMPLES IN LARGE-SCALE PRACTICAL ENGINEERING So far the new explicit and the new predictor-corrector integration algorithms have already been widely applied in China to solving large-scale dynamic problems in practical railway engineering. For example, it was successfully used as the calculation algorithm of dynamic simulation in China's Train Driving Simulators (TDS), which are used for training railway drivers. Another example of the applications in Reference 10 by Huanp and Zhan, in which the present predictor-corrector algorithm was used to solve the locomotive equations of motion in order to analyse non-linear lateral running stability of the new designed six-axle locomotive DF9. Here the application of present algorithms to railway vehicle-track coupling dynamics will be presented in detail. Railway vehicle-track coupling system consists of a moving vehicle and a continuous track structure. Dynamic model of the overall system is given in Figure 7. The considered vehicle model is supported on two double-axle bogies at each end and is described as a 10-degree-of-freedom lumped mass system comprising the carbody mass, M , , and its moment of inertia, J,, the two bogie masses, M , , and the associated moments of inertia, J , , and four wheelset unsprung masses, M , . The bogie sideframe mass is linked with the wheel unsprung mass through the primary suspension, K,, , CS1,and linked with the carbody mass through the secondary suspension, K , 2 , C,z. The vehicle is assumed to move on the track at a constant velocity, V . In the track model, the rail is modelled as an infinite Euler- Bernoulli beam (mrr EI) discretely supported at rail/sleeper junctions by a series of springs, dampers and masses. The three layers of discrete springs and TWO SIMPLE FAST INTEGRATION METHODS 421 1 Table IV. Vehicle and track parameters Vehicle Track (per railseat) M , = 77000 kg M , = 1100 kg M , = 1200 kg .I, = 1.2 x lo6 kgm2 J , = 760 kgm2 K q l = 2.14 x lo6 N/m C,, = 4.9 x 104Ns/m K.1 = 5.32 x lo6 N/m C.2 = 7 x IO'N s/m L, = 4.25 m (half of bogie distance) L, = 0.875 m (half of wheelset distance) R = 0.42 m (wheel radius) *C61 type car has no primary suspension 111, = 60.64 kg/m EI = 6 . 6 2 ~ 106Nm' L, = 0.545m (sleeper spacing) M , = 125.5 kg Mb = 682.6 kg K , = 7.8 x lo7N/m C,=5x104Ns/m Kb = 2.4 x lo*N/m C, = 5.88 x lo4 N s/m Kf = 6.5 x lo7N/m C, = 3.1 x 104Ns/m K , = 7.8 x lo7 N/m C, = 8.0 x lo4 N s/m dampers represent the elasticity and damping effects of the railpads (Kpi,C,,), the ballast (KhirChi), and the subgrade (Kfi. Cfi),respectively. The two layers of discrete masses below the rail represent the sleepers, M s i , and the ballast, Mbi,respectively. In order to account for the shear continuity of the interlocking ballast particles, shear springs, KWi,and dampers, CWi, are introduced between the ballast masses to model the shear coupling effects in the ballast. Physical parameter values of the vehicle and the track are listed in Table IV, which refer to the properties of the Chinese C61 type train car and 60 kg/m rail track on the Da-Qing railway line. The system interaction between the vehicle and the track is achieved at the wheel/rail interfaces through wheel/rail force compatibility. This is described by the non-linear Hertzian contact theory commonly used in wheel/rail interaction problems: where P ( t ) is the wheel/rail contact force, G is the Hertzian wheel/rail contact coefficient, and d Z ( t ) is the elastic compression of the wheel/rail contact point. The equations of motion of the vehicle-track coupling system consists of those of the vehicle subsystem. which can be entirely described as second-order ordinary differential equations in time domain by applying the D'Alembert principle (see Reference ll), and those of the track subsystem. Equation of motion of the rail is in form of a partial differential equation with fourth order, which must be transformed into a series of second-order ordinary equations by means of Ritz's method so as to be numerically integrated with present time-stepping methods. All the track dynamic equations and corresponding transformation have been given in References 11 and 12. The final equations of the overall system can be expressed in terms of the standard matrix form as equation (I), in which the mass matrix M is diagonal. To solve such a large complex system with about 500 degrees of freedom involving the non-linear coupling relationship, the proposed explicit or predictorsorrector integration method is of particular advantage. There is no need to solve a large-scale algebraic equation group during each time step and thus the computational efficiency is greatly enhanced. Numerical experiments show that the maximum effective time steps of both methods are almost the same, Ar = 1 x 10-4s, for the solution of present vehicle-track dynamics. With this time step, these two kinds of integration give similar results with somewhat difference, see Figure 8(a). When the time 4212 W.-M ZHAI (a) d t =1X10-'s 50 7 9 . 9 3 * - 7 . 1 , 1 (s) , 1 - , , .T 100 50 n hd W 2 0 -50 - 100 0.00 0.02 0.04 0.06 T(s) 0.08 0.10 0.00 0.02 0.04 0.06 0.08 0.10 T(9) Figure 9. Comparison of time histories of rail acceleration a, between (a) measurement and (b) calculation, obtained for a rail joint with height differcnce h = 0.5 mm step is reduced to half of the original, there will be no distinguishable difference between two integration results, see Figure 8(b). In the following calculation, the time step will be fixed at a very small value, At = 2 x s, so that high-frequency responses of the wheel/rail system may be caught. Of course, the calculated results at this time step are same for both new integration methods. The validity of new integration methods applied to the solution of present large-scale dynamic problem has been verified by a full-scale field experiment carried out by us at the Da-Qing railway line in north China. The properties of the tested vehicle and track are given in Table IV. A jagged rail joint with height difference h for travel in the step-up direction was set to induce the wheel/rail impact vibration. The initial impact velocity of the wheel for this kind of joint is calculated as13 TWO SIMPLE FAST INTEGRATION METHODS 190 I25 -- . 1 -a- Measured Measured 100 4213 Calculated 0 0 10 20 30 40 50 60 V(k/h) Figure 10. Comparison of peak values of rail accelerations a, between measurement (-A-) and calculation (-O--) versus vehicle running speed V . The height difference of the jagged rail joint is 1.0 mm in this case where R is the radius of the wheel, V is vehicle velocity, m, is the mass of wheel, and meq is the equivalent impact rail mass.I3 Accelerations of track components were measured by accelerometers mounted on the rail and the sleeper and placed in the ballast, respectively. As an example of measurement results, Figure 9(a) gives the rail acceleration time history measured at the vehicle running speed of 52 km/h. Corresponding calculated result is given in Figure 9(b). It is shown that the numerical integration gives close result to the field measured result. Peak values of calculated results and measured data are compared in Figure 10, which also shows good correlation between the calculated and the actually measured results. All the comparisons above indicate the feasibility of the newly developed integration methods in applications to large-scale dynamic problems in practical engineering. It should be noted that the computation time for simulating an entire process of the wheel/rail impact at the rail joint is only about 10 min with the explicit method and about 15 min with the predictor-corrector method on an IBM 486/66 PC, respectively. 6. CONCLUSIONS A new family of explicit two-step algorithms has been presented. The method is of second-order accuracy and with similar algorithmic damping characteristics to the Newmark method. The stability limit of the new explicit method is not lower than that of the central difference method, e.g. when cp = $ = f the critical stable time step is 2/ro. A new family of predictor-corrector procedures has been developed to improve accuracy and stability of the new explicit method. The basic predictor-corrector method usually has higher accuracy than the pure explicit method, and the accuracy can be further improved by the modified predictor-corrector method. The computation time, however, will be somewhat longer than that of the pure explicit method. Numerical results of examples have shown that good accuracy and stability of the explicit method will be obtained if cp = $ = f, and when these parameters are matched with 7 = f and /3 = 1% the predictor-corrector method can get its best accuracy. The integration procedures of both methods are rather simple and it is not required to solve any linear algebraic equations as long as the mass matrix of the solved system is diagonal. Therefore, it is convenient and economical to use these algorithms to analyse large-scale non-linear dynamic problems in engineering. 4214 W.-M. ZHAI ACKNOWLEDGMENTS This study was financially supported through a research grant by Chinese Science Foundation for Outstanding Youths and supported by Trans-Century Training Programme Foundation for the Talents by State Education Commission of China. The author thanks Professor D. P. Chen at Southwest Jiaotong University for his constant encouragement and constructive discussions. REFERENCES 1. T. Belytschko and T. J. R. Hughes, Computationul Methods JiJr Trunsirnt Analy.sis (Compututionul Methods in M~c.ha~iics, V o l . I ) , North-Holland Elsevier, Amsterdam, 1983. 2. N. A. Newmark. ‘A method of computation for structural dynamics’, .I. Eiig. Mech. Dio. ASCE, 85(EM3), 67 -94 (1959). 3. E. L. Wilson, 1. Farhoomand and K. J . Bathe. ‘Nonlinear dynamic analysis of complex structure’, Earthquake. rng. .WUL‘~. dTn., I, 241 -252 (1973). 4. C. HofTand R. L. Taylor, ‘Higher derivative explicit one step methods for non-linear dynamic problems. Part I: design and theory’. Int. j . numer. methods eng., 29, 275-290 (1990). 5. M. G . Katona and 0. C. Zienkiewicz. ‘A unified set of single step algorithms. Part 3: The beta-m method, a generalization of the Newmark scheme’, In!. ,j. numer. mcJrhorls errg., 21, 1345-1359 (1985). 6. T. J . R. Hughes and W. K. Liu, ‘Implicit explicit finite elements in transient analysis: stability theory’, 3. A p p l . Mc& A S M E , 45, 371 374 (1978). 7. H. M. Hilber, T. J. R. Hughes and R. L. Taylor, ’Improved numerical dissipation for time integration algorithms in structural dynamics‘. Earrhquake. eng. striict. dj,n.. 5, 283 292 (1977). 8. S. P. Chan, H. L. Cox and W. A. Benfield, ‘Transient analysis of forced vibrations of complex structural mechanical systems’, J . Roy. Aero. Soc., 66, 457-460 (1 962). 9. K. C . Park, ‘An improved stiffly stable method for dircct intcgration of nonlinear structural dynamic equations’. J . Appl. MiJCh. A S M E , 42, 464-470 (1975). 10. C. R. Huang and F. S. Zhan, ‘The numerical bifurcation method of nonlinear lateral stability analysis of a locomotivc’, in 2. Y. Shen (ed.), Proc. 13/h I A V S D S p p . on Dynamic-s of Vehic,les on Ronds arid on Trucks, Vehicle Systetn Dynufirics, Vol 23, (SUPPI.), 1994, pp. 234-245. 11. W. M. Zhai. ‘The vertical model of vehicle-track system and its coupling dynamics’. 3. Chinu R u h a y Soc.. 14(3), 10 21 (1992). 12. W. M. Zhai and X. Sun. ‘A detailed model for investigating vertical interaction between railway vehicle and track‘. in 2. Y . Shen (ed.), Proc. 13th I A V S D Symp. on Dynamics of Vehicles on Rouds and o n Tracks. Vehicle Systcw Dynumics, Vol. 23, (SUPPI.),1994, pp. 603.415. 13. I. L. Ver. C. S. Ventres and M. M. Myles, ‘Whecl/rail noise- Part 111: impact noise generation by wheel and rail discontinuitics’, .I. Sound Vih., 46, 395417 (1976).