Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Modeling Mechanical DAE Using Natural Coordinates

Modeling Mechanical DAE Using Natural Coordinates Mathematical and Computer Modelling of Dynamical Systems 1387-3954/01/0702-145$16.00 2001, Vol. 7, No. 2, pp. 145±158 Swets & Zeitlinger Modeling Mechanical DAE Using Natural Coordinates 1 1,2 1 C. KRAUS , M.WINCKLER and H.G. BOCK ABSTRACT An ef®cient modeling technique for multibody systems, which extends the concept of natural coordinates with closed kinematic loops is presented. By establishing a local coordinate system in each body a system with constant mass matrix is set up. The propagation of topological information into the model leads to the application of a block-oriented rational Cholesky decomposition of the system matrix. The overall algorithm shows linear complexity in the number of bodies for systems with a constant number of kinematic loops. To handle rank- de®cient constraint Jacobians arising from loop closing conditions the concept of constraint partitioning during decomposition is outlined. Restriction of the partitioning decision based on the topology information minimizes the monitoring effort and avoids disadvantages experi- enced in other partitioning methods. Numerical results for the 6-bar-mechanism proof the algorithm to run ef®ciently with projection stabilized index-1 integration methods. Keywords: constraint partitioning, dynamics, modeling, multibody, natural coordinates, nu- merical simulation; Models and Numerical Methods (65C20). 1 INTRODUCTION When modeling the dynamics of multibody systems (MBS) one of the key aspects is the choice of coordinates. Several approaches to this issue exist and the question of a suitable choice depends on the observed class of mechanisms and is an area of ongoing research. For tree-structured mechanisms without closed kinematic loops for example it is well known that a very ef®cient formulation for the equations of motion can be obtained using minimal Interdisciplinary Center for Scienti®c Computing (IWR), University of Heidelberg, Germany. Im Neuenheimer Feld 368, D-69120 Heidelberg, Germany. E-mail: Michel.Winckler@ IWR.Uni-Heidelberg.de 146 C. KRAUS ET AL. coordinates. Recursively traversing through the kinematic chains and using sets of variables that parameterize the degrees of freedom allowed by the joints between successive bodies introduces a set of minimal coordinates in a very intuitive way. This modeling technique leads to a set of ordinary differential equations (ODEs) for the equations of motion (see Featherstone [1] for example). With the advent of computer based modeling tools and due to the increased complexity of mechanical devices under surveillance, redundant coordinate systems for MBS on the basis of Lagrangian equations of the ®rst kind become an invaluable alternative to systems of minimal coordinates. The main reasons for using redundant coordinates are their ability to describe kinematically non- singular mechanical devices with mathematically non-singular sets of equa- tions and their capability to handle closed kinematic loops in a structured and well-de®ned way. When treating a general assembly of bodies in a MBS, a unifying idea for most modeling techniques is to establish a local coordinate system in each of these bodies. The coordinate choice is usually the choice of an ef®cient para- meterization of the af®ne transformation between the global (world) and each of the local (body) coordinate systems. Reference point coordinates e.g. model the Cartesian coordinates of the origin of the local coordinate system together with three angles (Euler angles). While this leads to a very compact set of coordinates, the drawback of singular coordinate con®gurations that are not introduced by a singularity in the mechanical device itself is the cause of a variety of numerical problems. Overcoming this drawback is only possible with redundant coordinate systems, since any parameterization of a free body in 3D with exactly 6 coordinates leads to at least two singular coordinate con®gurations. The equations of motion derived from Lagrangian equations of the ®rst kind lead to the well known set of differential algebraic equations of index 3 for MBS. Index reduction for numerical treatment of this DAE with index-1 semi-explicite integration methods with projection stabilization (such as MBSSIM by v. Schwerin, Winckler [6]) introduces the full descriptor form: p _ ˆ v v _ ˆ a MG a f ˆ …1† G 0 COORDINATES IN MBS - DAE 147 g … p†ˆ 0 g … p; v†ˆ 0 with positions p, velocity v ˆ dp=dt, acceleration a ˆ dv=dt, mass matrix M, forces f , kinematic constraints g , velocity constraints g , acceleration p v independent part of the acceleration constraints , Lagrange-multiplier  and constraint Jacobian G ˆ dg =dp. The goal of our modeling approach is to derive all system components of the full descriptor form (1) such as the mass matrix M and the constraint Jacobian G in an automatic fashion. For a fast and ef®cient simulation a low complexity count for the evaluation of these system components is vital. Contrary to this approach traditional modeling techniques try to reduce the number of dynamical variables instead, as especially the treatment of the matrix form components seems to depend quadratically or even cubicly on the number of variables. Our modeling method of choice uses an augmented and improved version   of the natural coordinates approach described by Garcõa de Jalon and Bayo in [2]. Although this method uses twelve coordinates per body, the exploitation of the resulting matrix structures and the low complexity count for the model component evaluation results in a very competitive formulation of the overall system. This leads directly to an algorithm of linear complexity that is fast and ef®cient if all structure information is treated properly by the simulation environment. 2 THE GENERAL AFFINE TRANSFORMATION The idea of a modeling technique using local, body ®xed coordinate systems and a global (reference) coordinate system is carried out in many different coordinate choices. A common subproblem when designing such a modeling technique is to parameterize the underlying af®ne transformation from the local to the global coordinates. Usually, this transformation is restricted to the special case of the combination of a translational part of the origin and a rotational part. An af®ne transformation of a body-®xed point r expressed in local coordinates to his global representation  r is given by r ˆ o‡ Xr …2† 148 C. KRAUS ET AL. 3 33 where o R is the origin of the local coordinate system and X R is the corresponding orientation-deformation matrix. This matrix does not have to be an element ofSO…3†. A more general view that allows scaling and shearing can be used to adjust the coordinate system to the local topology of the model. Our method does not parameterize X, but uses all nine entries of X as dynamic variables. Because of this relation between the dynamic variables and the matrix X every local point or vector can be transformed to its global representation by a linear equation of the dynamic variables. r ˆ o‡ Xr ˆ o‡…e je je †r 1 2 3 ˆ 1o‡ r e ‡ r e ‡ r e 1 1 2 2 3 3 0 1 0 1 1 r r r 1 2 3 B C B C B 1 C ˆ 1 r r r …3† @ A B C 1 2 3 @ A 1 r r r 1 2 3 |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚} |‚‚‚{z‚‚‚} ˆ C p 312 12 where C R is de®ned by the local af®ne coef®cients and p  R are the dynamic variables of each body. Since the transformation is linear, it holds for the ®rst and second derivative with the same matrix C . Remark The modeling process for natural coordinates not only uses af®ne points and their transformation as shown above, but also pure translation vectors. Replacing the leading 3 3 identity matrix block of C with a 3 3 zero block gives the correct transformation for such translation vectors while preserving the general matrix block structure of the transformation. 3 DERIVING THE MASS MATRIX The ®rst system component that we will investigate is the mass matrix M. Following the principles of theoretic mechanics, we can derive an expression for M from the general form of the kinetic energy of a system with generalized coordinates q: 1 1 T T T ˆ q _ Mq _ ˆ  q _ qd _ V: …4† 2 2 V COORDINATES IN MBS - DAE 149 Making use of the ®rst derivative of transformation (3) we can write the kinetic energy with respect to the dynamic variables p T T T ˆ  p _ C C pd _ V …5† While in local coordinates the mass matrix is dependent on the dynamical variables, the transformation to global coordinates removes this dependency allowing us to move p _ , p _ out of the integrand. This gives an analytical form for the mass matrix: 1 1 T T T T ˆ p _  C C dV p _ ˆ: p _ Mp _ …6† 2 2 We get a bi-linear form with respect to the velocity for the kinetic energy, and by this we have established the mass matrix for natural coordinates. The mass matrix is constant (as in many formalisms with a body-®xed coordinate system) and sparse, a fact that follows directly from the structure of the transformation matrix C . It can be computed directly or by using the center of mass and the matrix of inertia computed with respect to Euler angles. Remark This formulation does not introduce additional terms that corre- spond to coriolis or centrifugal forces because of the direct relation between position and velocity variables. External forces like gravitation and springs- damper elements are modeled directly as force terms in the righthand side. For a more detailed introduction see Garcõ Âa de Jalo  n and Bayo[2]. 4 MODELING RIGID BODIES As we have already mentioned, in natural coordinates the orientation- deformation matrix X must not be an element of SO…3†. However, when modeling rigid bodies a set of algebraic relations between the natural coordinates is imposed. These so-called rigid body constraints restrict the coordinate manifold to the suitable 6-dimensional subspace by ®xing lengths and angles between the local coordinate axis. This leads to a set of six simple constraint equations: All these equations are polynomial of second order, so they and their derivatives, which are used in the full descriptor form, are very cheap to 150 C. KRAUS ET AL. Fig. 1. The rigid body constraints. compute. The resulting constraint matrix is linear with respect to the dynamic variables and has a sparse vector structure. 0 1 0 e 00 B T C 00 e 0 B C B C RC @g B C 00 0 e RC B C G … p†ˆ … p†ˆ …8† B T T C @p 0 e e 0 B 2 1 C B C T T @ 0 e 0 e A 3 1 T T 00 e e 3 2 In many modeling techniques the computation of ± the acceleration independent residual of the acceleration constraint ± can be rather time consuming, due to the fact that the second derivative of a nonlinear kinematic constraint function in several variables usually is an expression with a high complexity count. In numerical experiments with different modeling approaches the computation of consumed up to 80% of the total time to evaluate the model. Comparing this to natural coordinates we ®nd that the computation of is about as cheap as the evaluation of the kinematic constraints. RC T ˆÿe_ e_ …iˆf1; 2; 3g† i i RC T ˆÿ2e_ e_ …iˆf4; 5; 6g† …9† j k RC The structure of G can be further improved by multiplication with the 4-block4-block block-diagonal 12 12 matrix 00 1 1 ÿ1 BB C C D… p†ˆ diag @@ e A A: …10† RC from the right. This reduces G to a matrix with only nine entries, all of them equal to one. Choosing the underlying local coordinate system to be orthonormal the matrix inverse of this coordinate transformation matrix is COORDINATES IN MBS - DAE 151 very simple to compute and can be stored for premultiplication of according blocks. 5 DEFINING JOINTS BY CONSTRAINTS SETS Introducing joints between free bodies can be understood as a further restric- tion of the manifold of consistent coordinate sets. So we impose additional joint constraints to model joints. The two key concepts for modeling standard joints are: restrict translational movement in the direction given by a direction vector and restrict rotational movement with respect to a rotational axis both of which can be derived using the concept of parallelism.   In contrast to the original approach by Garcõa de Jalon and Bayo that uses the vector product to model parallelism and always leads to a set of three linearly dependent equations, we are looking for joint constraints that achieve the restriction of the feasible manifold and have a Jacobian of full row rank. In the numerical treatment of the constraint matrix this result will be advantageous. So for both key concepts, translational and rotational restriction, we formulate the fact that akb using a basis for the orthogonal subspace: akb () 9 h ; h : h ? h ; a ? h ; a ? h ; b ? h ; b ? h : …11† 1 2 1 2 1 2 1 2 The vectors h and h are chosen in local coordinates of one of the two 1 2 bodies that the joint links together such that they satisfy their relations towards a and towards each other. Their relations to b are formulated as two constraints with a rank-2 Jacobian. A revolute joint (relative movement is restricted to an axis rotation) is realized in two steps. First, we ®x a common point with local representations b ; b . This special form of disallowing translational movement is modeled 1 2 by identifying b and b in a set of three equations of ®rst order. Using 1 2 transformation (3) to translate this relation to a relation between dynamical variables we get g … p ; p † :ˆ b … p †ÿ b … p †ˆ C p ÿ C p ˆ 0: …12† 1 2 1 1 2 2 b1 1 b2 2 p;1... 3 152 C. KRAUS ET AL. Fig. 2. Assembly of revolute joint. Furthermore, we have to forbid rotations other than about the rotation axis (local coordinates d ; d ). Using two vectors h ; h as in (11) and their 1 2 1 2 transformation to global coordinates C ; C we get h1 h2 T T g … p ; p † : ˆ d … p † h … p †ˆ…C p † C p ˆ 0 …13† 1 2 1 1 1 2 d1 1 h1 2 p;4 T T g … p ; p † : ˆ d … p † h … p †ˆ…C p † C p ˆ 0 …14† 1 2 1 1 2 2 d1 1 h2 2 p;5 where h and h are chosen perpendicular to d . 1 2 2 In the same fashion, cylindrical and prismatic joints can be formulated. In either case we restrict translational movement in speci®c directions. This is accomplished by de®ning a translational vector through the difference vector of two body-®xed points. With this idea we are able to reuse the method introduced for restricting rotation. More complex joints arising from industrial applications can be modeled by combining the formulation of these basic ones and extending the concept of parallelism. Remark Up to now we have not used the fact that the choice of a local coordinate system in each body is free. Using this additional modeling freedom we can adapt the local coordinate system in each body to one of the joints that connect this body to the overall mechanism. This technique leads to an even simpler inner structure of the transformation matrix C (see (3)), where some entries are zero. RC Remark While transformation (10) drastically simpli®es G , it also applies to the joint constraint matrix G and the mass matrix. To minimize computational effort, we formulate G in the transformed formulation and compute it directly instead of applying the transformation in a separate step. COORDINATES IN MBS - DAE 153 This operation can be implemented quite ef®ciently in computer code. Using it when formulating the theory distracts the view from the layout of the constraint equations. 6 EXPLOITING THE BLOCK-SPARSE STRUCTURE LEADS TO AN O…N† ALGORITHM The local topology of bodies and joints has lead to a block-sparse structure of the linear system (1). Exploiting this structure when solving the linear system during each integration step leads to a very ef®cient numerical treatment of the equation set. As quite a lot of equations in the set of joint constraints amount to identifying two points on different bodies with one another, it is quite tempting to reduce the set of dynamical variables by eliminating these point sharing conditions during the modeling step. In contrast to this our approach keeps the full structured matrix including the point sharing conditions. The reason for this is that the blocks representing the unconnected bodies before assembly (through joints) are kept untouched in our matrix system. They can be treated very ef®ciently, even preselecting the inner pattern of entries for faster decomposition. For systems with the topology of a linear chain there exists a block±ordering that results in a block-banded matrix. The augmented system of the mass matrix and the rigid-body constraints corresponding to each body are placed on the main block diagonal. The constraint matrix blocks corresponding to joints are placed on the ®rst sub-diagonals. This matrix can be decomposed in linear time using a block-oriented rational Cholesky algorithm. For systems with tree-structure (and a suitable ordering of the bodies) each branch from the tree root introduces a pair of off-diagonal blocks, while the chain itself is represented by a set of blocks in the ®rst sub-diagonals as in the example of the linear chain. However, the off-diagonal blocks do not introduce additional ®ll-in during the Cholesky decomposition, so this still results in an O…n†-Algorithm. Joints leading to closed kinematic loops introduce an additional pair of off- diagonal blocks. Due to the loop structure this blocks produces ®ll-in along the loop chain, leading to an additional effort of up to two blocks propagated along the set of bodies in the loop. Applying the Cholesky 154 C. KRAUS ET AL. decomposition properly the complexity is O…n† when comparing systems with an equal number of closed kinematic loops. This result is especially useful for the question whether a kinematic chain can be modeled in more detail. The matrix blocks that form the main diagonal correspond to the mass matrix and the rigid body constraints of each body. RC MG B ˆ …15† RC During decomposition of the block structured system, block-operations of the ÿ1 T type B J , where J are blocks corresponding to joint constraints, result in most of the computational effort. By using a null-space method and exploiting RC the very simple structure of the transformed and constant matrix…† G D ,we are able to further reduce the overall computation time by a factor of two, although additional transformations of the joint constraint matrix and the mass matrix are needed. 7 MODELING THE 6-BAR-MECHANISM The structure of the system matrix is illustrated by a model of the 6-bar- mechanism (described in [5]). The model consists of 6 bodies (one of them ®xed and not modeled) connected by six revolute joints in a special manner (see Fig. 4). The resulting system has one closed kinematic loop. It is well known that this model has one degree of freedom, so there exists one dependency in the set of constraint equations. In the structural diagram (see Fig. 3) numbers represent the augmented systems of each body consisting of the mass matrix and the rigid body constraint matrix. The joint constraint matrix and its transpose are represented by ÿ and j respectively. ‡ and  are blocks introduced by the decomposition algorithm, where  is additional ®ll-in propagated by the closed kinematic loop. Remark In the original formulation of Garcõ Âa de Jalo  n and Bayo the point sharing concept is used to shrink the system. A velocity projection algorithm is applied to this model, leading to quadratic complexity. That approach is superior only for a very small number of bodies. COORDINATES IN MBS - DAE 155 Fig. 3. Structural diagram. Fig. 4. 6-bar-mechanism. 8 TREATMENT OF REDUNDANT CONSTRAINTS Although the formulation of each joint avoids the problem of redundant constraints, closing kinematic loops can bring it in through the back door as we have seen in the 6-bar-mechanism example. So an automatic modeling of systems with closed kinematic loops could result in a system with rank- de®cient Jacobian. Direct treatment of this phenomenon can be accomplished by a rank analysis of the constraint matrix. Finding a set of constraints with maximum rank leads to the idea of constraint partitioning, a technique that is closely related to the coordinate partitioning approach, where a set of independent dynamical variables is selected. The active set of independent constraints in the constraint partitioning has to be kept under close surveillance by a monitor on the matrix condition of the active Jacobian. 156 C. KRAUS ET AL. To minimize the effort of such an approach it is convenient to keep the rigid body constraints in the active set all the time. This is justi®ed by their necessity in the non-redundant open loop case. Furthermore we can restrict ourselves to monitor each diagonal joint block only, applying the same argument. This reduces the extra work for the constraint partitioning to a neglectable additional effort. Again this result is achieved due to the fact that in the non-reduced system the blocks correspond to topological equivalents. In practice we use Cholesky decomposition for the joint-block. If a (near) zero diagonal element appears during decomposition, the corresponding constraint is dynamically eliminated by setting the corresponding solution vector element to zero during solution. This is completely equivalent to eliminating the constraint altogether, while still providing information about redundancy of the corresponding constraint. It is worth mentioning, that a restart of the method can be avoided, because the discontinuous Lagrange multipliers do not go into the differential part of the DAE. This drastically improves the performance of this method in comparison to coordinate partitioning. However, the method is only valid if one stays near enough to the consistent manifold. This is guaranteed by the invariant-projection methods supplied by MBSSIM. 9 CONCLUSION Natural coordinates are a powerful tool for modeling mechanical systems. The representation of the local topology in the system matrix and a simple and ¯exible underlying transformation over-compensate for the apparent drawback of a redundant coordinate approach. Investigations by v. Schwerin about computation time show that natural coordinates can compete with other modeling techniques [6]. The O…n†-algorithm outlined in this article and a sparse substructure make our approach especially competitive for systems with a large number of bodies. As a runtime example we have chosen a more handy model, namely the already mentioned 6-bar-mechanism. As a standard benchmark model with no global minimal coordinate formulation it is well suited to illustrate the treatment of a set of redundant constraints. In Figure 5 we compare runtime results for different projection accuracies. The Figure represents a series of numerical experiments where we have COORDINATES IN MBS - DAE 157 Fig. 5. Runtime (in secs) vs integration tolerance (± log(err)). plotted the computation time (in CPU seconds on y-axis) for a simulation of 10 s over the required absolute integration tolerance (ÿlog…err† on x±axis). The platform we used was a Pentium II 400 Mhz running Windows NT. Each graph represents the runtime results for a speci®c tolerance of the norm of the constraint residual. Obviously that with a more strict projection tolerance the runtime is reduced. This is due to the fact that the projection algorithm in MBSSIM is very cheap compared to the effort to perform additional steps due to a change in the active set of the constraint partitioning. Another example with many kinematic loops (hexapod) was published by the authors in [7]. The model under observation has thirteen bodies with six closed kinematic loops. It was shown that the assembly is as easy as in the open loop case. The hexapod can be simulated in real-time on a Pentium III 800 Mhz. The software package MBSNAT (Multi Body System simulation using NATural coordinates) by Kraus supplies an object-oriented implementation of model generator and linear algebra outlined in this paper, for use in any suited MBS integration package (e.g., MBSSIM by v. Schwerin and Winckler). 158 C. KRAUS ET AL. ACKNOWLEDGMENTS This paper extends the ideas outlined in a proceedings article of the Mathmod 2000 in Vienna [4]. This work is funded by WIR-Forschungsverbund Parallele Simulationswerkzeuge (principal investigators H.G. Bock and H. Ruder). This work is funded by the German Department of Science and Technology within the research grant Optimierung von Motorkomponen- ten (principal investigators H.G. Bock and J. Schlo È der). REFERENCES 1. Featherstone, R.: Robot Dynamics Algorithms. Kluwer Academic Publisher, Boston / Dordrecht / Lancaster, 1987.   2. Garcõa de Jalon, J. and Bayo, E.: Kinematic and Dynamic Simulation of Multibody Systems. Springer, 1994. 3. Kraus, C.: Modellierung und Rekursive Algorithmen fu È r Mehrko È rpersysteme in Natu È rlichen Koordinaten. IWR, Universita Èt Heidelberg, 1997. 4. Kraus, C. and Winckler, M.: Natural Coordinates and Mechanical DAE. In: Proceedings of 3rd MATHMOD VIENNA, 2000, pp. 201±204. 5. Lester, G. and Schiehlen, W.: Benchmark-Beispiele des DFG-Schwerpunktprogrammes Dynamik von Mehrko È rpersystemen. Zwischenbericht ZB-64. Institut B fu È r Mechanik, Universitat Stuttgart, 1991. 6. v. Schwerin, R.: Numerical Methods, Algorithms and Software for Higher Index Nonlinear Differential-Algebraic Equations in Multibody System SIMulation, PhD thesis. Universita Èt Heidelberg, 1997. 7. Winckler, M. and Kraus, C.: Simulation of Hexapod Machine Tools by Using Natural Coordinates. In: Proceedings of Year 2000 PKM International Conference, 2000, pp. 109±117. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical and Computer Modelling of Dynamical Systems Taylor & Francis

