Benjamin Boudon
email: benjamin.boudon@sigma-clermont.fr

Thu Thuy Dang

Rebecca Margetts

François Malburet

SIMULATION METHODS OF RIGID HOLONOMIC MULTI-BODY SYSTEMS WITH BOND GRAPHS

Keywords: bond graph, multibond graphs, multibody system, modelling methodology, simulation, closed kinematic chain

Bond graph software can simulate bond graph models without the user needing to manually derive equations. This offers the power to model larger and more complex systems than in the past. Multibond Graphs (those with vector bonds) offer a compact model which further eases handling multibody systems. Although multibond graphs can be simulated successfully, the use of vector bonds can present difficulties. In addition, most qualitative, bond graph-based exploitation relies on the use of scalar bonds. This paper discusses the main methods for simulating bond graphs of multibody systems, using a graphical software platform. The transformation between models with vector and scalar bonds is presented. The methods are then compared with respect to both time and accuracy, through simulation of two benchmark models. This paper is a tutorial on the existing methods for simulating 3D rigid and holonomic multibody systems using bond graphs, and discusses the difficulties encountered. It then proposes and adapts methods for simulating this type of system directly from its bond graph within a software package. The value of this study is in giving practical guidance to modellers, so that they can implement the adapted method in software.

INTRODUCTION

There has been a resurgence of interest in bond graphs in the last twenty years, as computer science has progressed [START_REF] Borutzky | Bond Graph Methodology -Development and Analysis of Multidisciplinary Dynamic System Models[END_REF][START_REF] Borutzky | Bond Graph Modelling of Engineering Systems -Theory, Applications and Software Support[END_REF] . Software environments with graphical user interfaces enable the user to enter, modify and interpret bond graphs, allowing the graphical aspect to be fully exploited. Equations can be automatically generated, reducing the scope for human error, and highperformance numerical solvers can be utilised to solve them.

As far as multibody systems are concerned, several bond graph methods have been developed. Karnopp and Rosenberg define a procedure for constructing bond graphs from Lagrange equations [START_REF] Karnopp | Power-conserving transformations: physical interpretations and applications using bond graphs[END_REF][START_REF] Rosenberg Rc | Multiport models in mechanics[END_REF][START_REF] Dc | Systems Dynamics : Modeling, Simulation, and Control of Mechatronics Systems[END_REF] . Tiernego and Bos offer a modular approach based on Newton-Euler equations [START_REF] Bos | Multibody Systems in Terms of Multibond Graphs with Application to a Motorcycle Multibody system[END_REF] and Hamiltonian formalism [START_REF] Breedveld | Multibond graph representation of Lagrangian mechanics: the elimination of the Euler Junction Structure[END_REF] . Different specific physical aspects of multibody systems have also been modelled: high pairs of joints [START_REF]A representation for planar point contact joints in the multibody mechanical bond graph library[END_REF] , complex friction models ( [START_REF] Bera | Evaluation of antilock braking system with an integrated model of full vehicle system dynamics[END_REF][START_REF] Nakhaeinejad | Dynamic modeling of rolling element bearings with surface contact defects using bond graphs[END_REF] ), and flexible bodies [START_REF] Damić | Modelling flexible body systems: a bond graph component model approach[END_REF][START_REF] Cohodar | Comparison of different formulations of 2D beam elements based on Bond Graph technique[END_REF] ). In addition, bond graph models have been developed for industrial applications: motorcycle [START_REF] Bos | Multibody Systems in Terms of Multibond Graphs with Application to a Motorcycle Multibody system[END_REF] , vehicle dynamics [START_REF] Pacejka Hb | Modeling Complex Vehicules using Bond Graphs[END_REF][START_REF] Pacejka Hb | Development and validation of a Bond Graph Handling model of an automobile[END_REF] , automotive powertrains [START_REF] Hrovat | Bond Graph Modeling and Computer Simulation of Automative Torque Converters[END_REF][START_REF] Hrovat | Bond graph modeling of automative power trains[END_REF] , robots [START_REF] Pathak Pm | Impendance Control of Space Robots Using Passive Degrees Freedom in Controller Domain[END_REF] , and biomechanics 18 . Some books which have some parts on the modelling of multibody systems [START_REF] Borutzky | Bond Graph Methodology -Development and Analysis of Multidisciplinary Dynamic System Models[END_REF][START_REF]Les bonds graphs[END_REF][START_REF] Das | Mechatronic modeling and simulation using bond graphs[END_REF][START_REF] Merzouki R | Intelligent mechatronic systems -Modeling, control and diagnostic[END_REF] . However, it is relatively unusual to simulate mechanisms directly from the bond graph without manually deriving the associated equations [START_REF] Bos | Multibody Systems in Terms of Multibond Graphs with Application to a Motorcycle Multibody system[END_REF][START_REF] Felez | Simulation of multi-body systems using multi-bond graph[END_REF][START_REF] Rideout | System partitioning and physical domain proper modeling through assessment of powerconserving model structure[END_REF][START_REF] Ersal | A Bond Graph Based Modular Modeling Approach Towards an Automated Modeling Environment for Reconfigurable Machine Tools[END_REF][START_REF] Rahman T | Evaluation of Dynamic Performance of Non-Spherical Parallel Orientation Manipulators through Bond Graph Multi-Body Simulation[END_REF][START_REF] Marquis-Favre | Mechatronic bond graph modelling of an automotive vehicle[END_REF] . This paper therefore focuses on the methods for simulating 3D rigid and holonomic multibody systems from bond graphs. The authors define the class of multibody system studied as one where the joints are lower pairs of joints, only a viscous friction model is considered, and the bodies are assumed to be rigid. The simulation of such systems in bond graphs is problematic, with little detailed literature on the subject.

The aim of this paper is to overview and expand on existing methods, in order to provide guidance on simulating multibody systems directly from a bond graph i.e. automatically within a software package, as opposed to the user manually deriving the dynamic equations. It is addressed to bond graph specialists as well as engineers. The software used here is 20-SIM: a simulation package for dynamic systems using physical components, block diagrams, bond graphs and equations of motion. The methodology recommended to modellers depends on the type of system to be modelled (open chain (OC) or closed kinematic chain (CKC) systems). Several methods for simulating multibody systems with bond graphs will be compared numerically, by observing parameters such as time (computational load) and accuracy.

An important part of this methodology is to show how a bond graph with vector bonds can be transformed into a bond graph with scalar bonds. This is because bond graph exploitation largely relies on the use of scalar bonds: simplification and reduction [START_REF] Ersal | Structural simplification of modular BG models based on junction inactivity[END_REF] , structural analysis 28 , and model inversion [START_REF] Gawthrop | Bicausal bond graphs[END_REF][START_REF] Gawthrop | Physical interpretation of inverse dynamics using bi-causal bond graphs[END_REF][START_REF] Ngwompo R | Bond graph methodology for the design of an actuating system : application to a two-link manipulator[END_REF][START_REF] Ngwompo R | Dimensioning problems in system design using bicausal bond graphs[END_REF][START_REF] Ngwompo R | Physical model-based inversion in control systems design using bond graph representation -application[END_REF][START_REF] Ngwompo R | Physical model-based inversion in control systems design using BG representation -theory[END_REF] .

The method is demonstrated on two classic multibody case studies: a planar pendulum and a slider crank. The bond graph approach used in this paper has specific features which allow a structured and modular development of complex mechatronic systems: it is multiphysics, graphical, object-oriented, and acausal. Embedded electronics could easily be inserted into the models (e.g. modelling a torque delivered by an electrical actuator in the slider crank system). The graphical nature of bond graphs facilitates a global view and comprehension of large complex mechatronic systems, such as helicopter anti-vibratory system [START_REF] Boudon | Simulation of a helicopter's main gearbox semiactive suspension with bond graphs[END_REF] . The orientedobject and acausal features also permit a modular approach, allowing knowledge to be capitalised upon and the modelling task to be automated.

