Performance of Composite Implicit Time Integration Scheme for Nonlinear Dynamic Analysis

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: 𝑀𝑢̈+𝐶𝑢̇+𝑓(𝑢,𝑡)=𝑝(𝑡), (2.1) 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 𝑡_𝑛+1, such that 𝑡_𝑛+1=𝑡_𝑛+Δ𝑡. Consider 𝑡_𝑛+𝛾=𝑡_𝑛+𝛾Δ𝑡 as a time instance between 𝑡_𝑛 and 𝑡_𝑛+1, with 𝛾∈(0,1) 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: 𝑢̇_𝑛+𝛾=𝑢̇_𝑛+𝑢̈_𝑛+𝑢̈_𝑛+𝛾2𝑢𝛾Δ𝑡,_𝑛+𝛾=𝑢_𝑛+𝑢̇_𝑛+𝑢̇_𝑛+𝛾2𝛾Δ𝑡. (3.1)(3.2) Substituting (3.1) into (3.2), we obtain 𝑢_𝑛+𝛾=𝑢^∗_𝑛+𝛾+𝛾²Δ𝑡²4𝑢̈_𝑛+𝛾 (3.3) with 𝑢^∗_𝑛+𝛾=𝑢_𝑛+𝛾Δ𝑡𝑢̇_𝑛+𝛾²Δ𝑡²4𝑢̈_𝑛. (3.4) On the other hand, (3.1) may be rewritten as 𝑢̇_𝑛+𝛾=𝑢̇^∗_𝑛+𝛾+𝛾Δ𝑡2𝑢̈_𝑛+𝛾 (3.5) with 𝑢̇^∗_𝑛+𝛾=𝑢̇_𝑛+𝛾Δ𝑡2𝑢̈_𝑛. (3.6) Therefore, using (3.3) and (3.5), the accelerations and velocities may be obtained as 𝑢̈_𝑛+𝛾=4𝛾²Δ𝑡²(𝑢_𝑛+𝛾−𝑢^∗_𝑛+𝛾𝑢̇),_𝑛+𝛾=𝑢̇^∗_𝑛+𝛾+2𝛾Δ𝑡(𝑢_𝑛+𝛾−𝑢^∗_𝑛+𝛾). (3.7)(3.8) The equation of motion (2.1) at time 𝑡+𝛾Δ𝑡 may be rewritten as 𝑀𝑢̈_𝑛+𝛾+𝐶𝑢̇_𝑛+𝛾+𝑓_𝑛+𝛾(𝑢_𝑛+𝛾)=𝑝_𝑛+𝛾. (3.9) With (3.7), (3.8), and (3.9), the residual force vector 𝑟_𝑛+𝛾 is defined as 𝑟_𝑛+𝛾=4𝛾²Δ𝑡²𝑀(𝑢_𝑛+𝛾−𝑢^∗_𝑛+𝛾)+𝐶[𝑢̇^∗_𝑛+𝛾+2𝛾Δ𝑡(𝑢_𝑛+𝛾−𝑢^∗_𝑛+𝛾)]+𝑓_𝑛+𝛾−𝑝_𝑛+𝛾. (3.10) 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 (𝐾^𝑖_𝑛+𝛾+4𝛾²Δ𝑡²2𝑀+𝛾Δ𝑡𝐶)Δ𝑢^𝑖+1_𝑛+𝛾=𝑃_𝑛+𝛾−{𝑓^𝑖_𝑛+𝛾4+𝑀[𝛾²Δ𝑡²(𝑢^𝑖_𝑛+𝛾−𝑢^∗_𝑛+𝛾)]+𝐶[𝑢̇^∗_𝑛+𝛾+2𝛾Δ𝑡(𝑢^𝑖_𝑛+𝛾−𝑢^∗_𝑛+𝛾)]} (3.11) with 𝑢^𝑖+1_𝑛+𝛾=𝑢^𝑖_𝑛+𝛾+Δ𝑢^𝑖+1_𝑛+𝛾 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 𝑔̇_𝑛+1=𝑐₁𝑔_𝑛+𝑐₂𝑔_𝑛+𝛾+𝑐₃𝑔_𝑛+1. (4.1) In this case, the constants [6] 𝑐₁,𝑐₂, and 𝑐₃ may be expressed as 𝑐₁=(1−𝛾)𝛾Δ𝑡;𝑐₂=−1(1−𝛾)𝛾Δ𝑡;𝑐₃=(2−𝛾)(1−𝛾)Δ𝑡. (4.2) 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: 𝑢̇_𝑛+1=𝑐₁𝑢_𝑛+𝑐₂𝑢_𝑛+𝛾+𝑐₃𝑢_𝑛+1,𝑢̈_𝑛+1=𝑐₁𝑢̇_𝑛+𝑐₂𝑢̇_𝑛+𝛾+𝑐₃𝑢̇_𝑛+1. (4.3) Figure 3 illustrates the three-point Backward Euler method in which the quantities at 𝑡_𝑛+1 are calculated from the values at 𝑡_𝑛 and 𝑡_𝑛+𝛾. These equations may be rewritten as 𝑢̇_𝑛+1=𝑢̇^∗_𝑛+𝛾+𝑐₃𝑢_𝑛+1,𝑢̈_𝑛+1=𝑢̈^∗_𝑛+𝛾+𝑐₃𝑢̇_𝑛+1 (4.4)(4.5) with 𝑢̇^∗_𝑛+1=𝑐₁𝑢_𝑛+𝑐₂𝑢_𝑛+𝛾,𝑢̈^∗_𝑛+1=𝑐₁𝑢̇_𝑛+𝑐₂𝑢̇_𝑛+𝛾. (4.6) Substituting (4.4) into (4.5), we obtain 𝑢̈_𝑛+1=𝑢̈^∗_𝑛+1+𝑐₃𝑢̇^∗_𝑛+1+𝑐²₃𝑢_𝑛+1. (4.7) Equation (2.1) at time 𝑡+Δ𝑡 may be rewritten as 𝑀𝑢̈_𝑛+1+𝐶𝑢̇_𝑛+1+𝑓_𝑛+1(𝑢_𝑛+1)=𝑃_𝑛+1. (4.8)Putting (4.4) and (4.5) into (4.8), the residual force vector is defined as 𝑟_𝑛+1=𝑀(𝑢̈^∗_𝑛+1+𝑐₃𝑢̇^∗_𝑛+1+𝑐²₃𝑢_𝑛+1)+𝐶(𝑢̇^∗_𝑛+1+𝑐₃𝑢_𝑛+1)+𝑓_𝑛+1−𝑝_𝑛+1. (4.9) Expanding the resulting equation (4.9) in a Taylor’s series up to the first-order terms, as a function of the displacements 𝑢_𝑛+1, we obtain (𝐾^𝑖_𝑛+1+𝑐²₃𝑀+𝑐₃𝐶)Δ𝑢^𝑖+1_𝑛+1=𝑃_𝑛+1−{𝑓^𝑖_𝑛+1+𝑀(𝑢̈^∗_𝑛+1+𝑐₃𝑢̇^∗_𝑛+1+𝑐²₃𝑢_𝑛+1)+𝐶(𝑢̇^∗_𝑛+1+𝑐₃𝑢_𝑛+1)} (4.10) with 𝑢^𝑖+1_𝑛+1=𝑢^𝑖_𝑛+1+Δ𝑢^𝑖+1_𝑛+1. The tangent stiffness matrix 𝐾^𝑖_𝑛+1=𝜕𝑓^𝑖_𝑛+1/𝜕𝑢^𝑖_𝑛+1 and the internal forces 𝑓^𝑖_𝑛+1 are obtained in a consistent way at the configuration corresponding to the displacements 𝑢^𝑖_𝑛+1. 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 𝑢_𝑛+1=𝑢_𝑛+𝛾+(1−𝛾)Δ𝑡𝑢̇_𝑛+𝛾+(1−𝛾)²Δ𝑡²4𝑢̈_𝑛+𝛾. (4.11) 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 10⁻⁵ on the norm of the residual forces. The mid-point time step varies according to different values 𝛾=0.4,0.45,0.5,0.55,0.6. 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 𝐸𝐴=10¹⁰N. 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, 𝑚𝑢̇²_𝑜/2=𝑚𝑔𝑙. Therefore, the initial velocity is given by 𝑢̇_𝑜=2𝑔𝑙=7.72m/s, and the initial acceleration is 𝑢̈_𝑜=2𝑔=19.6m/s². 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 𝜋_𝑜=𝑚𝑢̇²_𝑜/2=298Nm, and the angular momentum 𝐻_𝑜=𝑙𝑚𝑢̇₀. The period of this pendulum is give by 𝑇=𝜋2𝑙/𝑔=2.47seconds, which corresponds to an angle of 360^∘, that is, 1 cycle or a complete turnaround in 2.47 seconds.

