Decoupling Kinematic Loops for Real-Time Multibody Dynamic Simulations

Earth moving equipment are typically equipped with hydraulic cylinder actuators to perform the designated tasks. Multibody modelling of such systems results in models with kinematic loops that couples the motion variables of the loop bodies. Iterative solutions will be needed to satisfy the loop constraints and the applied constraints, which require evaluation of the constraint Jacobean matrix. The size of the Jacobean matrix and the associated projections depends on the number of motion variables in each kinematic loop. Consequently, the computational cost significantly increases as the number of variables in the kinematic loop increases. Real-time control and hybrid hardware-in-the-loop systems both require efficient and fast computational algorithms. Eliminating the kinematic loops can improve the computational efficiency and effectiveness of the control algorithms. This paper presents an efficient approach to eliminate the coupling due to the cylinder-rod connections and consequently the kinematic loops in the multibody models leading to efficient simulation. The proposed approach calculates the spatial accelerations and inertia forces of the actuator bodies and the interaction forces with other components. The actuator forces are then projected onto the connecting bodies leading to exact dynamics of the system.


Introduction
Most earth-moving machineries carry main attachments that perform the designated functions of the machine e.g.loader bucket, dozer blade, excavator bucket, dozer ripper, etc.As shown in Figure 1, the motions of these attachments are controlled by linkages and mechanisms.Depending on the machine function and design, the control linkages topology could contain four-bar linkages (FBL), four-bar slider linkages (FBSL), or compound linkage.In complex compound linkages, FBL and FBSL could be connected in series or parallel combinations.Some product lines could have different model capacities while performing the same function.Due to the workload capacity, significant difference in the design and the function requirements could lead to different design configurations and kinematic topology of the linkages.In many designs, the buckets are attached to the rocker link of a four-bar linkage that controls the bucket rotational motion.Also, blades and rippers could be attached to the coupler link of a four-bar linkage to control the translational motion.In both cases, the four-bar mechanism is actuated by an inverted FBSL mechanism with a translational joint to model the cylinder-rod motion.In most earth-moving equipment, the translational motion is controlled by a hydraulic actuator, as shown in Figure 1.The kinematic topology of the FBL and FBSL lead to a Corresponding author: momar@taibahu.edu.sa,momar23@gmail.comkinematic closed-loops [1,2,3].Real-time simulation for such systems poses a challenge for software developers and for analysts [4,5,6,7].This paper presents an approach for replacing the linear actuator and rigid link bodies in the multibody system model with passive force and constraint elements whose kinematic and dynamic equations are expressed in terms of the attached body state variables [8,9].The spatial displacement and velocity of any number of markers fixed in the element are computed, and this information is used to compute interaction between the element and other parts of a model.Spatial contact forces between the system bodies and element bodies are computed and applied to their respective markers.This process could significantly reduce model size and computer simulation time by eliminating excess bodies, state variables, and kinematic constraints from the equations of motion.The actuator is attached to the main driven bodies through primary and secondary ends through universal and spherical joints respectively.No other bodies may be kinematically coupled to the passive element, and no kinematic constraints, other than a distance constraint between the two attachment points, may be imposed upon it.In order to implement the passive element, the kinematic motion need to be expressed, the constraint equations need to be formulated, the inertia forces need to be calculated and the reaction forces need to be evaluated.
This paper is arranged as follows, in section 2 the passive element's spatial displacement, velocity, acceleration, and constraints are expressed.In section 3, the inertia and spatial forces acting on the element are derived and the reaction forces are evaluated.An example will be presented in Section 4. The algorithms in this paper refer to an actuator with two common attachment points that must not be allowed to coincide.The linear actuator element contains a single internal translational degree of freedom (DoF) parallel to a line passing through its two attachment points, and this DoF could be driven by force or kinematically constrained.A distance constraint between the primary and secondary attachment points may kinematically control this DoF.The constraint can be a function of time and if the constraint is constant, the passive actuator will act as a rigid link.The universal joint with intersecting orthogonal rotation axes is assumed.The universal joint connects two marker, the first fixed in a driven body and the second is fixed in the passive element, serving as its primary reference.The universal joint's first DoF is aligned with one of the body marker axes and the second with one of the primary reference marker axes, as shown in Figure 1.The spherical joint at the other end is considered the secondary attachment point, and it connects two marker frames, one fixed in a driven body and the other in the passive element.The translational joint between the element's two components constrains the primary and secondary markers to always have the same orientation.Thus there is no need to compute the spherical joint's rotation angles.Bodies of the passive elements are not included in the model's kinematic tree, so the interaction forces through the connecting joints are applied directly to the respective bodies.Thus an element's mass and inertia terms and the associated kinematic constraint equations aren't included in the equations of motion.Instead, all inertial and applied forces are added to the right hand side of the equations of motion.