The outline of the paper is as follows: Section 2 gives a brief review of bond graph modelling in the context of multibody systems. Section 3 is dedicated to the methodology used to simulate bond graphs for Multibody systems both with vector and scalar bonds. Section 4 presents the applications of the presented methods and a numerical comparison. Section 5 draws general conclusions. IGi Body i's inertia matrix about mass-centre of body i expressed in its frame[kg.m²]

EJSGi

Eulerian Junction Structure (EJS) matrix about mass-centre of body i expressed in its frame

BOND GRAPH MODELING DEDICATED TO MULTIBODY SYSTEMS

Brief review

This section reviews the main contributions on modelling dynamic behaviour of three-dimensional multibody systems. A more detailed review can be found in [START_REF] Borutzky | Bond Graph Methodology -Development and Analysis of Multidisciplinary Dynamic System Models[END_REF] .

The multibond graph formalism 36 37 is an extension of bond graph method, where the scalar power bonds become vector bonds and the elements multiports. It extended the application of the bond graph to the study of multibody systems with three dimensions.

The bond graph approach for multibody systems was introduced by Bos [START_REF] Mjl | Modelling the dynamics and kinematics of mechanical systems with multibond graphs[END_REF] [START_REF] Bos | Multibody Systems in Terms of Multibond Graphs with Application to a Motorcycle Multibody system[END_REF] . In his PhD thesis, Bos developed bond-graph models for three-dimensional multibody systems and discussed how to derive the equations of motion from the bond-graph in several different forms. He conducted simulations of a 3D motorcycle, although the equations were derived by hand.

Library models for a rigid body and for various types of joints have been provided by Zeid and Chung 39 so that bond graph models of rigid multibody systems can be assembled in a systematic manner.

Felez [START_REF] Felez | BONDYN : a bond graph based simulation program for multibody systems[END_REF] develops a software program that models multibody systems using bond graphs. To manage derivative causalities with this software, he proposes introducing Lagrange multipliers into the system to eliminate derivative causality.

Van Dijk and Breedveld present different methods for simulating bond graph models [START_REF] Van Dijk | Simulation of system models containing zero-order causal paths-II Numerical implications of class 1 zero-order causal paths[END_REF] . Simulations were conducted with a predecessor of the 20-sim software, and numerically compared on the basis of computing time and accuracy. The potential to use multibond graphs was mentioned, but the difficulty of implementing bond graphs with vector bonds was not detailed. This point will be discussed later in this paper.

Marquis-Favre and Scavarda [START_REF]Bond graph representation of multibody systems with kinematic loops[END_REF] propose a method to simplify bond graph models for multibody systems with kinematic loops. Nevertheless, few complex multibody systems with closed kinematic loops have been simulated directly from dedicated software.

More recently, a body of work has conducted simulations of complex kinematic closed systems with multi-bond graphs directly from 20-sim software: e.g. Rideout [START_REF] Rideout | System partitioning and physical domain proper modeling through assessment of powerconserving model structure[END_REF] , Ersal , Rahman 25 and Boudon [START_REF] Boudon | Simulation of a helicopter's main gearbox semiactive suspension with bond graphs[END_REF] .

Approach chosen

The authors selected the Bos and Tiernego method [START_REF] Mjl | Modelling the dynamics and kinematics of mechanical systems with multibond graphs[END_REF] for modelling multibody systems with bond graphs, because it allows a modular approach. This method enables a multibody system to be built as an assembly of bodies and joints, and is based on the use of absolute coordinate systems (Figure 1) and Newton-Euler equations. The dynamic equations of a rigid body therefore depend on its mass/inertia properties and on geometric parameters for the body under consideration. The kinematic joints constrain the effort and flow vectors associated with the articulation points in the assembly of two bodies, so that the desired relative motion can be achieved. Consequently, the dynamic equations of the complete system consist of the dynamic equations of each body, in terms of its own parameters and the constraint equations of each joint.

 

0 0 0 /0 0 j h j pes j j G h d F P V dt      m    (1)     /0 i/0 i/0 0 i i i i i i i i i G G i G h d M dt              I I     (2)
The upper part of the bond graph in Figure 2 represents the rotational dynamic part expressed in the body frame, while the lower part is for the translational dynamic part

O j M i G i z  i y  k M i x  y y          i z z 0             i R x OG y z z     0 z  i y    0 y    0 x x     x   i x  i G I i G EJS /0   i i /0  i i G V /0  k i M V /0  j i M V /0 i i i i k G M     /0 i i i i j G M       0 , ,     i P f i m 0    g i W
expressed in the inertial reference frame (or Galilean frame). Note that the power bonds corresponding to the rotational quantities are marked in purple while the power bonds corresponding to the translational quantities are in green. The two 1-junction arrays correspond to the angular velocity vector of body j with regards to the inertial frame

  i/ 0 i  
k i i j i i M R G R S i k M R G R S i j V V G M V V G M                 (3)
As the translational dynamic is expressed in the inertial reference frame, a modulated transformer 4)) between the body frame Ri and the inertial frame R0.

0 /0 0 /0 . i i i i G G V P V    (4)
In this paper, XYZ Cardan angles have been employed for the sake of simplicity. The rotation matrix 0 i P can be calculated from Cardan angles. The components of angular velocity for each body expressed in the body frame (the pseudo-velocities, or rates of change of Cardan angles) are used to determine the body's orientation and the corresponding coordinate transformation matrix. This classical process is used in the Cardan block detailed in Figure 3. It should be noted that the initial conditions for the integration of time derivatives must be consistent with regards to the kinematic constraints.

It can be noticed that the finite rotation transformation should be also defined with other coordinates systems since this transformation is powerful conservative : for example, Euler angles with angles which are compatible with the mechanism concerning the singularities aspects, the Rodrigues-Hamilton parameters or the Cayley-Klein parameters. The implementation of the finite rotation transformation with the Rodrigues-Hamilton parameters is given in [START_REF]Contribution à la représentation bond graph des systèmes mécaniques multicorps[END_REF] . Modelling Kinematic joints The joint models express the constraints that are introduced when rigid bodies are connected. As with the bond graph model of the rigid body, the joint models have been built in a modular way i.e. their models do not change when the whole model of the system is assembled. The idea of this section is to allow bond graph practitioners to use a library of all the common existing kinematic joint models.

General kinematic joint model

The modelling of kinematic joints determines the rotational or translation degree of freedom allowed by the joint.

Flow sources can be used to suppress the joint's degrees of freedom. However, in order to circumvent the causality constraints mentioned before, the joints models are presented with an additional L element which is either an R/C element or a modulated effort source MSe depending of the choice of simulation method (R/C element methods or the use of Lagrange multipliers).

For the unconstrained degrees of freedom, the choice of modelling assumptions can dictate additional elements. If the joints are assumed to be perfect (i.e. without any dissipation or energy storage), there are no additional elements. However, if dissipation or compliance are assumed in the joints, R or C elements are added at the corresponding 1-junctions. These R/C /0

ix i i iy iz                                       /0 / / 0 cos 0 cos tan sin tan sin cos tan sin tan 1 i i i i i i ix i i i i iy i i i i z i i i                                                          i i i i                     i i i i                 i i i ini du           0 cos( ) cos( )
cos( )sin( ) sin( ) cos( )sin( )sin( ) cos( )sin( ) cos( ) cos( ) sin( )sin( )sin( ) cos( )sin( ) sin( )sin( ) cos( )cos( )sin( ) cos( )sin( ) cos( )sin( )sin( ) c

i i i i i i i i i i i i i i i i i i i i i i i i i i i i P                                   os( ) cos( ) i i             i ini i ini i ini i ini                 0 i P i/0 i  
elements will be called functional elements and should not be confused with the R/C elements used as parasitic elements for purpose simulation in section 3.3.3. The joint model is built from the following kinematic relationships:

1 1 1 1/ 0 1/ 2 2 / 0         (5)
Where

1 1/ 0   is the absolute velocity of body 1 expressed in 1 R frame 1 1/ 2  
is the relative angular velocity of body 1 with regard to body 2 expressed in 1 R frame

