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

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

#### 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 bywhere is the mass matrix, is the vector of internal forces, and is the vector of external forces. , , and represent the displacement, velocity, and acceleration vectors, respectively. The initial conditions are given bywhere and are, respectively, initial displacement and velocity vectors required to be given before solving the equation. 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 , , and , respectively. The parameter satisfies the condition 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 arewhere the subscript ‘’ defines the value at time .

The formulas of the second substep arewhere , , .

The formulas of the third substep arewhere , , , .

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 , which is solved by Newton-Raphson iteration method. After simplification, we obtain ()where . Once the displacement vector (i.e., ) is obtained, the velocity and acceleration vectors can be solved with equations (3) and (4).

Likewise, for the second substep, we obtainwhere . After is obtained, and can be solved with equations (6) and (7).

As for the third substep, we obtainwhere .

The minimal iteration number for each substep is determined by the criterionwhere the tolerance is .

This method is a single parameter method, and all algorithmic properties are related to the algorithm parameter . 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 , , and cases), the Bathe method [16] (the and cases), and the OTTBIF method [24] (the , , and 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:where , , and are the damping ratio, the undamped angular frequency of the system, and the force excitation, respectively. satisfies , where is the period of the system. By setting condition and in equation (16), the recurrence formula of the Wen method for any undamped SDOF system can be obtained.where the expression of amplification matrix can be found in [21].

The spectral radii of four implicit methods with different parameters are plotted in Figure 1, where . 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-*α* () and case OTTBIF ().

##### 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 Figures 2 and 3, respectively. For the Wen method, it can achieve the lowest AD and PE in the case of . 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.

#### 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 parameters and of the Wen method for linear dynamics are adopted for analysis. For comparison, the single-step Generalized-*α* method (the case) [11], the two substeps’ Bathe method (the case) [16], and the three substeps’ OTTBIF method (the 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 aswhere 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.

The displacement, velocity, and acceleration results of the system are plotted in Figures 4–9, respectively. The quasi-exact solutions shown in the figures are obtained from the explicit Bathe method [27] by using very small time step . 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 case of the Wen method shows highest solution accuracy. The accuracy of the OTTBIF method is slightly lower than that of the Wen method (). The Bathe method and the Wen method () have desirable accuracy, but it is lower than that of the OTTBIF method and the Wen method (). 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 () is much higher than the Generalized-*α* method, the Bathe method, and the Wen method (). In Table 1, time cost of all methods in time history of 0–500 s is given. The Generalized-*α* method with and , respectively, has the longest and shortest elapsed time among all considered methods. Time costs of the Bathe method (), the OTTBIF method (), and the Wen method () are very close. Combined with the calculation accuracy, the Wen method () and the OTTBIF method have higher computational efficiency than other methods.

##### 4.2. A Bead Sliding on the Wire

In the fixed vertical 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 . The motion equation of bead in direction is given.where is the acceleration of gravity. The initial conditions of the system are (0) = −5*,* (0) = 0. The quasi-exact solutions are calculated by the explicit Bathe method with time step Δ = 1 × 10^{−5} s. In Figures 11–14, Δ = 0.01 s and 0.02 s are, respectively, selected for the Generalized-*α* and Bathe methods, while Δ = 0.03 s is adopted for the OTTBIF and Wen methods. In Figures 15–17, all methods adopt the same time step Δ = 0.03 s.

The 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. 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 achieves highest accuracy among all considered methods. In addition, accuracy of the Wen method () 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 . The results illustrate that the Wen method in case has the lowest energy decay and consumes less time cost than the OTTBIF method. The Bathe method () needs less computation time than other methods, but begets the largest proportion of energy dissipation. In Figures 15–17, when the same time step is used for all methods, the Wen method () displays far higher solution accuracy than other methods.

##### 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 () with very small time step . 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 Figure 18 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.

The motion equation of linear structural system under seismic load can be expressed as [30]where , , and are mass, damping, and stiffness matrices of the considered system. is the vector of ground motion acceleration. 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:where is the peak value of inelastic force and satisfies . , , and are the total force at the critical point of yield, the stiffness before yield, and stiffness after yield, respectively; , , and satisfy . The hysteretic restoring force of each time step is solved by Runge-Kutta method [31].

The hysteretic restoring force curves of the bottom story are shown in 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. 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 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 is slightly longer than the Bathe method and the Generalized-*α* method with the same time step , while the Bathe method () and the Generalized-*α* method () have longer elapsed time than the Wen method () and the OTTBIF method ().

**(a)**

**(b)**

**(a)**

**(b)**

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. 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 Figure 26, 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 × 10^{11} N/m^{2}. The 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. 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 and 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 is adopted as 4770 N which is far less than that in [32]. Likewise, the nonlinear parameters and 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 = 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.

The 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. 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 for various methods are, respectively, plotted in Figures 30 and 31. The displacement and velocity errors of different methods with the same time step 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.

**(a)**

**(b)**

**(a)**

**(b)**

Figures 30 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-*α* > Bathe > OTTBIF > Wen . 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 parameter *p* 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 (the 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.

#### Data Availability

The raw/processed data can be obtained by contacting the corresponding author by e-mail (De Zhou [email protected]).

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was substantially supported by the National Natural Science Foundation of China (no.12072375).