Abstract

This paper presents a new robust and efficient time integration algorithm suitable for various complex nonlinear structural dynamic finite element problems. Based on the idea of composition, the three-point backward difference formula and a generalized central difference formula are combined to constitute the implicit algorithm. Theoretical analysis indicates that the composite algorithm is a single-solver algorithm with satisfactory accuracy, unconditional stability, and second-order convergence rate. Moreover, without any additional parameters, the composite algorithm maintains a symmetric effective stiffness matrix and the computational cost is the same as that of the trapezoidal rule. And more merits of the proposed algorithm are revealed through several representative finite element examples by comparing with analytical solutions or solutions provided by other numerical techniques. Results show that not only the linear stiff problem but also the nonlinear problems involving nonlinearities of geometry, contact, and material can be solved efficiently and successfully by this composite algorithm. Thus the prospect of its implementation in existing finite element codes can be foreseen.

1. Introduction

Modeling of engineering problems leads in many cases to ordinary and partial differential equations which are often of nonlinear nature. A powerful tool to solve these differential equations is the finite element method which was developed 50 years ago and has achieved many important breakthroughs including various solution algorithms [14]. And it is believed that the use of finite element method will continue to increase in the future because it provides a more rapid and less expensive way to evaluate design concepts and design details [5, 6].

Dynamic transient analysis is necessary for many engineering problems, such as earthquake analysis or the vibration of structures. In particular when structures possess nonlinear features, the step-by-step time integration scheme may be the only tool available for dynamic transient analysis. The widely used step-by-step integration algorithms can be grouped into two main categories: explicit and implicit [7]. The explicit algorithm is appropriate to the simulation of wave propagation problems where high frequency responses are dominated, and the most crucial difficulty in explicit methods is the numerical instability which may be time consuming because a very small time step should be utilized [8, 9]; thus giant computers and parallel computing technologies are usually necessary. In the view of the fact that high frequency responses which are even totally spurious and are caused by spatial discretization are not the dominated response in long-term structural dynamic analysis, implicit algorithms which allow for larger time steps turn out to be the more efficient strategies [10].

Many researchers have been dedicated to the development of effective unconditionally stable implicit algorithms [11], and the research efforts can be classified, in essence, into three categories [12]. The first category introduces numerical dissipation by using various parameters to satisfy the energy criterion. These well-known methods, including the generalized-α method [13], Wilson-θ method [14], HHT-α [15], and WBZ-α [16], have one thing in common: one or more parameters must be chosen before the solution procedure. Obviously personal judgement and problem specifications will affect the accuracy and analysis will not be automatic [17, 18]. In particular, for nonlinear problems, it is difficult for an analyst to identify the appropriate parameters. In practical engineering applications, integration schemes with no parameter are more popular [10, 19]. In the second category, energy conservation criterion is enforced by employing Lagrange multipliers. It will definitely increase the unknowns and result in nonsymmetric tangent stiffness matrices, thus leading to costly computation [20]. Despite all these, failure in Newton-Raphson iteration of equilibrium had been found in some cases [21]. And the last category obtains conservation algorithmically, such as the energy-momentum method proposed by Simo and Tarnow [22]. However, the energy-momentum method lacks generality since different schemes have to be developed for each different system of equations to be solved [23] and it must be modified to damp out the spurious high frequencies at the price of accuracy loss [24, 25]. Besides, the tangential stiffness matrices of all the generalized energy-momentum methods are unsymmetric [12].

Based on the above summary here and more developments described elsewhere [26, 27], we can see that a time integration algorithm suitable for various complex nonlinear dynamic structural systems is lacking, which should inherit the desirable properties including unconditional stability, at least second-order accuracy, ease of use (no additional parameter), ease of implementation, and high computational efficiency. Some delightful latest achievements are indicating that the conception of composition can be the outlet [2830]. Generally speaking, there are two kinds of composition regimes.

