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

Authors: Jinyue Zhang; Lei Shi; Tianhao Liu; De Zhou; Weibin Wen
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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


---

## Body

## 1. 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]. Time-frequency analysis also plays a very important role in many areas such model updating, vibration control, dynamic optimization, etc. [2]. The 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–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. The 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]. The numerical dissipation of time integration method for filtering or suppression of the spurious high-frequency 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. Theoretical 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–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. The 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 Section2, 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.


## 2. The Wen Method for Nonlinear Calculation

The governing equation of motion for nonlinear structural system discretized by the finite element method can be represented by(1)MX¨+FX,X˙,t=Rt,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. The initial conditions are given by(2a)X0=X0,(2b)X˙0=X˙0,(2c)X¨0=M−1R0−FX0,X˙0,0,where X0 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. The parameter p satisfies the condition 0<p<0.406 to ensure the unconditional stability of the method. The 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. The formulas of the first substep are(3)X˙t+pΔt=X˙t+12pΔtX¨t+X¨t+pΔt,(4)Xt+pΔt=Xt+12pΔtX˙t+X˙t+pΔt,(5)MX¨t+pΔt+Ft+pΔtX,X˙,t=Rt+pΔt,where the subscript ‘t+pΔt’ defines the value at time t+pΔt.The formulas of the second substep are(6)X˙t+1−pΔt=m1Xt+m2Xt+pΔt+m3Xt+1−pΔt,(7)X¨t+1−pΔt=m1X˙t+m2X˙t+pΔt+m3X˙t+1−pΔt,(8)MX¨t+1−pΔt+Ft+1−pΔtX,X˙,t=Rt+1−pΔt,where m1=2p−1/pp−1Δt, m2=1−p/p2p−1Δt, m3=2−3p/p−12p−1Δt.The formulas of the third substep are(9)X˙t+Δt=n1Xt+n2Xt+pΔt+n3Xt+1−pΔt+n4Xt+Δt,(10)X¨t+Δt=n1X˙t+n2X˙t+pΔt+n3X˙t+1−pΔt+n4X˙t+Δt,(11)MX¨t+Δt+Ft+ΔtX,X˙,t=Rt+Δt,where n1=−1/Δt, n2=1/p−12p−1Δt, n3=1/p2p−1Δt, n4=p2−p−1/pp−1Δt.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 Xt+pΔt, which is solved by Newton-Raphson iteration method. After simplification, we obtain (j=1,2,3,⋯)(12)4pΔt2MΔXt+pΔtj=Rt+pΔt−FXt+pΔtj−1,X˙t+pΔtj−1,t−M4pΔt2Xt+pΔtj−1−Xt−4pΔtX˙t−X¨t,where ΔXt+pΔtj=Xt+pΔtj−Xt+pΔtj−1. Once the displacement vector Xt+pΔt (i.e., Xt+pΔtj) is obtained, the velocity and acceleration vectors can be solved with equations (3) and (4).Likewise, for the second substep, we obtain(13)m32MΔXt+1−pΔtj=Rt+1−pΔt−FXt+1−pΔtj−1,X˙t+1−pΔtj−1,t−Mm32Xt+1−pΔtj−1+m3m1Xt+m1X˙t+m3m2Xt+pΔt+m2X˙t+pΔt,where ΔXt+1−pΔtj=Xt+1−pΔtj−Xt+1−pΔtj−1. After Xt+1−pΔt is obtained, X˙t+1−pΔt and X¨t+1−pΔt can be solved with equations (6) and (7).As for the third substep, we obtain(14)n42MΔXt+Δtj=Rt+Δt−FXt+Δtj−1,X˙t+Δtj−1,t−Mn42Xt+Δtj−1+n4n1Xt+n1X˙t+n4n2Xt+pΔt+n2X˙t+pΔt+n4n3Xt+1−pΔt+n3X˙t+1−pΔt,where ΔXt+Δtj=Xt+Δtj−Xt+Δtj−1.The minimal iteration numberj for each substep is determined by the criterion(15)ΔXt+rΔtjΔΔXt+rΔtjXt+rΔtjΔXt+rΔtj≤ε,r=p,1−p,1,where the tolerance is ε=10−10.This method is a single parameter method, and all algorithmic properties are related to the algorithm parameterp. 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.