Three time steps are taken: Δ𝑡=0.01second, Δ𝑡=0.1second, and Δ𝑡=0.6seconds, corresponding to the following ratios to the period Δ𝑡/𝑇 = 0.004, 0.04, and 0.24; and also to the following angles: 1.45^∘, 14.5^∘ and 87.3^∘, 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 Δ𝑡=0.01second, the numerical dissipation detected is clearly negligible either for the total potential energy as well as for the angular momentum.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 5 

Rigid pendulum. Solution with the implicit-composed algorithm.

However, for Δ𝑡=0.1second, the numerical dissipations along the time are noticeable. On the other hand, for Δ𝑡=0.6seconds, 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 Δ𝑡=0.1second and Δ𝑡=0.6seconds 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 Δ𝑡=0.1second, and errors increase from the seventh cycle on. For Δ𝑡=0.6seconds, 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 𝜀≤2×10⁻⁸ 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 [𝑡_𝑛+𝛾;𝑡_𝑛+1] with 𝛾=0.5. Finally, it is worth mentioning that the algorithm presented here showed numerical stability even when too large time step is used, for example, for Δ𝑡=0.6seconds. 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 (1−𝛾)Δ𝑡 over the numerical dissipation generated by the method, the analysis of this problem is performed with 𝛾=0.4, 𝛾=0.45, 𝛾=0.5, 𝛾=0.55, and 𝛾=0.6. For time step Δ𝑡=0.1second, 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 𝑡=50seconds and a time step Δ𝑡=0.1second. In that table, one can observe how small such decays are. For time step Δ𝑡=0.6seconds, 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 𝛾<0.5, the numerical dissipation continues high. Therefore, one can conclude the following. (a) For 𝛾<0.5, the numerical dissipation is reduced. (b) For 𝛾>0.5, the numerical dissipation grows. (c) For Δ𝑡=0.1second, there are minor decreases for the potential energy and angular momentum. (d) For Δ𝑡=0.6seconds, there are strong decreases for the potential energy and angular momentum.