The first regime is assembling different independent integration schemes in several substeps into a composite algorithm in one whole time step, whose solution depends on both these independent schemes and the partition of these substeps. High accuracy of the composite algorithm relies on the equilibrium equation satisfaction at each moment of the substeps. The most representative one is the Bathe method, which combines the trapezoidal rule with a three-point backward difference method in two equal substeps [10, 28, 31]. Unfortunately, the computational cost of the Bathe algorithm is about twice as expensive as the trapezoidal rule due to more iterations in the substeps [28]. Moreover, generalizing the Bathe scheme to more substeps receives algorithms with less period elongation but amplitude decay and computational cost are increasing [28]. On the other hand, different segmentations of the substeps require additional parameters to be introduced [32].

The second regime is the composition of different difference formulae in one whole time step to inherit their advantages [30, 33, 34]. For example, Liu et al. have proposed an efficient time integration algorithm by compositing the two-point and three-point backward Euler formulae [30], which presents good stability for nonlinear structural dynamic problems, with no parameter and the same computational efficiency as the trapezoidal rule. However, the accuracy of this algorithm is less than satisfactory for its high numerical damping in the intermediate-frequency range.

In this paper, we devised a new efficient time integration algorithm based on the second composition regime. The three-point backward difference formula was retained due to its good performance on stiff problems [30], and at the same time, in order to achieve satisfactory accuracy, a generalized central difference formula was derived to calculate the acceleration term. It will be displayed that, maintaining all those above mentioned desired properties for nonlinear structural analysis, this new proposed composite method has gained significant improved accuracy comparing to the backward Euler algorithm. In the following, the composite time integration method is introduced in Section 2, the properties such as accuracy and stability are evaluated theoretically, and also a linear stiff system will be investigated deeply with the help of the new composite algorithm. In Section 3, to illustrate its potential and to assess its accuracy and efficiency on nonlinear structural dynamic analysis, three finite element examples, involving nonlinearities on geometry, contact, and material, respectively, are discussed. Finally, some concluding remarks are summarized in Section 4.

2. The Composite Time Integration Algorithm

In this section the detailed derivation process of the composite scheme is first presented, and then the accuracy, stability, and convergence rate of the scheme are evaluated theoretically. In addition, a typical stiff linear system is analyzed to probe valuable insights into the properties of the scheme.

2.1. The Basic Equations of the Integration Algorithm

For nonlinear transient structural dynamic problems, the governing equations of equilibrium discretized by finite elements take the common second-order ordinary differential form aswith the initial conditionswhere is the mass matrix, , , and are the nodal displacement, velocity, and acceleration, respectively, and are the nodal damping force and the elastic force corresponding to the element internal stresses, is the externally applied load, and and are the given initial displacement and velocity, respectively. In linear cases, and , where and are the damping matrix and stiffness matrix.

The step-by-step time integration solves (1) and (2) at a set of discrete time points. Assume that the solutions are completely known up to time , and then the solution at time is to be calculated. The equilibrium equation at time can be rewritten as

For the good performance on stiff problems, the second-order accuracy three-point Euler backward difference formula (4) is reserved here to combine with the generalized central difference formula (5), whose derivation will be given briefly below:

The acceleration can be obtained by the classical forward Euler integration formula:

To obtain an algorithm with higher accuracy, substitution of the common central difference formulae of the acceleration and its first-order derivative, as shown in (7), into (6) yields the generalized central difference formula (5). Consider

Since the nodal damping forces and the internal forces depend nonlinearly on the current values, the linearization of (3) yields the following matrices:

Substituting (4) and (5) into linearized equation (3) and employing the displacement incrementyield the following Newton-Raphson iteration:whereare the effective stiffness matrix and effective load vector, respectively, with iteration number .  Clearly, the effective stiffness matrix is symmetric as long as the damping matrix possesses symmetry.

Due to (4) and (5) being three time points, the composite algorithm requires a starting step algorithm, which is recommended as follows:

Equation (13) results from the Newmark algorithm with . Using this particular value of parameter leads to the same effective stiffness matrix of the starting step as the subsequent time steps, but the effective load of (12) should be replaced with