## 3. Basic Properties of the Wen Method

In this section, spectral stability, amplitude decay (AD), and period elongation (PE) of the Wen method are studied. For comparison, the Generalized-α method [11] (the ρ∞=0, ρ∞=0.5, and ρ∞=1.0 cases), the Bathe method [16] (the r=0.5 and r=2−2 cases), and the OTTBIF method [24] (the ρ∞=0, ρ∞=0.5, and ρ∞=1.0 cases) are also considered in this section. The parameters for these methods are selected according to the corresponding references.


### 3.1. Spectral Stability

The 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 high-frequency modes, and the numerical damping for low frequency modes should be as small as possible.In order to conduct stability analysis, a linear single-degree-of-freedom (SDOF) system is used, and the equation of the system is as follows:(16)x¨+2ξωx˙+ω2x=rt,where ξ, ω, and rt 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 rt=0 and ξ=0 in equation (16), the recurrence formula of the Wen method for any undamped SDOF system can be obtained.(17)xt+ΔtΔt⋅x˙t+ΔtΔt2⋅x¨t+Δt=AxtΔt⋅x˙tΔt2⋅x¨t,where the expression of amplification matrix A can be found in [21].The spectral radii of four implicit methods with different parameters are plotted in Figure1, 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).Figure 1
Spectral radii of various methods for undamped system (ξ=0).


### 3.2. 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. The AD and PE results of the considered time integration methods are plotted in Figures2 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. The Generalized-α method shows the lowest accuracy.Figure 2
The percentage AD results of various methods.Figure 3
The percentage PE results of various methods.


## 3.1. Spectral Stability

The 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 high-frequency modes, and the numerical damping for low frequency modes should be as small as possible.In order to conduct stability analysis, a linear single-degree-of-freedom (SDOF) system is used, and the equation of the system is as follows:(16)x¨+2ξωx˙+ω2x=rt,where ξ, ω, and rt 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 rt=0 and ξ=0 in equation (16), the recurrence formula of the Wen method for any undamped SDOF system can be obtained.(17)xt+ΔtΔt⋅x˙t+ΔtΔt2⋅x¨t+Δt=AxtΔt⋅x˙tΔt2⋅x¨t,where the expression of amplification matrix A can be found in [21].The spectral radii of four implicit methods with different parameters are plotted in Figure1, 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).Figure 1
Spectral radii of various methods for undamped system (ξ=0).


## 3.2. 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. The AD and PE results of the considered time integration methods are plotted in Figures2 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. The Generalized-α method shows the lowest accuracy.Figure 2
The percentage AD results of various methods.Figure 3
The percentage PE results of various methods.


## 4. 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. Then an eleven-story shear building and a space truss subjected to seismic load are considered. The suggested algorithm parametersp=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. The 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.


### 4.1. Van der Pol’s Equation

A nonlinear equation called Van der Pol’ equation is considered here [26]. The equation is written as(18)1+x2x¨+c+x2x˙+x=0,where the constant c represents the damping coefficient of the Van der Pol’s equation and is negative here. The value of c and the initial conditions of the system are given.(19)c=−0.2,x0=−0.02,x˙0=0.The displacement, velocity, and acceleration results of the system are plotted in Figures4–9, respectively. The quasi-exact solutions shown in the figures are obtained from the explicit Bathe method [27] by using very small time step Δt=1×10−5s. As shown in the figures, different time steps are adopted for different time integration methods to ensure roughly the same time cost.Figure 4
The displacement results of various methods with different time steps for example 4.1.Figure 5
The velocity results of various methods with different time steps for example 4.1.Figure 6
The acceleration results of various methods with different time steps for example 4.1.Figure 7
The displacement results of various methods with the same time step for example 4.1.Figure 8
The velocity results of various methods with the same time step for example 4.1.Figure 9
The acceleration results of various methods with the same time step for example 4.1.Figures4–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. The accuracy of the OTTBIF method is slightly lower than that of the Wen method (p=0.4). The 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. The 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/3or0.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.Table 1
Time cost of various methods in time history of 0–500 s, for example 4.1
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time step (s)0.120.360.240.360.360.360.36The elapsed time (s)120.344.397.165.199.495.495.4


