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

Authors: Xiangsheng Lei; Yanfeng Wang; Xinghua Wang; Gang Lin; Xinyan Qin
Journal: Mathematical Problems in Engineering
 (2022)

---

## 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.


---

## Body

## 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](1)Mu¨+Cu˙+Ku=Ft,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 stiffness matrix, respectively; and Ft denotes the load vector.Expressing the initial conditions of the displacement and velocity asu0 and u˙0, the initial acceleration can be given by (1)(2)Mu¨0=−Ku0−Cu˙0+F0,F0=F0.Introducingv=u˙, the second-order equations of motion in (1) can be transformed into a first-order form as(3)M00Mu˙v˙=−0−MKCuv+0Ft.


## 3. Two-Stage Fourth-Order Gauss Implicit Runge-Kutta Method

The first-order ODEs can be denoted as(4)y˙=fy,t,where y and f are unknown and known function vectors, respectively.Discretizing the time domain by time stepΔt and assuming that yn is known at time tn, yn+1 at time tn+1=tn+Δt can be determined by the two-stage fourth-order Gauss implicit Runge-Kutta (GIRK24) method [13]. The GIRK24 method solving (4) gives(5)yn+1=yn+12Δty˙α1+y˙α2,(6)y˙α1=fyn+14Δty˙α1+p1Δty˙α2,tn+p2Δt,(7)y˙α2=fyn+p3Δty˙α1+14Δty˙α2,tn+p4Δt,with(8)p1=14−36,p2=12−36,p3=14+36,p4=12+36,where y˙α1 and y˙α2 are auxiliary variables.Asy˙α1 and y˙α2 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 as(9)Ay˙=By+φt,or(10)y˙=A−1By+A−1φt,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(11)Ay˙α1=Byn+14Δty˙α1+p1Δty˙α2+φtn+p2Δt,(12)Ay˙α2=Byn+p3Δty˙α1+14Δty˙α2+φtn+p4Δt.Combining equations (11) and (12) gives(13)A−14ΔtB−p1ΔtB−p3ΔtBA−14ΔtBy˙α1y˙α2=IIByn+φtn+p2Δtφtn+p4Δt.Denotingy=uT,vTT, equation (3) can be expressed as the form of equation (9), in which(14)A=M00M,B=0M−K−C,φt=0Ft.Substituting equation (14) into equation (13) yields the GRK24 method for the equations of motion as(15)M^y˙^α=K^yn+F^,in which(16)y˙^α=u˙α1T,v˙α1T,u˙α2Tv˙α2TT,M^=M−14ΔtM0−p1ΔtM14ΔtKM+14ΔtCp1ΔtKp1ΔtC0−p3ΔtMM−14ΔtMp3ΔtKp3ΔtC14ΔtKM+14ΔtC,K^=I0I00I0IT0M−K−C,F^=0,FTtn+p2Δt,0,FTtn+p4ΔtT.Solvingy˙^α from equation (15) and substituting into (5), the displacement and velocity (u˙n=vn) can be updated by(17)un+1=un+12Δtu˙α1+u˙α2,vn+1=vn+12Δtv˙α1+v˙α2.To update the acceleration solutionu¨n+1, the equations of motion at time tn+1 are considered in the present study, which yields(18)Mu¨n+1=−Kun+1−Cvn+1+Ftn+1.As the solution ofy in the GIRK24 method is fourth-order accurate, the displacement and the velocity are fourth-order accurate as well. Denoting Eun+1¨ as the error of the acceleration solution, equation (16) gives(19)Mu¨tn+1+Eun+1¨=−Kutn+1+OΔt4−Cu˙tn+1+OΔt4+Ftn+1.Substituting (1) at time tn+1 into (19) yields(20)Eun+1¨=OΔt4,which 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 as(21)u¨+2ξωu˙+ω2u=0,where ξ is the physical damping ratio with ξ=0.2, ω is the physical frequency with ω=2π/T, and the natural period has T=1.0. The initial conditions are u0=1.0 and u˙0=1.0. The dimensionless parameters are considered in the SDOF model.Figure1 shows the solution errors of the displacement, velocity, and acceleration at time t=0.4 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.Figure 1
The solution errors of the single degree of freedom system.


## 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

Figure2 shows a vibration isolation system with m=10kg and k=40kN/m. The external force F0sinωt has parameters of F0=100N and ω/2π=50Hz. 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 Δ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.Figure 2
A vibration isolation system model.Figure 3
The displacement response results.Figure 4
The error of displacement response.


### 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 Δt=0.0001s.Figure 5
The pylon structure model.Figure 6
The vertical displacement response of the right mass point.


#### 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.Figure 7
The vertical velocity response of the right mass point.Figure 8
The vertical acceleration response of the right mass point.


#### 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.Figure 9
The horizontal acceleration response peak-value of the right mass point.


## 5.1. A Vibration Isolation System Model

Figure2 shows a vibration isolation system with m=10kg and k=40kN/m. The external force F0sinωt has parameters of F0=100N and ω/2π=50Hz. 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 Δ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.Figure 2
A vibration isolation system model.Figure 3
The displacement response results.Figure 4
The error of displacement response.


## 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 Δt=0.0001s.Figure 5
The pylon structure model.Figure 6
The vertical displacement response of the right mass point.


### 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.Figure 7
The vertical velocity response of the right mass point.Figure 8
The vertical acceleration response of the right mass point.


### 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.Figure 9
The horizontal acceleration response peak-value of the right mass point.


## 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.Figure 7
The vertical velocity response of the right mass point.Figure 8
The vertical acceleration response of the right mass point.


## 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.Figure 9
The horizontal acceleration response peak-value of the right mass point.


## 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.



---

*Source: 3563393-2022-11-02.xml*