`Mathematical Problems in EngineeringVolume 2008 (2008), Article ID 815029, 16 pageshttp://dx.doi.org/10.1155/2008/815029`
Research Article

Performance of Composite Implicit Time Integration Scheme for Nonlinear Dynamic Analysis

Department of Civil Engineering, University of Brasília, 70910-900 Brasília, Brazil

Received 26 February 2008; Accepted 12 May 2008

Copyright © 2008 William Taylor Matias Silva and Luciano Mendes Bezerra. 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 simple implicit time integration scheme for transient response solution of structures under large deformations and long-time durations. The authors focus on a practical method using implicit time integration scheme applied to structural dynamic analyses in which the widely used Newmark time integration procedure is unstable, and not energy-momentum conserving. In this integration scheme, the time step is divided in two substeps. For too large time steps, the method is stable but shows excessive numerical dissipation. The influence of different substep sizes on the numerical dissipation of the method is studied throughout three practical examples. The method shows good performance and may be considered good for nonlinear transient response of structures.

1. Introduction

In the last four decades, the computational mechanics community has accomplished many researches trying to propose effective methods for nonlinear dynamic analysis in the framework of the finite element method. For fast transient analyses, for example, impact problems, explicit methods are largely used. However, for methods conditionally stable, very small time steps are required to get reasonable solutions. For transient analyses of long-time duration, as in vibration problems of structural systems, the implicit methods are more effective. According to Bathe [1], the first implicit integration procedures used are Houbolt, Newmark, and Wilson-. Among these methods, the Newmark method and its particular case, the trapezoidal rule, became very popular and effective for linear dynamic analysis of practical problems. The trapezoidal rule scheme is the most effective one because it is a second-order method and uses single time step. However, in nonlinear dynamic analysis, the trapezoidal rule becomes considerably unstable. Such instability is due to the pathological growth of the total potential energy and angular momentum. The trapezoidal rule integration scheme does not guarantee the conservation of the energy-momentum along the time duration. To overcome this adverse characteristic, many implicit algorithms were additionally proposed based on the following ideas explained by Kuhl and Crisfield in [2]: (a) numeric dissipation [3]; (b) conservation of energy-momentum throughout the use of Lagrange multipliers [4]; and (c) imposition, in the algorithm, of energy-momentum conservation [5].

The present work analyzes the application of the trapezoidal rule to nonlinear dynamic analyses. To keep the conservation of the energy-momentum, the trapezoidal rule is combined to the three-point backward Euler method. This combination is very much employed in numerical procedures to solve ordinary and partial differential equations (o.d.e.s and p.d.e.s) [6]. Bank etal. [7] use this combination to solve first-order o.d.e.s to simulate the behavior of silicon devices and circuits. Recently, Bathe [8], Bathe and Baig [9] utilized these mixed algorithms to get solutions of second-order p.d.e.s describing the dynamic equilibrium of structural systems. They obtained transient responses for beams and plates discretizing them with solid 2D finite elements. The beams and plates studied by such authors were subjected to large translations and rotations due to rigid-body motions. In the next sections, the coupling of the trapezoidal rule and the three-point backward Euler method is explained in details, and the numerical dissipation is studied in front of different substeps.

2. The Implicit-Composed Algorithm

The equation of motion of a deformable body discretized by the finite element method may be expressed by the following matricial equation: where is the mass matrix, is the damping matrix, is the vector of internal forces, and is the vector of external forces. Moreover, , , and are, respectively, the vectors of acceleration, velocity, and displacement. We assume that and are constant matrices and we observe that (2.1) is a nonlinear equation because the internal force vector is a function of the displacement vector . Vectors and vector components are hereafter written with the same notation without loss of meaning. In general, time integration algorithms to solve (2.1) are formulated throughout the finite difference schemes and such schemes show numerical dissipation. In computational mechanics, numerical dissipation means an unexpected lost of energy in the numerical solution. This dissipation property may be good in getting better numerical stabilization for such integration schemes. The implicit-composed scheme divides the time step in two substeps. In the first substep, the trapezoidal rule is applied while in the second substep, we make use of the three-point backward Euler method. As the application of the algorithm aims to nonlinear analyses, it is necessary to establish an incremental-iterative strategy to get the final solution. In this work, the Newton-Raphson method, in the iterative phase, is used to dissipate the residual forces. The equation of motion may be written as a function of the displacements, developed in a Taylor’s series up to the first-order terms, and an incremental-iterative strategy is established for the time-step dynamic analysis.