### 4.2. A Bead Sliding on the Wire

In the fixed verticalx,y plane as shown in Figure 10, a bead of mass m slides on the smooth wire (the frictional force is neglected). The curve of wire can be expressed by the equation y=x2. The motion equation of bead in x direction is given.(20)1+4x2x¨+4xx˙2+2gx=0,where g is the acceleration of gravity. The initial conditions of the system are (0) = −5,x (0) = 0. The 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.Figure 10
A bead and the wire.Figure 11
The displacement results of bead for time duration 0–50s.Figure 12
The displacement results of bead for time duration 180–200 s.Figure 13
The velocity results of bead for time duration 180–200 s.Figure 14
The acceleration results of bead for time duration 180–200 s.Figure 15
The displacement results of various methods with the same time stepΔt=0.03s.Figure 16
The velocity results of various methods with the same time stepΔt=0.03s.Figure 17
The acceleration results of various methods with the same time stepΔt=0.03s.The displacement results of various methods for time duration 0–50 s are shown in Figure11. It is observed that the results of the Generalized-α method are divergent. The displacement, velocity, and acceleration responses of the bead for time duration 180–200 s are plotted in dslf-14, respectively. The 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. The 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. The 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.Table 2
The energy dissipation and time cost of various methods for example 4.2.
MethodBathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time step (s)0.020.030.030.030.03Elapsed time (s)837.7481.4737.51699.9699.9Energy dissipation percentage3.218.914.366.740.49


### 4.3. Nonlinear Structural Seismic Response Analysis

Here two nonlinear structural systems under seismic load are tested. The two considered structural systems are an eleven-story shear building and a space truss. The 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. The quasi-exact solution is obtained using the Bathe method (r=0.5) with very small time step Δt=1×10−5s. For composite time integration methods, the required loads at the substeps are obtained by employing linear interpolation of two adjacent seismic acquisition load values.


#### 4.3.1. An Eleven-Story Shear Building

The eleven-story shear building shown in Figure18 is subjected to the 1940 EI Centro North-South seismic load (see Figure 19). The system is idealized as a simple shear model and has 11 degrees of freedom. The lumped mass and shear stiffness of each story are illustrated in Figure 18. The 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]. The acquisition frequency of seismic load is 50 Hz, and we intercept 30 s of the seismic load for dynamic analysis. The curve of seismic acceleration load in terms of time is shown in Figure 19.Figure 18
Eleven-story shear building.Figure 19
The North-South direction of EI Centro seismic load curve.Figure 20
Bouc-Wen hysteresis model.The motion equation of linear structural system under seismic load can be expressed as [30](21)MX¨+CX˙+KX=−MRX¨g,where M, C, and K are mass, damping, and stiffness matrices of the considered system. Xg¨ 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.The motion equation of the nonlinear part of the system can be written as follows:(22)MX¨+CX˙+K+KpreX+αZ=−MRX¨g,(23)z˙i=Ax˙i−βu˙izi−γzix˙i,where α is the peak value of inelastic force and satisfies α=Qy1−kpost/kpre. Qy, kpre, and kpost are the total force at the critical point of yield, the stiffness before yield, and stiffness after yield, respectively; A, β, and γ satisfy A=2β=2γ=kpre/Qy. The hysteretic restoring force Z of each time step is solved by Runge-Kutta method [31].The hysteretic restoring force curves of the bottom story are shown in Figure21. 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. The 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. The elapsed time of the Wen method and the OTTBIF method with Δt=0.024s is slightly 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).Figure 21
The hysteretic force of the bottom element.Figure 22
The displacement results and relative errors of top story obtained from various methods. (a) Displacement results of 0–10 s. (b) The relative errors between 4 and 7s.
(a)(b)Figure 23
The velocity results and relative errors of top story obtained from various methods: (a) Velocity results of 0–10 s. (b) The relative errors between 4 and 7s.
(a)(b)Figure 24
The displacement errors of top story obtained from various methods with the same time stepΔt=0.024s.Figure 25
The velocity errors of top story obtained from various methods with the same time stepΔt=0.024s.Table 3
Time cost of various methods in time history of 0–30 s for example 4.3.1.
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time Step (s)0.0080.0240.0160.0240.0240.0240.024Elapsed time (s)0.4300.1520.2210.1600.1840.1800.180Different 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. The 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. Therefore, the composite method can greatly reduce time cost for nonlinear hysteretic force calculation.This example verifies effectiveness of the Wen method and other composite time integration methods for structural nonlinear seismic response analysis. The three-step Wen method and the OTTBIF method have higher efficiency than the Generalized-α and Bathe methods.


