Performance of a Three-Substep Time Integration Method on Structural Nonlinear Seismic Analysis

In this work, a study of a three substeps’ implicit time integration method called the Wen method for nonlinear finite element analysis is conducted. 'e calculation procedure of the Wen method for nonlinear analysis is proposed. 'e basic algorithmic property analysis shows that the Wen method has good performance on numerical dissipation, amplitude decay, and period elongation. 'ree nonlinear dynamic problems are analyzed by the Wen method and other competitive methods. 'e result comparison indicates that the Wen method is feasible and efficient in the calculation of nonlinear dynamic problems. 'eoretical analysis and numerical simulation illustrate that theWen method has desirable solution accuracy and can be a good candidate for nonlinear dynamic problems.


Introduction
In last several decades, time dependent hyperbolic equations coupled with their calculation methods are widely used in various field. Generally, it is difficult to obtain analytical responses of these equations. Time-frequency method and time integration method are often used to deal with dynamic problems. For instance, the time-frequency method is an important technology in dynamic response analysis, especially for seismic data processing and interpolation [1]. Timefrequency analysis also plays a very important role in many areas such model updating, vibration control, dynamic optimization, etc. [2]. e time integration method has been employed extensively in linear and nonlinear dynamic problems. In general, time integration methods can be categorized into two major categories: explicit methods [3][4][5] and implicit methods [3,6,7]. Almost all implicit methods are unconditionally stable, and thus a large time step can be adopted. By contrast, most of explicit methods need to adopt a small time step to ensure stability. However, in terms of time cost in one time step, explicit methods have advantage over implicit methods, because the implicit methods spend more computation time for the calculation of stiffness matrix inversion and iteration than explicit methods.
For nonlinear dynamics, both explicit and implicit methods are used to solve different types of problems. Explicit methods are widely used to solve high-frequency transient problems such as impact problems. However, for structural dynamics problems where the dynamic response is mainly composed of low frequency modes, implicit methods are more suitable [8]. For implicit methods, a series of single-step methods such as Newmark method [9], HHT-α method [10], and Generalized-α method [11] are firstly proposed. e trapezoidal rule is extensively employed in these methods, because it performs well in the dynamic analysis of linear systems [12]. However, for nonlinear problems, it is necessary to realize a dynamic balance at each time step in the calculation of matrix iteration [13]. But for trapezoidal rule, it become quite unstable when large deformation occurs in the dynamic analysis of structural system, because it cannot satisfy the conservation of energy and moment of inertia. Actually, Kuhl and Crisfield provide energy-conserving and decaying algorithms for predicting the performance of time integration method in nonlinear analysis [14]. e numerical dissipation of time integration method for filtering or suppression of the spurious highfrequency mode is still very important for nonlinear analysis [15]. However, it is difficult for single-step time integration method to achieve desirable numerical dissipation.
In order to overcome the shortcomings of single-step methods, Bathe and Baig proposed a composite time integration method [16] (called the Bathe method), which employs trapezoidal rule in the first substep and employs three-point backward Euler formula in the second substep. eoretical analysis and numerical calculation show that the Bathe method performs well in nonlinear problems [17] and has excellent numerical dissipation property [18]. Following the Bathe method, a lot of composite methods [19][20][21][22][23][24] have been proposed by adopting different combinations of trapezoidal rule and other formulas. Among these methods, the Wen method [21] has excellent numerical dissipative properties and achieves high accuracy in linear analysis. e performance of the Wen method on nonlinear analysis deserves to be expected.
In this paper, the performance of the Wen method on nonlinear analysis will be studied through theoretical analysis and numerical simulation. In Section 2, the calculation procedure of the Wen method for nonlinear system is given. In Section 3, the basic algorithmic properties of the Wen method and other competitive methods are analyzed. In Section 4, numerical simulations of two well-known nonlinear dynamic examples and two nonlinear seismic examples are conducted by using the Wen method and other competitive methods. Finally, some concluding remarks are summarized in Section 5.