Obviously, the new composite time integration algorithm is of second-order accuracy with no additional artificial parameter and only one set of nonlinear equations with symmetric effective stiffness matrix needs to be solved at each time step; hence the computational effort is the same as that of the trapezoidal rule.

2.2. Stability Analysis

In this subsection, the stability of the proposed algorithm for linear problems is performed theoretically while the stability performance for nonlinear cases will be investigated numerically in Section 3.

Now, we consider the following single degree-of-freedom system:where is the free vibration natural frequency and is the damping ratio.

Considering (4) and (5) at any three successive time points results inwhere and is the transfer matrix expressed aswhere .

Figure 1 shows the spectral radii of the transfer matrices of the proposed method without damping as well as those of the trapezoidal rule [35], the generalized-α method [13], the Bathe method [31], and backward Euler method [30]. It can be seen that all of these mentioned methods except the trapezoidal rule perform diminishing spectral radii with the increase of ratio value. The proposed method presents the spectral radius of 0.994 at the ratio , with almost no amplitude decrease. As the ratio increases, the spectral radius of the proposed method rapidly reduces to zero, indicating desired numerical damping in the high frequency range, which is an expected property for practical engineering analysis. In particular, the spectral radius of the proposed method is smaller than that of the Bathe method but larger than that of the backward Euler method as well as the generalized-α method. Therefore, the proposed composite method is unconditionally stable in linear case.

2.3. Accuracy Analysis

To evaluate the accuracy of the proposed composite method, the following spring-mass system without damping is considered:

The initial conditions are , , and .

The analytical solution of (19) is . Figures 2 and 3 show the period elongation and amplitude decay of the proposed composite scheme, respectively. For comparison, the corresponding results of the trapezoidal rule [35], the generalized-α method [13], the Bathe method [31], and backward Euler method [30] are also presented. Clearly, when the ratio is no more than 0.01, all methods present satisfactory accuracy. With the increase of the ratio , the new composite method maintains high accuracy, much better than the backward Euler difference method under the same efficiency, just a little worse than the Bathe method which costs more calculation effort.

2.4. Convergence Rate Analysis

In order to evaluate the convergence rate, we consider the following linear resonance problem as discussed in [36]:with the initial conditions , .

The exact solutions of this system for displacement, velocity, and acceleration, respectively, are

The circular frequency is assigned to  rad/s corresponding to the period  s. Figure 4 shows the relative error curves of displacement, velocity, and acceleration at the time of 1.0 s. Clearly, this new composite method is second-order and convergent in terms of the displacement, velocity, and acceleration.

2.5. Linear Stiff System Demonstration

In this subsection, a linear spring-mass model analyzed by Bathe and Noh [31], which seems simple but can represent complex practical structural system, is investigated to test the proposed composite algorithm.

The governing equations of the stiff linear system, with parameters shown in Figure 5, are

It can be seen that the two springs have quite different stiffness. The left stiff spring represents rigid connections in the complex models, while the right flexible one represents the flexible parts. Indeed, the high stiffness parts used in practice usually have little physical meaning other than to provide constraints, whose response would naturally not be included in the mode superposition solution. In order to verify the high frequency damping property of the new proposed composite method, the initial conditions , , , and  rad/s are supposed, and the time step used is 0.1 s. Hence we have , , and , where  s,  s are the two distinct natural periods of the system and  s is the period of the prescribed displacement boundary condition.

Comparing to the analytical solutions, the displacement, velocity, and acceleration responses shown in Figure 6 demonstrate that the proposed composite method performs very well on this stiff problem without adjustment of any parameter. However, with the same computational cost, the trapezoidal rule cannot even achieve stable result because of lacking appropriate numerical damping [31]; as shown in Figure 7, the acceleration response of node 1 displays distinct instability. As in practical finite element simulations, not only the rigid connection but also the spatial finite element discretization can bring in spurious high frequencies; thus a favorable algorithm should have appropriate numerical damping in the high frequency range [19].