#### 4.3.2. A Space Truss

As shown in Figure26, a space truss is considered. The total height of the structural system is 84 m, including 8 m for the first story and 76 m for the next 19 stories. The 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. The total degree of freedom of the structural system is 810. The elastic modulus of all linear elements is 2.1 × 1011 N/m2. The mass of all elements is 7850 kg/m2, and the cross-sectional area is 25 cm2. All elements are connected by rigid nodes. The 16 cross braces (marked in blue) in the first story are nonlinear elements. The Bouc-Wen hysteresis model [28] is adopted for these elements, and the expression of motion equation is given by equations (21)–(23). The nonlinear parameters kpre and kpost 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 Qy is adopted as 4770 N which is far less than that in [32]. Likewise, the nonlinear parameters kpre and kpost 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 Qy  = 3860 N is adopted to obtain strong nonlinearity. This structure system is also used in [32]; however, in this analysis, bar element is used and different parameters are selected.Figure 26
A space truss.The system is subjected to East-West (shown in Figure27), 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. The dynamic analysis in point A (seen in Figure 26) is considered. The 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. The displacement and velocity results (in time duration 8–10 s) of point A using different Δt for various methods are, respectively, plotted in Figures 30 and 31. The displacement and velocity errors of different methods with the same time step Δt=0.015s are, respectively, plotted in Figures 32 and 33. Time cost of various methods with different time steps in time history of 0–30s is listed in Table 4.Figure 27
The East-West direction of EI Centro seismic load curve.Figure 28
The Up-Down direction of EI Centro seismic load curve.Figure 29
The hysteretic force of all 16 nonlinear elements.Figure 30
The displacement results of point A and its relative error from various methods. (a) Displacement results of 8–10 s. (b) The relative error results of displacement of 8–10 s.
(a)(b)Figure 31
The velocity results of point A and the relative errors from various methods. (a) Velocity results of 8–10 s. (b) The relative error results of velocity of 8–10 s.
(a)(b)Figure 32
The displacement errors of various methods with the same time stepΔt=0.015s.Figure 33
The velocity errors of various methods with the same time stepΔt=0.015s.Table 4
Time cost of various methods in time history of 0–30s, for example 4.3.2.
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time Step (s)0.0050.0150.010.0150.0150.0150.015Elapsed time (s)180.158.491.058.988.485.185.1Figures30 and 31 illustrate that four methods using different time steps have acceptable accuracy. However, the Generalized-α method has the largest errors. The Bathe method shows slightly larger error than that of the Wen method and OTTBIF method. The Wen method has close calculation accuracy as the OTTBIF method. Time cost of four methods follows Generalized-αΔt=0.005 > Bathe Δt=0.01 > OTTBIF Δt=0.015 > Wen Δt=0.015. Therefore, 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. This example verifies effectiveness of the Wen method and other composite time integration methods for nonlinear seismic response analysis of space truss.