3. First Substep

At first, it is assumed that the solution of the equation of motion is known at time and we wish to get a solution at time , such that . Consider as a time instance between and , with according to Figure 1. Applying now the trapezoidal rule over the time step, , we can get the velocities and displacements for the time , by means of the following finite difference equations, respectively: Substituting (3.1) into (3.2), we obtain with On the other hand, (3.1) may be rewritten as with Therefore, using (3.3) and (3.5), the accelerations and velocities may be obtained as The equation of motion (2.1) at time may be rewritten as With (3.7), (3.8), and (3.9), the residual force vector is defined as Expanding the resulting equation (3.10) into a Taylor’s series as a function of the displacements and considering only the first-order terms, we can get with and being the consistent tangent stiffness matrix at the configuration corresponding to the displacements .

Figure 1: Generalized substep sizes.

Once the displacements are determined, the accelerations and velocities may be obtained by means of (3.7) and (3.8), respectively. For more details, see the incremental-iterative flow diagram described in Figure 2.

Figure 2: Incremental-iterative scheme of the implicit-composed algorithm.

4. Second Substep

Let the derivative [6] of a continuous function at time be written in terms of the derivatives of the function at times , and as In this case, the constants [6] and may be expressed as Thus, the velocities as functions of the displacements, and the accelerations as functions of velocities at time may be determined by the following equations in the same order: Figure 3 illustrates the three-point Backward Euler method in which the quantities at are calculated from the values at and . These equations may be rewritten as with Substituting (4.4) into (4.5), we obtain Equation (2.1) at time may be rewritten as Putting (4.4) and (4.5) into (4.8), the residual force vector is defined as Expanding the resulting equation (4.9) in a Taylor’s series up to the first-order terms, as a function of the displacements , we obtain with . The tangent stiffness matrix and the internal forces are obtained in a consistent way at the configuration corresponding to the displacements . Once the displacements are determined, the velocities and accelerations may be calculated according to (4.4) and (4.5), respectively. For more details, examine the incremental-iterative flow diagram represented in Figure 2. In [8, 9], the prediction displacement adopted in the first iteration of the second substep is not clear. In the present paper, the trapezoidal rule is employed to obtain the prediction displacement as In short, the method has the following main characteristics: (a) it has no additional variables, like Lagrange multipliers, are used; (b) it is suitable to elastic and inelastic analyses; (c) it has symmetry of the tangent stiffness matrix.

Figure 3: Three-point backward Euler method.

5. Numerical Examples

In the following examples, one finite element in 2D space representing a biarticulated bar is used. The internal forces’ vector and the stiffness matrix of such finite element are obtained from a total Lagrangian formulation. For more details of such formulation, see [10]. The mass matrix in the following examples considers the mass of the bar element as massless and lumped masses concentrated at the two end nodes.

To find the transient response, the incremental iterative scheme illustrated in Figure 2 is used with a convergence tolerance of on the norm of the residual forces. The mid-point time step varies according to different values . The objective here is to exam the performance of the implicit-composed algorithm described in Section 2 for varying and when different time step is adopted for long-time duration. With this in mind, it is important to analyze if the algorithm presents the following undesirable aspects: (1) excessive errors in the period and in the amplitude of the transient response, (2) strong growth of the total potential energy and the angular momentum, (3) strong decay of the total potential energy and the angular momentum, and (4) lack of convergence during the iterative process.

5.1. Rigid Pendulum

Among other authors, Crisfield and Shi [11], Kuhl and Crisfield [2] and Bathe [8] analyzed this example. The geometrical and physical characteristics of the rigid pendulum, the initial conditions, the boundary conditions and other data of the problem are in Figure 4. The rigid pendulum was discretized with one biarticulated finite element bar in 2D space with a very high-axial stiffness. The pendulum has two degrees of freedom restrained and two degrees of freedom released.