01001-p.2
Let m p and m s denote the primary and secondary markers each fixed in a respective driven body.The absolute position vector of these markers are represented by 0 0mp t and 0 0ms t .Since the primary and secondary attachment points coincide with the respective body marker origins, it follows that equate the global position of the mating marker pair origins.Let 0 0 0 0 0 ps s p t t t represent coordinates of a vector from the primary to secondary attachment point.Let d ps represent the distance between the attachment points and u ps represent coordinates of a unit vector parallel to the centerline and oriented from primary toward secondary attachment points as follows Equation ( 2) may be written in constraint form where ps d t is a function of time or constant as: Assuming that the axes of a right-handed Cartesian coordinate system are numbered (1,2,3), and let k denote the common centerline axis of the primary and secondary markers.The identifier k may be assigned a value of 1, 2, or 3, depending on desired definition.The (i, j) pair defines the universal joint's orthogonal rotation axes and the angles T and I represent the rotation angles about these axes.
Let (i, j, k) denote a triad of 3x1 vectors, each containing a 1 in the row associated with the numeric value assigned to the respective identifier in (i, j, k).Coordinates of the unit vector along the centerline, as defined in (2) and referenced to the primary frame, are p ps u k .This unit vector defined by k points from the primary toward the secondary attachment point.Since k contains a single unit entry in the kth row, the matrix product indicates that 0 ps u defines the kth column of  I Z (17) Now that the element's absolute spatial velocities are known, the spatial accelerations must be determined, so the actuator inertial forces can be computed and applied to the parent body attachment markers.Since high accuracy is not required for inertial forces and because an iterative method is not used in solving the equations of motion, the parent marker accelerations from the previous time step will be used in these calculations.Also needed is the relative rotational acceleration p mp Z for calculating the element's spatial accelerations.It could be obtained by differentiating the equation used to solve for 0 mp Z .Coordinates of the translational accelerations of markers m p and m s are 0   f as needed in the equations of motion.To determine these forces, we take advantage of additional information about the system by first expanding the spatial forces to show their rotational and translational components (32) The resulting spatial forces must be transformed and subtracted from the appropriate body marker frames to complete the equations.The passive actuator equations may appear to involve many multiplications, but the overall count when compared to implementing a regular actuator model is small.Also, eliminating the associated constraint Jacobian and mass matrix operations provides a substantial computational savings.01001-p.5