## 4.1. Van der Pol’s Equation

A nonlinear equation called Van der Pol’ equation is considered here [26]. The equation is written as(18)1+x2x¨+c+x2x˙+x=0,where the constant c represents the damping coefficient of the Van der Pol’s equation and is negative here. The value of c and the initial conditions of the system are given.(19)c=−0.2,x0=−0.02,x˙0=0.The displacement, velocity, and acceleration results of the system are plotted in Figures4–9, respectively. The quasi-exact solutions shown in the figures are obtained from the explicit Bathe method [27] by using very small time step Δt=1×10−5s. As shown in the figures, different time steps are adopted for different time integration methods to ensure roughly the same time cost.Figure 4
The displacement results of various methods with different time steps for example 4.1.Figure 5
The velocity results of various methods with different time steps for example 4.1.Figure 6
The acceleration results of various methods with different time steps for example 4.1.Figure 7
The displacement results of various methods with the same time step for example 4.1.Figure 8
The velocity results of various methods with the same time step for example 4.1.Figure 9
The acceleration results of various methods with the same time step for example 4.1.Figures4–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. The accuracy of the OTTBIF method is slightly lower than that of the Wen method (p=0.4). The 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. The 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/3or0.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.Table 1
Time cost of various methods in time history of 0–500 s, for example 4.1
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time step (s)0.120.360.240.360.360.360.36The elapsed time (s)120.344.397.165.199.495.495.4


## 4.2. A Bead Sliding on the Wire

In the fixed verticalx,y plane as shown in Figure 10, a bead of mass m slides on the smooth wire (the frictional force is neglected). The curve of wire can be expressed by the equation y=x2. The motion equation of bead in x direction is given.(20)1+4x2x¨+4xx˙2+2gx=0,where g is the acceleration of gravity. The initial conditions of the system are (0) = −5,x (0) = 0. The 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.Figure 10
A bead and the wire.Figure 11
The displacement results of bead for time duration 0–50s.Figure 12
The displacement results of bead for time duration 180–200 s.Figure 13
The velocity results of bead for time duration 180–200 s.Figure 14
The acceleration results of bead for time duration 180–200 s.Figure 15
The displacement results of various methods with the same time stepΔt=0.03s.Figure 16
The velocity results of various methods with the same time stepΔt=0.03s.Figure 17
The acceleration results of various methods with the same time stepΔt=0.03s.The displacement results of various methods for time duration 0–50 s are shown in Figure11. It is observed that the results of the Generalized-α method are divergent. The displacement, velocity, and acceleration responses of the bead for time duration 180–200 s are plotted in dslf-14, respectively. The 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. The 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. The 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.Table 2
The energy dissipation and time cost of various methods for example 4.2.
MethodBathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time step (s)0.020.030.030.030.03Elapsed time (s)837.7481.4737.51699.9699.9Energy dissipation percentage3.218.914.366.740.49


## 4.3. Nonlinear Structural Seismic Response Analysis

Here two nonlinear structural systems under seismic load are tested. The two considered structural systems are an eleven-story shear building and a space truss. The 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. The quasi-exact solution is obtained using the Bathe method (r=0.5) with very small time step Δt=1×10−5s. For composite time integration methods, the required loads at the substeps are obtained by employing linear interpolation of two adjacent seismic acquisition load values.


### 4.3.1. An Eleven-Story Shear Building