Figure 4: Rigid pendulum. Data and initial conditions.

The rigid pendulum has an axial stiffness of . The pendulum displacement is treated as a rigid body rotation around a fixed axis and with a constant angular velocity, where , and . Consider also the energy conservation, . Therefore, the initial velocity is given by , and the initial acceleration is . No external force is applied at the free end-node of the pendulum. Therefore, the total potential energy and the angular momentum are kept constants. The value of total potential energy is , and the angular momentum . The period of this pendulum is give by , which corresponds to an angle of , that is, 1 cycle or a complete turnaround in 2.47 seconds.

Three time steps are taken: , , and , corresponding to the following ratios to the period = , 0.04, and 0.24; and also to the following angles: , and , respectively. Correspondingly, these angles represent small, moderate, and large rotations. The transient analysis is carried out for a total time duration of 50 seconds which means 20 cycles. Figure 5(a) shows the mass trajectories for the three different time steps adopted; observe the coincidence between the trajectories. Examining Figure 5(b), for , the numerical dissipation detected is clearly negligible either for the total potential energy as well as for the angular momentum.

Figure 5: Rigid pendulum. Solution with the implicit-composed algorithm.

However, for , the numerical dissipations along the time are noticeable. On the other hand, for , an excessive numerical dissipation of the total potential energy and the angular momentum is observed. Consequently, errors of great magnitude in the period and in the amplitude response may be observed in Figures 5(c), 5(d), and 5(e), respectively, for displacements, velocities, and accelerations. Errors in the periods of the displacements may be noticed for time steps of and from the seventh cycle on. Those errors increase along the next cycles.

With respect to velocity and acceleration, it may be observed that there are errors in the period and in the amplitude for , and errors increase from the seventh cycle on. For , the errors are meaningful and the transient responses are short of precision to represent the physical model under analysis. In Figure 5(f), the magnitude of axial strains do not exceed due to the hypothesis of rigid-body motion. Figure 5(g) shows the evolution of the number of iterations along the time necessary to get convergence. It is important to point out that such figure deals with the sum of the iterations corresponding to the two substeps, that is, and with . Finally, it is worth mentioning that the algorithm presented here showed numerical stability even when too large time step is used, for example, for . In this case, no growth was observed for the energy momentum of the system, as can be seen in Figure 5(b).

To study the influence of the substeps’ sizes and over the numerical dissipation generated by the method, the analysis of this problem is performed with , , , and . For time step , the energy-momentum decays are shown in Figures 6(a) and 6(b), respectively. In both figures, it is noticed that the numerical dissipations grow proportional to the values. In Table 1, such decays are reported for and a time step . In that table, one can observe how small such decays are. For time step , the decays of the total potential energy and of the angular momentum, illustrated in Figures 6(c) and 6(d), also grow with . Table 2 shows that such decays are excessive which means that the solution for this case is inaccurate. Even for the numerical dissipation continues high. Therefore, one can conclude the following. (a) For the numerical dissipation is reduced. (b) For the numerical dissipation grows. (c) For , there are minor decreases for the potential energy and angular momentum. (d) For , there are strong decreases for the potential energy and angular momentum.

Table 1: Energy-momentum decaying at with .
Table 2: Energy-momentum decaying at with .
Figure 6: Rigid pendulum. Energy-momentum decaying with different substep sizes.
5.2. System with Five Spheres Connected with Massless Rigid Rods

Crisfield and Shi [11] analyzed this example. Figure 7(a) shows a chain of pinned bars (truss element) that is free to fly in the absence of gravity. Initially, the bars lie horizontally with no velocity in the -direction but a linear distribution of vertical velocity. Under such conditions, the chain should remain straight moving downwards and rotating at the same time. The system is made of 5 spherical masses connected by massless rigid rods.

Figure 7: Four-bar-chain. Solution with the implicit-composed algorithm.

The geometrical and physical characteristics of the five connected spheres, the initial conditions, the boundary conditions, and other applicable data are summarized in Figure 7(a). The initial conditions of the system are (a) an angular velocity of around the axis at pole B (node 5) parallel to the -axis which is equivalent to a linear distribution of vertical velocities, and (b) a zero-angular acceleration .

