An Efficient Newton-Type Iteration for the Numerical Solution of Highly Oscillatory Constrained Multibody Dynamic Systems.
In this paper we present a coordinate-split (CS) technique for the numerical solution of the equations of motion of constrained multibody dynamic systems.  We show how the CS technique can be implemented within the context of commonly used solution methods, for increased efficiency and reliability.A particularly challenging problem for multibody dynamics is the numerical solution of highly oscillatory nonlinear mechanical systems. Highly stable implicit integration methods with large stepsizes can be used to damp the oscillation, if it is of small amplitude.  However, the standard Newton iteration is known to experience severe convergence difficulties which force a restriction of the stepsize.  We introduce a modified coordinate-split (CM) iteration which overcomes these problems. Convergence analysis explains the improved convergence for nonlinear oscillatory systems, and numerical experiments illustrate the effectiveness of the new method.
Introduction
The equations of motion of a constrained multibody system can be written as [10]
are the generalized coordinates, - 1 , - 2 , ., -m ] are the
Lagrange multipliers, M(q) 2 IR n\Thetan is the mass-inertia matrix, f 2 IR n is the force
applied to the system, -
dt
is the velocity and -
is the acceleration vector.
The constraints are m smooth functions of q, whose Jacobian
is assumed of full row-rank. We assume that G(q)M(q)G T (q) is symmetric and positive
definite for every q 2 IR n to obtain a consistent physics represented by (1.1). The
degrees of freedom for the system (1.1) is Equation (1.1) is a well-known
index-3 DAE [3, 13].
Many methods have been proposed for modeling multibody systems. Direct numerical
integration of the index-3 DAE (1.1) suffers from the well-known difficulties
inherent in the solution of high-index DAEs [13]. One way to lower the index involves
introducing derivatives of the constraint g(q), along with additional Lagrange multipliers
-. This yields the stabilized index-2 or GGL formulation of the constrained
equations of motion [8]
which has been widely used in simulation. The Lagrange multiplier variables - and
- fulfill the role of projecting the solution onto the position (1.3d) and the velocity
(1.3c) constraints, respectively. Equations (1.3) and related systems have been solved
by a variety of methods. Here we will consider solution by implicit numerical methods
such as BDF or RADAU. A closely related approach is based on explicitly projecting
the numerical solution onto the constraints [18, 20, 22, 23] and involves many of the
same issues for the implementation that are considered here.
Many of the numerical methods for multibody systems solve the system (1.3)
directly. It is also possible to eliminate the Lagrange multipliers and reduce the size
of the system to the number of degrees of freedom. One way to accomplish this
begins with the stabilized index-2 system (1.3). Suppose that G(p) is full-rank on the
constraint manifold 0g. Then one can find an annihilation
matrix P (q) 2 IR (n\Gammam)\Thetan such that P (q)G T Premultiplying (1.3a)
and (1.3b) by P (q) yields an index-1 DAE
There is a potential gain in efficiency for this formulation due to the size-reduction
of the nonlinear system, compared to (1.3). An important practical consequence of
(1.4) is that (-) have been eliminated from the DAE, via multiplication of (1.3a,
1.3b) by the nonlinear P (q). Thus, the error test and Newton iteration convergence
test in a numerical implementation of (1.4) no longer need to include (-). These
higher-index variables can cause problems in the direct numerical solution of (1.3).
One could in principle also consider removing (-) from the test in the solution
of (1.3), however it is not usually possible to justify this action, particularly in the
case of the Newton convergence test. Elimination of these variables from the Newton
convergence test in the solution of (1.3) can lead to a code which sometimes produces
incorrect solutions. It is the fact that multiplying by the nonlinear P (q) eliminates
(-) from the nonlinear system, which allows these variables to be excluded from
the tests in the solution of (1.4).
Direct numerical solution of presents some challenges. First we must have
a means of generating P (q) which is reliable and cheap. Further, we note that the
Jacobian matrix for the Newton iteration involves complicated terms which arise
from the derivatives of P (q). We need a means of generating the Jacobian matrix.
Finally, practical issues such as the error test and Newton convergence test must be
considered.
Our approach for obtaining a cheap representation of P (q) is based on a coordinate-
splitting of the variables. A widely-used method which is related in the sense of also
making use of a splitting of the coordinates is the generalized coordinate partitioning
method [22]. In order to obtain the independent generalized coordinates, the Jacobian
G is calculated at an assembled configuration q and an LU-decomposition of G is
carried out to identify a nonsingular m \Theta m sub-matrix GY , where
such that X 2 IR n\Thetap and Y 2 IR n\Thetam . The columns of X and Y constitute the
Cartesian basis for IR n . The implicit function theorem assures that the dependent
variables y can be uniquely determined by the constraint g(q), at least locally, as a
function of the remaining generalized coordinates x, i.e.,
thus are declared to be independent coordinates. The twice differentiable function
arbitrary x in some neighborhood of
the assembled configuration. Differentiating twice the constraints (1.1c) with respect
to time, then multiplying (1.1b) by GM \Gamma1 and substituting for G -
v from the twice
differentiated constraints to solve for -, we obtain
dGv
dq
Substituting the solution -(q; v) into (1.1), and eliminating the dependent variables
y, -
y, and -
y, yields n-m differential equations
dh
dx
which is a state-space form representation of (1.1) [7, 14]. However, this differs
substantially from the approach we outline here because P (q) associated with this
method is not orthogonal to G T (q). Hence the index-reduction by differentiating the
constraints and projecting to the invariant space must be carried out explicitly. In
particular, this requires forming the derivative of the velocity constraints (i.e., the
acceleration constraints) explicitly. Another method for (1.4) has been proposed by
[9, 18, 19, 20], where P (q) is chosen to be an orthonormal basis of the local tangent
space of the constraint manifold. Choosing a smoothly varying P (q) is required and
may cause some practical difficulties.
Direct numerical solution of (1.4) via our coordinate-split approach yields an efficient
and reliable method for solving equations of motion for most multibody mechanical
systems. However, there is a class of multibody systems which present additional
computational challenges. These are the problems with high-frequency nonlinear os-
cillations. Highly oscillatory components are often used to model devices with strong
potential energy. Typical examples of such problems arise from modeling flexible
multibody mechanical, and molecular dynamic systems. For many problems, oscillations
of a sufficiently small amplitude are not important for the model, but they
severely restrict the stepsize for numerical methods. For these types of problems,
stiffly stable implicit numerical integration methods can be used to damp out the
oscillation [17]. However, the stepsize may still be severely restricted due to difficulties
in converging the Newton iteration for larger stepsizes. We have studied this
class of oscillating problems in [24]. The solutions are composed of a low-amplitude
high-frequency oscillation around the smooth solution [17, 21]. Along the smooth
solution, the eigenstructure of the local Jacobian matrix varies smoothly. However,
along the solutions which are nearby to the smooth solution, the local eigenstructure
oscillates with the high frequency, and is very badly behaved. The standard Newton
iteration inside a damping numerical method starts from a predictor which is on a
nearby solution, and attempts to find the smooth solution. It evaluates its Jacobian
matrix on the nearby solution, which determines the direction it takes toward the
smooth solution. Unfortunately, these Jacobian matrices do not yield good directions
for nonlinear oscillating problems as described above, unless the predictor is already
extremely close to the smooth solution. Thus, the standard Newton method must be
coupled with a severe reduction in the timestep to achieve an adequate predictor.
In Section 2, we outline the CS-iteration carried out by an efficient method for
its derivative. For the purpose of finding the smooth solution, we introduce
a modification to the Newton iteration, i.e., the CM-iteration. This iteration is easy
to implement, effective for non-oscillatory problems, and particularly effective for
nonlinear highly oscillatory problems. The basic idea of the CM-iteration is that
there are terms in the Jacobian which involve derivatives of the projection onto the
constraint manifold. These terms are large and complicated to compute, but small
on the slow manifold [21]. The CM-iteration sets these terms to zero, yielding a
reliable direction towards the smooth solution for the Newton-type iteration. We
prove the convergence of the CM-iteration and give an error estimate for the numerical
solution. In Section 3 we describe in more detail the structure of nonlinear oscillatory
mechanical systems, and derive estimates for the rates of convergence of the CS and
CM-iterations applied to these oscillatory systems. The difference in convergence
rate explains why the CM-iteration is highly effective for oscillatory systems, and
shows that its rate of convergence for non-oscillatory systems is similar to that of
the CS iteration. In Section 4, numerical experiments are given which demonstrate
the effectiveness of these methods, particularly for oscillatory nonlinear mechanical
systems.
2 The Coordinate-Split Technique
In this section we present the coordinate-split (CS) technique, which defines P (q)
in (1.4a) and (1.4b) via coordinate-splitting, and computes this matrix cheaply. Although
at first glance it would appear that implementation of this method would be
difficult due to complications in computing the derivatives of P (q), we show that the
special form of the pseudo-inverse can be used to give a much simpler derivation of
the Jacobian.
The construction of the annihilation matrix P (q) involves the solution of a class
of pseudo-inverses of the constraint Jacobian G(q). Given a smooth vector-valued
function r(q), a effective way to obtain the projected vector P (q)r is to use a splitting
of the original coordinates.
Definition 2.1 [Coordinate-Splitting Matrix] Let X and Y be the matrices whose
columns constitute the standard Cartesian basis of IR n\Thetan such that k(G(q)Y
bounded in a neighborhood U 0 of some q 0 , which is on the constraint manifold
0g. The p \Theta n coordinate-splitting matrix for (1.1) is defined by
Remark 2.1 Note that X and Y are piecewise constant with respect to q 0 2 M.
From the construction of the CS matrix P (q), one can easily see that P (q)G T
for all q 2 IR n , i.e., P (q) is orthogonal to range(G T ). Furthermore, the row vectors
of P (q) are orthonormal, i.e., P (q) T P is the identity matrix in IR p .
The computation of P (q) can be carried out using the LU-factorization or QR-
factorization of the constraint Jacobian matrix [12]. Then the projected vector P (q)r
can be computed relatively cheaply. In addition, directly applying the formulas given
in [11] (Theorem 4.3, pp. 420), we obtain
dq
dq
)s; with
where P (q) is defined by (2.1) and r(q) 2 IR n\Theta1 .
Using (2.1) and (2.2), the coordinate-split iteration for solving the nonlinear system
at each time step can be carried out. Applying, for example, a BDF formula to
yields the nonlinear system
ae h is the discretization operator, and h the stepsize of the time discretization.
Given an initial prediction (q (0)
applying Newton-type methods to (2.3) requires
the solution of a linear system
such that \Deltaq n and \Deltav n are the increments of q n and v n ,
@vn
and
We remark that the Coordinate-Split iteration leads
to a natural and reliable error estimator for the numerical integration method. In
particular, the local error estimator can be based on the independent coordinates and
velocities only.
To analyze the solution of (2.4), we rewrite the first two equations of (2.4), i.e.,
corresponding to the derivatives of (2.3a) and (2.3b),
\Deltaq n
\Deltav n
where the 2n \Theta 2n matrix J h is
h( dM(qn)
dvn
Since for small enough h, J h is invertible under the assumption of M(q n ) nonsingular,
the solution of (2.4) can be computed by
\Deltaq n
\Deltav n
where
Thus, (\Deltaq n ; \Deltav n ) can be obtained by the successive solution of two linear systems,
and
\Deltaq n
\Deltav n
The CS iteration is an efficient formulation of the Newton iteration for multibody
systems, so standard convergence theorems apply. We will now state conditions to
ensure the convergence of the CS iteration in the context of numerical integration
methods for (1.4). For simplicity we consider, instead of the second-order constrained
equations of motion (1.1), a first-order system
since the convergence of (2.9) can be trivially extended to (1.1). Applying stiffly
stable numerical methods, the convergence result is well-known, see [13] pp. 494-498.
Convergence of discretization methods for the index-1 system,
obtained by applying the coordinate-splitting matrix P (q) to (2.9a), is also well-
developed. By the construction of P (q), it is easy to see that the solution of the
CS iteration is equivalent to that of the local state-space ODE of the independent
coordinate x [7].
Applying a linear discretization operator ae h with stepsize h to (2.10) and scaling
(2.10a) by h yields a nonlinear system
where the residual function is
r h (q;
The convergence of the CS iteration can be carried out on a smooth constraint manifold
M. Assume that for any q 0 2 M, there exist h 0 and X 2 IR p\Thetan and Y 2 IR m\Thetan
such that
is
where J
Similarly, we assume that
for some C 2 , where q 1 and q 2 are in a neighborhood U(q 0 ) of q 0 .
Convergence of the CS iteration can be assured for a sufficiently accurate initial
guess. Under the assumptions (2.13) and (2.14), a sufficient condition for convergence
requires that the numerical integration satisfies
dr h
dq
dq
)s
for some C 3 ? 0 in a neighborhood of the true solution q   [6].
For highly oscillatory dynamic systems, the stepsize of a numerical method may
be restricted by the error test to follow the oscillation. Applying a stiffly stable
numerical method to (2.10), such as BDF of order - 2, one may be able to damp
the oscillation and take a larger stepsize to follow the trajectory of the equilibrium,
However, convergence of the Newton iteration requires a further
restriction on the stepsize. Depending on how close the predictor is to the equilibrium
of highly oscillatory components, the Newton direction imposed by the Jacobian can
excite the high-frequency oscillations. When applying the Newton method directly
to the discretization of the Lagrangian form (2.9), an even more severe problem in
Newton convergence is observed, and illustrated by the numerical experiments in
Section 4. The limitation on the stepsize due to the Newton convergence failures for
highly oscillatory nonlinear multibody systems can be overcome via a modification
to the CS iteration which we call the CM iteration.
The CM iteration
In large-scale multibody mechanical systems, most of the unwanted oscillations are
due to the noise of high-frequency forces, where the amplitude is well below the
solution tolerance. However, small perturbations in the position can cause drastic
changes in the Newton direction. This results in difficulties for convergence of Newton-type
methods. To remedy this problem in the CS iteration, we reduce the noise from
the oscillations by setting ( dP (q)
in the Newton iteration matrix, since it is
the main source contributing to the rapidly changing Newton direction. This term
is large away from the smooth solution but negligible on the smooth solution. The
CM iteration approximates the Newton direction of (2.10) via an oblique projection
to the unconstrained ODE
for a q 0 close to the solution q, e.g., G(q 0 )G T (q) invertible. When applying a stiffly
stable numerical integrator to highly oscillatory problems, this modification, for q 0
close enough to the smooth solution, overcomes the difficulties in the CS iteration.
Applying a stiffly stable discretization method to (2.16) coupled with constraint
(2.10b) leads to the nonlinear system
A convergence result for the modified CS iteration, denoted by CM , is given next.
We first give an upper bound on the difference between the derivative of the
projected vector P (q)r(q) and the projected derivative P (q)( dr(q)
Lemma 2.1 Suppose conditions (2.13) and (2.14) hold. Then for
is the disc in IR n with center q 0 and radius % 0 ,
dq
dq
Proof. From (2.2),
dq
dq
dq
)sk:
Since the row vectors of P (q) are p orthonormal vectors in IR n , applying the Cauchy
inequality gives
dq
dq
for all q 2 D(q Condition (2.13), with h ! 0, implies the result in
Convergence of the CM iteration and an estimate of the distance between the
solutions of (2.11) and (2.17) are given in the following.
Theorem 2.1 Suppose conditions (2.13)-(2.15) hold for some h 0 in a neighborhood
of q   , such that fq j g ! q   is generated by the CS iteration. Choosing -
sequence f-q k g generated by the CM iteration
dq , converges to -
q   . Furthermore, the distance between -
q   and q   is
bounded above by
for some moderate constant C.
Proof. Convergence of f-q k g follows immediately by noting that the CM iteration
is an inexact Newton method [5] with the property that the terms which the
CM-modification deletes from the Jacobian, hence the perturbation to the Newton
iteration, can be made as small as we want by taking the initial guess accurate enough.
Next, we show the estimate (2.20). Since the CS and CM iterations converge to
q   and -
respectively, we have
Subtracting
from both sides of the above equation we obtain
Expanding around q 0 yields
dq
dq
dq
dq r h (q). In matrix form, the above equations
become
O( -
By the assumption that -
J is invertible, we can write
):Note that since dP
when the constraints g(q) are linear, fq k g of the
CS iteration and f-q k g of the CM iteration are the same in that case. In general,
the rate of convergence of the CM iteration is superlinear, using the Dennis-Mor'e
Characterization Theorem [6].
3 Rate of convergence for highly oscillatory multi-body
systems
High frequency oscillatory forces often appear in the modeling of vehicle suspension
systems, modal analysis in structural dynamics, or modeling of oscillations in
computer-aided engineering etc. For simplicity, we consider the constrained dynamic
system of (1.1) with a dominant oscillatory force
may be, for example, the coefficients of stiff springs; i.e., In
practice, j(q) is usually oblique towards KerP (q), i.e., the oscillatory force(s) acts on
both the independent and the dependent coordinates. For the purpose of obtaining
a smooth solution with large stepsizes, we will show that the CM iteration can be
very effective for many classes of nonlinear oscillatory forces.
Deformation forces are the most common potential forces that can produce small
amplitude high-frequency oscillations. They are usually linear with respect to the
local coordinates [4, 25]. For these reasons, we consider the class of oscillatory forces
of the form
where the components of B and b 0 are slowly varying. In particular, B and b 0 may
be functions of some constraint-driven generalized coordinates. For example, B(') in
the 2D bushing problem in [24] has the form
\Gammak y sin ' k x cos ' k '54
are positive constants. When '(t) is smooth or constrained,
assumption (3.2) is valid.
Using a linear oscillatory force, the Lagrange equations of motion of the multibody
system can be written as
@(q;v) k. From assumption (2.14) on the constraint manifold, we can also
assume
dq
for all q.
In the context of the CS iteration, the problem of convergence of the Newton
iteration can be explained by analyzing the reduced potential function. The reduced
potential of (3.1) is
reduced potential force generated by (3.6) is
dq
At each iteration, the reduced potential force acts along the normal direction of the
constraint manifold. The gradient of the correction term yields
dq
Applying Y T to (3.8) gives
dq
and applying X T to (3.8) yields
dq
)s:
When high-frequency oscillations appear in the system, e.g., ffl ! 0, the reduced
potential force also becomes oscillatory if Y T r is nonzero. This is the general case
when the solution is not at an equilibrium position. Nevertheless, convergence of the
CS iteration can be achieved by using a small enough stepsize [17].
For small stepsizes, a convergence result for the CM iteration can be obtained
provided the assumptions of Theorem 2.1 are valid. In many applications, following
the oscillations is not of interest. Instead, one wants to use a large time step to damp
out the oscillations of small amplitude but high frequency. For this reason, we now
consider only the multistep numerical integration methods that are strictly stable at
infinity and A-stable, such as the lower order (i.e., - 2) BDF methods [13]. The
convergence of L-stable implicit Runge-Kutta methods to the smooth solution of the
highly oscillatory ODE of multibody mechanical systems can be found in [17]. Here
we focus on the convergence of the CM iteration for constrained multibody systems
with oscillatory forces when applying the above-mentioned linear multistep methods.
Numerical solutions on the slow manifold can be evaluated using the equilibrium
of (3.1), i.e., the slow solution [2, 15] satisfies
and the smooth solution is its asymptotic expansion to some order of ffl around the
manifold 0g. In the linear form, the smooth solution of (3.1) is not
far from
k. For the strongly damped numerical
solution 1. During the iterative solution onto the
slow manifold, the constraints may not be satisfied, which causes a large reaction
force in the form of (3.7). This may cause oscillations in the CS iteration, while
the CM iteration annihilates these nonlinear oscillations generated by the reduced
potential. This yields a superior performance of the CM iteration as compared to the
CS iteration for computing the smooth solution of (3.1). The result is explained in
the following.
Theorem 3.1 Consider the numerical solution of (3.1), with j(q) linear, and h the
stepsize of the stiffly stable multistep integration method. Assume the starting values
-neighborhood of the smooth solution. If both the CS and CM
iterations converge to q   and -
respectively, then the rate of convergence of the CS
iteration oe (CS) compared to that of the CM iteration oe (CM) is given by
k, and
k. In the case that the constraints are linear, the two
iterations yield the same result.
Proof. Recall that the CS iteration is a Newton iteration for the solution of
. Using the notation defined
in (2.8), and comparing the Jacobian matrix for the CS iteration at iteration k with
the Jacobian at the solution of the CS iteration, we obtain
For equation (3.1) solved by a multistep method with leading coefficient fi using the
CS iteration,
dq
and . Thus, from (3.9), the local rate of
convergence of the CS iteration is given, via the Contraction Mapping Theorem [6],
by
Now, recall that the CM iteration is a Newton iteration for the solution of
Comparing the Jacobian matrix for the CM iteration at iteration k with the Jacobian
at the solution of the CM iteration, we obtain
where for problem (3.1),
fiI \GammahI
Noting that in a -
-neighborhood of the smooth solution,
fiI \GammahI
equation (3.11) yields
Thus,
By a similar argument, this result holds for all linear numerical integration methods.4 Numerical Experiments
4.1 Point-mass with oscillatory force
The first example is a simple constrained multibody system under the influence of
a highly oscillatory force. Consider a unit point-mass constrained to the 2D unit
circle, using
and the constraint equation
1). The equations of motion are
is the applied force. Differentiating the constraint
with respect to time, an explicit form of the multiplier - is obtained by
For a highly oscillatory force f(q), one can see that -(t) is oscillating with the frequency
of f(q), and with amplitude proportional to the magnitude of kfk.
The numerical experiments are carried out using BDF of order - 2 in DASSL [3],
where the local error estimation has been modified. For the stabilized index-2 DAE
denoted by GGL, the local error is estimated using only the position, i.e., q.
Moreover, we have also included some experiments where the Newton convergence test
of GGL has been modified to exclude the multipliers. The corresponding numerical
solution is denoted by GGL   . For the coordinate-split and modified coordinate-split
iterations, denoted by CS and CM , respectively, the local error is estimated using
the independent variable X T q, as recommended in [24]. The CM iteration updates
the matrix P (-q) when a new Jacobian is required.
Linear oscillation
Let a unit gravitational force act along the negative y-direction, and apply a linear
oscillatory force
There is a stable equilibrium at 0]. The natural frequency of
the system is
ffl , and no dissipative force is present.
The numerical solution has been carried out with a moderate solution tolerance
. For a 0 to 0:25 second simulation, the results of
several combinations of the stiffness coefficient ffl are presented in Table 4.1, where
the initial values are 0]. The CM iterations show better
efficiency than those of CS, GGL and GGL   in all cases, i.e., comparing the numbers
of function and Jacobian evaluations in Table 4.1. In the Table, etfs and ctfs denote
the number of failures of the error test and Newton convergence test, respectively, in
DASSL. Comparing the results of GGL   with those of GGL, we observe an improved
Newton convergence. As ffl ! 0, i.e., for higher frequencies of the oscillation, the
CM iteration becomes even more efficient. In Figure 5.1, we plot the total energy of
each numerical solution. The CM iteration achieves the strongest damping because
DASSL is able to increase the stepsize faster with the CM iteration.
Linear spring force
In the next test, we replace the linear oscillatory force in the previous constrained
system by a spring forceffl 2
l
attachment point of the spring, l 0 the
natural length, and 1
the stiffness coefficient, as shown schematically in Figure 5.2.
Method TOL ffl no: steps no: fevals: no: jevals: no: etfs: no: ctfs:

Table

4.1: Results of the Constrained Point-Mass with a Linear Oscillatory Force
Time
(Logarithmic
Scale)
Total Energy (EPSILON=1.e-6, TOL=1.e-3, order=1,2)
CM
CS
GGL
GGL*

Figure

4.1: Total Energy Comparison of Linear Oscillatory Force Example;
and
Method TOL ffl no: steps no: fevals: no: jevals: no: etfs: no: ctfs:

Table

4.2: Results of the Constrained Point-Mass with an Oscillatory Linear Spring
Force
For unit mass and unit gravitational force, we set the spring attached to
(0; \Gamma0:5) and the natural length l 0:4, such that the equilibrium is at (0; \Gamma1; 0; 0).
Using the initial conditions [0:04471; \Gamma0:999; 0; 0], the results of the 0-0:05 second
simulation by the GGL, GGL   , CS, and CM iterations are shown in Table 4.2.
Because DASSL is able to increase the stepsize, and hence damp the solution faster
with the CM iteration, the CM method is quite effective in these tests. In Figure
5.3, we plot the total energy of each solution. The numerical solutions of x, w, and -
are presented in Figure 5.4.
4.2 Two-body pendulum with bushing
The second example is a two-body pendulum in 2D Cartesian coordinates. Six generalized
coordinates, locate the centers of mass and the
orientation of the bodies. The first body is grounded, and the second body is constrained
such that the distance between a point A of the first body and another point
B of the second body is fixed, and its orientation is held constant. This leads to five
constraint equations
Attachment Point
Spring
Y
Point-Mass
(0.0,

Figure

4.2: Constrained Point-Mass with a Linear Spring
Time
Total Energy v.s. Time
CM
GGL (dot)
GGL* (solid)

Figure

4.3: Total Energy Comparison of Oscillatory Spring Example;
x-velocity
Velocity
Time
lambda
Multiplier
x-coordinate
Position

Figure

4.4: Results of Oscillatory Spring Example;
where l AB and ' are constant, and
such that in the local reference coordinate systems of
body 1 and body 2, respectively. In this example, we use
l for the constraint equations (4.1).
We apply a nonlinear oscillatory force formulated using nonlinear beam theory [4].
This type of force arises commonly in flexible multibody dynamics [25]. As described
in (3.2), the deformation force between the ith and jth components is a function of
the relative displacement of the reference frames X 0
j , as shown
schematically in Figure 5.5. Typically, the relative displacement is measured by
are constant vectors with respect to the origins of the force reference
frames in their respective local coordinate systems, i.e.,
are the corresponding origins in a global coordinate system and A i and A j are
the transformation matrices from the global to the local coordinate system [10]. The
are calculated as
l) is the component of the kth row and lth column of A ij . The matrix
A ij is the relative orientation matrix of two force reference frames, i.e., B i and B j
are constant. The relative velocity is the time derivative of the relative displacement
and the relative angular velocity is are the
angular velocities of bodies i and j respectively.
Using the above defined notation, the force acting between the ith and jth components
due to the deformation can be written as
where K f is a 3 \Theta 3 structural stiffness matrix and C f is the 3 \Theta 3 damping coefficient
matrix. Similarly, the torque acting between the components is
Zi'
Zj'
Deformation force
reference frames
Component i Component j

Figure

4.5: Deformation Force of a Flexible Body
where K - and C - are analogous to K f and C f . Note that the force and torque in
this form are linear functions of the relative displacement (d ij and the relative
Here, the 2D bushing force has stiffness matrix
are O(1), and damping matrix
where c x , c y , and c ' are O(1). For this 2D bushing example, (3.3) becomes
ffl4
The attachment points of the force device are s 0= [0:5; 0] and s 0= [\Gamma0:5; 0] in the
body-fixed reference frames of bodies 1 and 2, respectively. The bushing force introduces
oscillatory applied forces, causing small oscillations of the numerical solution,
and yielding highly oscillatory multipliers in the index-2 DAE (1.3). The multipliers
associated with the highly oscillatory components exhibit high-frequency oscillations
with large amplitude. The standard convergence test of the Newton iteration depends
heavily on these multipliers. Therefore, we modified the convergence test in DASSL
to exclude the test for the multipliers. In addition, the multipliers are computed
in the GGL   by applying the pseudo-inverse (GY ) \GammaT Y T to r and to
is the predictor in DASSL, and ae h is the
discretization operator of BDF. The local error is estimated by the predictor-corrector
difference of (X T q; X T v) for CS, CM , and GGL   , and of (q; v) for GGL.
Using the initial values
numerical results with are shown in Table
4.3, where 10. For this moderate
stiffness all the methods perform well. Note that the constraint violation
with these initial values is O(10 \Gamma3 ). This implies that the difference of the constraint
reaction force at the initial value q 0 and that of a q   on the constraint manifold is
According to Theorem 3.1, the rate
of convergence of the CS iteration is proportional to ffi
ffl O(h). Increasing numbers of
convergence test failures in DASSL are expected as ffl ! 0. In this example, frequent
convergence test failures occured when We observe the
same difficulties in the Newton convergence of the GGL and GGL   iterations. On
the other hand, the CM iteration with its better Newton convergence, as explained
by Theorem 3.1, is able to take much larger time steps, and the nonlinear oscillation
is damped effectively. In Table 4.4, the results of
are shown. In Figure 5.6, we plot the stepsize taken by DASSL for
GGL, GGL   , CS, and CM using the stiffness coefficient
Method TOL no: steps no: fevals: no: jevals: no: etfs: no: ctfs:

Table

4.3: Results of Two-Body Pendulum with a Bushing Force,
Method TOL no: steps no: fevals: no: jevals: no: etfs: no: ctfs:

Table

4.4: Results of Bushing Problem,
time
stepsize
GGL .x.
CS -o-
GGL* -.

Figure

Time Steps Used in Solving the Bushing Problem,



--R

Projected implicit Runge-Kutta methods for differential-algebraic equations
Asymptotic Methods in the Theory of Nonlinear Oscillations
Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations
John Wiley

Numerical Methods for Unconstrained Optimization and Nonlinear Equations
On the error control for multistep methods applied to ODEs with invariants and DAEs in multibody dynamics
Automatic integration of the Euler-Lagrange equations with constraints
Index reduction methods for differential-algebraic equations

The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate
Matrix Computations
Solving Ordinary Differential Equations II: Stiff and Differential- Algebraic Problems
Computer Aided Kinematics and Dynamics of Mechanical Systems Volume
Perturbation Methods in Applied Mathematics
The Variational Principles of Mechanics
Integration of stiff mechanical systems by Runge-Kutta methods
On the numerical solution of the Euler-Lagrange equations
On the numerical solution of the Euler-Lagrange equations
On a geometric interpretation of differential-algebraic equations
approach to the torsion dynamics of macromolecules
Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems
Constrained equations of motion in multibody dynamics as ODEs on manifolds
Computational challenges in the solution of nonlinear oscillatory multibody dynamic systems
Dynamics of articulated structures
--TR

--CTR
B. Fox , L. S. Jennings , A. Y. Zomaya, On the modelling and computation of nonlinear multibody cable systems, International Journal of Computers and Applications, v.28 n.2, p.154-161, April 2006
highly oscillatory systems;error bound;differential-algebraic equations;initial value problems;multibody systems