1 2 / 0   is the absolute velocity of body 2 expressed in 1 R frame And   2 2 2 2 1 1 1 1 2/0 2/1 1/0 1 1 1 1 2/1 1/0 1/0 1 2 O O O O O V V V V V O O                     (6)
Which can be also written as

  2 2 1 1 1 1 1 1 /0 2/1 /0 1/0 1 2 O O O V V V O O             (7)
Where

2 1 /0 O V  is the absolute velocity of point O2 1 1 /0 O V  is the absolute velocity of point O1 2 1 2/1 O V  
is the relative velocity of point O2 with regard to the frame 1 R

1 1 1/0 1 2 O O     
is the velocity component due to the rotation of frame 1 R with regards to the inertial reference frame

The general kinematic joint model is then detailed in Figure 4.

Kinematic joint models

Based on the general kinematic joint model, models of common kinematic joints are given in Figure 5, Figure 6 and Figure 7.

1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 2 1 O O   2 1 2/1 O V   2 2/0   2 2 2/0 O V   2 1 P 1 0 P 2 0 P 1 2/1   2 1 P 2 1 P : : 1 1 1/0 2 1 O O      Figure 5 -Kinematic joint models -part 1 1 2/0   1 1/0   2 2/0   1 2/1   1 1/0 O V   2 2/0 O V   2 1 P 1 0 P 2 0 P 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 2 1 O O   2 1 2/1 O V   1 1/0   1 1/0 O V   2 2/0   2 2/0 O V   1 2 P 1 2/1   1 0 P 2 0 P 2 2/0   2 2 2/0 O V   2 1 P 1 0 P 2 0 P 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 1 1/0 2 1 O O      2 1 2/1 O V   2 2/0   2 2 2/0 O V   2 1 P 1 0 P 2 0 P 1 2 z z    1 y  1 x  2 x  2 y  1 x  2 x  2 y  2 z  1 z  1 y  1 2 z z    1 y  2 y  2 O 1 2/1   1 2 y y    1 2 z z    1 2 x x    1 O 2 O 1 P 1 P : : 1 2/1   1 1 1/0 2 1 O O      : : 2 1 P 2 1 P 2 P 2 P 1 2 1 O O   : : 1 x  2 x  1 O 1 2 P 1 2 P : :   1 2 O O O   1 2 O O O
0 0 ò 2 1 1 1 MTF 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 1 1/0 2 1 O O      2 1 2/1 O V   1 1 MTF MTF 2 2/0   2 2 2/0 O V   2 1 P 1 0 P P 1 1 L L 1 3 0 0 ò 3 2 1 1 1 MTF 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 2 1 O O   2 1 2/1 O V   1 0 MTF MTF 2/0   2 2/0 O V   2 1 P 1 0 P 2 0 P 1 1 L L 1 L 1 2/1   L 1 2 1 3 1 y  2 x  2 y  1 x  2 z  1 z  1 O 2 O 1 2 z z    1 y  1 x  2 x  2 y  1 O 2 O L TF:p 1 1 x  2 x  2 y  2 z  1 z  1 y  1 O 2 O 1 2 z z    1 y  1 x  1 O 2 O 1 2/1   1 0 0 ò 2 1 1 1 MTF 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 2 1 O O  1 1 MTF MTF 2 2/0   2 2 2/0 O V   2 1 P 1 0 P 2 P 1 1 L 1 1 3 1 2/1   1 0 0 ò 3 2 1 1 1 MTF 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 1 1/0 2 1 O O      2 1 2/1 O V   1 1 MTF MTF 2 2/0   2 2 2/0 O V   2 1 P 1 0 P P 1 1 L 1 1 1 1 2/1   L 1 2 1 3 1 dof 2 1 2/1 O V   4 dofs 5 dofs 4 dofs 1 2 1 P 2 1 P 1 1 1/0 2 1 O O      : 1 2 1 O O   : 2 1 P 2 1 P 2 1 P 2 1 P 2 1 P 2 1 P 1 1 1/0 2 1 O O      : 1 2 1 O O   :

METHODOLOGY OF MODELING FOR SIMULATION

General aspects

The first two steps of modelling a multibody system with bond graphs are choosing the dimension of the multi-body model (2D or 3D), and choosing the dimension of the bonds (scalars or vectors).

In keeping with the philosophy of a modular approach, the authors propose a 3D model and vector bonds. Depending on the modelling objectives, the user can transform the model to one with scalar bonds in order to exploit it later (as described in Figure 8).

Modelling with vector bonds

The first step in carrying out a simulation of a bond graph is the assignment of causalities. Two causality constraints appear with vector bonds, as noted by [START_REF] Ersal | A Bond Graph Based Modular Modeling Approach Towards an Automated Modeling Environment for Reconfigurable Machine Tools[END_REF] and 44 .

The first causality constraint (C1) is:

C1: each component of a vector bond must have the same causality.

Consequently, it is not possible to constrain the motion using Sf elements in only one or two dimensions without introducing some parasitic elements into the remaining unconstrained dimension(s).

The second causality constraint (C2) is: C2: the causality of transformers implied in the cross products, and the causality of gyrators in the bond graph model of the rigid body, are intrinsically fixed.  The transformers implied in the cross products must have flow-in-flow-out causality  The gyrators must have flow-in causality.

Indeed, contrary to scalar bond graphs, the ideal 2-port elements cannot propagate causality in both directions when the moduli are not invertible. In the context studied here (i.e. multibody system modelling with bond graphs and the Tiernego and Bos approach), the cases in which moduli (matrices associated to the elements) are not invertible are present for two elements: the transformer (TF) between the rotational and translational domain and the gyrator (GY) since both elements implement cross products. Consequently, transformers and gyrators have the mandatory fixed causality assignment specified above. The transformers and gyrators with the acceptable causal forms mentioned are given in Figure 9.

BG model Kinematic scheme

Plane joint

Fixed joint The causalities imposed on the transformers lead to some specific results when multibody systems with kinematic loops are considered. Figure 10 presents the bond graph model of a rigid body with vector bond and the imposed causalities. Consequently, in the bond graph of a rigid body, attaching flow sources (Sf) to more than one hinge point or the centre of mass at the same time is not possible and leads to causality conflicts. This situation typically occurs when the multibody system is composed of kinematic loops. In [START_REF] Ersal | A Bond Graph Based Modular Modeling Approach Towards an Automated Modeling Environment for Reconfigurable Machine Tools[END_REF] , the author gives the example of an oscillating bar which, with its two joints, is a closed kinematic chain system with one body. This case will appear in the slider crank model presented in this paper.

0 0 ò 3 2 1 1 1 MTF 1 1 1/0 O V   1 1/0   1 2/0   2 1 2/0 O V   1 1 1/0 2 1 O O      2 1 2/1 O V   1 1 MTF MTF 2 2/0   2 1 P 1 0 P 2 0 P 1 1 L 1 1 L 1 2/1   L 1 2 1 3 1 y  1 x  1 2 z z    2 x  2 y  1 O 2 O  1 z 1 y  1 x  2 x  2 y   1 2 O O  2 z 1 1/0   1 1/0 O V   2 2/0   2 2/0 O V   1 2 P 1 2/1   MTF MTF 0 0 1 1 1 1 1 0 P 2 0 P L L 2 2 2/0 O V   3 dofs 0 dof
Figure 10 -Bond graph model of rigid body with vector bonds and the imposed causalities

Modelling with scalar bonds

Bond graph models with scalar bonds do not have the same causality constraints as vector multi-bonds. However, in order to obtain the bond graph model with scalar bonds in a systematic way from the vector bond model, the same causality constraints (C1 and C2) on the gyrators and the transformers used for cross products are kept. Consequently, following the method chosen for the simulation of the bond graph model (presented in section 3.3), the scalar equivalent transformers representing coordinate transformation (here called coordinate transformation subsystems) should be adapted to respect the causality imposed by the rest of the bond graph. This is true of both the rigid bodies and the joints. The causality of these transformers depends on the elements used to constrain the motion.