3. Nonlinear Dynamic Finite Element Examples

In practical engineering, nonlinear behavior is more common. According to the different sources of nonlinearities, they can be divided into three categories: geometrical nonlinearity, contact nonlinearity, and material nonlinearity. In this section, three typical finite element examples involving each nonlinearity are solved to prove the precision and effectiveness of the proposed composite time integration algorithm.

3.1. Geometrical Nonlinear Example

The stiff single pendulum model, which involves large displacements and rotations, is a classical geometrical nonlinear example to demonstrate the ability of a new time integration algorithm in solving nonlinear problems, and the failure of the trapezoidal rule in this case has been illustrated by many scholars [10, 30, 37].

Figure 8 shows the geometry, material, boundary, and initial conditions of the single pendulum model, whose analytical rotational period is 2.4777 s and the total energy is 149 J. We solve this nonlinear problem using geometrical nonlinear finite elements with the time step of 0.01 s for as long as 2500 s (more than 1000 cycles). The solutions of the first 15 seconds and the last 15 seconds are shown in Figure 9. The total amplitude decay of the whole 1000 cycles is only about 0.0036% and the total period elongation is less than 0.0928%. We also ran this problem using a much larger time step of 0.5 s, and still a stable solution but of course with less accuracy was achieved, which evidently proved the robust stability of the new composite scheme for geometrical nonlinear problems.

This problem was also solved using the Bathe method under the same calculation conditions and convergence precisions. Results showed that at the time of 150 s, the energy decay of the Bathe method was 0.0248%, much larger than that of the composite method (0.0005%). And on average, the Bathe method needed 6 iterations while the composite method only needed 4 iterations in one time step to reach the convergence precisions. In addition, [30, 37] have proved that the trapezoidal rule cannot achieve stable results for this nonlinear problem. Thus the new composite method is of high accuracy and efficiency for geometrical nonlinear problems.

3.2. Contact Nonlinear Example

Another particular difficult nonlinear behavior is the contact between two or more bodies, which can be considered stemming from boundary nonlinearity. In this subsection, a benchmark test about center-to-center collision without friction between two identical elastic bars [38, 39] using Lagrange multiplier to impose contact constraints is performed with time step of 0.008 s. Each bar is 100 meters long and discretized by 50 equal eight-node 3D finite elements with the same size of 1 m × 1 m × 2 m. The detailed material parameters and initial conditions are shown in Figure 10.

The longitudinal velocity responses of the contacting pair nodes of each bar are illustrated in Figure 11, which clearly shows that impact occurs at the time of 0.1 s, and since then the two contacting pair nodes move at the same velocity of 0.5 m/s. And the contacting pair nodes keep moving together until the stress waves caused by the impact propagate back to the contacting ends at the time of 2.1 s. After that, the two bars separate and accomplish the exchange of their initial velocities. One distinct thing is that steep fluctuations with an interval of 2.0 s appear on the velocity responses after separation; the reason is that the stress waves are propagating without loss back and forth at the elastic wave velocity of 100 m/s. As shown in Figures 12 and 13, there is some little oscillation in the contact pressure and the total energy, which are attributed to the fact that the impenetrability constraints are only imposed on displacements in our procedure for simplicity, which should also be imposed on the rates of displacement for better accuracy [4042]. As a whole, the composite method displays robust stability for dynamic contact nonlinear problems.

3.3. Material Nonlinear Example

Material nonlinearity, characterized by a nonlinear response function between stresses and strains or by a set of evolution equations, often occurs in structural engineering. For example, in seismic design of buildings, the lead-rubber bearing (LRB for short) is an effective device designed with nonlinear constitutive relationships to reduce the destructive transmission energy of the group motion to the superstructure. The LRB is made of alternate layers of rubber and steel plates with one or more lead plugs inserted into the predrilled holes, and the lead core deforms in shear providing both the initial rigidity and the bilinear response, which is crucial for seismic isolation [43].