The eleven-story shear building shown in Figure18 is subjected to the 1940 EI Centro North-South seismic load (see Figure 19). The system is idealized as a simple shear model and has 11 degrees of freedom. The lumped mass and shear stiffness of each story are illustrated in Figure 18. The 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]. The acquisition frequency of seismic load is 50 Hz, and we intercept 30 s of the seismic load for dynamic analysis. The curve of seismic acceleration load in terms of time is shown in Figure 19.Figure 18
Eleven-story shear building.Figure 19
The North-South direction of EI Centro seismic load curve.Figure 20
Bouc-Wen hysteresis model.The motion equation of linear structural system under seismic load can be expressed as [30](21)MX¨+CX˙+KX=−MRX¨g,where M, C, and K are mass, damping, and stiffness matrices of the considered system. Xg¨ 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.The motion equation of the nonlinear part of the system can be written as follows:(22)MX¨+CX˙+K+KpreX+αZ=−MRX¨g,(23)z˙i=Ax˙i−βu˙izi−γzix˙i,where α is the peak value of inelastic force and satisfies α=Qy1−kpost/kpre. Qy, kpre, and kpost are the total force at the critical point of yield, the stiffness before yield, and stiffness after yield, respectively; A, β, and γ satisfy A=2β=2γ=kpre/Qy. The hysteretic restoring force Z of each time step is solved by Runge-Kutta method [31].The hysteretic restoring force curves of the bottom story are shown in Figure21. 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. The 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. The elapsed time of the Wen method and the OTTBIF method with Δt=0.024s is slightly 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).Figure 21
The hysteretic force of the bottom element.Figure 22
The displacement results and relative errors of top story obtained from various methods. (a) Displacement results of 0–10 s. (b) The relative errors between 4 and 7s.
(a)(b)Figure 23
The velocity results and relative errors of top story obtained from various methods: (a) Velocity results of 0–10 s. (b) The relative errors between 4 and 7s.
(a)(b)Figure 24
The displacement errors of top story obtained from various methods with the same time stepΔt=0.024s.Figure 25
The velocity errors of top story obtained from various methods with the same time stepΔt=0.024s.Table 3
Time cost of various methods in time history of 0–30 s for example 4.3.1.
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time Step (s)0.0080.0240.0160.0240.0240.0240.024Elapsed time (s)0.4300.1520.2210.1600.1840.1800.180Different 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. The 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. Therefore, the composite method can greatly reduce time cost for nonlinear hysteretic force calculation.This example verifies effectiveness of the Wen method and other composite time integration methods for structural nonlinear seismic response analysis. The three-step Wen method and the OTTBIF method have higher efficiency than the Generalized-α and Bathe methods.


### 4.3.2. A Space Truss

As shown in Figure26, a space truss is considered. The total height of the structural system is 84 m, including 8 m for the first story and 76 m for the next 19 stories. The 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. The total degree of freedom of the structural system is 810. The elastic modulus of all linear elements is 2.1 × 1011 N/m2. The mass of all elements is 7850 kg/m2, and the cross-sectional area is 25 cm2. All elements are connected by rigid nodes. The 16 cross braces (marked in blue) in the first story are nonlinear elements. The Bouc-Wen hysteresis model [28] is adopted for these elements, and the expression of motion equation is given by equations (21)–(23). The nonlinear parameters kpre and kpost 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 Qy is adopted as 4770 N which is far less than that in [32]. Likewise, the nonlinear parameters kpre and kpost 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 Qy  = 3860 N is adopted to obtain strong nonlinearity. This structure system is also used in [32]; however, in this analysis, bar element is used and different parameters are selected.Figure 26
A space truss.The system is subjected to East-West (shown in Figure27), 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. The dynamic analysis in point A (seen in Figure 26) is considered. The 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. The displacement and velocity results (in time duration 8–10 s) of point A using different Δt for various methods are, respectively, plotted in Figures 30 and 31. The displacement and velocity errors of different methods with the same time step Δt=0.015s are, respectively, plotted in Figures 32 and 33. Time cost of various methods with different time steps in time history of 0–30s is listed in Table 4.Figure 27
The East-West direction of EI Centro seismic load curve.Figure 28
The Up-Down direction of EI Centro seismic load curve.Figure 29
The hysteretic force of all 16 nonlinear elements.Figure 30
The displacement results of point A and its relative error from various methods. (a) Displacement results of 8–10 s. (b) The relative error results of displacement of 8–10 s.
(a)(b)Figure 31
The velocity results of point A and the relative errors from various methods. (a) Velocity results of 8–10 s. (b) The relative error results of velocity of 8–10 s.
(a)(b)Figure 32
The displacement errors of various methods with the same time stepΔt=0.015s.Figure 33
The velocity errors of various methods with the same time stepΔt=0.015s.Table 4
Time cost of various methods in time history of 0–30s, for example 4.3.2.
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time Step (s)0.0050.0150.010.0150.0150.0150.015Elapsed time (s)180.158.491.058.988.485.185.1Figures30 and 31 illustrate that four methods using different time steps have acceptable accuracy. However, the Generalized-α method has the largest errors. The Bathe method shows slightly larger error than that of the Wen method and OTTBIF method. The Wen method has close calculation accuracy as the OTTBIF method. Time cost of four methods follows Generalized-αΔt=0.005 > Bathe Δt=0.01 > OTTBIF Δt=0.015 > Wen Δt=0.015. Therefore, 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. This example verifies effectiveness of the Wen method and other composite time integration methods for nonlinear seismic response analysis of space truss.