Descriptive example
Figure 1 shows CAD model and the kinematic tree diagram of a partial hydraulic excavator.In this example, the machine's undercarriage and the lower structure is not taken into consideration.The kinematic model consists of 16 dynamic bodies with 6-DoF.This model will result in 96 second order differential equations in the Cartesian space.Since each dynamic body is connected to its parent with single DoF joint (translational or rotational), the kinematic tree has total of 16 joint variables.According to the system topology 6 closed kinematic loops are formed.
In the joint-based multibody dynamics formulation, chord-joints could be modeled by introducing a massless chord-body that is rigidly attached to a marker on the parent body of the chord-joint.The chord-body is then connected to the child body of the chord joint by kinematic joint, as shown in Figure 2, [1,2,9].In order to maintain the integrity of the chord joint, 6 kinematic constraints are used to represent the rigid connection between the chord-body and chord-joint parent body at the joint locations.For example, in the kinematic loop L1, the massless chord-body 17 is rigidly attached to the joint parent body 2 and connected to the joint child body 5 with a spherical joint.The other loops are modeled similarly as shown in Figure 2 Figure 2. using the loop description in Table 1.In joint-coordinates multibody dynamics formulation, the equation of motion of system with closed loop could be written as follows [2,9]: where M is the system inertia matrix, T a is the system topology matrix, H a is the influence coefficient matrix of the system kinematic joints, H c is the influence coefficient matrix of the system chord joints, a a is the vector of Cartesian accelerations of all bodies in the system, F a is the vector of spatial reaction forces at all the kinematic joints, q a is the vector of kinematic joint variables and the corresponding derivatives, Q c is the vector of chord-joint forces , G is the vector of externally applied forces to bodies in the system, T v P is the vector of quadratic inertia forces, J a is the vector of quadratic velocity terms associated with all kinematic joints, Q a is the vector of kinematic joint forces, and J c is the vector of quadratic velocity terms associated with all chord joints.
The reduced matrix in (33) is shown partitioned where the block 3 by 3 block matrix corresponds to the augmented open-loop system equations.To obtain the joint variable accelerations, the Cartesian variables a a and F a can be systematically eliminated.Eliminating the Cartesian accelerations and reaction forces will lead to the system equations of motion as follows: The equation of motion of the multibody system with appended constraints could be rearranged as: where M q is the generalized mass matrix as defined in (9), J is the Jacobian matrix, and O q is the vector of Lagrange multipliers associated with the constraints.The problem dimensionality of the system under consideration has total of 16 spatial bodies where each body has 6 spatial accelerations.The dimensions of the inertia matrix M is (96x96), the connectivity matrix T a is (96x96), and the joint influence matrix H a is (96x16).The amount of expensive calculations required to transform the equation of motion from the augmented form in (33) 01001-p.6 to the joint-subspace form in (34) should be clear to the reader.Any saving or reduction in this projection step will significantly impact the simulation speed.To this end, the resulting equation of motion in (34) will be highly condensed and requires the evaluation of the Jacobean matrix J at every time step.The Jacobean matrix could be calculated from the derivatives of the chord joint constraint equations and from the other motion constraints imposed on the system as follows: where i q ) is the Jacobean matrix associated with the kinematic constraints of the closed loop i.Each kinematic loop requires 6 constraint equations.Iteratively satisfying the loop constraints require computing the Jacobean matrices.Decoupling the loop variables requires additional transformation and projections that will consume extra processing time.The Jacobean matrix needs to be analyzed to determine the dependent and independent variables that will be used by the integrator.Using the proposed approach leads to eliminating the kinematic closed-loops and the associated passive bodies from the kinematic tree.As a result, the differential equations of the system will significantly decrease.Figure 2. shows the kinematic tree of the resulting system after eliminating the kinematic loops.The first four kinematic loops are replaced with primary and secondary bodies.A time-dependent constraint or a joint force could be used to drive the translational joint.While the last two loops, two constant-distance constraints are imposed between the primary and secondary markers to represent the length of the coupler links.The remaining number of kinematic bodies reduced to 5, the number of joint variables becomes 5, and the kinematic loops are eliminated as shown in Figure 2.This led to significant reduction in the computational time of the whole model.After implementation, the dynamics simulation of the system runs more than 10-times faster than real-time with no accuracy loss.This achieved reduction enables the integration of multibody solver for more advanced simulations applications: Virtual reality, Hardware in the loop, real-time control etc.
DOI: 10.1051/ C Owned by the authors, published by EDP Sciences

Figure 1 .
Figure 1.A) CAD model of a hydraulic excavator.B) Model kinematic tree.C) Kinematic representation

Z
is the primary body marker's known absolute rotational velocity, and 0 mp Z is the element's rotational velocity relative to its parent marker.Because a translational joint connects the element's two components, 0 0 p Z represents the absolute rotational velocity of both parts.Manipulating the kinematic equations, equation for p mp Z could be expressed as cos sin the scalar quantities (Z 1 , Z 2 , Z 3 ) represent the coordinates of p mp Z in (16), we can have

gff
are the net spatial forces applied to the respective primary and secondary parts. of the corresponding universal and spherical joint reaction forces applied to the elements, and the forces , are the internal reaction forces developed between the mating elements at their common translational joint interface.must be computed and subtracted from the respective parent body attachment marker forces to complete the equations of motion.To simplify the development let:

Figure 2 .
Figure 2. Schematic description of the passive element kinematic connectionUsing the proposed approach leads to eliminating the kinematic closed-loops and the associated passive bodies from the kinematic tree.As a result, the differential equations of the system will significantly decrease.Figure2.shows the kinematic tree of the resulting system after eliminating the kinematic loops.The first four kinematic loops are replaced with primary and secondary bodies.A time-dependent constraint or a joint force could be used to drive the translational joint.While the last two loops, two constant-distance constraints are imposed between the primary and secondary markers to represent the length of the coupler links.The remaining number of kinematic bodies reduced to 5, the number of joint variables becomes 5, and the kinematic loops are eliminated as shown in Figure2.This led to significant reduction in the computational time of the whole model.After implementation, the dynamics simulation of the system runs more than 10-times faster than real-time with no accuracy loss.This achieved reduction enables the integration of multibody solver for more advanced simulations applications: Virtual reality, Hardware in the loop, real-time control etc.
R i and R j define rotational transformations about the respective universal joint axes identified by i and j.To determine the angles T and I, express (6) as follows represent the element's translational velocity components at the respective connecting points.The secondary marker's velocity relative to the primary marker is obtained from the general equation 0 p R .Considering p kps t is the kth and only nonzero entry of p ps t and substituting 0 ps u into (3) yields a more compact local constraint equation as follows Once the element's spatial displacements, velocities, and accelerations are known, the spatial forces acting on the respective element bodies are computed and applied to their parent markers.Separate free-body diagrams of the passive element's components yield