The Wen Method for Nonlinear Calculation
e governing equation of motion for nonlinear structural system discretized by the finite element method can be represented by where M is the mass matrix, F is the vector of internal forces, and R is the vector of external forces. X, _ X, and X represent the displacement, velocity, and acceleration vectors, respectively. e initial conditions are given by where X 0 and _ X 0 are, respectively, initial displacement and velocity vectors required to be given before solving the equation. X 0 is the initial acceleration.
According to the calculation procedure of the Wen method for linear system [21], the time steps corresponding to the three substeps are pΔt, (1 − p)Δt, and Δt, respectively. e parameter p satisfies the condition 0 < p < 0.406 to ensure the unconditional stability of the method. e first step of the method uses the trapezoidal rule, the second step adopts the three-point Euler backward formula, and the third step employs the Houbolt formula. e formulas of the first substep are where the subscript 't + pΔt' defines the value at time t + pΔt. e formulas of the second substep are where where In following analysis, the detailed calculation procedure for nonlinear dynamics is deduced. Substituting equations (3) and (4) into equation (5) leads to a residual function in terms of unknown displacement vector X t+pΔt , which is solved by Newton-Raphson iteration method. After simplification, we obtain (j � 1, 2, 3, · · ·) where ΔX t+pΔt . Once the displacement vector X t+pΔt (i.e., X (j) t+pΔt ) is obtained, the velocity and acceleration vectors can be solved with equations (3) and (4).

Mathematical Problems in Engineering
As for the third substep, we obtain t+Δt + n 4 n 1 X t + n 1 _ X t + n 4 n 2 X t+pΔt where ΔX t+Δt . e minimal iteration number j for each substep is determined by the criterion �������������� where the tolerance is ε � 10 − 10 . is method is a single parameter method, and all algorithmic properties are related to the algorithm parameter p. In linear analysis, a parameter case of the Wen method is suggested for general linear dynamic problems to obtain desirable algorithmic properties including calculation accuracy, numerical dissipation, etc. [7,21]. In structural linear seismic response analysis, the suggested parameter case also has been proved to have higher solution accuracy than other parameter cases [25]. In Section 4, the influence of algorithmic parameters on solution accuracy of nonlinear problems will be analyzed through numerical examples.

Spectral Stability.
e spectral characteristics of the time integration method are the basis for analyzing the stability and numerical dissipation properties. In structural dynamics, especially for nonlinear dynamics or transient problems of large structures, an excellent time integration method should have strong numerical dissipation for highfrequency modes, and the numerical damping for low frequency modes should be as small as possible.
In order to conduct stability analysis, a linear singledegree-of-freedom (SDOF) system is used, and the equation of the system is as follows: where ξ, ω, and r(t) are the damping ratio, the undamped angular frequency of the system, and the force excitation, respectively. ω satisfies T � 2π/ω, where T is the period of the system. By setting condition r(t) � 0 and ξ � 0 in equation (16), the recurrence formula of the Wen method for any undamped SDOF system can be obtained.
x t+Δt where the expression of amplification matrix A can be found in [21]. e spectral radii of four implicit methods with different parameters are plotted in Figure 1, where Ω � ωΔt � 2πΔt/T. For the Wen method, the solutions of five parameter cases are considered. As shown in the figure, all cases are unconditionally stable, and all methods show numerical dissipation for high-frequency modes except the case G-α (ρ ∞ � 1.0) and case OTTBIF (ρ ∞ � 1.0).

Accuracy Analysis.
Generally speaking, the numerical dissipation and accuracy of time integration method are usually measured by the amplitude decay (AD) and the period elongation (PE), respectively. e AD and PE results of the considered time integration methods are plotted in Figures 2 and 3, respectively. For the Wen method, it can achieve the lowest AD and PE in the case of p � 0.4. Among the considered methods, the OTTBIF method has the lowest AD and PE results, namely, the highest accuracy for linear dynamics. And the AD and PE results of all Wen method cases are smaller than those of the Bathe method.
e Generalized-α method shows the lowest accuracy.