## 4.3.1. An Eleven-Story Shear Building

The eleven-story shear building shown in Figure18 is subjected to the 1940 EI Centro North-South seismic load (see Figure 19). The system is idealized as a simple shear model and has 11 degrees of freedom. The lumped mass and shear stiffness of each story are illustrated in Figure 18. The 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]. The acquisition frequency of seismic load is 50 Hz, and we intercept 30 s of the seismic load for dynamic analysis. The curve of seismic acceleration load in terms of time is shown in Figure 19.Figure 18
Eleven-story shear building.Figure 19
The North-South direction of EI Centro seismic load curve.Figure 20
Bouc-Wen hysteresis model.The motion equation of linear structural system under seismic load can be expressed as [30](21)MX¨+CX˙+KX=−MRX¨g,where M, C, and K are mass, damping, and stiffness matrices of the considered system. Xg¨ 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.The motion equation of the nonlinear part of the system can be written as follows:(22)MX¨+CX˙+K+KpreX+αZ=−MRX¨g,(23)z˙i=Ax˙i−βu˙izi−γzix˙i,where α is the peak value of inelastic force and satisfies α=Qy1−kpost/kpre. Qy, kpre, and kpost are the total force at the critical point of yield, the stiffness before yield, and stiffness after yield, respectively; A, β, and γ satisfy A=2β=2γ=kpre/Qy. The hysteretic restoring force Z of each time step is solved by Runge-Kutta method [31].The hysteretic restoring force curves of the bottom story are shown in Figure21. 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. The 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. The elapsed time of the Wen method and the OTTBIF method with Δt=0.024s is slightly 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).Figure 21
The hysteretic force of the bottom element.Figure 22
The displacement results and relative errors of top story obtained from various methods. (a) Displacement results of 0–10 s. (b) The relative errors between 4 and 7s.
(a)(b)Figure 23
The velocity results and relative errors of top story obtained from various methods: (a) Velocity results of 0–10 s. (b) The relative errors between 4 and 7s.
(a)(b)Figure 24
The displacement errors of top story obtained from various methods with the same time stepΔt=0.024s.Figure 25
The velocity errors of top story obtained from various methods with the same time stepΔt=0.024s.Table 3
Time cost of various methods in time history of 0–30 s for example 4.3.1.
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time Step (s)0.0080.0240.0160.0240.0240.0240.024Elapsed time (s)0.4300.1520.2210.1600.1840.1800.180Different 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. The 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. Therefore, the composite method can greatly reduce time cost for nonlinear hysteretic force calculation.This example verifies effectiveness of the Wen method and other composite time integration methods for structural nonlinear seismic response analysis. The three-step Wen method and the OTTBIF method have higher efficiency than the Generalized-α and Bathe methods.


## 4.3.2. A Space Truss

