Abstract
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. The present method possesses the features of unconditional stability and self-starting and can achieve fourth-order accuracy in displacement, velocity, and acceleration simultaneously. The algorithm is a matrix form and no need to iterate in the calculation. The convergent accuracy is verified by a numerical example, and the effectiveness is also verified by solving the dynamic responses of a vibration isolation system and the vibration responses of a pylon structure with cyclic loads and earthquake loads.
1. Introduction
The governing equations of structural dynamics after spatial discretized by the finite 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–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 first-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 differential equations into first-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 effective algorithms to solve the first-order ordinary differential equations (ODEs), such as the Runge-Kutta methods family [13], which can achieve an accuracy order higher than two with the unconditional stability. The 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.
The 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. Third, the effectiveness of the method in solving complicated problems with high accuracy will be verified by some numerical examples. Finally, the conclusion will be drawn.
2. Governing Equations of Motion
After spatial discretization, the second-order governing equations of motion can be denoted by [14]where denote the time; , , and denote displacement, velocity, and acceleration, respectively; , , and denote the mass matrix, damping matrix, and stiffness matrix, respectively; and denotes the load vector.
Expressing the initial conditions of the displacement and velocity as and , the initial acceleration can be given by (1)
Introducing , the second-order equations of motion in (1) can be transformed into a first-order form as
3. Two-Stage Fourth-Order Gauss Implicit Runge-Kutta Method
The first-order ODEs can be denoted aswhere and are unknown and known function vectors, respectively.
Discretizing the time domain by time step and assuming that is known at time , at time can be determined by the two-stage fourth-order Gauss implicit Runge-Kutta (GIRK24) method [13]. The GIRK24 method solving (4) giveswithwhere and are auxiliary variables.
As and are fully implicit in the algorithm, Attili [15] and González-Pinto and Rojas-Bello [16] have developed some effective solving procedure with iteration for the second-order initial value problems. The present study only considers the structural vibration responses in a linear range, and a matrix form-solving procedure is developed in the following.
The linear first-order ODEs can be expressed asorwhere and are constant matrices, is inversible, and is a function of .
Substituting the GIRK24 method of equation (6) and equation (7) into equation (9) or equation (10) yields
Combining equations (11) and (12) gives
Denoting , equation (3) can be expressed as the form of equation (9), in which
Substituting equation (14) into equation (13) yields the GRK24 method for the equations of motion asin which
Solving from equation (15) and substituting into (5), the displacement and velocity () can be updated by
To update the acceleration solution , the equations of motion at time are considered in the present study, which yields
As the solution of in the GIRK24 method is fourth-order accurate, the displacement and the velocity are fourth-order accurate as well. Denoting as the error of the acceleration solution, equation (16) gives
Substituting (1) at time into (19) yieldswhich indicates that the acceleration can also achieve the fourth-order accuracy.
4. Accuracy Validation of the Solving Procedure
The complete solving procedure is determined by equations (15), (17), and (18), which is self-started and can be solved step by step. The 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 first-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. The stability of the GIRK24 can be inherited from the first-order type ODEs to the second-order 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. The governing equation of the SDOF model is considered aswhere is the physical damping ratio with , is the physical frequency with , and the natural period has . The initial conditions are and . The dimensionless parameters are considered in the SDOF model.
Figure 1 shows the solution errors of the displacement, velocity, and acceleration at time with various time steps. The 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 fourth-order of accuracy simultaneously as indicated in the theoretical prediction. Therefore, the present developed method has a higher order of accuracy than the traditional second-order accurate time integration methods for structural analysis, especially in the acceleration solutions.
5. Numerical Examples
To verify the high-accurate performance and efficiency 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. The 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.
5.1. A Vibration Isolation System Model
Figure 2 shows a vibration isolation system with and . The external force has parameters of and . The 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. Thereafter, the dynamic response is dominated by the steady-state response with external frequencies. The 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 steady-state response period is stable where the magnitudes are not divergent even for the large time step of , 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 and Newmark-β method with . 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.
5.2. The Pylon Structure System
The 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. The material properties of the structure are considered as follows: the density to be 7850 kg/m³ and Young’s modulus to be 2 × 105 MPa. The truss members are constructed from two kinds of sections, where the areas of the section labeled by s are 600 mm2 and the rest of the areas are 3000 mm2, 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. The two-dimensional bar element is used to discretize the structure spatially, see references [17–19] for details, and the weights of the members are ignored. The first two circular frequencies are 43.995 rad/s and 117.769 rad/s, respectively. The Rayleigh damping is set as a damping ratio of 2% in the first 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 fine time step .
5.2.1. Two Joints with Cyclic Forces
The concentrated forces applied at the top-left and top-right joints are 120000 sin (0.5πt) N and 120000 cos (0.5πt) N, respectively. Figures 6–8 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.
5.2.2. Ground Earthquake Motion
The El Centro earthquake motion is applied to the system as in reference [19], and the right direction is set as the positive direction. The 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.
6. Conclusions
In this study, a high-accurate time integration method is proposed for solving the structural vibration responses. The method is derived from the two-stage fourth-order Gauss implicit Runge-Kutta (GIRK24) method for the first-order ODEs and is extended and developed for the equations of motion. The 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) self-starting; and (4) the method can achieve fourth-order accuracy in displacement, velocity, and acceleration simultaneously. The 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
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This project was supported by the Science and Technology Project of China Southern Power Grid Corporation (No. 037700KK52200002).