These coordinate transformation subsystems can have two different structures: the structure with flow-inflow-out causality and the structure with effort-in-effortout causality. The causalities of these coordinate transformation subsystems are imposed by the position of the 0-junctions or 1-junctions with regard to the scalar MTFs. The coordinate transformation subsystem with flow-in-flow-out causality is built with the 1-junctions at its input in order to propagate causality to the rest of the elements. Following the same logic, the coordinate transformation subsystem with effort-in-effort-out causality is built with the 0-junctions at the input.

Figure 11 presents the two possible structures of the coordinate transformation subsystem and the cross product model for the BG model of the rigid body. The use of these coordinate transformation subsystems in the joints is illustrated in section 4 of this paper (see Figure 19). For the sake of clarity, the modulation signals for the MTFs are not displayed. Figure 12 and Figure 13 presents respectively the translational dynamics model and rotational dynamics model in scalar bond graphs. To the authors' knowledge, these models have been presented for the first time in this form in [START_REF] Ersal | Structural simplification of modular BG models based on junction inactivity[END_REF][START_REF] Ersal | Orienting body coordinate frames using Karhunen-Loève expansion for more effective structural simplification[END_REF] .

A procedure for transforming a multibody model with vector bonds into one with scalar bonds is as follows: 1) Identify the causality of the transformers linked to coordinate transformation (in bodies and joints). 2) Replace each of these transformers with the appropriate transformer subsystems, and assign the same causality as that seen on the vectorial MTF in the model with vector bonds: a) If the vectorial MTF is flow-in-flow-out causality, then use the coordinate transformation subsystem with flow-in-flow out causality (with 1-junctions at the input and 0-junctions at the output). b) If the vectorial MTF is effort-in-effort-out causality, then use the coordinate transformation subsystem with effort-in-effort-out causality (with 0-junctions at the input and 1-junctions at the output).

   2 1 1 2 e P e 2 e  1 f     2 2 1 1 T f P f    2 1 1 2 e P f    2 2 1 1 T e P f 1 f  2 f  0 0 1 TF: TF : 1 1 1 1 MTF : i G I i G EJS /0 i i G V  /0 i i   /0 k i M V  /0 j i M V  /0 i i i i k G M     i i j G M  0 i P Inertial frame Body frame Mandatory causalities /0 i i i i j G M     I: MGY: i i k G M  Se: 0 1 g W   1 m I:
i P     0 11 i P     0 12 i P     0 13 i P     0 21 i P     0 22 i P     0 23 i P     0 31 i P     0 32 i P     0 33 i P     0 11 i P     0 12 i P     0 13 i P     0 21 i P     0 22 i P     0 23 i P     0 31 i P     0 32 i P     0 33 i P     0 0 / i x G R V 0 / i iz G R V 0 / i ix G R V /0 ix i  /0 iz i  0 / j ix M R V 0 / j iy M R V 0 / j iz M R V 0 0 / i y G R V 0 0 / i z G R V 0 / i iy G R V /0 iy i  0 / i iz G R V 0 / i ix G R V 0 / i ix G R V 0 0 / i x G R V 0 0 / i y G R V 0 0 / i z G R V 0 / i iz G R V  0 / i ix G R V  0 / i ix G R V  0 0 / i x G R V  0 0 / i y G R V  0 0 / i z G R V  ix k GM 0 / k ix M R V 0 / k iy M R V 0 / k iz M R V /0 ix i  /0 iz i  /0 iy i  ix k GM ix k GM ix k GM ix k GM ix k GM 0 / k ix M R V 0 / k iy M R V 0 / k iz M R V /0  iz iz i J /0  ix ix i J /0 ix i  /0 iz i  /0 iy i  /0 ix i  /0 iz i  /0 iy i  /0  iy iy i J iy J iz J ix J 0 0 / i x G R V  0 0 / i y G R V  0 0 / i z G R V  0 0 9.81 0 g              3.2.

DAE formulation

In this paper, absolute coordinates are selected in order to keep a modular approach. Consequently, due to the kinematic constraints, derivative causality appears on the inertial elements and leads to differential-algebraic equations. It is important to notice that both open chain (OC) and closed kinematic closed (CKC) systems therefore lead to a DAE formulation. One of the priorities of the simulation methods proposed in this paper will therefore be to handle DAEs.

Review of the simulation Methods

The methods presented in this section come from the references given in section 2.1. The value of this study is in giving practical guidelines to modellers, so that they can implement the adapted method in 20-sim software.

Minimal coordinates method Description

The minimal coordinates method is based on formulating the dynamics using minimal set of joint coordinates (as generalized coordinates) to describe the mechanism's degrees of freedom.

Tiernego and Bos method leads to a model with a mixed differential-integral causalities on its inertial elements. The number of variables of inertial elements with derivative causalities depends on the variables of the inertial elements with integral causalities. The principle of the method is to transform the dependent inertia storage element (using transfomers) and to combine them with to the independent elements. This idea was first described in the context of bond graphs by Karnopp [START_REF] Karnopp | Power-conserving transformations: physical interpretations and applications using bond graphs[END_REF] .

Allen [START_REF] Allen Rr | Multiport representation of inertia properties of kinematic mechanisms[END_REF] uses this transformation on multibody systems where the dynamics are described with Lagrange equations and generalized coordinates. Breedveld 7 utilises it on a rigid body dynamics model corresponding to Newton-Euler equations with absolute coordinates. From a bond graph point of view, the minimal coordinates method consists of a two-level transformation through transformer elements as detailed in Figure 14. One level of transformers converts the velocities associated with inertial elements in absolute coordinates to the generalized velocities chosen as the joint's coordinates (also called relative coordinates). A second level of transformers converts the generalized velocities to the independent generalized velocities vector. For a system of n degrees of freedom, there are n independent velocities grouped in vector i q  . The generalized velocity vector q  defined for a mechanism can be partitioned into independent and dependent generalized velocities as follows:

i d q q q           (8)
The dependent generalized velocities can be written as a function of the independent generalized velocities.

d i q Tq    ( 9
)
Where T is a displacement-dependent velocity transformation between dependent and independent generalized velocities.

From equations ( 8) and ( 9), the generalized velocities can be written as a function of the independent generalized coordinates.

i q M q     (10) With 1 M T        
. This notation is inspired by [START_REF] Allen Rr | Multiport representation of inertia properties of kinematic mechanisms[END_REF] .

The inertial velocities of each body can be written as a function of the generalized velocities.

i i V T q   (11)
Where

i i i V T T or T   are displacement-dependent
velocity transformations between inertial and generalized velocities.

Allen [START_REF] Allen Rr | Multiport representation of inertia properties of kinematic mechanisms[END_REF] shows that the transformation of multiport inertia elements over a modulated transformer transformation leads to a virtual inertia Ĩ and a modulated gyristor element MGR. The 'virtual' denomination comes from the fact that the terms of this inertia are not constant. The modulated gyristor MGR comes from the

ki q   kd q   T k q     1 1/0 i T  n T    i/0 i 1 G I 1 G EJS 1 m i G I i G EJS n m n G I n G EJS  1 0 /0 G V  0 /0 n G V   n/0 n 1 T  1 V T i m  0 /0 i G V i V T n V T
modulated transformer with a variable modulus. In 7 , Breedveld formulated the displacement of I-elements and the gyrator GY element over a transformation from the inertial velocities to the generalized velocities. This leads to the following equations.

. . . m . i i i i i T T G V i V i i I T I T T T        (12) 1 .
. m .

        i i i i i T T G V i V i i MGR T I T T T (13) 1 .
.

    i i i T G i EJS T EJS T (14)
The second transformation leads to the virtual inertia Ĩ2, gyristor MGR2 and junction structure EJS2 associated with the independent generalized velocities. These are both modulated by the independent generalized parameters qi: . . . . . m . .