As shown in Figure26, a space truss is considered. The total height of the structural system is 84 m, including 8 m for the first story and 76 m for the next 19 stories. The 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. The total degree of freedom of the structural system is 810. The elastic modulus of all linear elements is 2.1 × 1011 N/m2. The mass of all elements is 7850 kg/m2, and the cross-sectional area is 25 cm2. All elements are connected by rigid nodes. The 16 cross braces (marked in blue) in the first story are nonlinear elements. The Bouc-Wen hysteresis model [28] is adopted for these elements, and the expression of motion equation is given by equations (21)–(23). The nonlinear parameters kpre and kpost 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 Qy is adopted as 4770 N which is far less than that in [32]. Likewise, the nonlinear parameters kpre and kpost 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 Qy  = 3860 N is adopted to obtain strong nonlinearity. This structure system is also used in [32]; however, in this analysis, bar element is used and different parameters are selected.Figure 26
A space truss.The system is subjected to East-West (shown in Figure27), 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. The dynamic analysis in point A (seen in Figure 26) is considered. The 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. The displacement and velocity results (in time duration 8–10 s) of point A using different Δt for various methods are, respectively, plotted in Figures 30 and 31. The displacement and velocity errors of different methods with the same time step Δt=0.015s are, respectively, plotted in Figures 32 and 33. Time cost of various methods with different time steps in time history of 0–30s is listed in Table 4.Figure 27
The East-West direction of EI Centro seismic load curve.Figure 28
The Up-Down direction of EI Centro seismic load curve.Figure 29
The hysteretic force of all 16 nonlinear elements.Figure 30
The displacement results of point A and its relative error from various methods. (a) Displacement results of 8–10 s. (b) The relative error results of displacement of 8–10 s.
(a)(b)Figure 31
The velocity results of point A and the relative errors from various methods. (a) Velocity results of 8–10 s. (b) The relative error results of velocity of 8–10 s.
(a)(b)Figure 32
The displacement errors of various methods with the same time stepΔt=0.015s.Figure 33
The velocity errors of various methods with the same time stepΔt=0.015s.Table 4
Time cost of various methods in time history of 0–30s, for example 4.3.2.
MethodG-α (ρ∞=0)G-α (ρ∞=0)Bathe (r=0.5)Bathe (r=0.5)OTTBIF (ρ∞=0)Wen (p=1/3)Wen (p=0.4)Time Step (s)0.0050.0150.010.0150.0150.0150.015Elapsed time (s)180.158.491.058.988.485.185.1Figures30 and 31 illustrate that four methods using different time steps have acceptable accuracy. However, the Generalized-α method has the largest errors. The Bathe method shows slightly larger error than that of the Wen method and OTTBIF method. The Wen method has close calculation accuracy as the OTTBIF method. Time cost of four methods follows Generalized-αΔt=0.005 > Bathe Δt=0.01 > OTTBIF Δt=0.015 > Wen Δt=0.015. Therefore, 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. This example verifies effectiveness of the Wen method and other composite time integration methods for nonlinear seismic response analysis of space truss.


## 5. Conclusions

In this work, the Wen method is applied to the analysis of nonlinear dynamics. The calculation procedure of the Wen method for nonlinear analysis is proposed. The 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. The conclusion of this study can be summarized as follows.(1)
The Wen method shows concise calculation procedure for nonlinear dynamics analysis.(2)
It is observed that, for the Wen method, the parameterp values of desirable algorithmic properties (including AD, PE, and analytical accuracy) for linear dynamics are also suggested for nonlinear dynamic analysis.(3)
Numerical solutions of three nonlinear systems illustrate that the Wen method (thep=0.4 case) achieves desirable accuracy and efficiency for nonlinear dynamic examples and structural nonlinear seismic response analysis.(4)
In general, the theoretical analysis and numerical simulation demonstrate that the Wen method is a desirable candidate for nonlinear dynamic problems.



---

*Source: 6442260-2021-12-21.xml*