Numerical Examples
In order to demonstrate performance of the Wen method on nonlinear dynamic analysis, four representative nonlinear dynamic problems are considered in this section. First, two classical nonlinear systems are tested to illustrate solution accuracy of the Wen method. en an eleven-story shear building and a space truss subjected to seismic load are considered. e suggested algorithm parameters p � 1/3 and p � 0.4 of the Wen method for linear dynamics are adopted for analysis. For comparison, the single-step Generalized-α method (the ρ ∞ � 0 case) [11], the two substeps' Bathe method (the r � 0.5 case) [16], and the three substeps' OTTBIF method (the ρ ∞ � 0 case) [24] are considered. e parameters of three methods are selected according to the corresponding references. In order to achieve almost the same time cost for each method, different methods use different time steps. For comparison, the same time step for all methods is also considered.

Van der Pol's Equation.
A nonlinear equation called Van der Pol' equation is considered here [26]. e equation is written as where the constant c represents the damping coefficient of the Van der Pol's equation and is negative here. e value of c and the initial conditions of the system are given.
e displacement, velocity, and acceleration results of the system are plotted in Figures 4-9, respectively. e quasi-exact solutions shown in the figures are obtained from the explicit Bathe method [27] by using very small time step Δt � 1 × 10 − 5 s. As shown in the figures, different time steps are adopted for different time integration methods to ensure roughly the same time cost. Figures 4-6 show that, among all considered methods and cases with different time steps, the p � 0.4 case of the Wen method shows highest solution accuracy. e accuracy of the OTTBIF method is slightly lower than that of the Wen method (p � 0.4). e Bathe method and the Wen method (p � 1/3) have desirable accuracy, but it is lower than that of the OTTBIF method and the Wen method (p � 0.4). In addition, although the Generalized-α method adopts a relatively small time step, it shows lower accuracy than that of other composite time integration methods. In dsl-9, the same time step is adopted to compare different methods. It is observed that accuracy of the OTTBIF method and the Wen method (p � 0.4) is much higher than the Generalized-α method, the Bathe method, and the Wen method (p � 1/3). In Table 1, time cost of all methods in time history of 0-500 s is given.
e Generalized-α method with Δt � 0.12s and Δt � 0.36s, respectively, has the longest and shortest elapsed time among all considered methods. Time costs of the Bathe method (Δt � 0.24s), the OTTBIF method (Δt � 0.36s), and the Wen method (p � 1/3 or 0.4, Δt � 0.36s) are very close. Combined with the calculation accuracy, the Wen method (p � 0.4) and the OTTBIF method have higher computational efficiency than other methods.

A Bead
Sliding on the Wire. In the fixed vertical (x, y) plane as shown in Figure 10, a bead of mass m slides on the smooth wire (the frictional force is neglected). e curve of wire can be expressed by the equation y � x 2 . e motion equation of bead in x direction is given.
where g is the acceleration of gravity. e initial conditions of the system are (0) � −5, x (0) � 0. e quasi-exact solutions are calculated by the explicit Bathe method with time step Δ t � 1 × 10 −5 s. In Figures 11-14, Δ t � 0.01 s and 0.02 s are, respectively, selected for the Generalized-α and Bathe methods, while Δ t � 0.03 s is adopted for the OTTBIF and Wen methods. In Figures 15-17, all methods adopt the same time step Δ t � 0.03 s. e displacement results of various methods for time duration 0-50 s are shown in Figure 11. It is observed that the results of the Generalized-α method are divergent. e displacement, velocity, and acceleration responses of the bead for time duration 180-200 s are plotted in dslf-14, respectively. e results show that the Wen method in case p � 0.4 achieves highest accuracy among all considered methods. In addition, accuracy of the Wen method (p � 1/3) and the OTTBIF method is lower than that of the Bathe method, and the energy dissipation of the OTTBIF method is obviously larger than that of other methods. Table 2 illustrates the elapsed time of all methods at t � 200s. e results illustrate that the Wen method in case p � 0.4 has the lowest energy decay and consumes less time cost than the OTTBIF method. e Bathe method (Δt � 0.03s) needs less computation time than other methods, but begets the largest proportion of energy dissipation. In Figures 15-17, when the same time step Δt � 0.03s is used for all methods, the Wen method (p � 0.4) displays far higher solution accuracy than other methods.

