Performance of Composite Implicit Time Integration Scheme for Nonlinear Dynamic Analysis

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.


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 longtime 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 Therefore, using 3.3 and 3.5 , the accelerations and velocities may be obtained as The equation of motion 2.1 at time t γΔt may be rewritten as M ün γ C un γ f n γ u n γ p n γ .

3.9
Mathematical Problems in Engineering With 3.7 , 3.8 , and 3.9 , the residual force vector r n γ is defined as f n γ − p n γ .

3.10
Expanding the resulting equation 3.10 into a Taylor's series as a function of the displacements u n λ and considering only the first-order terms, we can get Δu i 1 n γ and K i n γ ∂f i n γ /∂u i n γ being the consistent tangent stiffness matrix at the configuration corresponding to the displacements u i n γ .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.

Second substep
Let the derivative 6 of a continuous function g at time t Δt be written in terms of the derivatives of the function g at times t, t γΔt and t Δt as 4.1 In this case, the constants 6 c 1 , c 2 , and c 3 may be expressed as Thus, the velocities as functions of the displacements, and the accelerations as functions of velocities at time t Δt 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 t n 1 are calculated from the values at t n and t n γ .These equations may be rewritten as Convergence ?
Calculate the iteration matrix Update the kinematics variables Residual vector evaluation Calculate the iteration matrix Update the kinematics variables  Putting 4.4 and 4.5 into 4.8 , the residual force vector is defined as 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 u n 1 , we obtain ∂f i n 1 /∂u i n 1 and the internal forces f i n 1 are obtained in a consistent way at the configuration corresponding to the displacements u i n 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 incrementaliterative 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.

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 −5 on the norm of the residual forces.The midpoint 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.

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.
The rigid pendulum has an axial stiffness of EA 10 10 N. The pendulum displacement is treated as a rigid body rotation around a fixed axis and with a constant angular velocity, where θ ω o , uo θl and üo θ 2 l u2 o /l.Consider also the energy conservation, m u2 o /2 mgl.Therefore, the initial velocity is given by uo 2gl 7.72 m/s, and the initial acceleration is üo 2g 19.6 m/s 2 .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 π o m u2 o /2 298 Nm, and the angular momentum H o lm u0 .The period of this pendulum is give by T π 2l/g 2.47 seconds, which corresponds to an angle of 360 • , that is, 1 cycle or a complete turnaround in 2.47 seconds.
Three time steps are taken: Δt 0.01 second, Δt 0.1 second, and Δt 0.6 seconds, corresponding to the following ratios to the period Δt/T 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 Δt 0.01 second, the numerical dissipation detected is clearly negligible either for the total potential energy as well as for the angular momentum.
However, for Δt 0.1 second, the numerical dissipations along the time are noticeable.On the other hand, for Δt 0.6 seconds, 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 Δt 0.1 second and Δt 0.6 seconds 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 Δt 0.1 second, and errors increase from the seventh cycle on.For Δt 0.6 seconds, 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 −8 due to the hypothesis of rigid-body motion.is important to point out that such figure deals with the sum of the iterations corresponding to the two substeps, that is, t n ; t n γ and t n γ ; t n 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 Δt 0.6 seconds.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 γΔt and 1 − γ Δt 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 Δt 0.1 second, 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 t 50 seconds and a time step Δt 0.1 second.In that table, one can observe how small such  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 Δt 0.1 second, there are minor decreases for the potential energy and angular momentum.d For Δt 0.6 seconds, there are strong decreases for the potential energy and angular momentum.

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 x-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.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 ω o θo 1.0 rad/s around the axis at pole B node 5 parallel to the z-axis which is equivalent to a linear distribution of vertical velocities, and b a zero-angular acceleration α o θo 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, t 50 seconds 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 xy-plane.The total potential energy is given by the expression π 22ml following expression:

5.1
The period of this system is given by T 2π/ω o 6.28 seconds, which corresponds to a turnaround of 360 • .Three time steps are used for this example: Δt 0.01 second, Δt 0.1 second, and Δt 1 second, which correspond to the following ratios to the period, that is, Δt/T 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 t 50 seconds or, approximately, 8 cycles.Figures 7 b , 7 c , and 7 d plot as functions of time the xdisplacement, the y-velocity, and the y-acceleration of the node-1, respectively.These transient responses are compared to the exact solution in 5.1 .For the time steps Δt 0.01 second and Δt 0.1 second, there is an excellent agreement to the exact solution.Conversely, for time step Δt 1 second, 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 Δt 1 second 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 −6 due to the hypothesis of rigid-body motion.Figure 7 f shows, for Δt 0.01 second and Δt 0.1 second, the numerical dissipation of the total potential energy and the angular momentum are insignificant.However, for Δt 1 second, 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, Δt 1 second.This can be seen in Figure 7 f noticing that no excessive increase of energy-momentum of the system is observed.For Δt 1 second, the energy-momentum decay, shown in Figures 8 a and 8 b , increases proportional to γ. Table 3 shows these decays for t 50 seconds.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.

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 EA 10 4 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 π o m u2 o /2 298 N•m.The angular momentum is H o lm uo 235 kg • cm 2 /s.In this example, the period is given by T π 2l/g 2.47 seconds, which corresponds to 1 cycle or a complete turnaround of the pendulum in 2.47 seconds.In addition, Due to the axial elastic behavior of the pendulum bar, other oscillation frequency exists, a high axial frequency corresponding to T 0.28 seconds.To capture this axial frequency, two time steps are adopted: Δt 0.01 second and Δt 0.05 seconds which correspond to the following ratios to the period; that is, Δt/T 0.036 and 0.18, respectively.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 Δt 0.01 second, the numerical dissipation detected is minimal either for the total potential energy as well as for the angular momentum.However, for Δt 0.05 seconds, there are small numerical dissipations increasing along the time.Furthermore, for Δt 0.05 seconds, Figure 10 f shows the axial oscillations and depicts the significant errors in the amplitudes due to numerical dissipation.With this large Δt, 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 Δt 0.05 seconds, 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 t 30 seconds and it is clear that the decays are very small and almost at the same amount.

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 rigidbody motions.d For time steps Δt with ratios to the period Δt/T ≤ 0.1, the algorithm presents an insignificant numerical dissipations, however, for Δt/T > 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 Mathematical Problems in Engineering 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.
time instance between t n and t n 1 , with γ ∈ 0, 1 according to Figure1.Applying now the trapezoidal rule over the time step, γΔt, we can get the velocities and displacements for the time t n γ , by means of the following finite difference equations, respectively: a
decays are.For time step Δt 0.6 seconds, the decays of the total potential energy and of the angular momentum, illustrated in Figures6 c and 6 d, also grow with γ.Table