i i i i i T T T T G V i V i i I M T I T M M T T M            (15)     2 . . . . . . .m . . . i i i i i i i T T G i T T V i V V i MGR M T I T M T M M T T M T M                    (16) 2 .
. . .

i i i T T G i EJS M T EJS T M       (17)
The bond graph of the two-level transformation is given in Figure 15. This method requires three steps of symbolic manipulations. The first one is the determination of the i T and M  matrices from the kinematic relations. The second is the derivation of these matrices, which is required for the construction of the MGR elements. The third is the multiplication of non-linear matrices for the final computation of the added elements. Allen 46 conducts manual symbolic manipulation so as to determine the velocity relations from the displacement equations obtained by geometric analysis. Bos 6 has improved this method by using the bond graph to derive the velocity's transformation. Felez [START_REF] Felez | Bond graph formulation of in termes of relative coordinates[END_REF] and Borutzky 1 later presented a similar method to systematically deduce the kinematic relations from a modified bond graph model by the addition of controlled flow sources.

The technique of transforming storage elements from dependent to independent has been applied on industrial mechanisms: a forming machine [START_REF] Allen Rr | Multiport representation of inertia properties of kinematic mechanisms[END_REF][START_REF] Rr | Mechanisms as components of dynamic systems: a bond graph approach[END_REF] , a web cutting machine 49 , a pneumatic welding robot [START_REF] Det | Modélisation en bond graphs, modélisation et simulation d'une pince à souder[END_REF] and manipulator robots [START_REF] Field | Minimum energy planning trajectory for robotic manipulators[END_REF] . More recently, bond graph examples have been published: a planar pendulum in [START_REF] Breedveld | Concept-Oriented Modeling of Dynamic Behaviour[END_REF] , and a simple slider crank modelled from Lagrange equations using Karnopp procedure in [START_REF] Borutzky | Bond Graph Methodology -Development and Analysis of Multidisciplinary Dynamic System Models[END_REF] . In this paper, it is applied on a slider crank modelled with the modular Bos and Tiernego method, and can be easily used for more complex systems. As Borutzky points out [START_REF] Borutzky | Bond Graph Methodology -Development and Analysis of Multidisciplinary Dynamic System Models[END_REF] , this method essentially applies the well-known joint coordinate formulation of Nikravesh [START_REF] Nikravesh | Systematic Construction of the Equation of Motion for Multibody Systems Containing Closed Kinematic Loops[END_REF] to the bond graph modelling of multibody systems.

Advantages

All derivative causalities are removed, and only integral causality remains on the inertia associated with the independent generalized velocity vector. The virtual inertia will always be in integral causality because it is grouped to an inertia element in integral causality during the transformation.Consequently, the number of equations is reduced and an ODE formulation is achieved, which is much more compact than a DAE system. The ODE model can be solved with a fixed step solver such as 4th Order Runge-Kutta.

Drawbacks

The technique can be very demanding in terms of computer memory, and requires consistent mathematical simplifications by an expert user.

The modularity (i.e. the property whereby the mechanism can be assembled as a set of subsystems according to its structure) disappears. Hence, simulation of a mechanism's subsystems can no longer be conducted individually to check the whole model one piece at a time.

The ZCP opening method Description

Once causality is applied, different types of causal paths between elements can be identified. A closed causal path without any integration operations is called a zero-order causal path (ZCP) [START_REF] Felez | Simulation of multi-body systems using multi-bond graph[END_REF][START_REF] Cacho | Deriving simulation models from bond graphs with algebraic loops.: The extension to multibond graph systems[END_REF][START_REF] Van Dijk | Simulation of system models containing zero-order causal paths-I. Classification of zero-order causal paths[END_REF][START_REF] Romero | A minimal set of dynamic equations in systems modelled with bond graphs[END_REF][START_REF] Romero | Efficient simulation of mechanism kinematics using bond graphs[END_REF] . Bond graphs with ZCPs generate mathematical models with DAEs. There is a direct link between the nature of the ZCPs in the bond graphs and the index of the DAEs. The definitions of ZCPs are given by Felez [START_REF] Felez | Simulation of multi-body systems using multi-bond graph[END_REF]

q   kd q   T k q     1 1/0 i T  n T    i/0 i  1 0 /0 G V  0 /0 n G V   n/0 n 1 T  1 V T  0 /0 i G V i V T n V T I  Gr EJS
 Class 4 ZCPs: Causal cycles whose topological loops are closed. Only Class 4 ZCPs lead to DAEs with an index of 2. When the DAE's index is inferior to 2, a Backward Differentiation Formula (BDF) solver (such as the one in 20-sim) can handle these equations. In some cases (scalar bonds or vector bonds with no kinematic loops), the simulation can thus be conducted without additional specific elements.

Advantages

First, this method enables simulations to be conducted directly from the model without requiring the addition of elements (R/C elements or controlled effort sources) at the correction location. The initial physical model is therefore not changed. Second, this method is faster than the two following methods (R/C elements and controlled effort sources), as demonstrated in section 4.

Drawbacks

When multibody systems have kinematic loops, the bond graph may contain Class-4 ZCPs which lead to a DAE formulation with an index superior to 1. When this happens, the BDF solver often encounters difficulties in the numerical computation of the model. The ZCP opening method consists of opening the Class-4 ZCP using 'break variables.' This technique can be used at an equation level or at a graphical level. The classical technique uses modulated sources (MSe) at 1-junctions of the Class-4 ZCPs to open them. Without a systematic approach to the detection of the Class-4 ZCPs, this method may be difficult to use when dealing with a case that involves complex multibody systems. Hence, the authors have not used this method here.

The R/C parasitic element

Description

The R/C parasitic element method first appears as the "Stiff-compliance" approach in [START_REF] Dc | Analysis and simulation of planar mechanisms systems using bond graphs[END_REF] . In the literature, others terminologies can be found: singular perturbation 59 , parasitic elements [START_REF] Rideout | System partitioning and physical domain proper modeling through assessment of powerconserving model structure[END_REF][START_REF] Ersal | Realization-preserving simplification and reduction of dynamic system models at the graph level[END_REF] , virtual springs [START_REF] Wang | Modeling and Simulation of Robotic Systems with Closed Kinematic Chains Using the Virtual Spring Approach[END_REF] and coupling or pads [START_REF] Merzouki R | Intelligent mechatronic systems -Modeling, control and diagnostic[END_REF][START_REF] Bera | Bond graph modeling of planar prismatic joints[END_REF] This method is based on the introduction of parasitic elements: compliances and resistances in the bond graph models of the joints. When used with multibody systems, this method has two objectives: eliminating the kinematic constraints, and eliminating the derivative causality (which yields algebraic loops) so that an explicit solver may be used.

As mentioned in section 3.1.1, vector bond graphs impose some supplementary causality constraints (C1 and C2). These constraints can be enforced using flow sources (Sf), but they may lead to causal conflicts because of their flow-out causality. Enforcing the constraints with parasitic elements (R and C elements) can circumvent these conflicts since they can take effortout causality. Firstly, they allow the preservation of the causality assignments for all bonds of a multibond (C1 respected). Secondly, they allow the suppression of causality conflicts which may appear due to the causality of the transformers implied in the cross product and gyrators (C2 respected).

Advantages

With the R/C parasitic element method, all derivative causalities are removed. This leads to an explicit ODE model, which can be solved with a fixed step solver such as 4th Order Runge-Kutta.

Unlike the Lagrange multipliers method, the consistency between the initial conditions for CKC system are not mandatory. In [START_REF] Rideout | System partitioning and physical domain proper modeling through assessment of powerconserving model structure[END_REF] , the parasitic element method is even used to determine initials conditions. In addition, over-constrained (also called hyperstatic) systems can be simulated without any difficulties, as demonstrated in section 4. Karnopp and Rosenberg 5 state that "The idea of using artificial C-and R-elements to enforce constraints and thus to avoid derivative causality or differential-algebraic equations may appear to be a "brute-force" approach. This may be true, but first, a brute-force approach that is effective should not be discounted and, second, it has been argued that this approach is in many cases superior to the alternatives."