Nonlinear Structural Seismic Response Analysis.
Here two nonlinear structural systems under seismic load are tested. e two considered structural systems are an elevenstory shear building and a space truss.
e Bouc-Wen hysteresis model [28] is employed for the nonlinear dynamic analysis of these two structural systems. We adopt the Rayleigh damping to form damping matrix, and the damping ratios are 3% and 5% for the first and second vibration modes, respectively. Since the first two modes are intrinsically selected by experience, the numerical dissipation properties of the considered methods cannot be reflected in these two examples.
Different time steps are adopted for all time integration methods to ensure roughly the same time cost for the dynamic analysis. e quasi-exact solution is obtained using the Bathe method (r � 0.5) with very small time step Δt � 1 × 10 − 5 s. For composite time integration methods, the required loads at the substeps are obtained by employing linear interpolation of two adjacent seismic acquisition load values.

An Eleven-Story Shear Building.
e eleven-story shear building shown in Figure 18 is subjected to the 1940 EI Centro North-South seismic load (see Figure 19). e system is idealized as a simple shear model and has 11 degrees of freedom. e lumped mass and shear stiffness of each story are illustrated in Figure 18. e first story (isolation layer) of the structural system adopts Bouc-Wen hysteresis model [28] (see Figure 20). More details of the structure can be found in [29]. e acquisition frequency of seismic load is 50 Hz, and we intercept 30 s of the seismic load for dynamic analysis. e curve of seismic acceleration load in terms of time is shown in Figure 19.
e motion equation of linear structural system under seismic load can be expressed as [30] Mathematical Problems in Engineering 5

