MSC Adams Solver PDF
MSC Adams Solver PDF
MSC Adams Solver PDF
Ann Arbor
August, 2004
Dan Negrut
Dan.Negrut@mscsoftware.com
Andrew Dyer
Andrew.Dyer@mscsoftware.com
1.
Introduction ................................................................................................................. 3
1.1.
Solving Nonlinear Equations. The Newton Method .......................................... 3
2. Definitions, Notations, Conventions ........................................................................... 9
2.1.
Generalized Coordinates used in ADAMS ......................................................... 9
2.2.
Joints in ADAMS.............................................................................................. 11
2.3.
Motions in ADAMS.......................................................................................... 12
3. Initial Condition Analysis ......................................................................................... 13
3.1.
Position Initial Condition Analysis ................................................................... 13
3.2.
Velocity Initial Condition Analysis .................................................................. 17
3.3.
Force and Acceleration Initial Condition Analysis ........................................... 18
4. Kinematic Analysis ................................................................................................... 20
4.1.
Position-Level Kinematic Analysis .................................................................. 20
4.2.
Velocity-Level Kinematic Analysis .................................................................. 20
4.3.
Acceleration Level Kinematic Analysis ........................................................... 21
5. Dynamic Analysis ..................................................................................................... 22
5.1.
Nomenclature, Conventions, Definitions. ......................................................... 22
5.2.
Formulation of Equation Of Motion in ADAMS ............................................. 23
5.3.
Numerical Solution for Dynamic Analysis. Jacobian Computation. ............... 24
6. Statics Analysis ......................................................................................................... 28
6.1. The STATIC approach .......................................................................................... 28
6.2. The DYNAMIC Approach ................................................................................... 29
6.3. The STATIC_HOLD Attribute ............................................................................. 29
7. A Simple Pendulum Example ................................................................................... 31
7.1.
The Simple Pendulum Model ........................................................................... 31
7.2.
Initial Condition Analysis ................................................................................. 32
7.3.
Dynamic Analysis ............................................................................................. 36
Conversion to First-order System ............................................................................. 41
7.4.
Kinematic Analysis ........................................................................................... 46
7.5.
Static Analysis .................................................................................................. 48
8. Ways and Means for Improving your MSC.ADAMS Simulation ............................ 50
8.1.
Settings in the INTEGRATOR Statement ........................................................ 50
8.2.
Settings in the EQUILIBRIUM Statement ....................................................... 51
8.3.
Rules of Thumb for Creating Robust Models in ADAMS ............................... 53
8.4.
Validating your Simulation Results .................................................................. 59
8.5.
Debugging Support in ADAMS/Solver ............................................................ 61
8.5.1.
DEBUG/EPRINT ...................................................................................... 61
8.5.2.
DEBUG/RHSDUMP ................................................................................ 64
8.5.3.
DEBUG/JMDUMP ................................................................................... 65
Appendix A. Integration Jacobian ................................................................................... 66
1. Introduction
The purpose of this document is to introduce the ADAMS user to the theory that
underlies the ADAMS/Solver. A basic understanding of the theory and the ways in
which the solver works can help the user make good modeling decisions, or choose a set
of simulation parameters that will lead to faster and more reliable simulations.
This document provides in section 1.1 a very brief introduction to the NewtonRaphson method for the solution of non-linear systems. Virtually all analysis modes in
ADAMS use this algorithm, and the user is strongly encouraged to understand the basics.
This is essential for understanding the informational and debugging messages that
ADAMS/Solver provides, why an equilibrium analysis is challenging, how changing the
Jacobian evaluation pattern is going to impact the convergence properties of the solver,
and list goes on.
In the first part of the document (sections 2 through 6) the focus is on equation
formulation and solution. The reader is introduced to the concept of generalized
coordinates, as well as the equations that govern the most representative types of analyses
that can be performed by the ADAMS/Solver; i.e., Initial Condition Analysis, Kinematic
Analysis, Dynamics Analysis, and Static/Quasi-static Analysis. These equation sets are
rather abstract, and therefore this discussion is followed in section 7 by a simple test case
that reinforces in a simplified two-dimensional framework the process of deriving the
equations for each analysis type.
Section 8 discusses ways in which an ADAMS/Solver simulation can be made to
run faster and more robustly. Subsections 8.1 and 8.2 are focused on dynamics and
statics simulations as they are typically more challenging. The discussion will cover
relevant parameters associated with the INTEGRATOR and EQUILIBRIUM statements.
Subsection 8.3 provides a set of recommendations and good practices for ADAMS
modeling. Subsection 8.4 discusses the issue of validating your results in ADAMS, while
subsection 8.5 explains how the user can get useful debug information from
ADAMS/Solver.
(1)
x (0) of the root is provided, a new configuration x (1) , hopefully closer to the root x is
computed as
x (1) = x (0)
f ( x (0) )
f ( x (0) )
(2)
where f () is the derivative of the function with respect to the variables on which it
depends.
This strategy of updating the value of x is obtained by linearizing the function f at the
point x (0) . Thus,
(3)
Rather than solving f ( x) = 0 , the root of the right hand side expression in Eq.(3) is
sought. This leads to a linear equation, and its root is denoted by x (1) :
f ( x (0) )
f ( x (0) )
which is the configuration update defined in Eq.(2). The algorithm continues by setting
x (0) x (1) and performing another iteration as indicated in Eq.(2), and illustrated in
Figure 1.
The method just described is called Newton-Raphson, or sometimes simply
Newton. In general, this approach leads to a quadratic rate of convergence to the root x
provided the iterative process is started close enough to the root. The drawback of the
method is that as indicated in Eq.(2), the new update x (1) requires the computation of the
derivative both the function f ( x (0) ) and the derivative f ( x (0) ) . At the cost of one
function and one derivative evaluation per iteration, this method is considered to be rather
expensive. An alternative is to only update the derivative once in a while, and to recycle
the derivative computed at an older configuration for several iterations. This method is
called Newton-like, and the convergence rate is linear.
Referring to the ADAMS/Solver documentation, consider the INTEGRATOR
statement. One of its attributes is the PATTERN setting, which enables the user to
dictate how often the derivative is to be updated. Computing the derivative (also called
the Jacobian), is in general an expensive operation and should be done as seldom as
possible. On the other hand, if its done too seldom, the speed of convergence will suffer
as there will be degradation from quadratic (for Newton-Raphson), to linear rate of
convergence (for Newton-like methods). Figure 1 shows a geometric interpretation of
the Newton-Raphson method: starting with the configuration x (0) , the nonlinear function
is locally approximated by a straight line. The new configuration x (1) is taken to be the
place where the straight line intersects the x -axis (the solution of the linear equation of
Eq.(3)). Moving on, the new tangent is determined; i.e., a new derivative f ( x (1) ) is
computed, and the new configuration x (2) is obtained. The Newton-Raphson iteration
sequence is thus carried out as:
(1)
(2)
=x
(0)
f ( x (0) )
f ( x (0) )
f ( x (1) )
=x
f ( x (1) )
(1)
x (3) = x (2)
f ( x (2) )
f ( x (2) )
f ( x( n ) )
, or the value of the of the function,
f ( x ( n ) )
f ( x ( n ) ) , is below a specified threshold (either exit criteria could be used, in general).
Note in Figure 1 that the computed configurations get closer and closer to the point where
the non-linear function f ( x) intersects the horizontal axis, that is it gets closer and closer
to the solution of f ( x) = 0 .
In Figure 2, the Newton-like method is illustrated by adding dotted lines to the
solid Newton-Raphson lines from Figure 1. The dotted line represents the slope of the
curve as computed at the point x (0) , which is used to find not only x (1) (as is the case
with the Newton-Raphson method), but also x (2) , x (3) , and x (4) . Note that the slope is reevaluated at x (4) and the iterative process continues with finding x (5) .
The interesting thing to notice is that the slope of the line stays the same for four
consecutive iterations; in ADAMS syntax, this would be equivalent to specifying a
PATTERN=T:F:F:F, (this is, by the way, the default setting for this attribute in the
ADAMS/Solver INTEGRATOR statement).
(1)
=x
(0)
f ( x (0) )
f ( x (0) )
x (4) = x (3)
f ( x (3) )
f ( x (0) )
x (5) = x (4)
f ( x (4) )
f ( x (4) )
f1 (q)
f (q)
2
f (q) = = 0
f n1 (q)
f n (q)
(4)
f
f
= i
q q=q( 0) q j q=q( 0)
(5)
x
p = y
z
(6)
The orientation of a rigid body is defined by a set of three Euler angles that correspond to
the 3-1-3 sequence rotation: , , and , respectively. These three angles are stored in
an array
(7)
The set of generalized coordinates associated with rigid body i in ADAMS is denoted in
what follows by
p
qi = i
i
(8)
(9)
= B B
(10)
where
sin sin
B = cos sin
cos
cos
0 sin
1
0
(11)
and is the body angular velocity expressed in the body-fixed coordinate system.
Equation (11) is important as it defines the relationship between the angular velocity of
the body (an intrinsic characteristic of the body), and the choice of generalized
coordinates (which, as mentioned earlier, can be chosen in a variety of ways).
Finally, note the relationship between the time derivative of the body orientation matrix
A and angular velocity :
= A
(12)
where the orientation matrix A is defined in terms of the 3-1-3 Euler rotation sequence
, , and as:
sin sin
sin sin
cos sin
cos
0
a = a3
a2
a3
0
a1
a3 ] ,
a2
a1
0
Also, a vector quantity marked with an over-bar indicates that the vector is expressed in a
local body reference frame.
10
q = q1T
T
= [ q1
q T2 q nb
q2 qn ]
(13)
with n = 6 nb will describe at a given time the position and orientation of each body in
the system.
(q ) = 0
(14)
( q ) = 1T ( q ) 2T ( q ) njT ( q ) = 1 ( q ) 2 ( q ) m ( q ) (15)
where nj is the number of joints in the model, and m is the sum of the number of
constraints induced by all joints. Note that q R n , while R m . Typically, m < n ;
i.e., the number of generalized coordinates is larger than the number of constraints they
must satisfy.
By taking one time derivative of the position kinematic constraint equations of Eq.(15),
the velocity kinematic constraint equations are obtained as
qq = 0
(16)
By taking yet another time derivative of Eq.(16), the acceleration kinematic constraint
equations are obtained as
= ( qq ) q
qq
q
(17)
Equations (15) through (17) can be seen as conditions that the generalized coordinates
array q along with its first and second time derivatives must satisfy. This is to ensure
11
that the evolution of the mechanical system makes sense; i.e., the mechanism is
assembled, and the parts move such that the constraints imposed by joints are obeyed at
every time.
( q, t ) = 0
(18)
Revisiting the definition of the position, velocity, and acceleration kinematic constraint
equations, both joints and motion constraints may, in general, be written as:
( q, t ) = 0
(19)
q ( q, t ) q = t ( q, t )
(20)
= ( qq ) q 2qt q tt ( q, t )
q ( q, t ) q
q
(21)
Equations (20) and (21) are obtained by taking one and respectively two time derivatives
of the position kinematic constraint equation of Eq.(19). A set of generalized coordinates
is said to be consistent, if it satisfies the position kinematic constraint equations.
Likewise, a set of generalized velocities is considered consistent if, for a consistent
position configuration, q satisfies the velocity kinematic constraint equations of Eq.(20).
12
( q, t0 ) = 0
(22)
while for the generalized velocities to be consistent, they must satisfy the velocity
kinematic constraint equation
q ( q, t ) q = t ( q, t )
(23)
q14 = 1.0 , q16 = 0 , the solver should assemble the mechanism and, at the same time, do
its best to satisfy the prescribed conditions. This IC analysis is solved in ADAMS via an
optimization approach. The constrained optimization problem solved minimizes the cost
function
f ( q1 , , qn ) =
2
2
1
1
w1 ( q1 q10 ) + + wn ( qn qn0 )
2
2
(24)
subject to the constraint equations ( q,t0 ) = 0 . In Eq.(24), the values wi are weight
factors, while qi0 can be regarded as the initial generalized coordinates inducing an initial
configuration of the system. Notice that this configuration q 0 = q10
not be consistent (i.e., it may not satisfy the constraints).
13
The user may prescribe that some of the entries in the q 0 array; i.e., qi01 , qi02 , etc., are to
be regarded "exact". As it will be justified shortly, the corresponding weights wi1 , wi2 ,
etc., will be given large values; i.e., values like 1010 . The remaining weights, namely the
ones corresponding to generalized coordinates qi0 that the user did not specify are given
values of 1 . With this approach, the constrained optimization problem solution will keep
the user-imposed IC values almost unchanged, while adjusting the "free-to-change"
generalized coordinates. Notice that the solution of the problem means minimizing the
cost function while satisfying the constraints.
In matrix notation, the constrained optimization problem reads:
For q R n , minimize
f (q ) =
T
1
q q0 ) W ( q q0 )
(
2
(25)
subject to
( q,t0 ) = 0
(26)
W = diag ( w1 , w2 , , wn )
(27)
while the constraints of Eq.(26) that must be satisfied in the optimization problem are
exactly the position kinematic constraints of Eq.(19).
ADAMS approximates the non-convex optimization problem of Eqs.(25) and (26) by a
succession of convex problems that are guaranteed to have a solution, which can
potentially be found in one iteration. Thus, the set of non-linear constraint equations of
Eq.(26) is linearized in the vicinity of q 0
( q, t0 ) = ( q 0 , t0 ) + q ( q 0 , t0 )( q q 0 )
(28)
With Equation (28) replacing Eq.(26) and the notation d q q 0 the now convex
optimization problem reads
Minimize
1
f ( d ) = d T Wd
2
14
(29)
( q 0 , t0 ) + q ( q 0 , t0 ) d = 0
Subject to
(30)
To solve the convex constrained optimization problem of Eqs.(29) and (30), define the
optimization Lagrangian:
F ( d, ) = f ( d ) + T ( q 0 ) + q ( q 0 ) d
(31)
=0
d
(32)
=0
which lead to the following linear system of equations:
W
qT ( q 0 ) d 0
=
0
q ( q 0 )
( q )
0
(33)
Based on Eq.(33), ADAMS computes the value of d , and given q 0 computes the
solution of the convex optimization problem as
q = q0 + d
(34)
Typically, after a couple of iterations the solution of the linear convex optimization
problem will satisfy the non-linear constraint equations of the mechanical system. The
approach fails when the linearization in Eq.(28) is a poor approximation of the non-linear
manifold solution of the original constraint equations. However, even when the manifold
is highly non-linear, the solution sequence will converge if the starting point q 0 is close
enough to the final solution. For this reason, it is essential for the algorithm to have a
good starting point.
15
To gain an understanding of how the weights wi help to keep the user defined initial
conditions at their prescribed values, consider a case with only one constraint. Likewise,
assume that the system has two generalized coordinates x and y , and that the constraint
equation they must satisfy is ( x, y ) = x 2 y = 0 . Obviously this constraint equation is
satisfied by an infinite number of pairs like (2,4), (2.5, 6.25), etc., but in this simple case
the user would like to specify that the value of x is x 0 = 1 , while the value y 0 is free to
change. As the value for y 0 is not prescribed, assume our initial guess takes y 0 = 6 .
Since the user prescribed a value for x , the weight associated with this generalized
coordinate is large; i.e., w1 = 1010 . Since there is no condition imposed on the second
generalized coordinate, w2 = 1 . Notice that prescribed in the discussion above refers to
coordinates specified to be exact in the ADAMS modeling language (as in the ADAMS
.adm file).
1010
0
2
2 d1 0
1 1 d 2 = 0
1 0 5
0
Then d1 = 10 9 , d 2 = 5 , = 5 , and
xIC = 1 + 109 1
yIC = 6 5 = 1
As can be easily verified, ( xIC , yIC ) = 2 109 + 1018 . The non-linear constraint
equation is very well satisfied, and the correction applied in the user prescribed initial
condition x is negligible (order 109 ). In this manner, the weights wi bias the solution
toward changes in one subset of generalized coordinates rather than another.
It is worth noting that during position IC analysis ADAMS checks the compatibility of
the constraint equations induced by the joints and motions present in the model. During
16
this constraint analysis, some of the constraint equations might turn out to be redundant.
The benign case is when a redundant constraint is consistent. An example of such a case
is when ADAMS is presented with one constraint equation that looks like:
x2 y = 0
(35)
(36)
The constraint in Eq.(36) does not add anything to the picture; when the first equation is
satisfied the second one is automatically satisfied too. Thus, the second equation is
redundant, but consistent, and throughout the simulation ADAMS will monitor this
equation to make sure that the redundant constraint continues to be consistent.
The solver will fail if the previous constraint equation is replaced by:
2x2 2 y = 1
(37)
When incompatible redundant constraint equations are found, ADAMS/Solver will carry
out an LU factorization with full pivoting and inform the user about encountering this
situation. ADAMS will stop the simulation upon finding incompatible redundant
constraints because from a modeling perspective, there is something qualitatively wrong
with the system being simulated.
Redundant constraints are usually encountered when too many joints are used to model
the mechanical system and the number of constraint equations generated by these joints
exceeds the number of generalized coordinates of the model. In what follows, two or
more constraint equations will be called independent if they are not redundant.
It is highly advised for the ADAMS modeler to eliminate redundant constraints. While
sometimes the redundant constraint(s) can be benign, they can often have detrimental
effects on the simulation as the reaction forces may not be calculated as intended, causing
such things as poor simulation performance, unexpected eigenvectors from
ADAMS/Linear, etc. More information on this can be found in Section 8.
position constraint equations (see Eq.(28)) and the solution is guaranteed to be found in
one iteration. Thus, the convex constrained optimization problem solved to retrieve the
initial velocities qi minimizes the cost function
f ( q1 , , qn ) =
1
T
( q q 0 ) W ( q q 0 )
2
(38)
q ( q, t0 ) q + t ( q, t0 ) = 0
(39)
As in the displacement IC analysis, the weight diagonal matrix W has several very large
positive entries that ensure that the user-prescribed initial velocities are not changed
significantly by the optimization algorithm. From here on, the same procedure
previously used for position IC analysis is employed. Note that, because of the linearity
of the velocity kinematic constraint equations the algorithm is guaranteed to converge in
one iteration.
An exception to the linearity of this problem is if the ADAMS user includes a GCON
(generalized constraint) statement which defines a non-linear velocity constraint. In this
case, obviously a non-linear solver is used to generate the consistent set of ICs.
M
qT ( q 0 ) q
F
=
q ( q 0 )
0
(40)
where M is the generalized mass matrix. Note that this is a linear system, and the
iterative process involved typically converges in one iteration. The user cannot directly
prescribe any initial acceleration for the force/acceleration IC analysis. The reaction
force and initial accelerations are evaluated based on the computed initial position, initial
velocity, and the applied force acting on the system at the initial time. For a more
detailed explanation of how the EOM are obtained (the first row in the Eq.(40)), see the
Section on Dynamic Analysis in ADAMS.
, the solution of the linear system above also provides the Lagrange multipliers
Besides q
. The constraint force and torque induced by joint j on body i are computed as:
T
( j) ( j)
F =
vi
C
18
(41)
( j) ( j)
T =
i
C
(42)
In Eqs.(41) and (42) the superscript C indicates that the quantities are expressed in a
Cartesian coordinate system; v i is the Cartesian velocity of body i ; i is the global
angular velocity; ( ) represents the set of constraint equations associated with joint j .
An explanation of why the reaction forces expressed in the global reference frame are
computed as indicated in Eqs.(41) and (42) is beyond the scope of this discussion, but it
suffices to say that they are obtained by projecting the Lagrange multipliers ( j ) along
the Cartesian translational and rotational velocities v i , and i of body i , respectively.
j
19
4.
Kinematic Analysis
( q1 , t1 ) = ( q 0 , t1 ) + q ( q 0 , t1 )( q1 q 0 )
(43)
Since the number of constraints is equal to the number of generalized coordinates, the
matrix q ( q 0 ,t1 ) is square. As the constraint equations were assumed to be independent,
this matrix is also invertible. Based on an explicit integrator (e.g., Forward-Euler) an
0
initial starting configuration q1( ) is determined, and the iterative algorithm proceeds at
each iteration j 0 by finding the correction (
j)
q ( q 0 , t1 ) ( ) = q1( ) , t1
j
Then, q1(
j +1)
(44)
j)
As with position IC analysis, ADAMS can fail to find q1 if the linearization at the initial
guess turns out to be a poor approximation of the non-linear manifold. In these
situations, the remedy lies in decreasing the simulation step-size, causing q1 to lie closer
on the manifold to the last consistent configuration, q 0 .
20
analysis, the non-singular matrix q ( q1 ,t1 ) is evaluated and the linear system of Eq.(20)
is solved for the new velocity q1 .
qT = F Mq
(45)
This equation is identical to the first row of the linear system of Eq.(40) and in fact
represents precisely the equations of motion. More information on how Eq.(45) is
obtained is provided in the next Section.
21
5.
Dynamic Analysis
1 T
1
u Mu + T J
2
2
(46)
f
F ( q, q , t ) = R 6 - the vector of applied forces; f is the vector of applied forces
n
expressed in the global reference frame, while n represent the applied torque expressed
in the local Cartesian reference frame
( P )T f
Q=
R T n
)
(
(47)
where with v P being the velocity of the point of application P of the external force F ,
the projection operators are computed like
P =
v P
u
(48)
R =
(49)
22
d K
dt q
K T
+ qT = Q
(50)
K T K T
P T
T
d u p p ( ) f
+
=
dt K T K T T R T n
)
(
(51)
It is worth pointing out that when dealing with a full system of rigid bodies connected
through joints, the system Equation Of Motions (EOM) are obtained by simply stacking
together the EOM for the bodies in the system.
Since
T
d K
= Mu
dt u
(52)
K
p = 0
(53)
K
= BT JB
(54)
Mu + pT = ( P ) f
T
K
+ T = ( R ) n
(55)
The first order differential equations above are called in what follows kinetic differential
equations, and they indicate how external forces determine the time variation of the
translational and angular momenta.
Finally, the time variation of the generalized coordinates is related to the translational and
angular momenta by means of the kinematic differential equations. By assembling the
kinetic and kinematic differential equations ADAMS generates a set of 15 equations for
23
each rigid body that provide the information necessary to find a numerical solution for
the dynamic analysis of a mechanical system. These equations are as follows:
T
Mu + pT ( P ) f = 0
(56)
BT JB = 0
(57)
T
K
+ T ( R ) n = 0
(58)
p u = 0
(59)
= 0
(60)
24
A second, more refined algorithm reduces the original index 3 problem to an analytically
equivalent yet numerically different index 2 DAE problem. Thus, instead of considering
the position, the velocity level kinematic constraint equations of Eq.(20) are solved for
along with the kinematic differential equations. In the ADAMS solver, this algorithm is
called SI2, and while typically slower than the index 3 approach it turns out to be more
accurate and robust.
In what follows the index 3 approach is presented in a reasonable amount of detail. In an
attempt to keep the presentation simple, the index 3 DAE will be integrated via an order 1
implicit integration formula, which converts the DAEs into a set of algebraic equations.
This formula is the backward Euler formula - a one-step, A-stable algorithm that
qualitatively captures all the relevant details characteristic to higher order methods.
Backward Euler integration formula replaces the derivative y 1 at time t1 with
y 1 =
1
1
y1 y 0
h
h
(61)
(62)
The system of equations in Eq.(62) is called a discretization system since the derivative
in the original IVP problem was discretized using the integration formula of Eq.(61).
Since almost always the function g is non-linear, a non-linear algebraic system needs to
be solved to retrieve y1 . This is done in ADAMS by using a Newton-Raphson type
iterative algorithm.
Based on the implicit Euler discretization formula introduced above, all the first order
time derivatives that appear in the equations of motion in Eqs.(56) through (60) are
discretized to produce a set of algebraic non-linear equations. In the index 3 approach,
the position kinematic constraint equations are appended to these equations, along with
the force function definition F and T . This appending of the force functions is done for
the sole purpose of increasing the number of unknowns and thus inducing a larger but yet
sparser Jacobian matrix. Thus, after the implicit Euler based discretization Eqs.(19), and
(56) through (60) along with the force/torque definition equations assume the following
form:
25
T
1
1
Mu Mu 0 + pT ( P ) f = 0
h
h
B T JB = 0
T
1
1
K
T
R T
0
+ ( ) n = 0
h
h
1
1
p p0 u = 0
h
h
1
1
0 = 0
h
h
( p, , t1 ) = 0
(63)
f F ( u, , p, , f , n, t1 ) = 0
n T ( u, , p, , f , n, t1 ) = 0
The unknowns in this non-linear system are u , , , p , , , f , n . The subscript 1
indicating the time step was dropped for convenience.
Introducing the array
u
p
y=
f
n
(64)
(y) = 0
(65)
A Newton-Raphson type algorithm finds the solution of this system. First a prediction
0
y ( ) of the solution is computed, typically by using a predictor constructed around an
explicit integrator. Once an initial guess of the solution is provided, iterations
( ))
y ( y 0 ) ( ) = y (
j
( j +1)
=y
( j)
( j)
(66)
( ))
are carried out until the correction ( ) are small enough (note that the residual y (
j
26
process might fail to converge. This is more likely to happen with dynamic analysis than
with other types of analysis, as it is clear that the system that needs to be solved at each
integration step is highly non-linear.
If the iterative process fails, the integration step-size is decreased and another step is
attempted. ADAMS users are familiar in this context with messages informing them that
the step-size was decreased too much, and yet the convergence was not attained. Getting
such a message is a bad omen, as typically the user will have to revisit the model, make
modifications in the simulation defining parameters, or to try a different integrator like
SI2, for example. More details on this can be found in Section 8.
Finally, although the discretization formula used to convey the dynamic analysis solution
message was backward Euler it conceptually captures the essence of the ADAMS
solution sequence. ADAMS typically uses higher order integrators whenever the signals
sent over by the problem being solved suggest that this would improve performance. The
expression of the integration Jacobian is qualitatively the same, with very minor and
insignificant changes for example the denominator of the fraction 1 h would become
27
6.
Statics Analysis
pT ( A P ) f = 0
B T JB = 0
T
K
T
R T
+ ( A ) n = 0
u=0
=0
(67)
( p, , t 0 ) = 0
f F ( u , , p , , f , n , x, t 0 ) = 0
n T ( u , , p , , f , n , x, t 0 ) = 0
d(u, , p, , f , n, x, t0 ) = 0
Note that unlike Eqs. (56) through (60), this set of non-linear equations was modified to
include the ADAMS DIFF elements, which for the purpose of this discussion were
assumed to have the form
x d(u, , p, , f , n, x, t ) = 0
Here x is the state associated with the model DIFFs, and for simplicity the DIFFs were
assumed in explicit form. Setting x = 0 in the definition of the DIFFs leads to the last
equation in Eq.(67). It is important to underline that the DIFFs should have been present
28
in the dynamic analysis discussion in section 5.2 as well, and they were omitted there
only to keep the presentation simple.
In the STATIC approach, an equilibrium configuration
y T = [u p f
n x]
is found as the solution of the algebraic non-linear system of Eq.(67), and to this end the
Newton-Raphson algorithm of section 1.1 is used. Recall that for this algorithm to work,
the initial starting point of the iteration sequence should be close enough to the solution.
This is precisely what makes an equilibrium analysis challenging. Unless the user
ensures that the initial mechanical system configuration is close to equilibrium, the
algorithm might fail.
NOTE: The numerical solution of the dynamic analysis discussed in detail in section 5.3
also requires the solution of a similar non-linear system during the integration corrector
stage. The situation is more tractable in the dynamic analysis case due to the predictor
employed by the integrator. It is one of the roles of the predictor to produce a good initial
starting point for the Newton-Raphson algorithm. The convenience of a predictor is not
available in the STATIC approach, and therefore the user has to either (a) provide a good
starting point by setting up the model to be in a configuration close to equilibrium, or (b)
change the EQUILIBRIUM parameters to control the convergence of the NewtonRaphson method (see sections 8.2 and 8.3 for a discussion in this sense)
29
In this example DX represents the distance along the X-axis between markers 23 and 11
(assumed to be defined somewhere else in the model). If the STATIC_HOLD is not
present in the above definition, during an equilibrium analysis the Solver will adjust the
value of DIFF(10), that is, the state associated with the DIFF such that the expression
DIFF(10)-DX(23,11) will become zero. This is because at equilibrium, the time
derivative of any quantity should be zero, and therefore DIFF1(10), which according to
the DIFF definition is equal to DIFF(10)-DX(23,11), should be zero. Thus, the value
IC=2.0 will be overridden by the solver, and the actual value of the state associated with
the DIFF after the statics analysis; i.e., the true IC value will be DX(23,11). Sometimes
this is what the user wants, but sometimes keeping the value of IC to what was originally
set is preferable. This is precisely what STATIC_HOLD does. Under these
circumstances, for all purposes the equilibrium analysis proceeds as though the DIFF
element is not present in the model. Furthermore, whenever the quantity DIFF(10) is
referenced during the static analysis by another modeling element, for instance a force
definition, the value returned would be that indicated by the IC (2.0 in our example),
while the quantity DIFF1(10) will be zero. The key observation when STATIC_HOLD
is present is that even if the Solver indicates that the static analysis converged, the
residual in satisfying the DIFF equation is not zero. In our case, the residual would be
2.0-DX(23,11), which would be zero only by chance.
Finally, if the STATIC_HOLD attribute is not present in the definition of a DIFF,
when the static analysis will have converged the value of the state associated with the
DIFF would not be 2.0, but precisely the value DX(23,11). The function DIFF1(10)
would evaluate to 0.0, and there would be a small residual in satisfying the DIFF, but it
would be at the most as large as the one specified in the IMBALANCE setting associated
with the EQUILIBRIUM statement/command.
30
7.
Figure 3
31
R
R/C
C
Figure 4
R = C + R/C
or, equivalently:
32
C + R/C R = 0
Where:
The constant R is a vector which locates the pin (revolute joint) in the ground
reference frame
The variable C is a vector which locates the center of mass of the pendulum in the
ground reference frame
The variable R/C is the vector from the center-of-mass of the pendulum to the pin
(revolute joint) attaching it to ground.
The condition imposed by the mechanism can be rewritten in vector form as the
constraint equation:
( x, y, ) = ( x + R / C1 R1 )i + ( y + R / C2 R2 ) j = 0
R / C1 = l cos
R / C2 = l sin
where i and j are unit vectors along the x and y coordinate axes in the global (ground)
reference frame. We can, alternatively, express this vector constraint in matrix notation:
x + R / C1 R1
( x, y , ) =
=0
y + R / C2 R2
The Jacobian matrix for is the 2 by 3 array:
1 0 l sin
q =
=0
0 1 l cos
Lets now formulate the constrained optimization problem for our constraint equations:
w1
0
1
0
0
w2
0
0
1
0
0
w3
l sin 0
l cos 0
1
0
l sin 0
0
0
0
0
d1
d
0
1
0
l cos 0 d3 =
0
0
0
1 x l cos R1
2 y 0 l sin 0 R2
0
For example, lets say we have a broken joint, and the pendulum is horizontal, length
2l=2, and we have decided to fix the orientation such that:
33
0 = 0 (horizontal )
w3 = 1010 ( fixed )
R1 = 4
= broken
x0 = 6
R2 = 0
= broken
y 0 = 2
w1 = w2 = 1( free)
l =1
In other words, our pendulum is horizontal, and the pendulum center of mass position
(x,y) is at (6,-2), while the revolute joint location is at (4,0). This puts the end of the
pendulum of length 2l at (5,-2), which is not coincident with the revolute location, as
required to make the constraint consistent.
Notice that we have put a large weight for the orientation of the pendulum, which is
achieved by fixing orientation like so:
34
1
0
1
0
0 0
1 0
0 1010
0 0
1 1
1 0 d1
0
0 1 d2
0
0 1 d 3 =
0
0 0 1 6 1 4 = 1
0 0 2 2 0 0 = 2
d1 = 1 = 1
d 2 = 2 = 2
d3 =
2
10
= 210
10
x = x 0 + d1 = 6 + (1) = 5
1
y1 = y 0 + d 2 = 2 + 2 = 0
1 = 0 + d3 = 0 + ( 210 ) 0
35
This puts our pendulum center of mass position (x,y), at (5,0) at 0 degrees orientation
(horizontal). And, since the half-length of the pendulum (l) is one, this is the correct
location given our revolute joint is at (4,0).
d L
dt q
L T
+ qT = Q
q
(68)
q is the column matrix of n generalized coordinates of the rigid bodies which describe the
configuration of the system at any given instant in time. L defines the Lagrangian, which
is the difference between the kinetic energy (T) of the mechanical system and the
potential energy (V).
The first expression represents the inertial forces; the second expression represents the
potential forces; the third expression represents the constraint (i.e., joints and motions)
forces, and Q represents the externally applied forces.
The general form of the this equation appears a bit difficult to understand as a whole, but
examining the components of the equations piece by piece will help in understanding
how it generates the equations of motion. The application of the mechanical system
consisting of a single moving part will illustrate the general case.
36
x
q = y
More generally, q contains all of the n coordinates of the bodies that make up a
mechanical system. For a two-dimensional system, n = 3N, where N = number of rigid
bodies.
The Euler-Lagrange equations above require the Lagrangian to generate the inertial and
potential forces:
L=TV
For our two-dimensional pendulum, the kinetic and potential energy are, respecitively:
T = (mx 2 + my 2 + I 2 )
V = mgy
Lets take the Euler-Lagrange equations one term at a time, starting with the left-most
term for the inertial forces.
First, notice that
L
is the momentum, and is an n x1 array:
q
L
x
mx
L L
=
= my
q y
I
L
The inertia forces on the body in the respective directions (q) are then:
37
L
x
mx mx
T
d L d L d
=
my = my
=
dt q dt y dt
I I
L
L
, which is an n x 1 array and indicates the sensitivity
q
of the Lagrangian for the mechanical system to one of the coordinates. These are the
potential forces; note that ADAMS/Solver includes only the potential force due to gravity
in this term all other potential forces are actually included in Q.
L
x
0
L L
=
= mg
q y
0
L
x 0
m 0 0
0 m 0
T
y + mg + q = Q
0 0 I 0
This is a system of second order differential equations. Note that this system will be
converted to a set of first order differential equations later by introducing velocity
variables.
The forces in the translational direction are determined by the first two equations. In the
third equation, the inertial load due to the angular acceleration is balanced by the net
torque about the center of mass of the pendulum. Since there are no applied forces (Q) in
this model, the last part of the Euler-Lagrange equations to determine is the constraint
forces.
38
R
R/C
C
Figure 5
R = C + R/C
or, equivalently:
39
C + R/C R = 0
Where:
The constant R is a vector which locates the pin (revolute joint) in the ground
reference frame
The variable C is a vector which locates the center of mass of the pendulum in the
ground reference frame
The variable R/C is the vector from the center-of-mass of the pendulum to the pin
(revolute joint) attaching it to ground.
The condition imposed by the mechanism can be rewritten in vector form as the
constraint equation:
( x, y, ) = ( x + R / C1 R1 )i + ( y + R / C2 R2 ) j = 0
R / C1 = l cos
R / C2 = l sin
where i and j are unit vectors along the x and y coordinate axes in the global (ground)
reference frame. We can, alternatively, express this vector constraint in matrix notation:
x + R / C1 R1
( x, y , ) =
=0
y + R / C2 R2
The Jacobian matrix for is the 2 by 3 array:
1 0 l sin
q =
=0
0 1 l cos
Associated with each constraint is a Lagrange multiplier that provides the constraint
(reaction) forces. Using the constraint equations and the Lagrange multipliers, final term
in the Euler-Lagrange equations can be formed:
= 0
l sin
T
q
2
1 =
40
The first two entries are constraint forces, and the third is the torque.
In summary, the power of the formulation process that uses the Euler-Lagrange equation
is illustrated by the conversion of the constraint equations into the coupling terms. The
dynamics of the mechanism are governed by:
Combining the force balance and constraint equations gives us our mechanical pendulum:
x 0
1
m 0 0
0 m 0
=0
2
y + mg +
0 0 I 0 1l sin 2l cos
Or, equivalently, by carrying out the matrix algebra:
mx +1
my + 2 + mg
I + 1l sin 2l cos = 0
x l cos R1
y l sin R2
This set of DAEs for our initial value problem is well-defined, with five equations and
five unknowns.
41
These are kinematic differential equations, also known as the velocity equations.
Substituting the velocity equations into our system of DAEs for the pendulum yields
mU x +1
mU y + 2 + mg
IU + 1l sin 2l cos
U x x
=0
U y y
x l cos R1
y l sin R2
Thus, the second-order system has been reduced to a first order system with the
unknowns:
U x
U
y
U
x
y
1
2
This system is well-posed at this point, but ADAMS/Solver adds set of equations to make
the solution more computationally simple and to increase the scarcity of the equations.
The angular momentum of the system is calculated explicitly and added to the system of
equations, while the translational momentum is not calculated explicitly since it is a
relatively easy quantity to compute. Thus, for our pendulum, the angular momentum, p,
will be introduced to the system of equations as:
L
=0
q
p I = p IU = 0
p
42
mU x +1
mU y + 2 + mg
p + 1l sin 2l cos
p IU
U x x
=0
U y y
x l cos R1
y l sin R2
with these unknowns:
U x
U y
U
p
x
y
1
2
Note the scarcity in the relationships between the variables in the system of DAEs. It is
of considerable importance in the algorithms to invert the matrices for the solution
process as well as in forming the analytical derivatives in the C++ Solver (as opposed to
the FORTRAN solver which numerically perturbs the system to get the derivatives).
To express this system of DAEs in a more compact form, we can rewrite the equations
as:
G (Y , Y , t ) = 0
..where G represents the column matrix of function of the unknowns and where
43
U x
U y
U
p
Y = x
y
1
2
is the column matrix of the unknowns. Again, the differential equations in the system
are first order, and the system as a whole is nonlinear.
Notice that the system is implicit as opposed to an explicit form where
Y = f (Y , t )
and a set of ODEs could be solved. For the implicit integration algorithm, the DAEs
are transformed into a system of non-linear algebraic equations by approximating each
derivative in the column matrixY with a backward differentiation formula (BDF) with
GSTIFF. For numerically stiff systems (characterized by over-damped high-frequency
eigenvalues and under-damped low-frequency eigenvalues) the BDF methods provide
the largest stability region among the multi-step integration methods.
The Jacobian Matrix
One of the most expensive and difficult steps of solving the dynamics problem is in
generating and inverting the Jacobian matrix. In this section, we will form the Jacobian
matrix to give you a sense of both the structure and individual entries. In later sections,
this can help you understand how the model is solved and how problems can occur in the
simulation.
We need to convert our DAEs into an algebraic set of equations so that we can use
Newton-Raphson to solve the problem. To accomplish this, the backward differentiation
formula(s) can be used. The first order BDF is:
q
q qn 1
h
which defines the numerical relationship between a generalized coordinate and its first
derivative (h is the time step) for the first order backward differentiation formula (BDF).
For the higher order approximation, the equation above becomes
44
q qn 1
h
m
h
1
m
h
1
h
I
1
1
h
1
h
1
h
l sin
l cos
l sin
l cos
where the entries that are not filled in are identically zero (here, = 1 for Backward
Euler).
Following the process outlined in Section 5, the Jacobian is evaluated and inverted
repeatedly (as specified by PATTERN) to solve the equations of motion. Notice that
several terms have the step size, h, in the denominator. When the step size is small, these
terms become very large, and this produces an ill-conditioned Jacobian. Many options
exist for dealing with this scenario and are described in more detail in Section 8.
45
x R1 = C sin t
..where the line of action is along the x-axis, and R is the location of the revolute joint
for the pendulum in the ground reference frame.
With the addition of this constraint equation (motion) to our system, the new set of
equations:
mU x +1
mU y + 2 + mg
p + 1l sin 2l cos
p IU
U x x
=0
U y y
x l cos R1
y l sin R2
x R1 = C sin t
46
The addition of another constraint equation to the DAEs has brought a fundamental
change to the nature of the problem. Again, notice that the number of generalized
coordinated equals the number of independent constraint equations, and thus the
dynamics problem is now a kinematic problem characterized by zero degrees of freedom.
Since the motion of the system is completely described throughout all time, adding forces
and torques to the system will not affect the motion of this system. Adding forces and
torques will change the required forces for the constraint equations to impose this motion,
however.
The Jacobian Matrix
The motion of a kinematics problem is completely defined by the constraint equations in
the system. Although it is still possible to form the Jacobian matrix and integrate the full
system of equations, it is no longer necessary. We only need to solve the three by three
system of non-linear algebraic constraint equations, that is, the final three equations in the
system above.
Beginning with an initial estimate of:
x
y
the Newton-Raphson algorithm will be used to solve for the position of the mechanism
at any given time, t. Note that the initial estimate for the first time step is the design
position, and the second time step is the solution from the first time step. However, for
every time step afterward, a linear prediction based on the last two solutions is used to
come up with the first estimate for Newton-Raphson.
The Jacobian matrix for our kinematic system is:
1 0 l sin
0 1 l cos
1 0
0
The kinematic solution generates values only for the displacement variables. To compute
values for the remaining elements of the state vector Y, the constraint equations are
differentiated once and evaluated with the computed values of the variables to obtain the
velocities. Next, the equations are differentiated and evaluated a second time to obtain
the accelerations. Finally, the force balance equations are used to compute values for the
constraint forces.
47
mU x +1
mU y + 2 + mg
p + 1l sin 2l cos
p IU
Ux x
=0
U y y
x l cos R1
y l sin R2
Remove all derivatives leaves:
mU x +1
mU y + 2 + mg
p + 1l sin 2l cos
p IU
U x x
=0
U y y
x l cos R1
y l sin R2
48
1
m
h
1
h
I
1
1
h
1
h
1
h
l sin
l cos
l sin
l cos
Keep in mind that there as static uses Newton-Raphson to find the equilibrium solution, it
makes no consideration to find a minimum energy configuration. Thus, it is possible to
find an equilibrium position where the pendulum is straight down (stable) or straight up
(unstable). Use of ADAMS/Linear to look at the eigenvalues of the system would reveal
positive eigenvalues when the pendulum is straight up, and negative eigenvalues when
down. One option is to use the EQUILIBRIUM/DYNAMIC, which allows an exit
criteria for equilibrium based on kinetic energy of the system.
49
8.
50
g)
h)
i)
j)
51
52
having a very large step-size in the integrator. Especially for challenging models
that fail to find equilibrium, setting a larger value for STABILITY effectively
makes the static algorithm be less aggressive. Recall the funnel approach
discussed above (see c) above). In the beginning take large STABILITY values,
and when you get closer to the solution be more aggressive and take smaller
values. A large STABILITY value is in the range 0.1-5. A small value is 1.E-5
(which is the default).
Note: when experimenting with the STABILITY setting always turn on
DEBUG/EPRINT to be able to judge how this attribute is influencing the
convergence of the static analysis for your model.
g) TLIMIT qualitatively plays the same role that ALIMIT does, except that it
monitors the corrections in the positions rather than orientations. It has been
noted that in general changing the default setting of this attribute does not have a
great impact on the robustness of a static analysis, or at least not as much as
changing the ALIMIT attribute.
53
In the end, if you can run a simulation at large step-sizes and high order it is guaranteed
that your simulation will finish quickly, and the quality of the results is going to be very
good.
Finally, please note that even when a function looks smooth, what is important is
how many derivatives of the function continue to be smooth. For instance, if a cubic
spline is used to represent a motion then the function will have two continuous
derivatives. The third derivative will be already discontinuous, and therefore GSTIFF
will be prevented to increase the integration order to high values. An easy fix to this
problem is to actually specify the motion as a velocity rather than a position level motion.
This simple change will potentially increase the GSTIFF integration order by one, a very
positive change.
2. Choose the modeling units wisely
A poor choice of modeling units leads in different equations to entries that are
vastly different in magnitude. The end result is that during the solution sequence the
Solver will encounter matrices that have very large condition numbers. If the situation is
not extremely bad, a matrix with large condition number will lead to very poor quality
corrections in the Newton-Raphson algorithm for example (see section 1.1) and therefore
large number of iterations for convergence. If the condition number is very large, the
corrector stage of the solution can fail altogether, or the matrix is identified singular by
the linear sparse solver.
To prevent situations like these:
a. Choose units so that model states (displacements and velocities) assume
reasonable values. For example, choosing "mm" for displacements of a rocket,
which travels thousands of kilometers, is a poor choice.
b. Choose time units appropriate to the phenomena being studied. If duration of
dynamic event time is short, consider using MILLISECONDS units.
3. If its possible, avoid dummy parts (any part with zero or very small mass)
Keep in mind these recommendations:
a) Sometimes dummy parts are useful, but if possible, avoid using them.
b) Avoid connecting dummy parts with compliant connections (BEAM, BUSHING,
etc.). If the mass is very small, qualitatively, acceleration = Force/mass = very
large number. Thus, under certain circumstances small masses/moments of
inertia introduce high frequencies and/or jolts into the system, which has
detrimental effects on the solver performance.
c) If you must use dummy parts, then constrain all its degrees of freedom, since with
no degrees of freedom for the dummy part, qualitatively, acceleration=Force/mass
is not an issue.
d) Dummy parts should be massless, 0.0 (or unassigned), and not 1e-20.
54
looking at pivots in a constraint jacobian, which are in no particular order. As a result the
physical meaning may be disregarded. In general, model that has redundant constraints is
an example of poor modeling practices, and it's an indication that the nature of the system
being modeled has not been understood. Redundant constraints can cause eigenvectors
from ADAMS/Linear to produce unwanted results, as the DOF kept may not be the ones
you are interested in for the linear analysis.
5. MOTIONS
Ideally, motions should be defined through functions that are as smooth as possible; i.e.,
that have a large number of continuous derivatives. A simulation in which all the inputs
are smooth will run with high integration orders.
a) Remember the advice number 1 above: keep the expression of the MOTION as
smooth as possible, and avoid discontinuities
b) Avoid using SPLINE functions in a MOTION definition (they only have two
continuous derivatives).
c) If you have to use a spline function in a motion though, you are better off defining
the spline at the velocity rather than position level. For more information see the
MOTION statement in ADAMS/Solver User Manual (see also Knowledge Base
Article #9752)
d) Avoid defining a MOTION as a function of variables (i.e., states).
e) A cubic spline (CUBSPL) may in general work better on motions than the Akima
spline. The derivatives of the Akima are not as nice as those of the cubic spline
(they are more useful in forces rather than motions, see Knowledge Base Article
#7534).
f) Consider filtering data to smooth out spline data for motions.
6. FORCES
a) Remember advice number 1 above: keep the expression of the force as smooth as
possible, and avoid discontinuities.
b) If using data, approximate forces with smooth, continuous spline functions.
c) Make sure velocities are correct in force expressions. For example, in this
damping function: c VX (12, 25, 25) , the fourth marker is missing. This marker
defines the reference frame in which the time derivatives are taken, and this may
be important.
7. CONTACTS
a) Contacts should penetrate before statics. Models with impacts should have slight
penetration in model position when doing statics.
55
b) All tires should penetrate the road. Models with tires should have slight
penetration in model position when doing statics. For example, if only rear tires
penetrate, the static position could be a "handstand".
c) Contact properties are model dependent. See the CONTACT statement
documentation and Knowledge Base Article #10170 for a starting point.
Adjustment of the properties to match experimental results is expected.
d) For contact models, set the INTEGRATOR attribute HMAX=0.01 (a good initial
guess). If problems are encountered, some times a smaller 0.001 values is better.
Note that many times problems come from the use of thin shells. If the input
geometry is very thin, there is a possibility that one geometry may completely
pass through another, resulting in invalid volume of intersection calculations. This
can result in missed contacts, pass-through, or generation of unusually high
contact force generation
e) It's not a bad idea to set INTEGRATOR/HINIT also, if there is a possibility of
sudden movement at the beginning of the simulation. Set HINIT to a small value,
e.g. 1.E-7, to give the integrator a chance to ease into the simulation.
f) Don't turn contact friction on (if possible) until everything else is working
g) The contact penetration depth parameter (CONTACT/DMAX) can force the
integrator to take extremely small step sizes if it is too low (it should not usualy
be less than 0.1 mm or 0.004 in).
8. SUBROUTINES
a) Always use an ADAMS function over a subroutine, if possible. The quality of the
partial derivatives is incomparably better when a function for a force for instance,
is defined via a statement rather than a user subroutine. Consequently, the quality
of the Jacobian in a Newton-Raphson algorithm is going to be better, which
typically results in better convergence properties
b) If using ADAMS/SolverC++, you can provide the partial derivatives from within
your subroutine in order to improve the overall quality of the Jacobian. The
function call SYSPAR is the means through which the user can feed the partials
back into the ADAMS/SolverC++. As indicated earlier, better partials result in
better Newton-Raphson convergence.
c) If you receive errors in your model, for the purpose of debugging eliminate usersubroutines so as they are not the source of error.
d) Make sure that your compiler is compatible with the current version of
MSC.ADAMS.
SIMULATION/INTEGRATORS
a) When applicable, perform an initial static analysis first. Note that a static solution
may be more difficult to find than a dynamic solution. If you care only about the
dynamic solution and cannot find static equilibrium, then either increase your
equilibrium error tolerance, or forget about the static simulation.
56
b) If GSTIFF won't start, it is most likely a problem with initial conditions. Try to
set HINIT to values that are conservative, and thus the integrator will get a chance
to slowly ease into the simulation.
c) Read the documentation about the new corrector implemented in ADAMS and
use it when having problems with discontinuities. Understand the benefits and
drawbacks associated with CORRECTOR=NEW seeting.
d) Don't let the integrator step over important events. Short duration events, such as
an impulse, can be captured by setting the maximum step-size HMAX to values
that are less than the width of the impulse.
e) Use a conservative value on HMAX combined with a lax ERROR setting if you
want to cause the integrator to acts like a fixed-step integrator. There is a
consensus that preventing frequent step-size changes improves the quality of the
solution.
f) Spikes in results output may arise from changes in step size. Reduce HMAX or
try setting HINIT=HMAX. Run with SI2 instead, or use the
INTERPOLATE=ON feature.
g) ADAMS/Solver uses a body 313 rotation sequence (psi, theta, phi). When the
second rotation angle (theta), becomes zero, the model is in a singular
configuration (there is an infinite number of configurations for the psi and phi
angles which lead to the same overall orientation; likewise, the time derivative of
the orientation angles run into problems in a singular configuration). When the
system finds itself into a configuration very close to a singular one, internally the
solver changes the orientation of the body principal axes. This operation always
requires an integrator restart. As a rule of thumb, if the z-axis of the part principal
inertia axis is parallel to z-axis of ground there will be an Euler singularity. If you
have prior knowledge about how the orientation of the bodies in the model is
going to change time, maybe you can define the axes of the IM marker such that
Euler singularities are avoided. Euler singularities lead to simulation slow-down,
and spikes in results (since the integrator is restarted at a small step-size and at
order one).
h) In the ADAMS/SolverC++, if your model already runs fine and you are after
gaining some extra speed-ups, try to set the Jacobian evaluation pattern to F, like
PATTERN=F. This will effectively change the strategy of evaluating the
Jacobian. Recall that if you dont specify anything the Jacobian evaluation
PATTERN is T:F:F:F:T:F:F:F:T:F. However, if you set PATTERN=F, the
integrator will take over the task of deciding how often the Jacobian needs to be
evaluated (adaptive Jacobian mode). It was noticed that for some models this
strategy leads to speed-ups of up to 15%. Note that (i) if you want to actually
never have the Jacobian re-evaluated unless a convergence failure occurs (helpful
for linear models), you can specify PATTERN=F:F; (ii) the adaptive Jacobian
mode works only for GSTIFF (both I3 and SI2) in the ADAMS/SolverC++.
i) In the ADAMS/SolverC++, if your model has high frequencies or discontinuities,
you might want to use the new HHT integrator. In our experiences quite often the
HHT integrator produced sizeable speed-ups in comparison with GSTIFF I3 and
SI2.
57
SIMULATION/STATICS
Remember that the static analysis if a very challenging analysis, particularly when
the initial configuration of the system if far from equilibrium. The equilibrium
configuration is typically obtained as a solution of a non-linear system, and therefore if
the initial guess is far from the solution the Newton-Raphson algorithm diverges (see
section 1.1)
a) Allow MAXIT to assume large values (200 for example). This is not going to
affect the speed of convergence, but it will give the Solver quite a few chances to
find an equilibrium configuration
b) Unless you start close to the equilibrium configuration use a STABILITY value
that is large (0.1 to 2.0, rather than the 1.E-5, which is the default). A large value
of the STABILITY parameter is an indication that you are not in a hurry to reach
the equilibrium, but you are content with the Solver taking small steps towards the
equilibrium configuration. Large values of the STABILITY attribute indicate that
you allow the Solver to be very aggressive in making corrections in the model
configuration. It is a good idea to experiment with the value of this settings, as in
our experience we found it to have quite an effect in the success/failure of the
static analysis.
c) Use the attribute ALIMIT to limit the value of the orientation update that the
Solver is allowed to take. The default is 30D, but for problematic models more
conservative values such as ALIMIT=10D can actually help with the overall
convergence. For more challenging models you might want to also change the
default value of the TLIMIT attribute, although in our experience this had a
smaller impact.
FRICTION
a) ADAMS has its own friction model implemented. However, if you are not
satisfied with the existing friction model, new friction models (such as LuGre,
elasto-plastic, etc) can be implemented using ADAMS function and usersubroutine support.
b) Note that unlike the ADAMS/Solver C++, the ADAMS/SolverF77 does not allow
friction in joints connecting flexible bodies.
DEBUGGING
a) While trying to debug or validate your simulations, always turn on
DEBUG/EPRINT. You should only turn off the DEBUG/EPRINT command
when you are comfortable with the results and satisfied with the CPU time
required to complete a simulation.
b) Understand the output sent out to the user, and see where the Solver has a hard
time advancing the simulation. Monitor and correct any warnings sent out by the
solver
58
c) Try to understand the system that you are simulating from a physical standpoint.
d) Use "building blocks" of concepts that worked in the past. Add enhancements to
the model using a crawl-walk-run approach.
e) Test with small models to isolate problems.
f) Have graphics for visualizing motion.
g) Look at damping terms as a source of errors. Incorrect sign, missing terms are
typical mistakes.
h) Models should have no warnings during simulation (e.g., redundant constraints,
spline functions, etc.)
i) Understand numerical methods, and particularly understand your integrator.
j) Look for results which become very large in magnitude; this could indicate a
discontinuity
MISCELLANEOUS
a) Avoid very large numbers and very small numbers. Be wary when your model
contains numbers like 1e+23 or 1e-20. This is most likely a symptom of a bad
units choice (see the UNITS section above), and it typically results in matrices
with large condition numbers.
b) Extend the range of spline data beyond the range of need.
c) Don't write function expressions that can potentially divide by zero. The quick
and dirty way to do this is to use the MAX function, as in FUNCTION =
8/MAX(0.01,your_function). However, there are more intelligent ways of doing
this, since the function MAX is discontinuous, which leads to problems of a
different nature (remember the advice number 1 above)
d) Add moderate damping in your model so frequencies can die off. This will speed
up your simulation.
e) Avoid very high damping rates. High damping is causing a rapid decay in
response, which is difficult for an integrator to follow.
f) Avoid toggles, dual solutions, or bifurcations since they lead to discontinuities.
g) Don't use 1.0 for the exponent in IMPACT or BISTOP functions. This creates a
"corner" (i.e. a non-smooth function). Try 2.2 for the exponent instead (note that
this parameter depends on units).
h) The ADAMS/SolverC++ has support for 2D parts. A 2D part adds fewer
equations to the overall problem and especially for such elements as chains and
belts this can results in more robust and faster simulations.
59
of thumb for creating robust ADAMS models (see section 8.3). With these tasks out of
the way, set DTOUT (END/STEPS in the output statement/command) to suit your timesampling simulation requirements. For a correlation study, i.e., correlation with test data,
this value is typically driven by the sampling rate of the test data. For a given DTOUT (or
END/STEPS combo) take the following annotated steps:
1. Start with an initial error tolerance value for the INTEGRATOR statement. The
starting point could be the Solver default value of 1E-3.
2. Build/Run simulation with DEBUG/EPRINT turned on, with the initial estimate
of error value.
3. Reduce error by a factor of your choosing - a factor of 10 is commonly used.
4. Build/Run simulation with DEBUG/EPRINT turned on, with the reduced error
value.
5. Compare the results between the simulations with original and reduced ERROR
values. Check results, particularly displacements first if you are using I3
integrator(s). You will be looking at one of the following two scenarios.
a. If results don't match, reduce error tolerance further from the current value
by a factor of your choosing, and start again from Step 2.
b. If results match closely, compare (gdiff) message files between the two
simulations. Ensure that Solver is "doing more work" (taking more steps)
because of the reduced error tolerance value.
If you are not seeing more integration steps being taken in the simulation
with the reduced error tolerance value, then most likely you need to
decouple the effects of step-size and integration error control. In this
situation, either your output step size DTOUT, or HMAX setting is
smaller than the integration step-size that would otherwise be required by
your ERROR setting, and is thus controlling the accuracy of your results.
To decouple DTOUT (or HMAX) and ERROR, make DTOUT (and/or
HMAX), much bigger (fewer steps, bigger steps) and use ERROR to do
the study as before until you return to case a) above and find that the
reduced error tolerance value did cause solver to do more work.
6. At this point, reset your error tolerance to the larger of the ERROR values from
the two simulations that proved that you had a converged solution (the
simulations that gave the same results). At this point you can be comfortable that
you have achieved a converged solution for displacements.
NOTE: Dealing with corrector/integrator problems (will most likely happen), is
handled in the ADAMS/Solver documentation, under the INTEGRATOR
statement/command section. More information is available in the MSC.ADAMS
Knowledge Base at http://support.adams.com/kb, or in section 8.1 of this
document.
7. Follow steps 1 through 5 if you want to converge the acceleration results (when
using an I3 integrator the accelerations are known to sometimes display spikes,
which typically are more difficult to eliminate). In our experience the
acceleration spikes can be eliminated or reduced by slowly decreasing the value
60
8.5.1. DEBUG/EPRINT
This statement/command is used to obtain a synopsis of the processes that take
place during a solution sequence. A solution sequence should be regarded as any action
or algorithm invoked by the Solver to advance the simulation or to execute a user issued
command. For example, issuing a simulate/statics command will cause the Solver to a)
set up the set of equations that needs to be solved to find an equilibrium configuration;
and b) employ and algorithm to find the solution of the assembled equation set.
DEBUG/EPRINT will send out informative messages and warnings during both stages a)
61
and b). Acting on the warnings sent out by the solver is going to make for better models
and more robust/fast simulations. As an example of a warning message, you might get
---- START: WARNING ---An Out-of-range X value has been used in the spline calculation.
Extrapolation is required for curve 1 of SPLINE/123.
The X value (-100.000000) should be between -50.000000 and 55.000000
---- END: WARNING ----
In this case, the user defined the spline curve anticipating that the free parameter will
always assume values anywhere between 50 and 55. Nevertheless, the model
configuration is such that the spline is being evaluated for a value of 100. This is an
indication that either the spline range initial assumption was wrong and the bounds
should be altered, or that somewhere else there could possibly be problem in the model
that led to unreasonable values for the spline parameter.
There are some differences in the amount and type of information that is put out
by the ADAMS/SolverC++ and ADAMS/SolverF77. Either one though sends out
information summarizing the progress towards finding a solution. The algorithms
employed are specific to each analysis type (position IC, statics, kinematics, dynamics,
etc.), but the common denominator is that the solution is always computed as a set of
successive iterations that each corrects the previous system configuration to bring it
closer to the actual solution. This iterative process can and should be monitored for each
model that is problematic. Below is shown an actual portion of the output that you get by
having DEBUG/EPRINT and running a dynamics analysis with the GSTIFF integrator.
Integration step No. 7,
Order = 2
Cumulative iterations of
From 5.000000000E-04, Step = 8.33333E-05
the corrector: 315
Iteration
Residual (or equation error)
Change in the Variable
Jacobian
count
Maximum Element/ID Equation
Maximum Element/ID Variable
new?
___
________ ______________________ ________ ______________________ ___
1
-4.7E+05 Part/10 Phi Torque
6.2E+05 RevJnt/13 LAM
yes
2
2.5E+00 Part/11 Psi Torque
5.0E+02 RtnMot/2 LAM
no
3
-2.0E-03 Part/11 Psi Torque
2.6E-01 RevJnt/15 LAM
no
4
1.3E-04 Part/10 Theta Torque
1.4E-04 RevJnt/13 LAM
no
62
-8.6E-06
2.0E-03
TrnJnt/2 LAM
no
VMNPO:SING_MTX
-----------------------------------------------------------------------The system matrix has become numerically singular with this pivot order!
-----------------------------------------------------------------------Equation for the singular row ...... MOTION 270_1032.MOT10320301
...... (row=MOTION/10320301 CONSTR.)
Variable for the singular column ... PART 270_1032.PAR10323
...... (col=PART/10323 Theta)
Pivot magnitude .................... 4.05476E-17
------------------------------------------------------------------------
SOLVE:REGEN
Symbolic refactorization attempted at time = 1.7274.
When solving for the Newton-Raphson corrections (see section 1.1) the Jacobian needs to
be factored to determine its LU factorization. The pivots in the LU factorization; the
diagonal entries in U, are chosen at the beginning of the simulation and they are re-used
time and again upon a new call for the solution of this system. Especially for long
simulations after a while the original pivot sequence results in a singular matrix. This
requires a refactorization to compute a new pivot sequence in the LU factorization, and
possibly a decrease of the integration step-size. The user is flagged when the solver runs
into such a scenario, as indicated in the ADAMS/SolverF77 message above.
DEBUG/EPRINT is an extremely useful tool in diagnosing simulation problems.
By inspecting eprint output for problematic models, you can find the root cause of your
problem and thus be able to fix the problem yourself (see the discussion in sections 8.1
through 8.3), or call the ADAMS hotline and provide the support engineer with
additional information that can help him/her diagnose the issue. To conclude, always
start the process of tinkering with a model simulation by having the DEBUG/EPRINT on.
You are encouraged to have eprint on and experiment with different settings of the
INTEGRATOR, EQUILIBRIUM, etc., statements/commands. This is the way to gain
understanding of what parameters influence what aspect of the simulation, and eventually
become a power ADAMS/Solver user.
8.5.2. DEBUG/RHSDUMP
The DEBUG/RHSDUMP statement/command is intended for the power user who has a
basic understanding about the concept of ADAMS state, residual, Newton-Raphson
correction. This is the type of information that DEBUG/RHSDUMP will output and it
will help the user to quickly
a) Inspect all the entries in the array of force-balance residuals, constraint violation,
etc. to check that the values are as expected. For instance, a modeling error in an
obscure element might be easily noticed in a large value in a very large initial
residual in satisfying the force-balance equation. Note that each equation that is
formulated in ADAMS is given a name, and this name along with the violation
(residual) associated with that equation is sent out to the user
64
b) Inspect all the values of the state to find any state that has an abnormally large (or
small value). ADAMS/Solver has names associated with every variable, such that
the user can monitor name-values state information.
c) Inspect the values of the corrections computed at the end of each Newton
iteration. This information is of limited use, but if an extremely large value is
noticed, or the opposite when a zero value is obtained where a non-zero
correction is expected is an indication that either (i) some residual is not correctly
computed (see a) above), or (ii) some Jacobian entries are not computed correctly.
Note that the amount of information that is sent out to the user in DEBUG/RHSDUMP
mode is large. Keep this in mind when running long simulations with large models.
Only turn this feature on at the beginning of the simulation for a short while, or right
before the model runs into a problem that you try to diagnose.
8.5.3. DEBUG/JMDUMP
The information obtained with the DEBUG/JMDUMP command/statement is
primarily used for debugging of the code by ADAMS developers. It is basically a dump
of the Jacobian matrix assembled by the Solver for the purpose of carrying out NewtonRaphson iterations. The information sent out indicates each partial derivative that is
computed in the process of assembling the Jacobian; that is, a value is provided, along
with information regarding what equation this partial derivative is of, and with respect to
what variable the partial is taken.
The information output with DEBUG/JMDUMP is typically useful for small
models, and even then analyzing the output obtained with DEBUG/EPRINT and
DEBUG/RHSDUMP is more useful and should be considered first. Note that for large
models the amount of information generated can become very quickly overwhelming.
65
1
h M
y = I
Fu
Tu
T ( P )T f
p
p
T ( P )T f
p
pT
(P )
B T JB
T
( K )
T ( R )T n
p
T ( R )T n ( K ) T
0
F
T
0
I Ff
Tf
1
T
I ( K )
h
0
1
I
h
0
0
0
0
F
T
p
Fp
Tp
1
I
h
F
T
66
T
(R )
Fn
I Tn
0