Drawbacks

This method introduces new elements (R and C elements) to the initial bond graph, whose parameters and initial conditions (in the case of C-elements) need to be specified. The values of the compliant elements must reflect the compliances which exist in all mechanical joints. The stiffness' that are used should be high enough so as to approximate the constraints and not change the dynamics of the system. However, high stiffness' can introduce high-frequency dynamics. This forces the solver to take very small integration steps to meet the tolerance criteria and, consequently, the simulation is slowed. This method therefore leads to a compromise between the accuracy of the results and the simulation time: the stiffer the system is, the more numerical errors are reduced but the higher the simulation time.

Implementation

Even theoretically where explicit solvers can be used, the authors recommend the use of the implicit MBDF (Modified Backward Differentiation Formula) solver for three reasons. The first is because the MBDF solver is faster than classic explicit solver when a stiff ODE is concerned. The second is because, as an implicit method, the stability is guaranteed. The third is that controlled effort sources necessitate use of the MBDF solver. Using the MBDF solver with R/C Elements too allows comparison between the use of R/C elements and controlled effort sources. In this paper, all simulations have been conducted with MBDF solvers.

The Lagrange multipliers method

Description

In the Lagrange-multiplier method, the constraint forces are modelled by controlled effort sources (MSe) rather than by parasitic elements. This method originates from Bos and Felez [START_REF] Bos | Multibody Systems in Terms of Multibond Graphs with Application to a Motorcycle Multibody system[END_REF][START_REF] Felez | BONDYN : a bond graph based simulation program for multibody systems[END_REF] and has been implemented in 20-Sim by Borutzky and Van Dijk [START_REF] Borutzky | Bond Graph Methodology -Development and Analysis of Multidisciplinary Dynamic System Models[END_REF][START_REF] Van Dijk | Simulation of system models containing zero-order causal paths-II Numerical implications of class 1 zero-order causal paths[END_REF][START_REF] Borutzky | Representing discontinuities by means of sinks of fixed causality Proc of the[END_REF] . The same concept is also implemented by Ersal [START_REF] Ersal | A Bond Graph Based Modular Modeling Approach Towards an Automated Modeling Environment for Reconfigurable Machine Tools[END_REF] with pseudo flow source (PSf) elements.

This method is similar to the parasitic elements in that the causality conflict is suppressed by removing the flow sources (Sf) and enforcing the constraint by other elements which do not have flow causality (here the controlled effort sources MSe)

The controlled effort sources (MSe) are computed such that the difference between velocities for the constraints is zero. The principle is to apply an effort equivalent to the one that the system would impose on a flow source. It would have a practical operation as a flow source but with effort-out causality instead of a flow-out causality.

The constitutive laws of the controlled effort sources are as follows:

-the usual constitutive law for an effort source, -the constraint: effort = e such that flow(e)=f.

Figure 16 illustrates the implementation of this constraint. The effort e(t) applied iteratively to the system is determined so that the difference ε(t) between the flow measurement f (t) and its set point fc(t) tends to zero. This implementation uses the constraint() function in 20-Sim. At every simulation step, this function induces an iterative procedure to find the force that keeps the velocity offset at zero within a given error margin. This iterative procedure is only supported by the Modified Backward Differentiation Formulation in 20-Sim software.

The effort-out causality of the controlled effort sources ensures that all the inertial elements (I elements) receive integral causality. The dependent states are therefore not visible as derivative causality. However, the controlled effort source establishes within itself algebraic dependencies and thus indicates an implicit form of equations. As stated in 1 , the DAE system has an index of 2 due to the fact that the constraint forces do not appear in the algebraic constraints but is in a semiexplicit form which can be solved by the MBDF solver.

Advantages

Contrary to the parasitic elements method, the modulated sources (MSe) do not need additional tuning because no additional parameters are added to the system. In other words, the order of the system is not modified because no new states (with their associated parameters and initial conditions) are introduced. The absence of supplementary parameters creates a bondgraph that describes the system ideally, within the limit of the numerical tolerance on the constraints equations.

Although some iterations are required during the simulation to satisfy the constraint equations at each time step, the computational load is comparable or better than the parasitic elements method, where the differential equations are truly explicit but very stiff.

Drawbacks

The bond graph obtained with this approach leads to implicit differential equations. Explicit integration algorithms therefore cannot be used, and the constraints can only be met within some numerical tolerances during the simulation. Consequently, the implementation of this method requires care.

Implementation

The Lagrange multipliers method is more challenging to implement than RC elements, because the kinematic constraints imposed by the joints may produce computational conflicts. This issue occurs on closed kinematic chains with over constrained multibody systems. Due to the topology of the system, more than one joint could impose the same constraint on the system, in which case the simulation will no longer be possible. The number of the controlled effort sources must not exceed the degrees of freedom that need to be eliminated.

In order to solve this problem, the redundant constraints must be removed. For the slider crank example in the last section of this paper, a practical detection method will be applied to determine the redundant constraints.

The system's initial conditions should be consistent with constraints. If a controlled effort source is attached to a 1-junction to keep the velocity at zero, and a nonzero initial velocity is imposed on an I-element connected to the same 1-junction, then a simulation cannot be conducted. Note at this point that 20-Sim does not offer an automatic correction of inconsistent initial conditions, but the user can always code the calculation of consistent initial conditions. Figure 16 -Implementation of the controlled effort source

CASE STUDIES

The systems chosen here are intentionally kept simple, in order to demonstrate the methods for simulating bond graph models of multibody systems.

System   flow t     0  c flow t    t   effort t Constraint()   p.effort t   p. flow t

System

Two classical systems are chosen: the planar pendulum as an example of an open-chain (OC) system, and the slider crank as an example of a closed kinematic chain (CKC) system. The physical models parameterized with absolute coordinate systems are given in Figure 17 and Figure 22. For ease of representation, the physical models are given with a planar representation. However, the bond graph models have been realised on 3D physical models. This will impact the modelling of the crankshaft, which only contains redundant constraints when modelled in 3D.

. Modelling with vector bonds

The models of the planar pendulum with the vector bonds are presented in Figure 18. Note that, when the Tiernego/Bos method is applied (ZCP method), inertial elements with differential causalities are present. When R/C or controlled effort sources methods are used, integral causalities are conserved in the inertial elements.

Figure 18 -Bond graph models with vector bonds of a planar pendulum