Mathematical Problems in Engineering
where M, C, and K are mass, damping, and stiffness matrices of the considered system. X g is the vector of ground motion acceleration. R is the influence vector, which defines the displacement of the system when the foundation produces static unit displacement in the direction of load excitation. e motion equation of the nonlinear part of the system can be written as follows: where α is the peak value of inelastic force and satisfies α � Q y (1 − k post /k pre ). Q y , k pre , and k post are the total force at the critical point of yield, the stiffness before yield, and stiffness after yield, respectively; A, β, and c satisfy A � 2β � 2c � k pre /Q y . e hysteretic restoring force Z of each time step is solved by Runge-Kutta method [31].      Figure 21. It can be observed that the internal force exceeds the bottom story yield limit. In Figure 22, the displacement results and errors of top story obtained from various methods are illustrated. e velocity results and errors are plotted in Figure 23. For clarity, the displacement and velocity during 0-10 s are analyzed. It is noticed that all four methods can achieve close accuracy. In Figures 24 and   25, the displacement and velocity errors of different methods with the same time step Δt � 0.024s are, respectively, plotted. It is observed that the errors of the Wen method and the OTTBIF method are much smaller than those of the Bathe method and the Generalized-α method. In Table 3, time cost of various methods with different time step in time history of 0-30 s is given. e elapsed time of the Wen method and the OTTBIF method with Δt � 0.024s is slightly Step ( longer than the Bathe method and the Generalized-α method with the same time step Δt � 0.024s, while the Bathe method (Δt � 0.016s) and the Generalized-α method (Δt � 0.008s) have longer elapsed time than the Wen method (Δt � 0.024s) and the OTTBIF method (Δt � 0.024s). Different from the linear seismic response analysis, in each time step, the nonlinear seismic response analysis entails a large amount of computation time for the nonlinear hysteretic force. e composite time integration method adopts the linear interpolation scheme to calculate the hysteretic force required for each substep, which reduces time cost for nonlinear problems. erefore, the composite method can greatly reduce time cost for nonlinear hysteretic force calculation.
is example verifies effectiveness of the Wen method and other composite time integration methods for structural nonlinear seismic response analysis. e three-step Wen method and the OTTBIF method have higher efficiency than the Generalized-α and Bathe methods. Figure 26, a space truss is considered. e total height of the structural system is 84 m, including 8 m for the first story and 76 m for the next 19 stories. e horizontal size of the structure in the first five stories is 30 × 40 m, which is divided into 5 × 3 grids; the horizontal size of the next 5 stories is 20 × 24 meters, which is divided into 3 × 2 grids; and the horizontal size of 11 to 20 stories is 20 × 16 meters, which is divided into 2 × 2 grids. e total degree of freedom of the structural system is 810. e elastic modulus of all linear elements is 2.1 × 10 11 N/m 2 . e mass of all elements is 7850 kg/m 2 , and the cross-sectional area is 25 cm 2 . All elements are connected by rigid nodes. e 16 cross braces (marked in blue) in the first story are nonlinear elements.

A Space Truss. As shown in
e Bouc-Wen hysteresis model [28] is adopted for these elements, and the expression of motion equation is given by equations (21)- (23). e nonlinear parameters k pre and k post of the 6 cross braces (perpendicular to the y-axis) are, respectively, 29.8 MN/m and 14.9 MN/m [32]. To obtain strong nonlinearity, the critical point of yield Q y is adopted as 4770 N which is far less than that in [32]. Likewise, the nonlinear parameters k pre and k post of the 10 cross braces (perpendicular to the x-axis) are, respectively, 27.3 MN/m and 13.65 MN/m [32], while a small critical point of yield Q y � 3860 N is adopted to obtain strong nonlinearity. is structure system is also used in [32]; however, in this analysis, bar element is used and different parameters are selected.
e system is subjected to East-West (shown in Figure 27), North-South (shown in Figure 19), and Up-Down (shown in Figure 28) directions of 1940 EI Centro seismic loads in the X-axis, Y-axis, and Z-axis, respectively. e dynamic analysis in point A (seen in Figure 26) is considered. e calculated hysteretic restoring force curves of the 16 nonlinear elements are shown in Figure 29 where the internal force of all nonlinear elements exceeds yield limit. e displacement and velocity results (in time duration    erefore, computation efficiency of four methods follows Wen > OTTBIF > Bathe > Generalized-α. As for the same time step (see Figures 32 and 33), the OTTBIF method and the Wen method show much higher accuracy than the Bathe method and Generalized-α method. is example verifies effectiveness of the Wen method and other composite time integration methods for nonlinear seismic response analysis of space truss.

Conclusions
In this work, the Wen method is applied to the analysis of nonlinear dynamics. e calculation procedure of the Wen method for nonlinear analysis is proposed. e basic properties of the Wen method are provided. Numerical simulations of some nonlinear dynamic problems, especially for nonlinear seismic response analysis, are studied by the Wen method and other competitive methods. e conclusion of this study can be summarized as follows.
(1) e Wen method shows concise calculation procedure for nonlinear dynamics analysis. (2) It is observed that, for the Wen method, the parameter p values of desirable algorithmic properties (including AD, PE, and analytical accuracy) for linear dynamics are also suggested for nonlinear dynamic analysis. Data Availability e raw/processed data can be obtained by contacting the corresponding author by e-mail (De Zhou 210026@ csu.edu.cn).

Conflicts of Interest
e authors declare that they have no conflicts of interest.