Modeling Mechanical DAE Using Natural Coordinates

Download PDF
14 pages

Loading next page...
 
/lp/taylor-francis/modeling-mechanical-dae-using-natural-coordinates-y0UuRh8HrZ

References (7)

Publisher
Taylor & Francis
Copyright
Copyright Taylor & Francis Group, LLC
ISSN
1744-5051
eISSN
1387-3954
DOI
10.1076/mcmd.7.2.145.3645
Publisher site
See Article on Publisher Site

Abstract

Mathematical and Computer Modelling of Dynamical Systems 1387-3954/01/0702-145$16.00 2001, Vol. 7, No. 2, pp. 145±158 Swets & Zeitlinger Modeling Mechanical DAE Using Natural Coordinates 1 1,2 1 C. KRAUS , M.WINCKLER and H.G. BOCK ABSTRACT An ef®cient modeling technique for multibody systems, which extends the concept of natural coordinates with closed kinematic loops is presented. By establishing a local coordinate system in each body a system with constant mass matrix is set up. The propagation of topological information into the model leads to the application of a block-oriented rational Cholesky decomposition of the system matrix. The overall algorithm shows linear complexity in the number of bodies for systems with a constant number of kinematic loops. To handle rank- de®cient constraint Jacobians arising from loop closing conditions the concept of constraint partitioning during decomposition is outlined. Restriction of the partitioning decision based on the topology information minimizes the monitoring effort and avoids disadvantages experi- enced in other partitioning methods. Numerical results for the 6-bar-mechanism proof the algorithm to run ef®ciently with projection stabilized index-1 integration methods. Keywords: constraint partitioning, dynamics, modeling, multibody, natural coordinates, nu- merical simulation; Models and Numerical Methods (65C20). 1 INTRODUCTION When modeling the dynamics of multibody systems (MBS) one of the key aspects is the choice of coordinates. Several approaches to this issue exist and the question of a suitable choice depends on the observed class of mechanisms and is an area of ongoing research. For tree-structured mechanisms without closed kinematic loops for example it is well known that a very ef®cient formulation for the equations of motion can be obtained using minimal Interdisciplinary Center for Scienti®c Computing (IWR), University of Heidelberg, Germany. Im Neuenheimer Feld 368, D-69120 Heidelberg, Germany. E-mail: Michel.Winckler@ IWR.Uni-Heidelberg.de 146 C. KRAUS ET AL. coordinates. Recursively traversing through the kinematic chains and using sets of variables that parameterize the degrees of freedom allowed by the joints between successive bodies introduces a set of minimal coordinates in a very intuitive way. This modeling technique leads to a set of ordinary differential equations (ODEs) for the equations of motion (see Featherstone [1] for example). With the advent of computer based modeling tools and due to the increased complexity of mechanical devices under surveillance, redundant coordinate systems for MBS on the basis of Lagrangian equations of the ®rst kind become an invaluable alternative to systems of minimal coordinates. The main reasons for using redundant coordinates are their ability to describe kinematically non- singular mechanical devices with mathematically non-singular sets of equa- tions and their capability to handle closed kinematic loops in a structured and well-de®ned way. When treating a general assembly of bodies in a MBS, a unifying idea for most modeling techniques is to establish a local coordinate system in each of these bodies. The coordinate choice is usually the choice of an ef®cient para- meterization of the af®ne transformation between the global (world) and each of the local (body) coordinate systems. Reference point coordinates e.g. model the Cartesian coordinates of the origin of the local coordinate system together with three angles (Euler angles). While this leads to a very compact set of coordinates, the drawback of singular coordinate con®gurations that are not introduced by a singularity in the mechanical device itself is the cause of a variety of numerical problems. Overcoming this drawback is only possible with redundant coordinate systems, since any parameterization of a free body in 3D with exactly 6 coordinates leads to at least two singular coordinate con®gurations. The equations of motion derived from Lagrangian equations of the ®rst kind lead to the well known set of differential algebraic equations of index 3 for MBS. Index reduction for numerical treatment of this DAE with index-1 semi-explicite integration methods with projection stabilization (such as MBSSIM by v. Schwerin, Winckler [6]) introduces the full descriptor form: p _ ˆ v v _ ˆ a MG a f ˆ …1† G 0 COORDINATES IN MBS - DAE 147 g … p†ˆ 0 g … p; v†ˆ 0 with positions p, velocity v ˆ dp=dt, acceleration a ˆ dv=dt, mass matrix M, forces f , kinematic constraints g , velocity constraints g , acceleration p v independent part of the acceleration constraints , Lagrange-multiplier  and constraint Jacobian G ˆ dg =dp. The goal of our modeling approach is to derive all system components of the full descriptor form (1) such as the mass matrix M and the constraint Jacobian G in an automatic fashion. For a fast and ef®cient simulation a low complexity count for the evaluation of these system components is vital. Contrary to this approach traditional modeling techniques try to reduce the number of dynamical variables instead, as especially the treatment of the matrix form components seems to depend quadratically or even cubicly on the number of variables. Our modeling method of choice uses an augmented and improved version   of the natural coordinates approach described by Garcõa de Jalon and Bayo in [2]. Although this method uses twelve coordinates per body, the exploitation of the resulting matrix structures and the low complexity count for the model component evaluation results in a very competitive formulation of the overall system. This leads directly to an algorithm of linear complexity that is fast and ef®cient if all structure information is treated properly by the simulation environment. 2 THE GENERAL AFFINE TRANSFORMATION The idea of a modeling technique using local, body ®xed coordinate systems and a global (reference) coordinate system is carried out in many different coordinate choices. A common subproblem when designing such a modeling technique is to parameterize the underlying af®ne transformation from the local to the global coordinates. Usually, this transformation is restricted to the special case of the combination of a translational part of the origin and a rotational part. An af®ne transformation of a body-®xed point r expressed in local coordinates to his global representation  r is given by r ˆ o‡ Xr …2† 148 C. KRAUS ET AL. 3 33 where o R is the origin of the local coordinate system and X R is the corresponding orientation-deformation matrix. This matrix does not have to be an element ofSO…3†. A more general view that allows scaling and shearing can be used to adjust the coordinate system to the local topology of the model. Our method does not parameterize X, but uses all nine entries of X as dynamic variables. Because of this relation between the dynamic variables and the matrix X every local point or vector can be transformed to its global representation by a linear equation of the dynamic variables. r ˆ o‡ Xr ˆ o‡…e je je †r 1 2 3 ˆ 1o‡ r e ‡ r e ‡ r e 1 1 2 2 3 3 0 1 0 1 1 r r r 1 2 3 B C B C B 1 C ˆ 1 r r r …3† @ A B C 1 2 3 @ A 1 r r r 1 2 3 |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚} |‚‚‚{z‚‚‚} ˆ C p 312 12 where C R is de®ned by the local af®ne coef®cients and p  R are the dynamic variables of each body. Since the transformation is linear, it holds for the ®rst and second derivative with the same matrix C . Remark The modeling process for natural coordinates not only uses af®ne points and their transformation as shown above, but also pure translation vectors. Replacing the leading 3 3 identity matrix block of C with a 3 3 zero block gives the correct transformation for such translation vectors while preserving the general matrix block structure of the transformation. 3 DERIVING THE MASS MATRIX The ®rst system component that we will investigate is the mass matrix M. Following the principles of theoretic mechanics, we can derive an expression for M from the general form of the kinetic energy of a system with generalized coordinates q: 1 1 T T T ˆ q _ Mq _ ˆ  q _ qd _ V: …4† 2 2 V COORDINATES IN MBS - DAE 149 Making use of the ®rst derivative of transformation (3) we can write the kinetic energy with respect to the dynamic variables p T T T ˆ  p _ C C pd _ V …5† While in local coordinates the mass matrix is dependent on the dynamical variables, the transformation to global coordinates removes this dependency allowing us to move p _ , p _ out of the integrand. This gives an analytical form for the mass matrix: 1 1 T T T T ˆ p _  C C dV p _ ˆ: p _ Mp _ …6† 2 2 We get a bi-linear form with respect to the velocity for the kinetic energy, and by this we have established the mass matrix for natural coordinates. The mass matrix is constant (as in many formalisms with a body-®xed coordinate system) and sparse, a fact that follows directly from the structure of the transformation matrix C . It can be computed directly or by using the center of mass and the matrix of inertia computed with respect to Euler angles. Remark This formulation does not introduce additional terms that corre- spond to coriolis or centrifugal forces because of the direct relation between position and velocity variables. External forces like gravitation and springs- damper elements are modeled directly as force terms in the righthand side. For a more detailed introduction see Garcõ Âa de Jalo  n and Bayo[2]. 4 MODELING RIGID BODIES As we have already mentioned, in natural coordinates the orientation- deformation matrix X must not be an element of SO…3†. However, when modeling rigid bodies a set of algebraic relations between the natural coordinates is imposed. These so-called rigid body constraints restrict the coordinate manifold to the suitable 6-dimensional subspace by ®xing lengths and angles between the local coordinate axis. This leads to a set of six simple constraint equations: All these equations are polynomial of second order, so they and their derivatives, which are used in the full descriptor form, are very cheap to 150 C. KRAUS ET AL. Fig. 1. The rigid body constraints. compute. The resulting constraint matrix is linear with respect to the dynamic variables and has a sparse vector structure. 0 1 0 e 00 B T C 00 e 0 B C B C RC @g B C 00 0 e RC B C G … p†ˆ … p†ˆ …8† B T T C @p 0 e e 0 B 2 1 C B C T T @ 0 e 0 e A 3 1 T T 00 e e 3 2 In many modeling techniques the computation of ± the acceleration independent residual of the acceleration constraint ± can be rather time consuming, due to the fact that the second derivative of a nonlinear kinematic constraint function in several variables usually is an expression with a high complexity count. In numerical experiments with different modeling approaches the computation of consumed up to 80% of the total time to evaluate the model. Comparing this to natural coordinates we ®nd that the computation of is about as cheap as the evaluation of the kinematic constraints. RC T ˆÿe_ e_ …iˆf1; 2; 3g† i i RC T ˆÿ2e_ e_ …iˆf4; 5; 6g† …9† j k RC The structure of G can be further improved by multiplication with the 4-block4-block block-diagonal 12 12 matrix 00 1 1 ÿ1 BB C C D… p†ˆ diag @@ e A A: …10† RC from the right. This reduces G to a matrix with only nine entries, all of them equal to one. Choosing the underlying local coordinate system to be orthonormal the matrix inverse of this coordinate transformation matrix is COORDINATES IN MBS - DAE 151 very simple to compute and can be stored for premultiplication of according blocks. 5 DEFINING JOINTS BY CONSTRAINTS SETS Introducing joints between free bodies can be understood as a further restric- tion of the manifold of consistent coordinate sets. So we impose additional joint constraints to model joints. The two key concepts for modeling standard joints are: restrict translational movement in the direction given by a direction vector and restrict rotational movement with respect to a rotational axis both of which can be derived using the concept of parallelism.   In contrast to the original approach by Garcõa de Jalon and Bayo that uses the vector product to model parallelism and always leads to a set of three linearly dependent equations, we are looking for joint constraints that achieve the restriction of the feasible manifold and have a Jacobian of full row rank. In the numerical treatment of the constraint matrix this result will be advantageous. So for both key concepts, translational and rotational restriction, we formulate the fact that akb using a basis for the orthogonal subspace: akb () 9 h ; h : h ? h ; a ? h ; a ? h ; b ? h ; b ? h : …11† 1 2 1 2 1 2 1 2 The vectors h and h are chosen in local coordinates of one of the two 1 2 bodies that the joint links together such that they satisfy their relations towards a and towards each other. Their relations to b are formulated as two constraints with a rank-2 Jacobian. A revolute joint (relative movement is restricted to an axis rotation) is realized in two steps. First, we ®x a common point with local representations b ; b . This special form of disallowing translational movement is modeled 1 2 by identifying b and b in a set of three equations of ®rst order. Using 1 2 transformation (3) to translate this relation to a relation between dynamical variables we get g … p ; p † :ˆ b … p †ÿ b … p †ˆ C p ÿ C p ˆ 0: …12† 1 2 1 1 2 2 b1 1 b2 2 p;1... 3 152 C. KRAUS ET AL. Fig. 2. Assembly of revolute joint. Furthermore, we have to forbid rotations other than about the rotation axis (local coordinates d ; d ). Using two vectors h ; h as in (11) and their 1 2 1 2 transformation to global coordinates C ; C we get h1 h2 T T g … p ; p † : ˆ d … p † h … p †ˆ…C p † C p ˆ 0 …13† 1 2 1 1 1 2 d1 1 h1 2 p;4 T T g … p ; p † : ˆ d … p † h … p †ˆ…C p † C p ˆ 0 …14† 1 2 1 1 2 2 d1 1 h2 2 p;5 where h and h are chosen perpendicular to d . 1 2 2 In the same fashion, cylindrical and prismatic joints can be formulated. In either case we restrict translational movement in speci®c directions. This is accomplished by de®ning a translational vector through the difference vector of two body-®xed points. With this idea we are able to reuse the method introduced for restricting rotation. More complex joints arising from industrial applications can be modeled by combining the formulation of these basic ones and extending the concept of parallelism. Remark Up to now we have not used the fact that the choice of a local coordinate system in each body is free. Using this additional modeling freedom we can adapt the local coordinate system in each body to one of the joints that connect this body to the overall mechanism. This technique leads to an even simpler inner structure of the transformation matrix C (see (3)), where some entries are zero. RC Remark While transformation (10) drastically simpli®es G , it also applies to the joint constraint matrix G and the mass matrix. To minimize computational effort, we formulate G in the transformed formulation and compute it directly instead of applying the transformation in a separate step. COORDINATES IN MBS - DAE 153 This operation can be implemented quite ef®ciently in computer code. Using it when formulating the theory distracts the view from the layout of the constraint equations. 6 EXPLOITING THE BLOCK-SPARSE STRUCTURE LEADS TO AN O…N† ALGORITHM The local topology of bodies and joints has lead to a block-sparse structure of the linear system (1). Exploiting this structure when solving the linear system during each integration step leads to a very ef®cient numerical treatment of the equation set. As quite a lot of equations in the set of joint constraints amount to identifying two points on different bodies with one another, it is quite tempting to reduce the set of dynamical variables by eliminating these point sharing conditions during the modeling step. In contrast to this our approach keeps the full structured matrix including the point sharing conditions. The reason for this is that the blocks representing the unconnected bodies before assembly (through joints) are kept untouched in our matrix system. They can be treated very ef®ciently, even preselecting the inner pattern of entries for faster decomposition. For systems with the topology of a linear chain there exists a block±ordering that results in a block-banded matrix. The augmented system of the mass matrix and the rigid-body constraints corresponding to each body are placed on the main block diagonal. The constraint matrix blocks corresponding to joints are placed on the ®rst sub-diagonals. This matrix can be decomposed in linear time using a block-oriented rational Cholesky algorithm. For systems with tree-structure (and a suitable ordering of the bodies) each branch from the tree root introduces a pair of off-diagonal blocks, while the chain itself is represented by a set of blocks in the ®rst sub-diagonals as in the example of the linear chain. However, the off-diagonal blocks do not introduce additional ®ll-in during the Cholesky decomposition, so this still results in an O…n†-Algorithm. Joints leading to closed kinematic loops introduce an additional pair of off- diagonal blocks. Due to the loop structure this blocks produces ®ll-in along the loop chain, leading to an additional effort of up to two blocks propagated along the set of bodies in the loop. Applying the Cholesky 154 C. KRAUS ET AL. decomposition properly the complexity is O…n† when comparing systems with an equal number of closed kinematic loops. This result is especially useful for the question whether a kinematic chain can be modeled in more detail. The matrix blocks that form the main diagonal correspond to the mass matrix and the rigid body constraints of each body. RC MG B ˆ …15† RC During decomposition of the block structured system, block-operations of the ÿ1 T type B J , where J are blocks corresponding to joint constraints, result in most of the computational effort. By using a null-space method and exploiting RC the very simple structure of the transformed and constant matrix…† G D ,we are able to further reduce the overall computation time by a factor of two, although additional transformations of the joint constraint matrix and the mass matrix are needed. 7 MODELING THE 6-BAR-MECHANISM The structure of the system matrix is illustrated by a model of the 6-bar- mechanism (described in [5]). The model consists of 6 bodies (one of them ®xed and not modeled) connected by six revolute joints in a special manner (see Fig. 4). The resulting system has one closed kinematic loop. It is well known that this model has one degree of freedom, so there exists one dependency in the set of constraint equations. In the structural diagram (see Fig. 3) numbers represent the augmented systems of each body consisting of the mass matrix and the rigid body constraint matrix. The joint constraint matrix and its transpose are represented by ÿ and j respectively. ‡ and  are blocks introduced by the decomposition algorithm, where  is additional ®ll-in propagated by the closed kinematic loop. Remark In the original formulation of Garcõ Âa de Jalo  n and Bayo the point sharing concept is used to shrink the system. A velocity projection algorithm is applied to this model, leading to quadratic complexity. That approach is superior only for a very small number of bodies. COORDINATES IN MBS - DAE 155 Fig. 3. Structural diagram. Fig. 4. 6-bar-mechanism. 8 TREATMENT OF REDUNDANT CONSTRAINTS Although the formulation of each joint avoids the problem of redundant constraints, closing kinematic loops can bring it in through the back door as we have seen in the 6-bar-mechanism example. So an automatic modeling of systems with closed kinematic loops could result in a system with rank- de®cient Jacobian. Direct treatment of this phenomenon can be accomplished by a rank analysis of the constraint matrix. Finding a set of constraints with maximum rank leads to the idea of constraint partitioning, a technique that is closely related to the coordinate partitioning approach, where a set of independent dynamical variables is selected. The active set of independent constraints in the constraint partitioning has to be kept under close surveillance by a monitor on the matrix condition of the active Jacobian. 156 C. KRAUS ET AL. To minimize the effort of such an approach it is convenient to keep the rigid body constraints in the active set all the time. This is justi®ed by their necessity in the non-redundant open loop case. Furthermore we can restrict ourselves to monitor each diagonal joint block only, applying the same argument. This reduces the extra work for the constraint partitioning to a neglectable additional effort. Again this result is achieved due to the fact that in the non-reduced system the blocks correspond to topological equivalents. In practice we use Cholesky decomposition for the joint-block. If a (near) zero diagonal element appears during decomposition, the corresponding constraint is dynamically eliminated by setting the corresponding solution vector element to zero during solution. This is completely equivalent to eliminating the constraint altogether, while still providing information about redundancy of the corresponding constraint. It is worth mentioning, that a restart of the method can be avoided, because the discontinuous Lagrange multipliers do not go into the differential part of the DAE. This drastically improves the performance of this method in comparison to coordinate partitioning. However, the method is only valid if one stays near enough to the consistent manifold. This is guaranteed by the invariant-projection methods supplied by MBSSIM. 9 CONCLUSION Natural coordinates are a powerful tool for modeling mechanical systems. The representation of the local topology in the system matrix and a simple and ¯exible underlying transformation over-compensate for the apparent drawback of a redundant coordinate approach. Investigations by v. Schwerin about computation time show that natural coordinates can compete with other modeling techniques [6]. The O…n†-algorithm outlined in this article and a sparse substructure make our approach especially competitive for systems with a large number of bodies. As a runtime example we have chosen a more handy model, namely the already mentioned 6-bar-mechanism. As a standard benchmark model with no global minimal coordinate formulation it is well suited to illustrate the treatment of a set of redundant constraints. In Figure 5 we compare runtime results for different projection accuracies. The Figure represents a series of numerical experiments where we have COORDINATES IN MBS - DAE 157 Fig. 5. Runtime (in secs) vs integration tolerance (± log(err)). plotted the computation time (in CPU seconds on y-axis) for a simulation of 10 s over the required absolute integration tolerance (ÿlog…err† on x±axis). The platform we used was a Pentium II 400 Mhz running Windows NT. Each graph represents the runtime results for a speci®c tolerance of the norm of the constraint residual. Obviously that with a more strict projection tolerance the runtime is reduced. This is due to the fact that the projection algorithm in MBSSIM is very cheap compared to the effort to perform additional steps due to a change in the active set of the constraint partitioning. Another example with many kinematic loops (hexapod) was published by the authors in [7]. The model under observation has thirteen bodies with six closed kinematic loops. It was shown that the assembly is as easy as in the open loop case. The hexapod can be simulated in real-time on a Pentium III 800 Mhz. The software package MBSNAT (Multi Body System simulation using NATural coordinates) by Kraus supplies an object-oriented implementation of model generator and linear algebra outlined in this paper, for use in any suited MBS integration package (e.g., MBSSIM by v. Schwerin and Winckler). 158 C. KRAUS ET AL. ACKNOWLEDGMENTS This paper extends the ideas outlined in a proceedings article of the Mathmod 2000 in Vienna [4]. This work is funded by WIR-Forschungsverbund Parallele Simulationswerkzeuge (principal investigators H.G. Bock and H. Ruder). This work is funded by the German Department of Science and Technology within the research grant Optimierung von Motorkomponen- ten (principal investigators H.G. Bock and J. Schlo È der). REFERENCES 1. Featherstone, R.: Robot Dynamics Algorithms. Kluwer Academic Publisher, Boston / Dordrecht / Lancaster, 1987.   2. Garcõa de Jalon, J. and Bayo, E.: Kinematic and Dynamic Simulation of Multibody Systems. Springer, 1994. 3. Kraus, C.: Modellierung und Rekursive Algorithmen fu È r Mehrko È rpersysteme in Natu È rlichen Koordinaten. IWR, Universita Èt Heidelberg, 1997. 4. Kraus, C. and Winckler, M.: Natural Coordinates and Mechanical DAE. In: Proceedings of 3rd MATHMOD VIENNA, 2000, pp. 201±204. 5. Lester, G. and Schiehlen, W.: Benchmark-Beispiele des DFG-Schwerpunktprogrammes Dynamik von Mehrko È rpersystemen. Zwischenbericht ZB-64. Institut B fu È r Mechanik, Universitat Stuttgart, 1991. 6. v. Schwerin, R.: Numerical Methods, Algorithms and Software for Higher Index Nonlinear Differential-Algebraic Equations in Multibody System SIMulation, PhD thesis. Universita Èt Heidelberg, 1997. 7. Winckler, M. and Kraus, C.: Simulation of Hexapod Machine Tools by Using Natural Coordinates. In: Proceedings of Year 2000 PKM International Conference, 2000, pp. 109±117.

Journal

Mathematical and Computer Modelling of Dynamical Systems – Taylor & Francis

Published: Jun 1, 2001

There are no references for this article.