A High-Accurate Time Integration Method for Solving Structural Vibration Responses

A time integration method for the equations of motion is developed based on the Gauss implicit Runge-Kutta method to high-accurate solving the responses in structural vibration. Te present method possesses the features of unconditional stability and self-starting and can achieve fourth-order accuracy in displacement, velocity, and acceleration simultaneously. Te algorithm is a matrix form and no need to iterate in the calculation. Te convergent accuracy is verifed by a numerical example, and the efectiveness is also verifed by solving the dynamic responses of a vibration isolation system and the vibration responses of a pylon structure with cyclic loads and earthquake loads.


Introduction
Te governing equations of structural dynamics after spatial discretized by the fnite element method or the boundary element method lead to a set of second-order equations of motion in the time domain. As the implicit time integration methods generally can use a large time step with the unconditional stability, they have been developed to solve the equations of motion, especially for the structural vibration response. Typical implicit time integration methods include the Newmark method [1], the Houbolt method [2], the HHT-α method [3,4], the three parameters optimal (TPO) method [5,6] (the generalized-α method [7]), the composite methods [8][9][10]. Although these methods can solve the second-order equations of motion directly, they can only achieve second-order accuracy with the unconditional stability. Moreover, the acceleration response of some methods can only achieve frst-order accuracy [11], such as the HHT-α method and the TPO method.
To increase the solution accuracy of the equations of motion, an alternative way is to transform the second-order diferential equations into frst-order ones. After that, higher accurate solution method can be developed, such as the precise time step integration method [12]. Moreover, there are already some efective algorithms to solve the frst-order ordinary diferential equations (ODEs), such as the Runge-Kutta methods family [13], which can achieve an accuracy order higher than two with the unconditional stability. Te obstacle to use the Runge-Kutta methods for the equations of motion is that the complicated algorithm structure increases the computational cost. Nevertheless, if the computational cost is acceptable, a high-order Runge-Kutta method can be considered.
Te present study aims to develop a high-accurate time integration method for the structural vibration analysis based on a Gauss-type implicit Runge-Kutta method. First, the solution procedure for the equations of motion will be derived, where the displacement, velocity, and acceleration to be the same order of accuracy are considered. Second, the high-order accuracy of the method will be validated. Tird, the efectiveness of the method in solving complicated problems with high accuracy will be verifed by some numerical examples. Finally, the conclusion will be drawn.

Governing Equations of Motion
After spatial discretization, the second-order governing equations of motion can be denoted by [14] where t denote the time; u, _ u, and € u denote displacement, velocity, and acceleration, respectively; M, C, and K denote the mass matrix, damping matrix, and stifness matrix, respectively; and F(t) denotes the load vector.
Expressing the initial conditions of the displacement and velocity as u 0 and _ u 0 , the initial acceleration can be given by Introducing v � _ u, the second-order equations of motion in (1) can be transformed into a frst-order form as

Two-Stage Fourth-Order Gauss Implicit Runge-Kutta Method
Te frst-order ODEs can be denoted as where y and f are unknown and known function vectors, respectively. Discretizing the time domain by time step Δt and assuming that y n is known at time t n , y n+1 at time t n+1 � t n + Δt can be determined by the two-stage fourth-order Gauss implicit Runge-Kutta (GIRK24) method [13]. Te GIRK24 method solving (4) gives _ y α2 � f y n + p 3 Δt _ y α1 + 1 4 Δt _ y α2 , t n + p 4 Δt , with where _ y α1 and _ y α2 are auxiliary variables. As _ y α1 and _ y α2 are fully implicit in the algorithm, Attili [15] and González-Pinto and Rojas-Bello [16] have developed some efective solving procedure with iteration for the second-order initial value problems. Te present study only considers the structural vibration responses in a linear range, and a matrix form-solving procedure is developed in the following.
Te linear frst-order ODEs can be expressed as where A and B are constant matrices, A is inversible, and φ(t) is a function of t.
Substituting the GIRK24 method of equation (6) and equation (7) into equation (9) or equation (10) yields Combining equations (11) and (12) gives Denoting y � [u T , v T ] T , equation (3) can be expressed as the form of equation (9), in which 2 Mathematical Problems in Engineering Substituting equation (14) into equation (13) yields the GRK24 method for the equations of motion as in which Solving _ y α from equation (15) and substituting into (5), the displacement and velocity ( _ u n � v n ) can be updated by To update the acceleration solution € u n+1 , the equations of motion at time t n+1 are considered in the present study, which yields As the solution of y in the GIRK24 method is fourthorder accurate, the displacement and the velocity are fourthorder accurate as well. Denoting E € un+1 as the error of the acceleration solution, equation (16) gives which indicates that the acceleration can also achieve the fourth-order accuracy.

Accuracy Validation of the Solving Procedure
Te complete solving procedure is determined by equations (15), (17), and (18), which is self-started and can be solved step by step. Te displacement, velocity, and acceleration can be obtained successively. If the acceleration is no need to be solved, equation (18) can be omitted, and the method degenerates into the original GIRK24 method for the equations of motion in the frst-order form. It is indicated that the Gauss kind Runge-Kutta method possesses desirable properties to the Hamiltonian system, such as simplecticness [16], and the unconditional stability to linear systems. Te stability of the GIRK24 can be inherited from the frst-order type ODEs to the secondorder type equations of motion in linear structural dynamics analysis. A typical single degree-of-freedom (SDOF) model is adopted to verify the convergence accuracy of the developed method in the displacement, velocity, and acceleration. Te governing equation of the SDOF model is considered as where ξ is the physical damping ratio with ξ � 0.2, ω is the physical frequency with ω � 2π/T, and the natural period has T � 1.0. Te initial conditions are u 0 � 1.0 and _ u 0 � 1.0. Te dimensionless parameters are considered in the SDOF model. Figure 1 shows the solution errors of the displacement, velocity, and acceleration at time t � 0.4 with various time steps. Te slope of the logarithmic curves equals to the order of accuracy of the solutions. It can be validated that the displacement, velocity, and acceleration can achieve fourthorder of accuracy simultaneously as indicated in the theoretical prediction. Terefore, the present developed method has a higher order of accuracy than the traditional secondorder accurate time integration methods for structural analysis, especially in the acceleration solutions.