Table 1 

Energy-momentum decaying at 𝑡=50seconds with Δ𝑡=0.1second.

Table 2 

Energy-momentum decaying at 𝑡=50seconds with Δ𝑡=0.6seconds.

(a)

(b)

(c)

(d)

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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

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 𝜔_𝑜=𝜃̇_𝑜=1.0rad/s 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 𝛼_𝑜=𝜃̈_𝑜=0.

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, 𝑡=50seconds 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 𝜋=22𝑚𝑙²𝜔²_𝑜=0.11×10⁸N⋅cm, and the angular momentum, with respect to pole B, is 𝐻_𝑜=44𝑚𝑙²𝜔_𝑜=0.22×10⁸kg⋅cm²/s. The components of the displacement, velocity, and acceleration vectors of node 1 (pole A) may be obtained by the following expression: 𝑢_𝑥1=2𝑙[1−cos(𝜔_𝑜𝑡)];𝑢_𝑦1=−2𝑙[𝑤_𝑜𝑡+sin(𝜔_𝑜𝑢̇𝑡)];_𝑥1=2𝑙𝜔_𝑜[sin(𝜔_𝑜𝑡)];𝑢̇_𝑦1=−2𝑙𝜔_𝑜[1+cos(𝜔_𝑜𝑢̈𝑡)];_𝑥1=2𝑙𝜔²_𝑜[cos(𝜔_𝑜𝑡)];𝑢̈_𝑦1=2𝑙𝜔²_𝑜[sin(𝜔_𝑜𝑡)]. (5.1) The period of this system is given by 𝑇=2𝜋/𝜔_𝑜=6.28seconds, which corresponds to a turnaround of 360^∘. Three time steps are used for this example: Δ𝑡=0.01second, Δ𝑡=0.1second, and Δ𝑡=1second, which correspond to the following ratios to the period, that is, Δ𝑡/𝑇=0.0016, 0.016, and 0.16 in the same order, to the following angles 0.57^∘, 5.73^∘, and 57.3^∘. These angles represent small, moderate, and large rotations, respectively.

