# Performance of Composite Implicit Time Integration Scheme for Nonlinear Dynamic Analysis

Authors: William Taylor Matias Silva; Luciano Mendes Bezerra
Journal: Mathematical Problems in Engineering
 (2008)

---

## 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.


---

## Body

## 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)Mu¨+Cu˙+f(u,t)=p(t),
where M is the mass matrix, C is the damping matrix, f is the vector of internal forces, and p(t) is the vector of external forces. Moreover, u¨, u˙,
and u are, respectively, the vectors of acceleration,
velocity, and displacement. We assume that M and C are constant matrices and we observe that
(2.1) is a nonlinear equation because the internal force vector f is a function of the displacement vector u.
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
timetn and we wish to get a solution at time tn+1,
such that tn+1=tn+Δt.
Consider tn+γ=tn+γΔt as a time instance between tn and tn+1,
with γ∈(0,1) according to Figure 1. Applying now
the trapezoidal rule over the time step, γΔt,
we can get the velocities and displacements for the time tn+γ,
by means of the following finite difference equations, respectively:
(3.1)u˙n+γ=u˙n+u¨n+u¨n+γ2γΔt,(3.2)un+γ=un+u˙n+u˙n+γ2γΔt.
Substituting (3.1) into (3.2), we obtain
(3.3)un+γ=un+γ*+γ2Δt24u¨n+γ
with
(3.4)un+γ*=un+γΔtu˙n+γ2Δt24u¨n.
On the other hand, (3.1) may be rewritten as
(3.5)u˙n+γ=u˙n+γ*+γΔt2u¨n+γ
with
(3.6)u˙n+γ*=u˙n+γΔt2u¨n.
Therefore, using (3.3) and (3.5), the accelerations and velocities may be obtained as
(3.7)u¨n+γ=4γ2Δt2(un+γ−un+γ*),(3.8)u˙n+γ=u˙n+γ*+2γΔt(un+γ−un+γ*).
The equation of motion (2.1) at
time t+γΔt may be rewritten as
(3.9)Mu¨n+γ+Cu˙n+γ+fn+γ(un+γ)=pn+γ.
With (3.7), (3.8), and (3.9),
the residual force vector rn+γ is defined as
(3.10)rn+γ=4γ2Δt2M(un+γ−un+γ*)+C[u˙n+γ*+2γΔt(un+γ−un+γ*)]+fn+γ−pn+γ.
Expanding the resulting equation (3.10) into a Taylor’s series as a function of the displacements un+λ and considering only the first-order terms, we
can get
(3.11)(Kn+γi+4γ2Δt2M+2γΔtC)Δun+γi+1=Pn+γ−{fn+γi+M[4γ2Δt2(un+γi−un+γ*)]+C[u˙n+γ*+2γΔt(un+γi−un+γ*)]}
with un+γi+1=un+γi+Δun+γi+1 and Kn+γi=∂fn+γi/∂un+γi being the consistent tangent stiffness matrix
at the configuration corresponding to the displacements un+γi.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 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)g˙n+1=c1gn+c2gn+γ+c3gn+1.
In this case, the constants [6] c1,c2, and c3 may be expressed as
(4.2)c1=(1−γ)γΔt;c2=−1(1−γ)γΔt;c3=(2−γ)(1−γ)Δt.
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:
(4.3)u˙n+1=c1un+c2un+γ+c3un+1,u¨n+1=c1u˙n+c2u˙n+γ+c3u˙n+1.
Figure 3 illustrates the
three-point Backward Euler method in which the quantities at tn+1 are calculated from the values at tn and tn+γ.
These equations may be rewritten as
(4.4)u˙n+1=u˙n+γ*+c3un+1,(4.5)u¨n+1=u¨n+γ*+c3u˙n+1
with
(4.6)u˙n+1*=c1un+c2un+γ,u¨n+1*=c1u˙n+c2u˙n+γ.
Substituting (4.4) into (4.5), we obtain
(4.7)u¨n+1=u¨n+1*+c3u˙n+1*+c32un+1.
Equation (2.1) at time t+Δt may be rewritten as
(4.8)Mu¨n+1+Cu˙n+1+fn+1(un+1)=Pn+1.Putting (4.4) and (4.5) into (4.8), the residual force vector is defined as
(4.9)rn+1=M(u¨n+1*+c3u˙n+1*+c32un+1)+C(u˙n+1*+c3un+1)+fn+1−pn+1.
Expanding the resulting equation (4.9)
in a Taylor’s series up to the first-order terms, as a function of the displacements un+1, we obtain
(4.10)(Kn+1i+c32M+c3C)Δun+1i+1=Pn+1−{fn+1i+M(u¨n+1*+c3u˙n+1*+c32un+1)+C(u˙n+1*+c3un+1)}
with  un+1i+1=un+1i+Δun+1i+1.
The tangent stiffness matrix Kn+1i=∂fn+1i/∂un+1i and the internal forces fn+1i are obtained in a consistent way at the
configuration corresponding to the displacements un+1i.
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
(4.11)un+1=un+γ+(1−γ)Δtu˙n+γ+(1−γ)2Δt24u¨n+γ.
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
Figure2 is used with a convergence tolerance of 10−5 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 ofEA=1010N.
The pendulum displacement is treated as a rigid body rotation around a fixed
axis and with a constant angular velocity, where θ˙=ωo, u˙o=θ˙l and u¨o=θ˙2l=u˙o2/l.
Consider also the energy conservation, mu˙o2/2=mgl.
Therefore, the initial velocity is given by u˙o=2gl=7.72m/s,
and the initial acceleration is u¨o=2g=19.6m/s2.
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=mu˙o2/2=298Nm,
and the angular momentum Ho=lmu˙0.
The period of this pendulum is give by T=π2l/g=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:Δt=0.01second,  Δt=0.1second, and Δt=0.6seconds,
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.01second,
the numerical dissipation detected is clearly negligible either for the total
potential energy as well as for the angular momentum.Rigid pendulum. Solution with the implicit-composed algorithm.
(a)(b)(c)(d)(e)(f)(g)However, forΔt=0.1second,
the numerical dissipations along the time are noticeable. On the other hand, for Δt=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 Δt=0.1second and Δt=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Δt=0.1second,
and errors increase from the seventh cycle on. For Δt=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−8 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,
[tn;tn+γ] and [tn+γ;tn+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.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γΔ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.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 t=50seconds and a time step Δt=0.1second.
In that table, one can observe how small such decays are. For time step Δt=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 Δt=0.1second,
there are minor decreases for the potential energy and angular momentum. (d) For Δt=0.6seconds, there are strong decreases for the potential
energy and angular momentum.Table 1
Energy-momentum decaying att=50seconds with Δt=0.1second.
γTotal energyAngular momentum0.401.59%0.79%0.453.27%1.64%0.505.16%2.60%0.557.22%3.67%0.609.43%4.82%Table 2
Energy-momentum decaying att=50seconds with Δt=0.6seconds.
γTotal energyAngular momentum0.4055.86%33.56%0.4571.95%47.03%0.5078.15%53.25%0.5581.96%57.52%0.6084.66%60.83%Rigid pendulum. Energy-momentum decaying with different substep sizes.
(a)(b)(c)(d)


### 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 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.Four-bar-chain. Solution with the implicit-composed algorithm.
(a)(b)(c)(d)(e)(f)(g)The
geometrical and physical characteristics of the five connected spheres, the
initial conditions, the boundary conditions, and other applicable data are
summarized in Figure7(a). The initial conditions of the system are (a) an
angular velocity of ωo=θ˙o=1.0rad/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=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 xy-plane. The
total potential energy is given by the expression π=22ml2ωo2=0.11×108N⋅cm,
and the angular momentum, with respect to pole B, is Ho=44ml2ωo=0.22×108kg⋅cm2/s.
The components of the displacement, velocity, and acceleration vectors of node
1 (pole A) may be obtained by the following expression:
(5.1)ux1=2l[1−cos(ωot)];uy1=−2l[wot+sin(ωot)];u˙x1=2lωo[sin(ωot)];u˙y1=−2lωo[1+cos(ωot)];u¨x1=2lωo2[cos(ωot)];u¨y1=2lωo2[sin(ωot)].
The period of this system is given by T=2π/ωo=6.28seconds,
which corresponds to a turnaround of  360°.
Three time steps are used for this example: Δt=0.01second, Δt=0.1second, and Δt=1second,
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 oft=50seconds or, approximately, 8 cycles. Figures 7(b),
7(c), and 7(d) plot as functions of time the x-displacement, 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.01second and Δt=0.1second, there is an excellent agreement to the exact
solution. Conversely, for time step Δt=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 Δt=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−6 due to the hypothesis of rigid-body motion.
Figure 7(f) shows, for Δt=0.01second and Δt=0.1second,
the numerical dissipation of the total potential energy and the angular
momentum are insignificant. However, for Δt=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, Δt=1second.
This can be seen in Figure 7(f) noticing that no excessive increase of energy-momentum
of the system is observed.ForΔt=1second,
the energy-momentum decay, shown in Figures 8(a) and 8(b), increases
proportional to γ.
Table 3 shows these decays for t=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 att=50seconds with Δt=1second.
γTotal energyAngular momentum0.408.24%4.49%0.4511.25%6.37%0.5013.45%7.86%0.5515.19%9.12%0.6016.60%10.21%Four-bar-chain. Energy-momentum decaying with different sub-step sizes.
(a)(b)


### 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 EA=104N 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=mu˙o2/2=298N⋅m.
The angular momentum is Ho=lmu˙o=235kg⋅cm2/s.
In this example, the period is given by T=π2l/g=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 T=0.28seconds.
To capture this axial frequency, two time steps are adopted: Δt=0.01second and Δt=0.05seconds which correspond to the following ratios to the period; that is, Δt/T=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 Figures10(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.01second,
the numerical dissipation detected is minimal either for the total potential
energy as well as for the angular momentum. However, for Δt=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 y-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.Elastic pendulum. Solution with the implicit-composed algorithm.
(a)(b)(c)(d)(e)(f)(g)Furthermore, forΔt=0.05seconds, 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.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 t=30seconds and it is clear that the decays are very small
and almost at the same amount.Table 4
Energy-momentum decaying att=30seconds with Δt=0.05seconds.
γTotal energyAngular momentum0.402.03%0.21%0.452.10%0.22%0.502.19%0.25%0.552.25%0.27%0.602.26%0.27%Elastic pendulum. Energy-momentum decaying with different substep sizes.
(a)(b)


## 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 ofEA=1010N.
The pendulum displacement is treated as a rigid body rotation around a fixed
axis and with a constant angular velocity, where θ˙=ωo, u˙o=θ˙l and u¨o=θ˙2l=u˙o2/l.
Consider also the energy conservation, mu˙o2/2=mgl.
Therefore, the initial velocity is given by u˙o=2gl=7.72m/s,
and the initial acceleration is u¨o=2g=19.6m/s2.
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=mu˙o2/2=298Nm,
and the angular momentum Ho=lmu˙0.
The period of this pendulum is give by T=π2l/g=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:Δt=0.01second,  Δt=0.1second, and Δt=0.6seconds,
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.01second,
the numerical dissipation detected is clearly negligible either for the total
potential energy as well as for the angular momentum.Rigid pendulum. Solution with the implicit-composed algorithm.
(a)(b)(c)(d)(e)(f)(g)However, forΔt=0.1second,
the numerical dissipations along the time are noticeable. On the other hand, for Δt=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 Δt=0.1second and Δt=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Δt=0.1second,
and errors increase from the seventh cycle on. For Δt=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−8 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,
[tn;tn+γ] and [tn+γ;tn+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.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γΔ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.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 t=50seconds and a time step Δt=0.1second.
In that table, one can observe how small such decays are. For time step Δt=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 Δt=0.1second,
there are minor decreases for the potential energy and angular momentum. (d) For Δt=0.6seconds, there are strong decreases for the potential
energy and angular momentum.Table 1
Energy-momentum decaying att=50seconds with Δt=0.1second.
γTotal energyAngular momentum0.401.59%0.79%0.453.27%1.64%0.505.16%2.60%0.557.22%3.67%0.609.43%4.82%Table 2
Energy-momentum decaying att=50seconds with Δt=0.6seconds.
γTotal energyAngular momentum0.4055.86%33.56%0.4571.95%47.03%0.5078.15%53.25%0.5581.96%57.52%0.6084.66%60.83%Rigid pendulum. Energy-momentum decaying with different substep sizes.
(a)(b)(c)(d)


## 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 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.Four-bar-chain. Solution with the implicit-composed algorithm.
(a)(b)(c)(d)(e)(f)(g)The
geometrical and physical characteristics of the five connected spheres, the
initial conditions, the boundary conditions, and other applicable data are
summarized in Figure7(a). The initial conditions of the system are (a) an
angular velocity of ωo=θ˙o=1.0rad/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=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 xy-plane. The
total potential energy is given by the expression π=22ml2ωo2=0.11×108N⋅cm,
and the angular momentum, with respect to pole B, is Ho=44ml2ωo=0.22×108kg⋅cm2/s.
The components of the displacement, velocity, and acceleration vectors of node
1 (pole A) may be obtained by the following expression:
(5.1)ux1=2l[1−cos(ωot)];uy1=−2l[wot+sin(ωot)];u˙x1=2lωo[sin(ωot)];u˙y1=−2lωo[1+cos(ωot)];u¨x1=2lωo2[cos(ωot)];u¨y1=2lωo2[sin(ωot)].
The period of this system is given by T=2π/ωo=6.28seconds,
which corresponds to a turnaround of  360°.
Three time steps are used for this example: Δt=0.01second, Δt=0.1second, and Δt=1second,
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 oft=50seconds or, approximately, 8 cycles. Figures 7(b),
7(c), and 7(d) plot as functions of time the x-displacement, 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.01second and Δt=0.1second, there is an excellent agreement to the exact
solution. Conversely, for time step Δt=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 Δt=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−6 due to the hypothesis of rigid-body motion.
Figure 7(f) shows, for Δt=0.01second and Δt=0.1second,
the numerical dissipation of the total potential energy and the angular
momentum are insignificant. However, for Δt=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, Δt=1second.
This can be seen in Figure 7(f) noticing that no excessive increase of energy-momentum
of the system is observed.ForΔt=1second,
the energy-momentum decay, shown in Figures 8(a) and 8(b), increases
proportional to γ.
Table 3 shows these decays for t=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 att=50seconds with Δt=1second.
γTotal energyAngular momentum0.408.24%4.49%0.4511.25%6.37%0.5013.45%7.86%0.5515.19%9.12%0.6016.60%10.21%Four-bar-chain. Energy-momentum decaying with different sub-step sizes.
(a)(b)


## 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 EA=104N 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=mu˙o2/2=298N⋅m.
The angular momentum is Ho=lmu˙o=235kg⋅cm2/s.
In this example, the period is given by T=π2l/g=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 T=0.28seconds.
To capture this axial frequency, two time steps are adopted: Δt=0.01second and Δt=0.05seconds which correspond to the following ratios to the period; that is, Δt/T=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 Figures10(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.01second,
the numerical dissipation detected is minimal either for the total potential
energy as well as for the angular momentum. However, for Δt=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 y-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.Elastic pendulum. Solution with the implicit-composed algorithm.
(a)(b)(c)(d)(e)(f)(g)Furthermore, forΔt=0.05seconds, 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.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 t=30seconds and it is clear that the decays are very small
and almost at the same amount.Table 4
Energy-momentum decaying att=30seconds with Δt=0.05seconds.
γTotal energyAngular momentum0.402.03%0.21%0.452.10%0.22%0.502.19%0.25%0.552.25%0.27%0.602.26%0.27%Elastic pendulum. Energy-momentum decaying with different substep sizes.
(a)(b)


## 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Δ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 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.



---

*Source: 815029-2008-07-14.xml*