Numerical Examples
To verify the high-accurate performance and efciency of the developed method, numerical examples are performed in this section. A vibration isolation system and a pylon structure system are adopted as the examples of practical applications. Te dynamic response is also solved by the Newmark-β method as a comparison, where the parameters are chosen to make the method to be second-order accuracy. Figure 2 shows a vibration isolation system with m � 10kg and k � 40kN/m. Te external force F 0 sin ωt has parameters of F 0 � 100N      and ω/2π � 50Hz. Te physical damping ratio of 0.1 is considered, and the initial static condition is adopted. In the initial period, the dynamic response is dominated by the transient state response with natural frequencies. Tereafter, the dynamic response is dominated by the steady-state response with external frequencies. Te displacement solutions and their errors are shown in Figures 3 and 4, respectively, where the solution error is obtained by comparing them with the exact solution. It can be seen that the solution for the present developed method in the steadystate response period is stable where the magnitudes are not divergent even for the large time step of Δt � 0.008s, which is the unconditional stability feature inherited from the GIRK24 method. Moreover, the solution error is also bounded in the steady-state response period and is not higher than the amplitude of the exact solution. For the same time step, the present developed method shows much lower solution error than the Newmark-β method. Even if the former uses four times larger time step, it still shows a more accurate solution than the latter, see the results of the present developed method with Δt � 0.004s and Newmark-β method with Δt � 0.001s. For the steady-state response, the present method can still have considerably small solution error at the time points when only three time steps are used in a period.

Te Pylon Structure
System. Te two-dimensional pylon structure system [17] shown in Figure 5 is adopted to verify the stability and accuracy of the developed method in solving complicated structures. Te material properties of the structure are considered as follows: the density to be 7850 kg/m³ and Young's modulus to be 2 × 10 5 MPa. Te truss members are constructed from two kinds of sections, where the areas of the section labeled by s are 600 mm 2 and the rest of the areas are 3000 mm 2 , and their lengths are shown in Figure 5. Two lumped masses of 150 kg are applied at the top-left and top-right joints, respectively. Te twodimensional bar element is used to discretize the structure spatially, see references [17][18][19] for details, and the weights of the members are ignored. Te frst two circular frequencies are 43.995 rad/s and 117.769 rad/s, respectively. Te Rayleigh damping is set as a damping ratio of 2% in the frst two modes. Two kinds of loading conditions are considered as follows: (1) two cyclic concentrated forces applied at the top-left and top-right joints and (2) a ground earthquake load applied at the bottom of the structure; and the initial static condition is applied. As the exact solution is hard to be obtained, the reference solutions are adopted by the Newmark-β method with a fne time step Δt � 0.0001s.  show the solutions of the displacement, velocity, and acceleration at the top-right joint, respectively. It can be seen that the present developed method shows more accurate solutions than the Newmark-β method under the same time step. Even if the former uses four times larger time step, it still shows a more accurate solution than the latter, which is similar to the result in the previous example. Nevertheless, the present developed method and the Newmark-β method both show a feature that the error increases with the time step, and a large time step could lead to an unacceptable solution error which indicates that a proper time step should be chosen to compromise an acceptable computational cost.

Ground Earthquake Motion.
Te El Centro earthquake motion is applied to the system as in reference [19], and the right direction is set as the positive direction. Te peak response of the top-right joint happens at about 3.3 s as indicated in reference [19]. Figure 9 shows the acceleration solution near the peak response point. It can be seen again that the present developed method shows more accurate solutions than the Newmark-β method under the same time step. Even if the former uses four or much times larger time step, it still shows a more accurate solution than the latter at the time points. It should be noted that the time step cannot be too large as the time point of the peak would be skipped. As an accurate acceleration solution could be very important in the structure analysis under the earthquake load [19], the present developed method can provide a remarkably accurate solution and can be recommended in such practical applications.

Conclusions
In this study, a high-accurate time integration method is proposed for solving the structural vibration responses. Te method is derived from the two-stage fourth-order Gauss implicit Runge-Kutta (GIRK24) method for the frst-order ODEs and is extended and developed for the equations of motion. Te merits of the presented method are as follows: (1) the solving and updating formulations for the displacement, velocity, and acceleration are derived explicitly; (2) the method possesses the unconditional stability; (3) selfstarting; and (4) the method can achieve fourth-order accuracy in displacement, velocity, and acceleration simultaneously. Te numerical properties are validated by the theoretical analysis and numerical examples. Compared with the Newmark-β method, the present method needs to solve a higher-order matrix equation. Nevertheless, the present method is desirable and has advantages in solving high-accuracy demanded problems, such as the structure responses analysis under earthquake.

Data Availability
Te data used to support the fndings of this study are available from the corresponding author upon request.

Conflicts of Interest
Te authors declare that there are no conficts of interest regarding the publication of this paper.  Mathematical Problems in Engineering 7