This system was discretized with four finite elements, biarticulated bar elements in the 2D space. The finite element model has five nodes making a total of 10 degrees of freedom. There are no constraint nodes. Gravitational forces are not considered, and therefore the total potential energy and the angular momentum are kept constant along the time considered for the analysis of this problem, that is, Due to the mass symmetry, the center of mass of the system is in middle of the bar length. Therefore, the system of five concentrated masses is subjected to large translations and rotations in the -plane. The total potential energy is given by the expression , and the angular momentum, with respect to pole B, is . The components of the displacement, velocity, and acceleration vectors of node 1 (pole A) may be obtained by the following expression: The period of this system is given by , which corresponds to a turnaround of . Three time steps are used for this example: , , and , which correspond to the following ratios to the period, that is, , 0.016, and 0.16 in the same order, to the following angles , , and . These angles represent small, moderate, and large rotations, respectively.

The transient analysis is carried out for a total time duration of or, approximately, 8 cycles. Figures 7(b), 7(c), and 7(d) plot as functions of time the -displacement, the -velocity, and the -acceleration of the node-1, respectively. These transient responses are compared to the exact solution in (5.1). For the time steps and , there is an excellent agreement to the exact solution. Conversely, for time step , significant errors in the period and in the amplitude are observed. Note these errors increase in the subsequent cycles. Consequently, the displacement, velocity, and acceleration obtained for the time step are not suitable to represent the physical problem studied in this example. In Figure 7(e), the magnitude of element-1 axial strains do not exceed due to the hypothesis of rigid-body motion. Figure 7(f) shows, for and , the numerical dissipation of the total potential energy and the angular momentum are insignificant. However, for , an excessive numerical dissipation is observed and grows along the time. Figure 7(g) shows the number of iterations necessary to get convergence in the solution. As a final point, it is remarkable that the algorithm shows numerical stability even for too large time step, for example, . This can be seen in Figure 7(f) noticing that no excessive increase of energy-momentum of the system is observed.

For , the energy-momentum decay, shown in Figures 8(a) and 8(b), increases proportional to . Table 3 shows these decays for . In that table, it can be observed that the values of these decays are less than 17% for the potential energy and less than 11% for the angular momentum.

Table 3: Energy-momentum decaying at with .
Figure 8: Four-bar-chain. Energy-momentum decaying with different sub-step sizes.
5.3. Elastic Pendulum

This example was analyzed by Kuhl and Crisfield [2] and Bathe [8], among other researchers. The geometrical and physical characteristics of the elastic pendulum, the initial conditions, the boundary conditions, and other data of the problem are in Figure 9. The pendulum was discretized with one biarticulated 2D finite element bar which has two degrees of freedom restrained and two degrees of freedom released. An axial stiffness is assumed. A nonzero initial velocity is considered. No gravitational force is assumed to act, and therefore no external force is applied at the free end-node of the pendulum. Therefore, the total potential energy and the angular momentum are kept constants along the time. The potential energy is . The angular momentum is . In this example, the period is given by , which corresponds to 1 cycle or a complete turnaround of the pendulum in . In addition, Due to the axial elastic behavior of the pendulum bar, other oscillation frequency exists, a high axial frequency corresponding to . To capture this axial frequency, two time steps are adopted: and which correspond to the following ratios to the period; that is, and 0.18, respectively.

Figure 9: Elastic pendulum. Data and initial conditions.

Although, in this case, there are oscillations in high frequencies, no sudden growth is observed in the amplitude of the axial oscillations and in the energy-momentum of the pendulum system. These can be demonstrated in Figures 10(b) and 10(f), respectively. Figure 10(a) shows the pendulum trajectories. The complete agreement of the trajectories is clear. Examining Figure 10(b), for , the numerical dissipation detected is minimal either for the total potential energy as well as for the angular momentum. However, for , there are small numerical dissipations increasing along the time. Figures 10(c) and 10(d) show the displacement and velocity of node 2 in the -direction, respectively. In these figures, for both time steps used, the transient responses are almost coincident. However, in Figure 10(e), significant errors in the amplitude and in the acceleration period are detected.