1 x  1 y  0 x  0 y  A B  1 1 1 , x y 1 G 1 (z 0)  1 1 ( 0, 0)     1 1 1 1 0 0 Sf Sf Sf R 1 1 1 1 MTF: MTF: Sf Sf 1 1 0 0 TF TF 1 1 1 1 ò Se: 1 MTF: Cardan L L L 0 0 R 1 1 1 1 MTF: MTF: 0 0 TF TF 1 1 1 1 ò Se: 1 MTF: Fixed frame Revolute joint Body 1 1 1/0   1 1 /0 G V  1 A/0 V  1 0 0 / G R V  1 0 P  3 1 0 1 0 0* 0 I P P P  1 1/0   1 0 /0 G V  1 0 P 1 1 /0 G V  1 A/0 V  0 1 g W   0 1 g W   With ZCP method With L element (R/C

Flow in/ Flow out causality

Effort in/ Effort out causality

1 1/0   0 1/0   1 1/0   1 1/0   1 /0 B V  GA   GB  0 0/0*   1 A/0 V  0 A/0* V  0 0/0*   0 A/0* V  0 A 1/0 V   1 0 P 0 0* P 0 0* 3 P I  1 0 P 0 0* P 0 0/0*   0 A/0* V  1 /0 B V  0 A 1/0 V   1 G I 1 G EJS 1 m 1 m I:
MGY:

I P P P  1 1/0   1 A/0 V  1 0 P 1 0 P 0 1/0   Cardan 0 1 1 1 1 x OG y z             0 1 1 1 1 x OG y z             0 0/0*   0 A/0* V  1 1/0   1 A/0 V  1 A/0 V  1 A/0 V  4.1.

Simulation results

For the ZCP method, the equations are DAEs because of the inertial elements with differential causalities. However, since the ZCPs are class 2 or 3 (not 4), 20sim's BDF solver can handle the DAE. For the R/C method, the inertial elements with integral causalities lead to explicit differential equations, which could be easily integrated using explicit algorithms such as 4th Order Runge-Kutta. In order to compare the three methods (ZCP, R/C and MSe), the BDF solver in 20-Sim is used. For the Lagrange multipliers method, the dependencies are not visible in form of derivative causality, but the equations still take an implicit form. Consequently, as previously discussed, explicit integration algorithms cannot be used.

From the results obtained (Table 1), similar conclusions to 41 can be made. The ZCP method is the fastest of the three methods. The computational time of the Lagrange multipliers method is still comparable to the ZCP method but needs the addition of elements: the control effort sources. It can also be seen that there is a slight difference in computational time between the models with scalar bonds and vector bonds. The models with vector bonds take a little longer to simulate than the ones with scalar bonds. In order to test the accuracy of the solution, these methods are compared to a bond graph model with minimal coordinates. This is the simplest possible model, presented in [START_REF] Breedveld | Concept-Oriented Modeling of Dynamic Behaviour[END_REF] and recalled in Figure 20. From the results obtained (Figure 21), it can be seen that the ZCP method is almost identical to the bond graph with minimal coordinates. The R/C method yields a short 'peak,' which is due to the excitation of the modes introduced by these elements. Classically, the slider crank is composed of three bodies: the crank, the rod, the piston and four joints. Depending on the choice of the joints on CKC system, a hyperstatic system (also called an over-constrained system) can appear. That is the case when the slider crank comprises three revolute joints and one prismatic joint (Figure 22). In this case, the system is over-constrained because the number of constraints is higher than the relative degree of freedom after connection (shown in Table 2). The last line of the table corresponds to the remaining DOFs after having closed the kinematic loop. This mechanical model can be modified into an isostatic one (Figure 23). The method is based on the reduction of DOFs from the reduction of constraints, as shown in Table 3. The redundant coordinates can be determined automatically using a structural multibody tool such as MapleSIM [START_REF] Renaud R | Modeling and Design of An Active Suspension System with Maple and MapleSim[END_REF] . The process of identifying and removing constraints is done numerically. MapleSIM achieves this by obtaining the Jacobian and residual of all the constraints, and then performing numerical projections to determine which rows/constraints have the most effect on the condition of the matrix.

A summary of the two models studied is presented in Table 4.

Hyperstatic model Modelling with vector bonds

The bond graph model of the hyperstatic slider crank with vector bonds and R/C elements method is presented in Figure 24. The R/C method relaxes the kinematic joint constraints. The dynamic equations are therefore transformed into an ODE form with no geometric constraints to deal with. Consequently, this method easily permits the simulation of hyperstatic CKC systems. In fact, the added springs have a very strong stiffness to keep a physical behaviour close to a rigid joint but, from a computational point of view, the hyperstatic nature of the model has disappeared thanks to the addition of the R/C elements.

When over-constraints are present, the Lagrange multipliers create more difficulties than the use of R/C elements methods. In this case (CKC system), more than one joint could impose the same constraint on the system. In order to solve the model with Lagrange multipliers, the redundant constraints must be eliminated. For that purpose, the model has to be changed to an isostatic one by joint altercation.

Modelling with scalar bonds

The model with scalar bonds and R/C elements is built with the procedure given in Section 3.1.2.

Isostatic model Modelling with vector bonds

When the definition of joints is modified by removing the redundant constraints, the simulation of the isostatic model can be conducted with both the R/C elements and Lagrange multipliers methods. The bond graph model of the isostatic slider crank with vector bonds and MSe elements is presented in Figure 25. Gr EJS ò reduce the DAE index. There is no automatic method for locating ZCPs on large bond graphs within platforms such as 20-sim at the time of writing. However, an automatic ZCP-location method would not significantly improve computation time and accuracy.

x OG y z             ò 0 2 2 2 2 x OG y z             ò 0 3 3 3 3 x OG y z             ò  3
Figure 31 shows that the R/C method and the Lagrange multipliers method perform comparably to the minimal coordinates method, with errors of less than 4e-3 m/s. The errors are dependent of some parameters: the value of the R/C method or the value of the integration error accepted for the solver. The proposed values of R/C elements (parasitic stiffness=10 -7 N/m and parasitic damping=200 N.s/m) are therefore considered satisfactory for this application.

CONCLUSIONS

This paper presents three methods for simulating bond graph models of multibody systems: ZCP, R/C elements and Lagrange Multipliers. For each method, the authors suggest both conditions and practical rules for application implementing. Additional considerations include the nature of the chain of Multibody systems, the nature of the system towards its joint constraints, and the nature of the bonds.

A method for transforming vector bond graphs to scalar bond graphs is also provided. Future work will automate this.

Numerical comparisons complement the results given by Van Dijk 41 and Felez [START_REF] Felez | Simulation of multi-body systems using multi-bond graph[END_REF] . There is no best unique solution for conducting bond graph simulations of multibody systems: all of these methods are correct, and the best will depend on the application. The authors suggest the following guidelines.  The R/C method is perhaps the most convenient. It permits the simulation of both iso and over constrained multibody systems with kinematic closed loops. It is tolerant of small inconsistencies in the initial conditions, and allows the use of a classical explicit solver. The computational time and accuracy are comparable to the other methods.  Lagrange multipliers must be implemented with care in the case of the over-constrained multibody systems, but give an ideal description of the system within the limit of the numerical tolerance on the constraints equations.  The ZCP method can be easy and quick to implement on simple multibody systems. When applying it to large systems, a systematic way of detecting the class-4 ZCPs would be advantageous.

APPENDIX A: SIMULATION PARAMETERS

The parameters used in the two systems (simple pendulum and slider crank) are the following.

COORDINATES Relations between inertial kinematic quantities and generalized velocities

The rotational inertial velocity of body 1 expressed as a function of the generalized coordinates and the associated transformation matrix is given in equations line (18).

T            
, 1 0 0 0 0 0 0 0 0 0 0 0 0

T              (18)
The translational inertial velocity of body 1 expressed as a function of the generalized coordinates and the associated transformation matrix is given in equations line (19).

    1 0 /0 sin 2 cos 2 0 G r r V                          ,     1 0 /0 sin 0 0 0 2 cos 0 0 0 2 0 0 0 0 G r r V x                                           ,    
V r r T                             (19)
The rotational inertial velocity of body 2 expressed as a function of the generalized coordinates and the associated transformation matrix is given in equations line (20).

T            
, 1 0 0 0 0 0 0 0 0 0 0 0 0

T              (20)
The translational inertial velocity of body 2 expressed as a function of the generalized coordinates and the associated transformation matrix is given in equations line (21) and line (22).

                                             ,            
                                                               (22)
The kinematic quantities of body 3 expressed as a function of the generalized coordinates and the associated transformation matrix is given in equations line (23).

Relations between dependent generalized velocities and independent generalized velocities

The flow relations between independent  and dependent velocities       , , x used to build the BG model of the slider crank are the following:

  2 2 ( ) ( ) ( ( ) 1 and
) rcos l rsin

T T               (23)                              (23)                   . T x and                                                    2 2
                                1 . . M T x and        

WiVV

  Weight vector of body i expressed in the inertial reference frame 0 R i/0   Angular velocity vector of body i with regard to the inertial reference frame expressed in the frame i R Linear absolute velocity (in regards to the inertial frame) of point M expressed in the frame i Linear relative velocity vector of point M of body i with regards to the body k expressed in the frame i R Angular parameters of the body i with regard to the inertial reference frame [rad] Model parameters mi Body i's mass expressed the inertia frame [kg]

Figure 1 -

 1 Figure 1 -Parametrization of the free rigid body 2.3. Modelling Rigid bodies Consider the architecture of a rigid body multibond graph model based on 1 , 6 , 42 , 43 .

Figure 2 -

 2 Figure 2 -Bond graph model of the rigid body

  and the translational velocity vector of body i's centre of mass in regards to the inertial frame of Figure2describes the kinematic relations (equation (3)) between the velocities of the two points of the body i (

Figure 3 -

 3 Figure 3 -Calculation of the Cardan angles and rotation matrix from the angular velocity

Figure 4 -

 4 Figure 4 -General kinematic joint model In the general kinematic joint model, MTF elements modulated by coordinate transformation matrix are used to express the kinematic quantities in the frame associated with the body they are connected to.

Figure 6 -

 6 Figure 6 -Kinematic joint models -part 2

Figure 7 -

 7 Figure 7 -Kinematic joint models -part 3

Figure 8 -

 8 Figure 8 -Steps for the BG exploitation of MBS

Figure 9 -

 9 Figure 9 -Causalities imposed with [Modulated] TF and GY elements

Figure 11 -

 11 Figure 11 -Bond graph model of rigid body with scalar bonds

  0

Figure 14 -

 14 Figure 14 -Equivalent mechanism BG model in relative coordinatesAs mentioned in section 2.3, the dynamics of a mechanism are represented by two elements. The rotational dynamics of a mechanism are described by Ielements, characterised by the inertia tensor I of the body and the corresponding gyrator MGY. The translational dynamics are also described by an I-element, characterised by the mass tensor mi.For a system of n degrees of freedom, there are n independent velocities grouped in vector i q  . The generalized velocity vector q  defined for a mechanism

1

 1

2 .

 2

Figure 15 -

 15 Figure 15 -Displacement of I and GY elements through two-levels of transformers

  as follows:  Class 1 ZCPs: The causal path is set between storage ports with integral causality and storage ports with differential causality. The associated topological loops are flat.  Class 2 ZCPs: The causal path is set between elements whose constitutive relations are algebraic (resistors are the most typical case). The topological loops are flat.  Class 3 ZCPs: A causal cycle whose topological loops are open. The causal path starts and ends in the same port of an element.

  ki

Figure 17 -

 17 Figure 17 -Physical models of a planar pendulum

Figure 20 -

 20 Figure 20 -Planar pendulum with minimal coordinates

Figure 22 -

 22 Figure 22 -Physical model of a hyperstatic slider crank

Figure 23 -

 23 Figure 23 -Physical model of an isostatic slider crank The different structure of the mechanical model (hyperstatic or isostatic) plays an important role in the simulation as shown in the next section.

Figure 24 -

 24 Figure 24 -Bond graph model with vector bonds and R/C elements method of the hyperstatic slider crank

Figure 29 -Figure 30 -

 2930 Figure 29 -Evolution of the centre of gravities in the inertial frame of the different bodies

Figure 31 -

 31 Figure 31-Difference between the different methods (Lagrange multipliers and R/C methods) and the model with minimal coordinates on the piston velocity

  Lagrange multipliers (integration error=1e-8) vs Minimal coordinates R/C method vs Minimal coordinates

Table 1 -

 1 Numerical comparison of the simulation methods for the planar pendulum model

				Computing
	System Constraints Bonds Methods	time
				(s)
			ZCPs	0.045
		Vector	R/C	0.145
	Planar pendulum (OC)	Isostatic	MSe Minimal coordinates	0.050 0.015
		Scalar	ZCPs	0.031
			R/C	0.6
			MSe	0.040

Table 2 -

 2 Analysis of the hyperstatic slider crank model

						Relative DOF after
	Stage	Joint concerned	Relative DOF	Joint type	Number of constraints	connection
	1	fixed frame and crank	6	revolute joint	5	1
	2	crank and rod	6	revolute joint	5	1

Table 3 -

 3 Analysis of the isostatic slider crank model

						Relative DOF after
	Stage	Joint concerned	Relative DOF	Joint type	No of constraints	connection
	1	fixed frame and crank	6	revolute joint	5	1
	2	crank and rod	6	rotational joint	4	2
	3	rod and piston	6	spherical joint	3	3
	4	piston and Fixed frame	1+2+3=6	prismatic joint	5	1

Table 4 -

 4 Main characteristics of the two slider crank models

		Hyperstatic slider crank model	Isostatic slider crank model
	Number of coordinates	18 (3 bodies)	18 (3 bodies)
	Number of constraints	20
		3 revolute joints (5 DOFs)
		+ 1 prismatic joint (5 DOFs)

Table 5 -

 5 Numerical comparison of the simulation methods for the Slider crank models

	System	Constraints Bonds	Methods	Computing time (s)
		Hyperstatic	Vector Scalar	R/C R/C		0.5 0.45
				ZCPs	Hard to conduct
	Slider Crank (CKC)	Isostatic	Vector	R/C MSe Minimal coordinates	0.19	0.51 0.330
			Scalar	ZCPs R/C	Hard to conduct 0.4
				MSe		0.300

Table 6 -

 6 Model parameters of the simple pendulum

	Parameter	Description	Value Units
	L	Bar length	1	m
	M	Bar mass	1	kg
	I	Bar moment of inertia 0.083 kg.m²
	B	Rotational damping in the revolute joint	10	N.s/rad
		Parasitic stiffness'	10 6	N/m
		Parasitic damping	200	N.s/m

Table 7 -

 7 Model parameters of the slider crank

	Parameter	Description	Value Units
	r	Crank length	0.2	m
	l	Rod length	1	m
	m1	Crank mass	1	kg
	m2		Rod mass	1	kg
	m3	Piston mass	1	kg
	I1	Crank moment of inertia	1	kg.m²
	I2	Crank moment of inertia	1	kg.m²
		Parasitic stiffness'	10 7	N/m
		Parasitic damping	200	N.s/m
	T	Input torque	10	N.m
	b	Rotational damping in the C revolute joint	2	N.s/rad
	c	Viscous damping in the prismatic joint	5	N.s/m
	APPENDIX B: FLOW RELATIONS BETWEEN
	DEPENDANT	AND	INDEPENDANT

1/0

Modelling with scalar bonds

The models of the planar pendulum with scalar bonds are presented in Figure 19.

Figure 19 -Bong graph models with scalar bonds of a planar pendulum 1 0

Modelling with scalar bonds

The models with scalar bonds are built in a similar way for both the R/C method and the Lagrange multipliers method. The bond graph model of the isostatic slider crank with scalar bonds and MSe elements is presented in Figure 26.

0 1 0 0* 0

Relative dependent coordinates

Absolute coordinates

Figure 28 -Slider crank model with minimal coordinates

Simulation results

A simulation of the dynamics of these models has been conducted. Identical results to [START_REF] Rideout | Teaching Multi-Body Dynamics in an Undergraduate Curriculum: An Intuitive and Explicit Formalism Based on Parasitic Elements[END_REF] have been obtained (Figure 29 and Figure 30).

A comparison of results in terms of computing time of the simulation is presented in Table 5. The simulation of hyperstatic kinematic closed systems (such as the slider crank with three revolute joints) is possible with R/C elements as we have previously mentioned, and the simulation time of the hyperstatic system is comparable to the isostatic system. Both the R/C method and the Lagrange multipliers method have comparable computation times. Their computational loads are bigger than that of the minimal coordinates, but remain within acceptable limits. Finally, as in the planar pendulum, models with vector bonds always take a little longer to simulate than the ones with scalar bonds.

For the ZCP method with both scalar and vector bonds, CKC systems induce several zero-order causal paths of class 1 and 4. The Class1-ZCPs are associated with dependences between energy storage elements, whereas the Class4-ZCPs are associated with causal cycles along the junction structures. The DAE index will typically be greater than one. Consequently, the model with ZCPs does not permit direct simulation with the MBDF solver. In order to enable simulation, a specific causality assignment or break variables should be used to

MTF: