Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 907023, 11 pages

http://dx.doi.org/10.1155/2015/907023

## A Robust and Efficient Composite Time Integration Algorithm for Nonlinear Structural Dynamic Analysis

State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

Received 28 August 2014; Accepted 30 September 2014

Academic Editor: Song Cen

Copyright © 2015 Lihong Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 [1–4]. 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 [28–30]. 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.*

*α*