In this example, a seven-story building isolated by the LRB is analyzed using the composite algorithm under EL-Centro, 1940 earthquake motion (North-South component with peak acceleration equal to 3.417 m/s2). Rayleigh damping consisting of mass matrix and stiffness matrix with coefficients both equal to 0.01 is also considered. For simplicity, the mass and the stiffness of each floor of the superstructure are kept constant, with detailed parameters shown in Figure 14. And some additional assumptions are made: floors of each story of the superstructure are assumed to be rigid; the force-deformation behavior of the LRB is considered bilinear defined by three parameters , , and (as shown in Figure 14(b)), which denote the yield strength, postyield stiffness, and yield displacement, respectively.

Since it is difficult to obtain the analytical solution, results of the composite time integration algorithm are compared and validated with that of the commercial finite element software Ansys (version of 12.1) using Newmark integration method (, ). The absolute acceleration response of the top floor under time step  s is plotted in Figure 15 and the nonlinear force-displacement response of the base floor is depicted in Figure 16. It can be seen that the results of the composite algorithm are very close to the Ansys results where Combin40 and Combin14 elements [44] are utilized to simulate the bilinear and the linear floors, respectively. Meanwhile, we can see that, with the help of the LRB, there is prominent reduction in the top floor acceleration of the building, and the fundamental frequency of the superstructure is shifted away from the dominant frequencies of the earthquake ground motion; thus the seismic safety of the building is guaranteed.

To verify the accuracy, another much smaller time step of 0.002 s was also used in the solution procedure of both the Ansys Newmark method and the new composite algorithm. Focusing on the top acceleration response, the accuracy of the time integration algorithm was quantificationally evaluated using the differences between the big time step (0.02 s) and the small time step (0.002 s). For this material nonlinear dynamic system, it is acceptable that result of small time step is much closer to the true solution. And the 2-norm of the difference between the two methods both using small time step is about 0.7, indicating that the new composite integration and the Ansys Newmark integration are both accurate when time step is small enough. As shown in Figure 17, the difference between the solutions of the new composite algorithm using different two time steps, whose 2-norm is about 1.4, is much smaller than that of the Ansys Newmark solutions, with 2-norm equal to 25.9. Obviously the composite algorithm did an excellent job for its high accuracy. Moreover, considering its computational efficiency and simple implementation, the composite algorithm has great potential to be embedded into finite element codes for nonlinear structural dynamic engineering analysis.

4. Conclusions

In this paper, an excellent implicit time integration algorithm was proposed through composition of the three-point backward difference formula and a generalized central difference formula. Theoretical analysis and investigation on a simple frequently used stiff linear system which can represent complicated practical structures have been carried out to give considerable insights into the composite time integration algorithm. Besides, three finite element examples involving nonlinearities on geometry, contact, and material were analyzed in detail to verify the excellent performance on nonlinear problems.

Serving as a robust and efficient time integration method valuable for both linear and nonlinear dynamic engineering analysis, the attractive properties of the composite algorithm are emphasized herein: (1) it is of second-order accuracy and is unconditionally stable; (2) with appropriate numerical damping in high frequency range, it is suitable for not only linear stiff problems but also nonlinear problems; (3) it is second-order and convergent in terms of the displacement, velocity, and acceleration; (4) it only solves a symmetric matrix equation once within each time step, resulting in the same calculation efficiency as the trapezoidal rule but can deal with the nonlinear problems where the trapezoidal rule fails; (5) without any additional parameters to be determined, it is simple and convenient to be implemented into existing finite element codes.

It can be concluded that composition is a potential idea for constructing robust and effective algorithms which deserves intensive studies in the future study of the finite element method fields, especially in nonlinear cases.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This study was financially supported by the National Natural Science Foundation of China (11372156), the Research Project of State Key Laboratory of Hydroscience and Engineering of Tsinghua University (no. 2012-KY-4), and the National Basic Research Program of China (2013CB035902).