Figure 10: Elastic pendulum. Solution with the implicit-composed algorithm.

Furthermore, for , Figure 10(f) shows the axial oscillations and depicts the significant errors in the amplitudes due to numerical dissipation. With this large , it is impossible to have a more precise response of the system under high frequency. Finally, Figure 10(g) represents the number of iterations to get convergence in the solutions. For the time step , the decays of the total potential energy and the angular momentum, shown in Figures 11(a) and 11(b), respectively, are practically the same for different substeps used. In Table 4, these decays are reported for and it is clear that the decays are very small and almost at the same amount.

Table 4: Energy-momentum decaying at with .
Figure 11: Elastic pendulum. Energy-momentum decaying with different substep sizes.

6. Concluding Remarks

Concerning the performance of the implicit-composed algorithm applied to nonlinear dynamic analysis, the following conclusions may be taken. (a) The algorithm is easy to implement in a computer program. (b) The mathematical formulation of the algorithm is very simple. (c) The algorithm is effective to deal with large translations and rotations due to rigid-body motions. (d) For time steps with ratios to the period , the algorithm presents an insignificant numerical dissipations, however, for an increasing numerical dissipation is observed. (e) The computational cost of the algorithm is twice greater than the computational cost of the trapezoidal rule due to the two iterative cycles needed in each time step. (f) The algorithm preserves the energy-momentum without the need of Lagrange multipliers or without any imposition in the algorithm for energy-momentum conservation. (g) The algorithm allows the user to work with symmetric matrices. (h) The method is applicable to elastic and inelastic analyses. (i) No additional variables like Lagrange multipliers are used. (j) For too large time step, even with inaccurate solution, the method is stable.

From the view of the authors, the excessive numerical dissipation when using a too large time-step is the major drawback of the present scheme, for applications in nonlinear analyses in practical problems.

References

1. K.-J. Bathe, Finite Element Procedures, Prentice-Hall, Englewood Cliffs, NJ, USA, 1996.
2. D. Kuhl and M. A. Crisfield, “Energy-conserving and decaying algorithms in non-linear structural dynamics,” International Journal for Numerical Methods in Engineering, vol. 45, no. 5, pp. 569–599, 1999.
3. J. Chung and G. M. Hulbert, “A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method,” Journal of Applied Mechanics, vol. 60, no. 2, pp. 371–375, 1993.
4. D. Kuhl and E. Ramm, “Constraint energy momentum algorithm and its application to non-linear dynamics of shells,” Computer Methods in Applied Mechanics and Engineering, vol. 136, no. 3-4, pp. 293–315, 1996.
5. J. C. Simo and N. Tarnow, “The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics,” Zeitschrift für Angewandte Mathematik und Physik, vol. 43, no. 5, pp. 757–792, 1992.
6. L. Collatz, The Numerical Treatment of Differential Equations, Springer, New York, NY, USA, 3rd edition, 1966.
7. R. E. Bank, W. M. Coughran, Jr., W. Fichtner, E. H. Grosse, D. J. Rose, and R. K. Smith, “Transient simulation of silicon devices and circuits,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 4, no. 4, pp. 436–451, 1985.
8. K.-J. Bathe, “Conserving energy and momentum in nonlinear dynamics: a simple implicit time integration scheme,” Computers & Structures, vol. 85, no. 7-8, pp. 437–445, 2007.
9. K.-J. Bathe and M. M. I. Baig, “On a composite implicit time integration procedure for nonlinear dynamics,” Computers & Structures, vol. 83, no. 31-32, pp. 2513–2524, 2005.
10. W. T. Matias Silva, “Aplicaciones de algoritmos que conservan la energia-momentum en dinámica no-lineal,” in Congreso de Métodos Numéricos en Ingeniería, pp. 1–20, Granada, Spain, July 2005.
11. M. A. Crisfield and J. Shi, “A co-rotational element/time-integration strategy for nonlinear dynamics,” International Journal for Numerical Methods in Engineering, vol. 37, no. 11, pp. 1897–1913, 1994.