The transient analysis is carried out for a total time duration of 𝑡=50seconds 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 Δ𝑡=0.01second and Δ𝑡=0.1second, there is an excellent agreement to the exact solution. Conversely, for time step Δ𝑡=1second, 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 Δ𝑡=1second 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 𝜀≤1.1×10⁻⁶ due to the hypothesis of rigid-body motion. Figure 7(f) shows, for Δ𝑡=0.01second and Δ𝑡=0.1second, the numerical dissipation of the total potential energy and the angular momentum are insignificant. However, for Δ𝑡=1second, 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, Δ𝑡=1second. This can be seen in Figure 7(f) noticing that no excessive increase of energy-momentum of the system is observed.

For Δ𝑡=1second, the energy-momentum decay, shown in Figures 8(a) and 8(b), increases proportional to 𝛾. Table 3 shows these decays for 𝑡=50seconds. 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 𝑡=50seconds with Δ𝑡=1second.

(a)

(b)

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 𝐸𝐴=10⁴N 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 𝜋_𝑜=𝑚𝑢̇²_𝑜/2=298N⋅m. The angular momentum is 𝐻_𝑜=𝑙𝑚𝑢̇_𝑜=235kg⋅cm²/s. In this example, the period is given by 𝑇=𝜋2𝑙/𝑔=2.47seconds, which corresponds to 1 cycle or a complete turnaround of the pendulum in 2.47seconds. In addition, Due to the axial elastic behavior of the pendulum bar, other oscillation frequency exists, a high axial frequency corresponding to 𝑇=0.28seconds. To capture this axial frequency, two time steps are adopted: Δ𝑡=0.01second and Δ𝑡=0.05seconds which correspond to the following ratios to the period; that is, Δ𝑡/𝑇=0.036 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 Δ𝑡=0.01second, the numerical dissipation detected is minimal either for the total potential energy as well as for the angular momentum. However, for Δ𝑡=0.05seconds, 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 10 

Elastic pendulum. Solution with the implicit-composed algorithm.

Furthermore, for Δ𝑡=0.05seconds, 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 Δ𝑡=0.05seconds, 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 𝑡=30seconds and it is clear that the decays are very small and almost at the same amount.

Table 4 

Energy-momentum decaying at 𝑡=30seconds with Δ𝑡=0.05seconds.

(a)

(b)

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 Δ𝑡/𝑇≤0.1, the algorithm presents an insignificant numerical dissipations, however, for Δ𝑡/𝑇>0.1, 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.