Kinematic and dynamic analysis of spatial multibody systems based on a formulation with fully Cartesian coordinates and a generic rigid body
Abstract
This work introduces the Fully Cartesian Coordinates Formulation with a Generic Rigid Body (FCC-GRB), a novel global multibody formulation for three-dimensional mechanical system analysis. The formulation's intrinsic characteristics are thoroughly detailed and compared with other widely-used global formulations, enabling its application in both kinematic and dynamic analysis of complex mechanical systems and as a teaching tool in advanced multibody dynamics courses. FCC-GRB formulation is founded on two main premises: multibody systems are described using only Cartesian coordinates, and the rigid bodies are modeled with a fixed and predetermined structure. Consequently, the kinematic constraints are described by lower-degree equations and the system mass matrix is highly sparse. Additionally, the introduction of the generic rigid body simplifies the modeling process by making the definition of the bodies independent of system topology. To reduce the number of generalized coordinates, a reduced modeling approach using less coordinates for describing the generic rigid body is also introduced and compared with the fully-defined alternative. The formulation's accuracy was validated through forward dynamic analysis of benchmark problems. Simulations demonstrated excellent agreement with reference data, with both modeling approaches yielding comparable kinematic results. The reduced approach offered faster computational performance, particularly in more complex models.
Full text
Research paper Kinematic and dynamic analysis of spatial multibody systems based on a formulation with fully Cartesian coordinates and a generic rigid body S´ ergio B. Gonçalves a , Ivo Roupa b , Paulo Flores c , Miguel Tavares da Silva a,* a IDMEC, Instituto Superior T´ ecnico, Universidade de Lisboa, Lisboa, Portugal b ITI/LARSyS, Instituto Superior T´ ecnico, Universidade de Lisboa, Lisboa, Portugal c CMEMS‑UMinho, Department of Mechanical Engineering, University of Minho, Guimar˜ aes, Portugal ARTICLE INFO Keywords: Multibody dynamics formulations Fully Cartesian coordinates Generic rigid body Kinematic analysis Dynamic analysis ABSTRACT This work introduces the Fully Cartesian Coordinates Formulation with a Generic Rigid Body (FCC-GRB), a novel global multibody formulation for three-dimensional mechanical system analysis. The formulation’s intrinsic characteristics are thoroughly detailed and compared with other widely-used global formulations, enabling its application in both kinematic and dynamic analysis of complex mechanical systems and as a teaching tool in advanced multibody dynamics courses. FCC-GRB formulation is founded on two main premises: multibody systems are described using only Cartesian coordinates, and the rigid bodies are modeled with a fixed and predetermined structure. Consequently, the kinematic constraints are described by lower-degree equations and the system mass matrix is highly sparse. Additionally, the introduction of the generic rigid body simplifies the modeling process by making the definition of the bodies independent of system topology. To reduce the number of generalized coordinates, a reduced modeling approach using less coordinates for describing the generic rigid body is also introduced and compared with the fully-defined alternative. The formulation’s accuracy was validated through forward dynamic analysis of benchmark problems. Simulations demonstrated excellent agreement with reference data, with both modeling approaches yielding comparable kinematic results. The reduced approach offered faster computational performance, particularly in more complex models. Abbreviations Description 2D Two-dimensional 3D Three-dimensional (continued on next page) * Corresponding author. E-mail addresses: [email protected] (S.B. Gonçalves), [email protected] (I. Roupa), [email protected] (P. Flores), [email protected] (M.T. Silva). Contents lists available at ScienceDirect Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmt https://doi.org/10.1016/j.mechmachtheory.2025.105955 Received 23 December 2024; Received in revised form 3 February 2025; Accepted 10 February 2025 Mechanism and Machine Theory 209 (2025) 105955 Available online 8 March 2025 0094-114X/© 2025 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ).
(continued) AR Angular Relation CAR Constant Angular Relation CP Coincident Point CV Coincident Vector DAE Differential Algebraic Equations DC Driver Constraint DoF Degree-of-Freedom EoM Equations-of-Motion FCC Fully Cartesian Coordinates FCC-GRB Fully Cartesian Coordinates with a Generic Rigid Body FD Forward Dynamics GRB Generic Rigid Body GSJ Grounded Spherical Joint ICC Intraclass Correlation KC Kinematic Constraint KJ Kinematic Joint ODE Ordinary Differential Equations OV Orientation of a Generic Vector PJ Prismatic Joint PRJ Pinned Revolute Joint RB Rigid Body RJ Revolute Joint RoM Range-of-Motion RGRB Reduced Generic Rigid Body SJ Spherical Joint TP Translation of a Generic Point UJ Universal Joint Symbol (Latin) Description 03, I3Null and identity matrices b τ , b− τ Moment arms of the equivalent force pair f τ and f− τ cP 1, cP 2, cP 3Scaling factors representing the local coordinates of point P in the local reference frame of the rigid body CP iConstant transformation matrix that describes the kinematics of a generic point P with respect to body i Cs iConstant transformation matrix that describes the kinematics of a generic vector s with respect to body i Cs1s2Constant transformation matrix that relates the generic vectors s1 and s2 e x , e y , e z Basis vectors that define the global reference frame of the multibody system fGeneric external force f fk RInternal reaction forces associated to the kinematic constraint of type k f τ , f− τ Equivalent force pair of the external moment of force τ gGeneralized forces of the system g2viVelocity-dependent inertial force of the rigid body i defined in its reduced form gf iGeneralized equivalent force of a concentrated external force f applied in body i g τ iGeneralized equivalent force of an external moment of force τ applied in body i gΦGeneralized internal forces gΦk iGeneralized internal forces for the kinematic constraint of type k with respect to body i Iii Components of the inertia tensor of the rigid body i miMass of the rigid body i MMass matrix of the system MiMass matrix of the rigid body i M2viEquivalent mass matrix of the rigid body i defined in its reduced form nb, nc, nq, nv Number of bodies, kinematic constraints, generalized coordinates and vectors of the system PGeneric point P P*Reference point P* OiOrigin of the local reference frame of body i q, ˙ q, ¨ qGeneralized positions, velocities and accelerations of the system q3vi, ˙ q3vi, ¨ q3viPosition, velocity and acceleration of the fully-defined equivalent vector of the reduced rigid body i ˙ q* iGeneralized virtual velocities of the rigid body i ˙ q* 3viGeneralized virtual velocities of the rigid body i described in terms of the fully-defined equivalent form rP, ˙ rP, ¨ rPPosition, velocity and acceleration vectors of the generic point P in the global reference frame rP*, ˙ rP*, ¨ rP*Position, velocity and acceleration vectors of a reference point P* in the global reference frame ˙ r* PVirtual velocities of the generic point P s, ˙ s, ¨ sPosition, velocity and acceleration vectors of a generic unit vector s s*, ˙ s*, ¨ s*Position, velocity and acceleration vectors of a reference unit vector s* SConstant transformation matrix that converts the reduced rigid body i in its fully-defined equivalent form tTime TKinetic energy ui, ˙ ui, ¨ ui vi, ˙ vi, ¨ vi wi, ˙ wi, ¨ wiPosition, velocity and acceleration of the generic rigid body vectors u, v and w of body i in the global reference frame uʹ i, vʹ i, wʹ iPosition of the generic rigid body vectors u, v and w in the local reference frame of body i uSkew-symmetric matrix of vector u used to compute the cross product (continued on next page) S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 2
(continued) VPotential energy Vi,Vi, ˙ ViTransformation matrices that convert the reduced rigid body i in its fully-defined equivalent form W* iVirtual Power generated by the internal forces of body i ˙ yState vector including the generalized velocities and accelerations utilized in the direct integration method Symbol (Greek) Description α , βBaumgarte stabilization coefficients γRight-hand-side vector of the acceleration equations of the system γkContributions to the right-hand-side vector of the acceleration equations for the kinematic constraint of type k θs1s2Angle between two generic unit vectors s1 and s2 θ* s1s2, ˙ θ* s1s2, ¨ θ* s1s2Angular position, velocity and acceleration of a reference angle θ* between two generic unit vectors s1 and s2 λLagrange multipliers of the system λkLagrange multipliers for the kinematic constraint of type k ν Right-hand-side vector of the velocity equations of the system ν tVector of the partial derivatives of ν with respect to t ν kContributions to the right-hand-side vector of the velocity equations for the kinematic constraint of type k ξ, η , ζBasis vectors that define the local reference frame of a generic rigid body ξP i, η P i, ζP iCoordinate of the point P in the ξ, η and ζ axes of the local reference frame of body i ρ iMass density of rigid body i τ Generic external moment of force τ τ int kInternal joint moment of force associated to joint k τ int kmInternal moment of force of joint k associated to the m-th angular driving constraint Φ, ˙ Φ, ¨ ΦVector of the kinematic constraints and respective first derivative and second derivative with respect to time ΦqJacobian matrix of the system ΦtVector of the partial derivatives of Φ with respect to time ΦkKinematic constraint equation of type k Φk qContributions to the Jacobian matrix from the kinematic constraint of type k ΩiGeometric domain of the rigid body i 1. Introduction Multibody dynamics methodologies have been successfully applied in the kinematic and dynamic analysis of complex mechanical systems, as they allow for efficient modeling and simulation of numerous problems, providing results with a high physical meaning [1]. In fact, multibody methodologies are frequently applied in many different areas, namely vehicle dynamics [2–4], mechatronics and robotics [5–7], mechanisms and machines [8–12], biomechanics and medical devices [13–18], co-simulation with finite elements method [19–21], media and videogames [22], and in real-time problems [23–27]. Over the last decades, several multibody systems formulations have been proposed, varying in the way the bodies are defined, the nature of the coordinates utilized, and the type of constraint and governing equations needed to describe the system kinematics and dynamics [28–31]. In general, multibody systems formulations can be categorized into two main types: recursive and global. Recursive formulations, also referred to as topological since the model definition depends on the system topology, can be further split into fullyand semi-recursive approaches. It should be noticed that each type of formulation comes with its own set of advantages and drawbacks, which significantly impact the simplicity of the modeling procedure, the systematization of the equations of motion (EoM), as well as computational implementation and efficiency [32]. For a more in-depth discussion on the differences, benefits, and practical applications of each type of multibody system formulation, the interested reader is referred to the works of Jal´ on et al. [33–35], Cuadrado et al. [36], Bae et al. [37], Roupa et al. [30] or Yu et al. [38]. In global formulations, the position and orientation of each body is usually described in relation to the global reference frame [28–30]. This type of formulation tends to generate dependent generalized coordinates, requiring the incorporation of constraint equations to describe such dependencies and to model the mechanical system. Despite increasing the size of the problem to solve, this approach greatly reduces the dependency of the EoM from the system topology, allowing for an easy systematization of the modeling procedure and definition of the governing equations. By and large, two main families of global formulations are available in literature to model multibody systems, chiefly the reference point coordinates [28,39,40], also referred to as Cartesian coordinates, and the natural coordinates [29,41]. Essentially, these two formulations vary on the type of coordinates adopted to model the bodies, which therefore affect the structure and complexity of the kinematic equations these generate. In Cartesian coordinates formulation, rigid bodies are typically described by a fixed kinematic structure composed of one point and a set of angular coordinates to describe respectively their position and orientation [28,42]. This modeling approach generates less generalized coordinates per rigid body than the other global formulations and makes the modeling of the rigid bodies independent of the system topology. Moreover, this formulation produces diagonal and constant mass matrices, and the joint reaction forces and moments of force can be obtained directly from the EoM. However, the use of angular variables implies that the kinematic constraint equations are described by nonlinear terms, which eventually penalize the computational efficiency [29]. In natural coordinates, the rigid bodies are defined resorting only on the use of rectangular coordinates of points and vectors, and, therefore, the constraints equations for the most common kinematic pairs present a linear or quadratic dependency on the generalized coordinates [29,43]. Due to the type of the elements that constitute the rigid body, the formulation with natural coordinates utilizes S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 3
more coordinates per body than the Cartesian coordinates. However, a benefit associated with the natural coordinates approach is the possibility of sharing points and vectors between bodies, defining directly kinematic pairs. This modeling approach allows for the reduction of the number of generalized coordinates and kinematic constraint equations needed to completely describe the configuration of multibody mechanical systems. Therefore, and depending on the topology of the system, the natural coordinates formulation with an implicit definition of the joints tends to produce a smaller number of coordinates than the Cartesian coordinates [29,43]. The description of the mechanical system by means of points located in relevant positions of each model segment implies that the mathematical definition of the rigid bodies and their mass matrices is dependent on the system topology, making the modeling procedure less systematic. Moreover, the use of shared elements also implies that the bodies are not completely independent, meaning that the system mass matrix is coupled, and the joint reaction forces cannot be directly obtained when the EoM are solved. To obtain the reaction forces directly, an explicit definition of the joints can be adopted, increasing the number of generalized coordinates and the computational burden [43]. Based on the natural coordinates, Gameiro et al. [44] proposed a novel global formulation that preserves the main characteristics of the natural and Cartesian formulations, that is, the bodies are modeled with a predetermined kinematic structure and are algebraically defined using only Cartesian coordinates. Similar approaches have been proposed by Uhlar et al. [45], and Pappalardo et al. [46–48]. More recently, Roupa et al. [30] presented a detailed description of the theoretical basis of the formulation for a planar case, formulating the constraint equations and contributions to the Jacobian matrix and right-hand side vectors of velocity and acceleration for the most common kinematic constraints. The formulation was subsequently applied to the inverse and forward dynamic analysis of different classical and biomechanical models, showing an excellent agreement with the benchmark data available in the literature [30, 49–51]. Since the formulation was built upon the two aforementioned premises, it was named as fully Cartesian coordinates with a generic rigid body (FCC-GRB). It is important to note that the term fully Cartesian coordinates were initially proposed by Jal´ on and his co-authors to describe the natural coordinates formulation, as this formulation considers only the use of Cartesian coordinates to describe the kinematics of the system. However, as the kinematic pairs arise naturally from the sharing of points and vectors, this formulation later became known as natural coordinates formulation [52]. According to Roupa et al. [30], the use of a fixed structure to describe the rigid bodies allows for an easier systematization of the modeling procedure, as their definition becomes independent of the topology of the system. It must be highlighted that this approach differs from the traditional modeling with natural coordinates, where the definition of the system varies with the topology and structure of the system, approximating the formulation to the one with Cartesian coordinates. In turn, the exclusive use of rectangular coordinates implies that, as in the natural coordinates approach, the constraint equations present, at most, a quadratic dependency on the generalized coordinates for the most common kinematic constraints [30]. Moreover, due to the type of the elements that compose the rigid body, it is possible to determine the kinematics of a generic point or vector directly from the generalized coordinates of the system using a set of constant transformation matrices. This procedure, similar to the one utilized in the natural coordinates formulation [29,53], systematizes the evaluation of the coordinates, velocities and acceleration of any point or vector of the system, simplifying the kinematic definition of the constraint equations, the application of external forces and moments or the mathematical derivation of the system mass matrix. Since these transformation matrices are constant in time, the contributions to the Jacobian matrix are mostly constant or linear on the generalized coordinates, meaning that the contributions to the right-hand side vectors of velocity and acceleration can be null, linear, or quadratic [30]. Thus, the present work extends authors’ previous developments [30] to formulate spatial multibody mechanical systems using fully Cartesian coordinates (FCC) together with a generic rigid body (GRB) formulation. In the sequel of this process, the main ingredients related to kinematics and dynamics of spatial multibody systems are presented in a detailed manner. Moreover, the main features of the proposed formulation are compared with natural and Cartesian coordinates, which allows for a discussion about the main benefits and limitations of the proposed approach. Finally, the kinematic and dynamic differences between modeling the generic rigid body with a fully-defined and reduced approach are also analyzed and their computational performance examined. The validity and accuracy of the formulation using both the fully-defined and reduced modeling approaches are assessed by performing a forward dynamic (FD) analysis of five classical benchmark problems. The trajectory of relevant points of the model, the violation of the kinematic constraints, the variation of the mechanical energy and the computational performance are compared with the equivalent models available in the Library of Computational Benchmark Problems available at International Federation for the Promotion of Mechanism and Machine Science website (IFToMM) [54,55]. In order to evaluate the differences in the computational time between a fully-defined or a reduced modeling approach, the FCC-GRB formulation is applied in the simulation of a n-four-bar mechanism with a variable number of segments. One of the applications previously identified for the planar FCC-GRB formulation is its use in teaching multibody dynamics in higher education and advanced courses [30]. Since the theoretical principles are maintained in the spatial version, this formulation could have an even greater impact in the educational field. Therefore, the main steps required for the complete implementation of the formulation are presented in detail throughout this work, so that, along with the planar approach [30], they can be directly used as a learning tool. 2. Fully Cartesian coordinates formulation in spatial multibody systems 2.1. Fundamental equations One of the main differences among global, semi-recursive, and recursive formulations lies in the type of governing equations that describe the system kinematics and dynamics. This distinction arises from the type of coordinates used to describe the system S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 4
configuration and the strategy employed to characterize the degrees-of-freedom (DoF) of the model under analysis. The formulation of the kinematic constraints and EoM for the FCC-GRB formulation is similar to other global approaches. Therefore, only the fundamental equations needed to describe the system kinematics and dynamics are presented in this section. The interested reader is referred to the planar version of the formulation [30] or the works of Nikravesh [28], Jal´ on and Bayo [29], Haug [40], Shabana [39] and Amirouche [56] for a detailed deduction of the equations. 2.1.1. Kinematic analysis The configuration of a multibody model is fully defined by a set of independent or dependent coordinates, usually referred to as generalized coordinates (q). Considering a system composed of nb bodies, the corresponding set of coordinates can mathematically be expressed as q={qT 1⋯qT i…qT nb }T(1) where qi represents the vector of the generalized coordinates of the i-th body. Depending on the type of coordinates utilized, the vector q can describe the position of a set of points and vectors of the system, the global orientation of the bodies, or the relative angular displacements between adjacent elements. In the presence of dependent coordinates, such as in the global formulations, the number of generalized coordinates (nq) surpasses the number of DoFs of the system, requiring the definition of a set of algebraic equations to describe their dependencies. These are usually introduced in the form of kinematic constraints (KC), and gathered in the vector of the kinematic constraints of the system as Φ(q,t) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Φ1 ⋮ Φi ⋮ Φnc ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ =0 (2) in which Φi is the i-th set of kinematic constraint equations, and nc the number of kinematic constraints. Besides describing the dependencies between the generalized coordinates, the vector of kinematic constraints also defines the structure and topology of the system, and the kinematics of the motion under analysis. This desideratum is achieved by mathematically describing the geometric relations existing between the generalized coordinates, thus defining the structure of the bodies (e.g., KC of rigid body type), the kinematic pairs and their corresponding DoFs (e.g., KC of joint type), among other topological features. Vector Φ also includes equations to express the geometric relations between the generalized coordinates and other external variables, allowing for the definition of topologic relations with external elements (e.g., KC with respect to global reference frame), the implementation of passive and active actuators (e.g., KC of actuator type), or the prescription of the DoFs of the system (e.g., KC of driver type). Mathematically, vector Φ is composed by a set of algebraic equations in their holonomic form that, when properly solved, allows for the computation of the kinematic consistent positions of the system. As some kinematic constraints exhibit a nonlinear behavior, it should be noted that Eq. (2) represents a system of nonlinear equations. Different methods have been proposed to solve the kinematic constraint equations, varying in their complexity and rate of convergence [57,58]. A common approach to compute Eq. (2) deals with the use of the iterative Newton-Raphson method or, when in the presence of redundant constraints, the iterative Newton-Raphson method with the least square approach [29,30]. The full characterization of the system kinematics requires the calculation of the consistent velocities and accelerations for the motion under analysis. For that purpose, the approach available in [29] can be applied, meaning that the generalized velocities of a multibody system (˙ q) defined with an FCC-GRB formulation can be directly obtained by solving the velocity constraint equations ( ˙ Φ) in order to ˙ q as ˙ Φ(q, ˙ q,t) = 0⇔Φq ˙ q= ν (3) with ν (t) = − Φt(4) where Φq represents the Jacobian matrix of the system, containing the derivatives of the kinematic constraints with respect to the generalized coordinates, ν is the right-hand side vector of velocities, and Φt denotes the vector containing the partial derivatives of the kinematic constraint equations with respect to time. Similarly, the generalized accelerations of the system ( ¨ q) can be obtained solving the acceleration constraint equations as ¨ Φ(q, ˙ q, ¨ q,t) = 0⇔Φq ¨ q=γ(5) with γ(q, ˙ q,t) = ν t−(Φq ˙ q)q ˙ q(6) S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 5
in which γ and ν t represent respectively the vector of the right-hand side vector of accelerations and the vector containing the partial derivatives of vector ν with respect to time. 2.1.2. Dynamic analysis The dynamic analysis of a multibody system requires the proper establishment of the EoM, as these equations express the relations between the internal, external, inertial, and velocity-dependent forces and moments that act on the system and its kinematics [32]. In the case of global formulations, the EoM are generic and their derivation independent of the system topology. In fact, this approach simplifies the dynamic analysis of complex multibody systems, since it does not require the deduction of the EoM for each system under analysis, as in the case of recursive formulations [30]. Thus, the EoM for the FCC-GRB formulation with a generic rigid body follows the same approach applied in the natural or Cartesian coordinates and can be expressed as [29,40,56] M¨ q+ΦT qλ=g(7) where M represents the mass matrix of the system, λ is the vector of Lagrange multipliers, and g denotes the vector of the generalized forces. Equation (7) describes the relation between the inertial (M¨ q), internal (ΦT qλ) and external forces (g) that act on the system. When the objective of the analysis is to determine the dynamic response of the system to a set of know external forces (forward dynamics), Eq. (7) represents a second order ordinary differential equation (ODE) that needs to be solved for accelerations, and subsequently integrated over time to obtain the generalized velocities and positions. However, since both the generalized accelerations and Lagrange multipliers are unknown variables, the system becomes underdetermined, and a new set of equations needs to be introduced to allow for its solution. For this purpose, the position constraint equations, given by Eq. (2), can be considered, yielding the follow system ⎧ ⎨ ⎩ M¨ q+ΦT qλ=g Φ(q,t) = 0(8) Equation (8) represents a set of index-3 differential algebraic equations (DAEs), where the first expression describes the dynamics of the system, and the second one ensures the kinematic consistency of the multibody system. To overcome the numerical problems of solving higher index DAEs, Eq. (8) is usually replaced by a DAE of lower order [59,60]. A common method to reduce the order of the equations is to force the kinematic consistency of the system using the acceleration constraint equations, expressed by Eq. (5), instead of the position constraint equations. Thus, the EoM for a constrained multibody mechanical system can be expressed as a set of linear algebraic equations, which, when solved from a forward dynamics perspective, represent a set of index-1 DAEs as ⎡ ⎣MΦT q Φq0⎤ ⎦{¨ q λ}={g γ}(9) It must be noticed that the reduction of the index of the DAEs can result in the violation of the kinematic constrains during the integration procedure, leading to the appearance of numerical instabilities and drift problems [35,60,61]. Over the decades, different methods have been proposed to deal with these numerical difficulties, such as stabilization procedures [61–64], augmented Lagrangian formulations [59,65,66], penalty formulations [25,67–69], velocity projections methods [70–73] or other partition methodologies [74–81]. 2.2. Fully Cartesian coordinates with a generic rigid body One of the key characteristics that distinguishes the different global multibody formulations is the approach utilized to represent the bodies that compose the system under study. The different methodologies available influence the number and type of the coordinates and kinematic constraints needed to properly define the system, the methods required to apply the external forces and moments, the definition of the system mass matrix, among other features. When using Cartesian coordinates formulations, the position and orientation of a given rigid body are described by the Cartesian coordinates of one point and a set of angular-related variables. If applied in the study of spatial systems, six coordinates are required to properly define the six DoFs associated with a free rigid body, namely three Cartesian coordinates of one reference point, and three Euler angles or Rodrigues parameters to describe respectively the translation of the body and its orientation. These representations are characterized by the existence of singular configurations, which can lead to the appearance of numerical issues during the dynamic analysis of the system [39]. A common approach to handle this issue is to replace this type of representation with another (e.g. Euler parameters or the Rodrigues formula), increasing the total number of generalized coordinates per body to seven. However, the incorporation of an extra coordinate leads to the appearance of dependent coordinates, requiring the introduction of one rigid body kinematic constraint [28,39,82]. An alternative solution to describe the kinematics of a multibody system is to utilize the natural coordinates formulation, originally proposed by Jal´ on and his co-authors [29,41,83]. With this formulation, the position and orientation of each rigid body is defined resorting solely to the use of Cartesian coordinates of points of interest and unit vectors, reason why they are referred to as fully cartesian coordinates [29]. An important feature related to the natural coordinates formulation is the possibility of sharing points and S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 6
vectors between different bodies, resulting in an implicit definition of several kinematic joints (e.g., spherical joint and revolute joint). As a result, the number of generalized coordinates and joint kinematic constraints needed to describe the configuration of a system is significantly reduced when compared with a procedure based on explicit definition of the joints. This particular characteristic is, in fact, the reason for the name “natural” associated with this formulation, as the kinematic joints arise naturally from the sharing of points and vectors. This modeling procedure is achieved by adopting a variable configuration for the definition of a given rigid body. Depending on the topology of the system to be discretized, different types of bodies with a variable number of points and vectors can be employed, such as two points and two vectors, three non-collinear points and one vector, four or more non-collinear points. As an example, a typical modeling approach to describe a link as a rigid body in a spatial mechanism involves using two points, usually located in its extremities, and two unit vectors, generating a total of 12 generalized coordinates. Consequently, a set of kinematic constraints equations needs to be introduced, at least six, to express the relations between the dependent generalized coordinates. These kinematic restrictions are typically introduced in the form of rigid body kinematic constraints, which include algebraic relations to ensure: (i) the constant distance between the points that compose the body, that is constant distance condition; (ii) the unitary nature of the unit vectors, that is unitary module condition; (iii) the constant angle between unit vectors, that is constant angle condition; (iv) the geometric relation between multiple points and vectors with respect to a reference frame composed by three selected vectors, that is linear combination conditions [29,43]. An alternative modeling approach, which uses a smaller number of vectors per rigid body, can be utilized to reduce the dimension of the system to be solved. Taking as example the case of the link considered above, this can be described by the Cartesian coordinates of two points and one unit vector, resulting in a total of nine generalized coordinates. Despite requiring less generalized coordinates and rigid body constraints, this procedure leads to the existence of variable mass matrices, which depend on the generalized velocities of the system. A common methodology to handle this particular aspect consists of decomposing the rigid body mass matrix in two terms, one dependent of the generalized coordinates, which is assembled in system mass matrix, and a second variable term that is treated as a velocity-dependent inertial force [29,53]. 2.2.1. Fully-defined rigid body In this work, a new approach to model spatial multibody systems with fully Cartesian coordinates is proposed. Contrarily to the natural coordinates formulation, in which the definition of each rigid body is dependent on the topology of the system, the presented formulation considers the concept of a generic rigid body. Thus, each body can be described with a fixed and predetermined kinematic structure composed of one vector that defines the reference point (rOi) and three unit vectors (u i , v i , w i ) to represent, respectively, its position and orientation in space (see Fig. 1a). Thus, the generalized coordinates vector for a given rigid body i (q i ) defined according to the proposed FCC-GRB formulation can be written as qi={rT OiuT ivT iwT i}T(10) With the purpose to keep the analysis simple, in what follows, the rigid body reference point Oi is located at the center of mass (CoM) of each body, and the rigid body vectors are aligned with the principal axes of inertia of the corresponding rigid bodies. An important feature associated with the FCC-GRB formulation is the definition of the rigid bodies resorting solely to Cartesian coordinates of points and vectors. This modeling approach has the advantage of not requiring the use of explicit angular-related variables, which eventually produce kinematic constraints that have a linear or quadratic dependency on the system generalized coordinates. The rigid body definition described in Eq. (10) generates 12 dependent generalized coordinates per body, requiring, at least, six constraint equations to describe the topologic relations between them. As in the case of formulation with natural coordinates, the topologic relations are introduced in the form of rigid body constraint equations, in particular the unitary module condition and the Fig. 1. Representation of the kinematic structure of a generic rigid body defined with FCC: a) Fully-defined variation (12 generalized coordinates); b) Reduced variation (9 generalized coordinates). S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 7
constant angle condition. However, in the FCC-GRB formulation, the kinematic constraints required to properly define the bodies are independent of the topology of the segments being modeled. For the case of a rigid body defined with one point and three vectors, three unit module conditions need to be added to the vector of the kinematic constraints (Ф) to guarantee the unitary nature of the rigid body vectors and three constant angle conditions to ensure that their relative orientation remains constant along the period of analysis, as it will be discussed in Sections 2.3.1 and 2.3.2. When compared to other common global formulations, defining a multibody system with FCC-GRB tends to generate a larger number of generalized coordinates. In fact, the FCC formulation generates the same number of coordinates per rigid body as the natural coordinates formulation. However, because the natural coordinates formulation allows for the sharing of points and vectors between bodies, it ultimately produces fewer total coordinates. To address this, the next section explores an alternative modeling approach where the base generic rigid body is substituted with a reduced generic rigid body, which comprises fewer generalized coordinates. 2.2.2. Reduced rigid body Knowing that from two non-collinear vectors it is possible to obtain a vector basis in 3D using orthogonalization methods [84], the kinematic structure of the generic rigid body can be further simplified to reduce the total number of generalized coordinates that each body generates. This procedure is performed by considering one vector describing the position of the reference point (rOi) and only two directional unit vectors (u i , v i ) to describe the reduced generic rigid body (RGRB) (see Fig. 1b), resulting in a total of nine generalized coordinates as qi={rT OiuT ivT i}T(11) It is important to notice that although the reduced modeling approach decreases the total number of generalized coordinates and kinematic constraints of rigid body type (eliminating two constant angle conditions and one unit module condition), it leads to the appearance of variable mass matrices. Similar to the natural coordinates formulation, this issue can be addressed by dividing the mass matrix into two components: one that is incorporated into the system mass matrix and another that is treated as an external force (see Section 4.2). 2.2.3. Relation between a reduced and fully-defined rigid body It should be noted that, depending on the modeling strategy adopted to describe a generic rigid body (fully-defined or reduced), the vector containing its generalized coordinates (qi) differs. Consequently, the kinematic constraint equations that define both the topology and the guiding of the system are not identical. To address this issue, this section introduces the kinematic relations required to convert the reduced system into its fully-defined equivalent. These relations provide a consistent framework for defining the multibody system and the necessary mathematical tools to derive expressions specific to the reduced approach. Despite allowing for the decrease of the number of generalized coordinates, the use of a reduced rigid body implies that the vector basis, composed of vectors ui, vi and wi, and which defines the inertial frame of the body (ξ, η , ζ), is not fully determined. Therefore, a procedure to transform the set of generalized coordinates of the reduced rigid body i (qi) into the equivalent fully-defined vector (q3vi) is first required. This desideratum can be achieved by determining a vector basis using the rigid body vectors (ui and vi) and orthogonalization methods. Accordingly, let one consider an auxiliary vector wi perpendicular to the plane defined by the two rigid body vectors of body i, such that wi=ui×vi(12) Thus, it is possible to define a set of algebraic equations that relate the set of reduced generalized coordinates of the rigid body i (qi) and the fully-defined equivalent vector (q3vi) considering a transformation matrix Vi, such that q3vi={rT OiuT ivT iwT i}T=Viqi(13) with Vi=⎡ ⎢ ⎢ ⎣ I30303 03I303 0303I3 0303 ui ⎤ ⎥ ⎥ ⎦(12×9) (14) where I3 and 03 are a 3 ×3 identity and null matrices, respectively, and ui represents a skew-symmetric matrix that allows for the computation of the cross product presented on Eq. (12), as wi= uivi=⎡ ⎣ 0−uizuiy uiz0−uix −uiyuix0⎤ ⎦⎡ ⎣vix viy viz ⎤ ⎦(15) By algebraically manipulating Eq. (14), it is possible to rewrite it in terms of a constant matrix S and the transformation matrix Vi, such that S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 8
q3vi=SViqi(16) with S= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ I3030303 03I30303 0303I303 030303 1 2I3 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦(12×12) (17) and Vi=⎡ ⎢ ⎢ ⎣ I30303 03I303 0303I3 03− vi ui ⎤ ⎥ ⎥ ⎦(12×9) (18) The expression for the velocity ( ˙ q3vi) of the fully-defined equivalent vector can be obtained by differentiating Eq. (16) with respect to time, as ˙ q3vi={˙ rT Oi ˙ uT i ˙ vT i ˙ wT i}=S˙ Viqi+SVi ˙ qi(19) By employing the following cross product relation viui= uT ivi= − uivi(20) the expression presented in Eq. (19) can be reduced to ˙ q3vi= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 030303 030303 030303 03−1 2˙ vi 1 2˙ ui ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ rOi ui vi ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ + ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ I30303 03I303 0303I3 03−1 2 vi 1 2 ui ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ˙ rOi ˙ ui ˙ vi ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ =⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ˙ rOi ˙ ui ˙ vi ui ˙ vi− vi ˙ ui ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ =Vi ˙ qi(21) In turn, the acceleration (¨ q3vi) for the fully-defined equivalent vector can be obtained by differentiating the previous expression with respect to time, yielding ¨ q3vi={¨ rT Oi ¨ uT i ¨ vT i ¨ wT i}T=˙ Vi ˙ qi+Vi ¨ qi(22) with ˙ Vi=⎡ ⎢ ⎢ ⎣ 030303 030303 030303 03−˙ vi˙ ui ⎤ ⎥ ⎥ ⎦(12×9) (23) It should be noted that by utilizing the two transformation matrices V and ˙ V along with Eqs. (16) to (23), the kinematics of the reduced rigid bodies can be fully determined, providing the necessary mathematical relations for conducting the kinematic and dynamic analysis of the multibody system. For evaluating the kinematic constraints, Eq. (13) and matrix V can be used directly instead of Eq. (16), thereby reducing the computational effort of this step. 2.3. Kinematic constraints of rigid body The establishment of the kinematic constraints of rigid body in the vector of kinematic constraints is required to ensure the rigid nature of the body. These constraints are defined in the form of algebraic equations that relate the dependent generalized coordinates that describe the rigid body. It is worth noting that aligning the rigid body vectors with the principal axes of inertia of the segment implies the existence of a direct relation between the generalized coordinates of the rigid body and its rotation matrix. However, since the formulation does not include angular variables, the geometric properties intrinsic to the definition of the body topology, namely the orthogonality of the basis vectors defining the local reference frame and their unitary nature, must be explicitly defined as kinematic constraint equations. As described previously, with the definition of a generic rigid body, the modeling procedure is significantly simplified, requiring the S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 9
3.2.5. Constant angular relation between two vectors belonging to fully-defined rigid bodies The constant angular relation (CAR) condition ensures that two different generic vectors maintain their relative direction along the time. Besides ensuring the geometric relation between the two unit vectors, the CAR condition is also useful to define joints characterized by two rotational DoFs. By constraining the relative angle between the two vectors, the CAR condition allows for the free rotation of the bodies in the planes orthogonal to the axes defined by each vector. This condition, together with the CP condition, can be used to model universal-type joints. CAR condition can also be applied to define geometric relations that allow for the definition of prismatic joints, sliders, among other mechanical elements. As in the case of the constant angle condition, the CAR relation can be expressed by one algebraic equation that forces the dot product of the two vectors to be equal throughout the period of analysis. Thus, let one consider two generic unit vectors s1 and s2 belonging respectively to bodies i and j (see Fig. 5), it is possible to define the CAR constraint (ΦCAR) in its homogenous form as ΦCAR(qi,qj)=sT 1s2−cos(θs1s2) = (Cs1 iqi)T(Cs2 jqj)−cos(θs1s2) = 0 ΦCAR =qT iCs1s2qj−cos(θs1s2) (64) with Cs1s2=[Cs1 i TCs2 j](12×12)(65) Equation (64) represents one constraint equation that presents a quadratic relation between the generalized coordinates of bodies i and j, and, consequently, the contributions to the Jacobian matrix have a linear structure in the form, such that ΦCAR q={qT jCs1s2T ⏞⏟⏟⏞ ΦCAR qi qT iCs1s2 ⏞⏟⏟⏞ ΦCAR qj}(1×24) (66) As the matrices Cs1 i and Cs2 j are constant and not time-dependent, matrix Cs1s2 has also the same properties. Hence, Eq. (64) represents a scleronomic relation, meaning that its contribution to the right-hand side vector of velocity ( ν CAR) is null. In turn, the quadratic nature of the CAR condition implies that the contributions to the right-hand side vector of acceleration (γCAR) present a quadratic relation on the generalized velocities of the system, as presented below ν CAR =0 (67) γCAR = − 2˙ qT iCs1s2˙ qj= − 2˙ sT 1 ˙ s2(68) Although Eq. (64) presents a non-linear term, this only needs to be computed once during the analysis. Since matrix Cs1s2 is constant, it must be evaluated one time at the beginning of the analysis. Thus, during the analysis of the mechanism, the CAR constraint condition is described only by algebraic equations that are null, linear, or quadratic in type, having a strong influence on the computational performance of the formulation. Moreover, Eqs. (64) to (68) can be further simplified if the vectors being constrained are rigid body vectors. Accordingly, let one consider two rigid body vectors (e.g., ui and vj) belonging respectively to bodies i and j, the CAR constraint (ΦCAR) in its homogenous form can be given by ΦCAR(qi,qj)=uiTvj−cos(θuivj)=0 (69) ΦCAR q={vT j ⏞⏟⏟⏞ ΦCAR ui uT i ⏞⏟⏟⏞ ΦCAR vj}(1×6) (70) γCAR = − 2˙ uT i ˙ vj(71) In the case of rigid body vectors, the kinematic constraint equation, and the contributions to Φq, ν and γ for the CAR condition maintain the same degree of linearity as in the case with the generic vector, and hence it offers the same advantages. However, obtaining all the required values entails fewer mathematical operations, improving the computational efficiency of the formulation. 3.2.6. Constant angular relation between two vectors belonging to reduced rigid bodies The definition of the CAR condition for the reduced rigid bodies can also use the dot product constraint explored in the previous section and the transformation in the fully-defined equivalent vector presented in Section 2.2.3. Accordingly, let one consider two generic unit vectors s1 and s2 belonging respectively to bodies i and j defined in its reduced form, then, the CAR constraint (ΦCAR) in its homogenous form can be expressed as: ΦCAR(qi,qj)=sT 1s2−cos(θs1s2) = (Cs1 iSViqi)T(Cs2 jSVjqj)−cos(θs1s2) = 0 ΦCAR =qT iSViCs1s2SVjqj−cos(θs1s2) (72) Due to the dependency of the matrices V on the generalized coordinates of the system, Eq. (58) represents one constraint equation S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 16
of higher degree. It is important to note that its evaluation resorts only on algebraic multiplications involving linear matrices and vectors dependent on the generalized coordinates of the system. This higher complexity also results in more complex expressions for the contributions to the Jacobian matrix and vector γ, in the form ΦCAR q={qT jVT jSTCs1s2TVi ⏞⏟⏟⏞ ΦCAR qi qT iVT iSTCs1s2Vj ⏞⏟⏟⏞ ΦCAR qj}(1×18) (73) ν CAR =0 (74) γCAR = − (qT iVT iSTCs1s2˙ Vj ˙ qj+qT jVT jSTCs1s2T˙ Vi ˙ qi+2˙ qT iVT iCs1s2Vj ˙ qj) γCAR = − (sT 1Cs2 j ˙ Vj ˙ qj+sT 2Cs1 iV ⋅ iq ⋅ i+2˙ sT 1 ˙ s2)(75) 3.3. Kinematic constraints of driving type Often, the kinematic analysis of a multibody system requires the prescription of all the DoFs associated with the mechanical model. This information is included in the form of kinematic constraint equations of driving type, being usually divided into two major groups: (i) Linear driving constraints – comprise the expressions needed to prescribe the translations of specific points or the orientation of vectors; (ii) Angular driving constraints – include the equations required to depict the angular rotations between vectors/bodies. The driving constraints are usually formulated using algebraic equations that relate the generalized coordinates and the kinematics of the prescribed movement. Consequently, this type of constraints presents a rheonomic nature, meaning that their expressions have an explicit dependency on the time vector. As for the topological constraints, several conditions can be applied to fully describe the kinematics of the mechanical elements to drive. In this work, particular focus is given to the formulation of the linear translation of one point, the orientation of one vector and the angular relation between two vectors. It is noteworthy that, similar to the topological constraints, kinematic constraints of driving type can involve algebraic expressions to prescribe the relative movement between two elements within the system. Alternatively, they can encompass specific equations to describe the motion of one element of the system in relation to an exterior element. In the first group, one can consider the angular relation constraint to guide the angle between two system vectors, and in the second group the translation of one point, the orientation of one vector in space, or the relative angle between a system vector and an external vector. 3.3.1. Translation of one generic point belonging to a fully-defined body The translation of one generic point (TP) condition ensures that one point of interest belonging to a given rigid body of the model follows the trajectory of one reference point along the period of analysis. Mathematically, this condition can be defined by stating that the point of the model shares the same position as the reference point in the global reference frame. Hence, the TP condition can be defined in a similar way to the approach used for modeling the CP condition. Accordingly, let one consider a point P belonging to body i and the respective reference point P*, in which the position (rP*), velocity (˙ rP*) and acceleration (¨ rP*) are known (see Fig. 6a), the TP constraint (ΦTP) in its homogenous form can be stated as ΦTP(qi,t) = rPi−rP*(t) = CP iqi−rP*(t) = 0(76) Fig. 6. Schematic representation of the linear kinematic constraints of driver type: a) TP condition; b) OV condition. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 17
Equation (76) represents three linear equations that enforces the coordinates of the point P to be equal to those of the reference point P*. Hence, this constraint drives three DoFs, which are associated with the translation of the body in space. For that reason, the TP condition is usually applied to either prescribe the position of the entire model in space by guiding the position of one point in the model parent body, or to guide both the position and orientation of the rigid bodies of the system by minimizing the distance between the prescribed position of several reference points and their corresponding interest points [16]. As the second term of the Eq. (76) does not depend on the generalized coordinates of the system, the contributions of the TP constraint to the Jacobian matrix are constant and equal to the matrix C ΦTP q=[CP i ⏞⏟⏟⏞ ΦTP qi](3×12) (77) Since the TP condition presents a rheonomic nature, the contributions to the right-hand side vectors of velocity ( ν TP) and acceleration (γTP) are not null, being their values equal to the velocity (˙ rP*) and acceleration (¨ rP*) of the reference point P*, respectively ν TP =˙ rP*(t)(78) γTP =¨ rP*(t)(79) 3.3.2. Translation of one generic point belonging to a reduced rigid body The formulation of the TP condition for a reduced rigid body is similar to that presented for a fully-defined body. Since this approach requires calculating the fully-defined equivalent vector by means of Eq. (16), both the kinematic constraint equations and the contributions to the Jacobian matrix and γTP reflect the quadratic nature of the constraint as a result of the multiplication of the transformation matrix V by the generalized coordinates of the body. Accordingly, the expression for the TP condition in its homogenous form and the contributions for the Φq, ν and γ are defined as ΦTP(qi,t) = CP iq3vi−r* P(t) = CP iSViqi−r* P(t) = 0 (80) ΦTP q=[CP iVi ⏞⏟⏟⏞ ΦTP qi](3×9) (81) ν TP =˙ r* P(t)(82) γTP =¨ r* P(t) − (CP i ˙ Vi ˙ qi)(83) 3.3.3. Orientation of one generic vector belonging to a fully-defined body The orientation of a generic vector (OV) constraint is a kinematic driver condition that guides the Cartesian components of one vector in space. This constraint can be utilized to prescribe the direction of any generic vector, enabling the definition of the global orientation of the model segments or the rotation axes of the kinematic joints in space. From a mathematical point of view, the OV constraint is similar to the TP condition. However, the expression uses the matrix C for a generic vector s instead of the one used for a generic point P. Hence, the OV kinematic constraint (ΦOV) in its homogenous form can be stated as ΦOV(qi,t) = si−s*(t) = Cs iqi−s*(t) = 0(84) where si is a generic vector belonging to body i, and s* represents the prescribed Cartesian components of the reference vector to follow along the period of analysis (see Fig. 6b). Consequently, the respective contributions to the Jacobian matrix (ΦOV q) and right-hand side vectors of velocity ( ν OV) and acceleration (γOV) are similar to the ones presented in Eqs. (77) to (79), but considering instead the matrix substitution referred before and the velocity (˙ s*) and acceleration (¨ s*) of the vector s*. If the vector to drive is a rigid body vector, then the OV kinematic constraint can be simplified. Accordingly, let one consider the vector u from the rigid body i, the OV algebraic condition and the respective contributions to the Jacobian matrix are given by ΦOV(qi,t) = ui−s*(t) = 0(85) ΦOV q=[I3 ⏞⏟⏟⏞ ΦOV ui](3×3) (86) In turn, the contributions to vectors ν and γ are equal to ones presented for a generic vector s. 3.3.4. Orientation of one generic vector belonging to a reduced rigid body As in the case of the fully-defined body, the OV constraint equation for a reduced rigid body is similar to the one presented in the TP condition, considering the substitution of the constant matrix C for a generic point P by the equivalent matrix for a generic vector s. As S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 18
so, the OV kinematic constraint (ΦOV) for a reduced rigid body in its homogenous form can be stated as ΦOV(qi,t) = Cs iq3vi−s*(t) = Cs iSViqi−s*(t) = 0(87) Similarly, the contributions to the Jacobian matrix and vectors ν and γ will be analogous to the ones presented in Eqs. (81) to (83), using the equivalent matrix Cs i and the velocity and acceleration vectors of the reference vector s*. 3.3.5. Angular relation between two vectors belonging to fully-defined rigid bodies The angular relation (AR) kinematic driver enforces that the angular displacement between two vectors is equal to a prescribed angle θ*. Consequently, this type of constraint can be used to guide the rotational DoFs associated with the kinematic joints, as it is schematically represented in Fig. 7. Mathematically, the angular driver constraint can be defined using the dot product between two vectors as it was described in Section 3.2.5. However, since the angular displacement between the two vectors varies according to the prescribed angle, the AR condition becomes a rheonomic constraint, meaning that additional terms need to be included to the contributions to the right-hand side vectors of velocity and acceleration. Hence, let one consider the DoF associated with the relative angle between two generic unit vectors s1 and s2 belonging respectively to bodies i and j, the AR kinematic driver constraint can be defined as ΦAR(qi,qj,t)=sT 1s2−cos(θ* s1s2(t))=(Cs1 iqi)T(Cs2 jqj)−cos(θ* s1s2(t))=0 ΦAR =qT iCs1s2qj−cos(θ* s1s2(t))(88) where θ* s1s2 is the prescribed angle between vectors s1 and s2 along the time. As in the CAR case, Eq. (88) exhibit a quadratic relation between the generalized coordinates of the system, meaning that the computational advantages identified for the CAR condition are maintained. Moreover, as the second term of Eq. (88) is independent of the generalized coordinates, the contributions to the Jacobian matrix are, in fact, equal to the ones presented in Eq. (66). In contrast, the contributions to the vectors ν and γ become dependent on the angular velocity (˙ θ* s1s2) and angular acceleration (¨ θ* s1s2) of the prescribed angle θ* s1s2, as presented hereafter ν AR = − sin(θ* s1s2(t))˙ θ* s1s2(t)(89) γAR = − 2˙ qT iCs1s2˙ qj−(cos(θ* s1s2(t))(˙ θ* s1s2(t))2+sin(θ* s1s2(t))¨ θ* s1s2(t))(90) Equation (88) describes one constraint equation, which implies that each AR kinematic driver guides one DoF of the multibody system. Hence, one AR constraint equation needs to be added for each rotational DoF of the joint being guided to allow for its full kinematic description. It is important to note that, due to the use of the dot product, the AR condition presents limitations when driving angles near 0 and π rad. For these particular cases, the vectors become aligned, resulting in the appearance of linear dependent lines in the Jacobian matrix. Moreover, since the codomain of the cosine function is positive to values in the 1 st and 4 th quadrants and negative in the 2 nd and 3 rd quadrants, the AR constraint can generate kinematically consistent positions, which do not match the prescribed movement. In order to overcome the addressed issues, the AR kinematic driver should be defined in such a way that it only drives angles between 0 and π rad or between π and 2 π . An alternative approach to avoid the issues related to the domain of the dot product equation is to introduce additional driving constraints. Instead of relying on a single constraint equation, this method employs a set of AR equations with different vector combinations to describe the same DoF. This approach allows for expanding the range of applicability of the AR kinematic driver constraint to cover the entire RoM of the joint. However, this modeling strategy introduces redundant constraints, requiring numerical methods specifically designed to handle such systems. For inverse dynamics analysis, the kinematic problem can be efficiently solved using the Newton-Raphson method with a least squares approach, without significantly Fig. 7. Schematic representation of the AR kinematic driver defined between rigid body i and j. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 19
impacting computational performance [35,51]. In contrast, forward dynamic analyses of guided systems may experience reduced computational efficiency due to the increased problem dimensions. Additionally, the presence of redundant constraints requires specialized methodologies to find a solution, such as those discussed in Section 2.1.2. A potential alternative to implementing these methodologies is to divide the simulation into multiple subproblems with shorter time intervals, ensuring that each subproblem employs a set of vectors that satisfies the AR condition domain. Equations (88) to (90) can be further simplified, if the AR driver constraint is defined between two rigid body vectors. Accordingly, let one consider two rigid body vectors u and v belonging to bodies i and j, the AR driving constraint can be described using the dot product as follows ΦAR(qi,qj)=uiTvj−cos(θ* uivj(t))=0 (91) As in the generic case, the second term of the Eq. (91) does not have an explicit dependency on the generalized coordinates of the system. Hence, the contributions to the Jacobian matrix are equal to the ones presented in Eq. (70). In turn, the contributions to the right-hand side vector of velocity are the same as the ones presented in the generic case (see Eq. (89)), while the contributions to the vector γ include the dot product of the velocity of vectors ui and vj, as expressed by Eq. (71), and the rheonomic terms presented on the generic case (see Eq. (90)). 3.3.6. Angular relation between two vectors belonging to reduced rigid bodies The AR kinematic drive for a reduced body uses the dot product definition to describe the angular relation between two generic vectors. However, the method to compute the fully-defined equivalent vectors needs to be applied first to allow for the calculation of the components of the generic vectors. Consequently, the AR kinematic driver in its homogenous form can be stated as ΦAR(qi,qj)=(Cs1 iSViqi)T(Cs2 jSVjqj)−cos(θ* uivj(t))=0 ΦAR =qT iSViCs1s2SVjqj−cos(θ* uivj(t))(92) As the rheonomic term is independent of the generalized coordinates, the contributions to the Jacobian matrix are equal to those obtained for the CAR condition and given by Eq. (73). Since this term is equal to the fully-defined case, the contributions to vector ν are the same as the ones presented in Eq. (89). Finally, the contributions to the right-hand side vector of acceleration include the quadratic terms dependent on the velocity and position of the system as presented in Eq. (75) and the rheonomic terms given by Eq. (90). 3.3.7. Angular relation between one vector belonging to a fully-defined rigid body and one external vector An additional angular relation driver to prescribe the relative motion of one vector belonging to a rigid body and an external vector is formulated. This condition can be used to guide the angular DoFs associated to the kinematic joints defined between a given body and an external element, such as the pinned joints introduced in Section 3.4. In terms of mathematical definition, this type of constraint can be formulated considering the dot product expression as presented for the generic case with two vectors (see Sections 3.3.5 and 3.3.6). Hence, let one consider a generic unit vector s1 belonging to a fully-defined rigid body i and an external unit vector s* 2, the AR* constraint driver in its homogeneous form can be established as ΦAR*(qi,t) = sT 1s* 2−cos(θ* s1s* 2(t))=(Cs1 iqi)Ts* 2−cos(θ* s1s* 2(t))=0 (93) where θ* s1s* 2 is the prescribed angle between vector s1 and s* 2 along the period of analysis. Since only one of the vectors presents a dependency on the generalized coordinates of the system, the contributions to the Jacobian matrix include only the terms relative to its derivatives with respect to the q, yielding ΦAR* q(qi,t) = {s* 2 TCs1 i ⏞⏟⏟⏞ ΦAR* qi}(1×12) (94) When the vector to guide is a rigid body vector, Eqs. (93) and (94) simplifies, resulting in ΦAR*(qi,t) = uT is* 2−cos(θ* uis* 2(t))=0 (95) ΦAR* q(qi,t) = {s* 2 T ⏞⏟⏟⏞ ΦAR* qi}(1×3) (96) where ui represents the vector u of rigid body i. The Jacobian contributions associated with Eqs. (94) and (96) do not present an explicit dependency on the generalized coordinates, which means that its value can be determined a priori if the time vector for the analysis is known. This independency of the q vector implies that the contributions to the right-hand side vector of velocity ( ν AR*) and acceleration (γAR*) have only the terms relative to the rheonomic term, as given by Eqs. (89) and (90). S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 20
3.3.8. Angular relation between one vector belonging to a reduced rigid body and one external vector In the presence of a reduced rigid body definition, the AR* constraint driver can be established in a similar manner described above using the fully-defined equivalent vector, that is ΦAR*(qi,t) = (Cs1 iSViqi)Ts* 2−cos(θ* s1s* 2(t))=0 (97) with ΦAR* q(qi,t) = ⎧ ⎪ ⎨ ⎪ ⎩s*T 2Cs1 iVi ⏞⏟⏟⏞ ΦAR* qi⎫ ⎪ ⎬ ⎪ ⎭(1×9) (98) As matrix Vi depends on the generalized coordinates, the contribution to the Jacobian matrix cannot be determined a priori as in the fully-defined case. This linear relation on the generalized coordinates implies that additional terms dependent on the generalized velocities of body i need to be included in the right-hand side vector of accelerations besides the rheonomic terms presented above, that is γAR*= − (s* 2 TCs1 i ˙ Vi ˙ qi)−(cos(θ* s1s* 2(t))(˙ θ* s1s* 2(t))2+sin(θ* s1s* 2(t))¨ θ* s1s* 2(t))(99) 3.4. Kinematic constraints with respect to the global reference frame The full definition of the multibody mechanical system may require additional constraints to define the topologic relations between the bodies and the global reference frame. These constraints allow for the establishment of specific joints, such as fixed or pinned joints, that relate the generalized coordinate of the system to external elements. Two main modeling approaches can be adopted to define these topologic relations. The first one utilizes the definition of ground bodies, that is, virtual massless bodies that are used to constrain the system bodies. This approach has the main advantage of using the same kinematic constraint equations that are utilized to relate two rigid bodies (see Section 3.2), simplifying the modeling procedure. However, the definition of the ground bodies requires the incorporation of extra coordinates to the vector of generalized coordinates (q), increasing the complexity of the analysis. The second approach employs specific kinematic constraints to describe the geometric relations between points and vectors of the model and external elements. Consequently, no additional coordinates need to be added to the vector of generalized coordinates, being these relations fully-described by the use of algebraic equations. This approach requires defining specific equations based on the mechanical elements being modeled. Mathematically, these constraints are similar to those used to constrain or drive two rigid bodies (see Sections 3.2 and 3.3). As the constraint equations depend only on the generalized coordinates of one body, the contributions to the Jacobian matrix and right-hand side vector of velocity and acceleration include only the terms relative to that body. Therefore, the topological constraints with respect to global reference frame are not detailed in this work, however, the interested reader is referred to [30] for a derivation of this type of kinematic constraints. 4. Dynamics of 3D multibody systems with FCC-GRB formulation The inverse and forward dynamic analysis of multibody systems requires the assembly of the EoM as described by Eq. (7) or Eq. (9), respectively. In both cases, beyond describing the system topology and prescribing the guided DoFs via the kinematic constraint equations presented in Section 3, the implementation of a multibody formulation also requires the definition of the system mass matrix and the development of methods for applying external forces and moments of force. Hence, this section explores the mathematical differences between using a fully-defined versus a reduced definition of a rigid body in the dynamic analysis of mechanical systems with FCC, presenting the main equations needed to evaluate inertial and external forces. 4.1. Mass matrix for a fully-defined generic rigid body One of the distinctive features, which differentiates the FCC-GRB formulation from the natural coordinates formulation, is the simplicity in the definition of the system mass matrix. This straightforwardness is a direct result of the kinematic structure adopted for the generic rigid body, which generates sparse and uncoupled mass matrices [30]. The rigid body mass matrix can be defined using the principle of the virtual power, as described in Jal´ on and Bayo [29]. In fact, the derivation of the mass matrix for a spatial rigid body follows the same approach applied in the 2D formulation [30] and, for that reason, this work focuses on presenting only the main differences. The interested reader is referred to these two works for a detailed derivation of the rigid body mass matrices. By applying the principle of the virtual power and the method for describing the kinematics of a generic point P (see Section 3.1.1), it is possible to obtain an expression that relates the mass matrix for a given rigid body i (Mi) and the constant matrix CP i for a generic point P [29,30,53] S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 21
Mi= ρ i∫ Ωi CP i TCP idΩ(100) where ρ i and Ωi are respectively the mass density and the geometric domain of the rigid body i. By solving the integral presented in Eq. (100), and considering that the reference point is strategically located at the CoM and the rigid bodies vectors are aligned with the inertia principal axes of the rigid body being modeled, the final form of the mass matrix for a generic rigid body is obtained as [29,53] Mi=⎡ ⎢ ⎢ ⎣ miI3030303 03IuiI30303 0303IviI303 030303IwiI3 ⎤ ⎥ ⎥ ⎦(12×12) (101) with Iui=1 2(I ηη i+Iζζi−Iξξi) Ivi=1 2(Iξξi+Iζζi−I ηη i) Iwi=1 2(Iξξi+I ηη i−Iζζi) (102) where mi is the mass of rigid body i and Iii the components of its inertia tensor. It is important to note that by locating the reference point at the CoM and aligning the vectors with the principal axes of inertia, the mass matrix of a generic rigid body becomes diagonal, with its entries equal to the mass and the inertia relations along the rigid body vectors (Iui, Ivi, Iwi). If the reference point is located at another point or the vectors are not aligned, the off-diagonal entries of the mass matrix will present a dependency on the products of inertia. From a modeling approach, the use of a generic rigid body approximates the definition of the mass matrix from the one utilized in the Cartesian coordinates, i.e., the matrix becomes constant and independent of the topology of the system. Moreover, its entries present a higher physical meaning than the ones obtained in the natural coordinates formulation, as they directly represent the inertial properties of the body. As mentioned above, the possibility of sharing vectors between bodies implies that the FCC formulation with a generic rigid body allows for both an explicit and implicit definition of some kinematic joints. Considering a full explicit model, the assembly of the system mass matrix is a straightforward process, and its entries are directly the rigid body mass matrices. Considering the vector of the generalized coordinates presented in Eq. (1), the global mass matrix of the system will be constant and diagonal, as shown below M=⎡ ⎢ ⎢ ⎢ ⎢ ⎣ M1 ⋱ Mi ⋱ Mnb ⎤ ⎥ ⎥ ⎥ ⎥ ⎦(nq×nq) (103) Despite generating more coordinates, this modeling approach has the main advantage of generating sparser and uncoupled mass matrices, which is a computational advantage if the proper numerical methods are used. In turn, the use of an implicit modeling approach results in coupled mass matrices, and consequently, less sparse matrices. Thus, let one consider a vector u shared by the generic rigid bodies i and j, the contributions of each body to the global mass matrix are given by M= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ M1 ⋱ MrOiMui+MujMviMwiMrOjMvjMwj ⋱ Mnb ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦(nq×nq) (104) It is important to note that in both cases, as the rigid body mass matrices are independent of the generalized coordinates, the mass matrix of the system is constant. This fact implies that matrix M only needs to be evaluated once at the beginning of the analysis, reducing the number of calculations required for each evaluation of the EoM. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 22
4.2. Mass matrix for a reduced generic rigid body The mass matrix given by Eq. (101) describes mathematically the inertial properties of the generic rigid body i, namely the mass and the mass moments of inertia around the three axes that constitute the local reference frame, which, accordingly with the kinematic structure adopted for the generic rigid body, are defined by the three rigid body vectors. Therefore, when in the presence of a reduced rigid body, the mass matrix needs to also include the inertia moment around a third virtual vector w to fully describe the inertial properties of the body. Once more, the procedure presented in Section 2.2.3 can be utilized to determine the respective fully-defined equivalent vector, and then define the mass matrix for the reduced rigid body i. As for the fully-defined case, the principle of virtual power can be used to determine the equivalent mass matrix. Thus, let one consider the fully-defined equivalent vector of the reduced rigid body i, the virtual power generated within this body by the inertial forces (W* i) is given by W* i= − ρ ∫˙ q*T 3viCP i TCP i ¨ q3vidΩ= − ρ ∫˙ q*T iVT iCP i TCP i(˙ Vi ˙ qi+Vi ¨ qi)dΩ(105) where ˙ q* i represents the generalized virtual velocities of the rigid body i and ˙ q* 3vi denotes the virtual velocities described in terms of the fully-defined equivalent form. As the transformations matrices Vi and ˙ Vi are independent of the volume of body i, these can be moved outside the integral. In this way, the integral in Eq. (105) is equal to the one in Eq. (100), which in turn is equal to the mass matrix of a fully-defined rigid body. Hence, Eq. (105) can be further simplified as W* i= − ˙ q*T iVT i ρ ⎛ ⎝∫Ω CP i TCP idΩ⎞ ⎠(˙ Vi ˙ qi+Vi ¨ qi) W* i= − ˙ q*T i(VT iMiVi ¨ qi+VT iMi ˙ Vi ˙ qi)(106) Equation (106) can be divided into two main terms. The first one, dependent on the generalized accelerations of the system, corresponds to the equivalent mass matrix of the reduced rigid body (M2vi) as M2vi=[VT iMiVi](9×9)(107) Matrix M2vi can be assembled in the global mass matrix of the system, following the same approach presented in Eqs. (103) or (104). However, since matrix Vi is dependent on the generalized coordinates, the global mass matrix becomes dependent on the state of the system. This fact implies that for each evaluation of the EoM, the global mass matrix needs to be updated with the contributions relative to the reduced rigid bodies. The second term of Eq. (106) can be treated as a velocity-dependent generalized inertial force (g2vi), which needs to be added to the vector of the generalized forces, such that g2vi=VT iMi ˙ Vi ˙ qi(108) The velocity-dependent inertial force is dependent on the state of the system, meaning that the vector of the generalized force vectors needs to be also updated for each time step. 4.3. Application of external forces and moments of force Since the external forces and external moments of force are not necessarily applied in the generalized coordinates of the system, those must be transformed into equivalent generalized forces. It is important to note that the methodology proposed hereafter is generic and can be applied to any external force or moment applied to the system. This method is also valid for state-dependent forces or moments that can be represented by an equivalent external force or moment of force. However, in this particular case, the state of the system needs to be firstly determined to allow for the calculation of the magnitude, orientation and application point of these forces. 4.3.1. External forces for a fully-defined rigid body The principle of virtual power can be utilized to derive the equivalent generalized force of any external force applied in the system. This principle states that the virtual power produced by an external force f applied in point P belonging to body i is equal to the virtual power produced by the generalized equivalent force, i.e., the product of the generalized equivalent force gf i by the vector of virtual velocities of rigid body i [29,53], as presented below ˙ r* P Tf=˙ q* i Tgf i(109) where ˙ r* P represents the vector of the virtual velocities of point P. By applying the methodology presented in Section 3.1.1, it is possible to express the algebraic relation presented in Eq. (109) as a function of the generalized coordinates of the system, such that S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 23
˙ q* i TCP i Tf=˙ q* i Tgf i⇒gf i=CP i Tf(110) Equation (110) shows that the equivalent generalized force can be expressed as the product of the vector of the external force f and the constant transformation matrix C of the point P described in relation to body i (see Fig. 8a). This relation implies that if both the external force f and its application point are constant, as in the case of the conservative gravitational force, the equivalent generalized force is also constant, meaning that its contribution to the generalized force vector only needs to be computed once at the beginning of the analysis. Similarly, if only the application point of force f is constant in the local reference of the body (e.g., forces generated by linear springs, dampers, or actuator forces with fixed attachment points), matrix C becomes constant along the period of analysis, implying that it can also be defined a priori. 4.3.2. External moment of force for a fully-defined rigid body As the FCC-GRB formulation does not utilize angular coordinates to describe the rigid bodies orientations, the external moments of force cannot be directly applied to the system as a generalized moment. Thus, each external moment of force needs to be first converted into the equivalent forces couple, which is then applied to the system using the same approach presented in the previous section. Accordingly, let one consider the external moment of force τ applied to body i, and the equivalent force couple f τ and f− τ , such that { τ =f τ ×b τ +f− τ ×b− τ f τ +f− τ =0(111) where b τ and b− τ are respectively the moment arms of force vectors f τ and f− τ . As the virtual power produced by the generalized equivalent moment of force (g τ i) is equal to the sum of the virtual power produced by the equivalent force couple [29,53], the following expression yields ˙ q* i Tg τ i=˙ r* Pf τ Tf τ −˙ r* Pf− τ Tf− τ (112) Defining the position of point Pf− τ at the origin of the local reference frame of body i and of point Pf τ at the tip of vector s (see Fig. 8b) and recalling Eqs. (36) and (40), Eq. (112) can be written in terms of the generalized coordinates of the system as ˙ q* i Tg τ i=˙ q* i T(CP i T−COi i T)f τ ⇒g τ i=Cs i Tf τ (113) By locating the point Pf− τ at the origin of the reference, the contribution of the force f− τ to the moment is null, being only responsible by counteracting the translation produced by the force f τ . Therefore, the vector s should be defined in such a way that it is contained in the plane normal to vector τ and the norm of the product Cs iTf is equal to the magnitude of the moment of force τ . It must be noted that this linear dependency on the matrix C implies that if the direction of the external moment of force τ does not change along the time, the matrix Cs i is constant and consequently only the product presented in Eq. (113) needs to be computed each time the EoM are evaluated. 4.3.3. External forces for a reduced rigid body The application of the external forces for a rigid body defined in the reduced form follows the same approach presented for the fully-defined rigid body. However, the expressions presented in Section 4.3.1 need to be revised considering the fully-defined equivalent vector (q3v) introduced in Section 2.2.3. Hence, and recalling the principle of virtual power for an external force f (see Fig. 8. a) Application of an external force f at point P belonging to the generic rigid body i; b) Application of an external moment of force τ to the generic rigid body i and respective transformation into the equivalent force couple f τ and f− τ . S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 24
Eq. (109)), the generalized equivalent force for a reduced rigid body can be described as ˙ q* i TVT iCP i Tf=˙ q* i Tgf i⇒gf i=VT iCP i Tf(114) In contrast to the fully-defined case, the generalized equivalent force for the reduced rigid body depends on the generalized coordinates of the system, meaning that each time the EoM are evaluated, the generalized force vector needs to be updated, even when constant forces are applied. 4.3.4. External moment of force for a reduced rigid body The same principles presented in Section 4.3.2 can be applied to determine the generalized moment of force for an external moment τ . However, the fully-defined equivalent vector needs to be used to allow for the systematic calculation of the contributions to the generalized force vector (g τ i). ˙ q* i Tg τ i=˙ q* i TVT i(CP i T−COi i T)f τ ⇒g τ i=VT iCs i Tf τ (115) It should be noticed that Eq. (115) presents an explicit dependency on the generalized coordinates of the system, which implies that, even for constant moments of force, the contributions to vector g need to be calculated for each instant of time. Nevertheless, as in this particular case, the vector s is constant in the local reference frame of the body i, matrix Cs i is also constant and can be defined a priori during the pre-analysis procedure. 4.4. Internal reaction forces and moments of force for a fully-defined rigid body One of the major objectives, when performing a dynamic analysis of mechanical systems, is the determination of the internal forces and moments of force generated during the motion of the system. From a mechanical standpoint, these represent the forces and moments that need to be generated within the system to comply with the imposed kinematic constraints. Their physical meaning is related to the type of kinematic constraint they are obtained. Algebraically, the internal forces of the system (gΦ) can be calculated by solving Eqs. (7) or (9) in order to the internal forces, such that gΦ={gΦ 1 T⋯gΦ i T⋯gΦ nb T}T=ΦT qλ(116) with gΦ i=gΦRB i+gΦKJ i+gΦDC i(117) where gΦ i represents the vector of the internal forces acting on body i, which, in turn, is equal to the sum of all contributions from kinematic constraints and drivers related to this body, including rigid body topological constraints (RB), kinematic joint topological constraints (KJ), and driving constraints (DC). The physical meaning of Eqs. (116) and (117) can be better interpreted if the contributions of each kinematic driver are evaluated separately. In fact, the Lagrangian multipliers associated with each kinematic constraint represent the magnitude of the internal forces, while the respective lines of the Jacobian matrix their direction. This relevant characteristic of the Lagrange multipliers is of particular relevance, as it allows for the determination of the joint reaction forces and moments directly from the EoM. Taking as an example the case of a spherical joint, defined by the CP condition (see Sections 3.2.1 and 3.2.2), the generalized internal forces generated by this joint are given by gΦCP =ΦCPT qλCP =⎡ ⎣CP i T −CP j T⎤ ⎦λCP =⎧ ⎨ ⎩ CP i TλCP CP j T(−λCP)⎫ ⎬ ⎭(24×1) (118) By comparing Eq. (118) with Eq. (110) and having in mind that λCP is a 3D vector, it can be observed that it expresses the application of two concentrated forces of equal magnitude and opposite direction (λCP ,−λCP) at point P. These two forces represent the action-reaction force couple generated at the joint. Hence, the joint reaction forces (fCP R) are obtained directly from the EoM, as they are equal to the Lagrangian multipliers associated with the CP condition as fCP R=fCP Ri= − fCP Rj={λCP}(3×1)(119) In turn, the internal forces generated by the AR kinematic drivers describe the equivalent forces that produce the internal driving moments of force. Therefore, the determination of the moments of force at the joints in an FCC formulation is not a direct procedure, requiring additional calculations. Accordingly, let one consider an angular driving constraint between vectors s1 and s2 belonging to bodies i and j, and recalling the contributions to the Jacobian matrix presented in Sections 3.3.5 and 3.3.6, the generalized internal forces for an AR driver constraint (gΦAR ) can be given by S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 25
rigid body vectors with joint axes, simplifying the definition of system topology and enabling the generation of more efficient mathematical expressions [33]. Additionally, the constant nature of matrix C for points and vectors fixed in the local reference frame allows it to be computed in a pre-kinematic step and assembled into the Jacobian matrix without requiring continuous updates during EoM evaluations. Nevertheless, the FCC-GRB formulation is not without limitations. Compared to both the Cartesian coordinates and natural coordinates formulations, the FCC-GRB typically generates a larger number of generalized coordinates and kinematic constraint equations. Specifically, even when using the reduced approach, the number of generalized coordinates is larger than in the classical angular-based global formulations. For instance, when compared to a Cartesian coordinates formulation with Euler parameters, the FCC-GRB formulation requires approximately two additional generalized coordinates, and, consequently, two additional kinematic constraints, per rigid body [28]. On the other hand, while the FCC-GRB generates the same number of generalized coordinates per rigid body as the natural coordinates, the ability of this formulation to share points and vectors between adjacent bodies can reduce the total number of coordinates and constraints in the system [29,53]. However, the FCC-GRB also supports vector sharing, which implicitly defines shared vector SV conditions. Although not explored in this study, this capability could further reduce the total number of generalized coordinates and kinematic constraints. For highly constrained systems with a large number of revolute joints, the FCC-GRB may result in a lower total coordinate count compared to the Cartesian coordinates formulation. Finally, one potential application identified for the planar FCC-GRB formulation is its use in teaching multibody dynamics topics in advanced courses [30]. The main features supporting this idea were its greater intuitiveness and simplicity in implementing and modeling multibody systems compared to other common global formulations. These premises are still valid for the spatial formulation with the fully-defined approach, namely: i) the introduction of the generic rigid body simplifies the modeling of the multibody system, allowing for easier systematization of this procedure. This is particularly more noticeable when compared with the natural coordinates formulation, where the model is strongly dependent on the topology of the system; ii) the rigid body mass matrix remains diagonal and constant, with its structure independent of the system topology. Moreover, its entries are directly the inertial properties of the body being modeled, resulting in a mass matrix with a higher physical meaning; iii) the rigid bodies are defined resorting only to Cartesian coordinates, meaning that no background knowledge on the parametrization of rotations in spatial models is required, simplifying the implementation of the formulation. This last point was precisely one of the major advantages attributed to the natural coordinates formulation, supporting its use in educational applications [33]. Considering that the FCC-GRB formulation is even more intuitive, it can be more easily applied in the teaching of multibody topics in advanced courses. For that reason, the authors opted for presenting in detail the specificities of the formulation in 3D, so that, together with the planar work [30], the undergraduate and graduate students can easily implement the formulation and understand the results it provides. Hence, these two works can be directly used as an education tool to support courses in the STEM areas, providing insights of how to model and analyze simple and complex mechanical systems in 3D. CRediT authorship contribution statement S´ ergio B. Gonçalves: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization. Ivo Roupa: Writing – review & editing, Validation, Conceptualization. Paulo Flores: Writing – review & editing, Supervision. Miguel Tavares da Silva: Writing – review & editing, Validation, Supervision, Project administration, Conceptualization. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: S´ ergio B. Goncalves, Ivo Roupa, Paulo Flores, and Miguel Tavares da Silva reports financial support was provided by Fundaç˜ ao para a Ciencia e a Tecnologia (FCT). If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors acknowledge Fundaç˜ ao para a Ciˆ encia e a Tecnologia (FCT) for its financial support via the projects LAETA Base Funding (DOI: 10.54499/UIDB/50022/2020), LAETA Programmatic Funding (DOI: 10.54499/UIDP/50022/2020), UIDB/04436/ 2020 and UIDP/04436/2020, and Portuguese Recovery and Resilience Program (PRR) for its financial support via IAPMEI/ANI/FCT under Agenda C645022399-00000057 (eGamesLab). Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.mechmachtheory.2025. 105955. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 32
Data availability No data was used for the research described in the article. References [1] W. Schiehlen, Research trends in multibody system dynamics, Multibody Syst. Dyn. 18 (2007) 3–13, https://doi.org/10.1007/S11044-007-9064-4. [2] S. Bruni, J.P. Meijaard, G. Rill, A.L. Schwab, State-of-the-art and challenges of railway and road vehicle dynamics with multibody dynamics approaches, Author (s) (2020), https://doi.org/10.1007/s11044-020-09735-z. [3] J.N. Costa, P. Antunes, H. Magalh˜ aes, J. Pombo, J. Ambr´ osio, A novel methodology to automatically include general track flexibility in railway vehicle dynamic analyses, Proc. Inst. Mech. Eng. F. J. Rail. Rapid. Transit. 235 (2021) 478–493, https://doi.org/10.1177/0954409720945420. [4] D. Apeng, L. Shu, Z. Wenguo, Multi-body coupling dynamic research of carrier-based aircraft catapult launch based on natural coordinate method, Proc. Inst. Mech. Eng., Part K: J. Multi-Body Dyn. 233 (2019) 195–207, https://doi.org/10.1177/1464419318785978. [5] J. Coelho, B. Dias, G. Lopes, F. Ribeiro, P. Flores, Development and implementation of a new approach for posture control of a hexapod robot to walk in irregular terrains, Robotica 42 (2024) 792–816, https://doi.org/10.1017/S0263574723001765. [6] J.C. Samin, O. Brüls, J.F. Collard, L. Sass, P. Fisette, Multiphysics modeling and optimization of mechatronic multibody systems, Multibody Syst. Dyn. 18 (2007) 345–373, https://doi.org/10.1007/S11044-007-9076-0. [7] N. Li, F. Li, H. Yang, H. Peng, Real-time control of a soft manipulator based on reduced order extended position-based dynamics, Mech. Mach. Theory. 202 (2024) 105774, https://doi.org/10.1016/J.MECHMACHTHEORY.2024.105774. [8] V. Pozhbelko, A unified structure theory of multibody open-, closed-, and mixed-loop mechanical systems with simple and multiple joint kinematic chains, Mech. Mach. Theory. 100 (2016) 1–16, https://doi.org/10.1016/J.MECHMACHTHEORY.2016.01.001. [9] D. Yuan, X. Sun, L. Hu, Q. Peng, X. Chen, Y. Li, S. Huang, L. Zhao, B. Li, Coupled dynamics modeling and analysis of a coring drilling equipment for hard-rock tunnel boring, Lecture Notes Mech. Eng. (2024) 3015–3034, https://doi.org/10.1007/978-981-99-8048-2_206. [10] S.B. Gonçalves, I. Roupa, M. Tavares da Silva, On the analysis of fully cartesian coordinates – a comparison between a reduced and full-defined modelling approach in spatial mechanisms, in: ECCOMAS Thematic Conference on Multibody Dynamics 2023, Lisbon, Portugal, 2023. [11] P. Flores, J. Ambr´ osio, J.P. Claro, Dynamic analysis for planar multibody mechanical systems with lubricated joints, Multibody Syst. Dyn. 12 (2004) 47–74, https://doi.org/10.1023/B:MUBO.0000042901.74498.3A. [12] M. Da Lio, V. Cossalter, R. Lot, On the use of natural coordinates in optimal synthesis of mechanisms, Mech. Mach. Theory. 35 (2000) 1367–1389, https://doi. org/10.1016/S0094-114X(00)00006-9. [13] I. Roupa, M.R. da Silva, F. Marques, S.B. Gonçalves, P. Flores, M.T. da Silva, On the modeling of biomechanical systems for human movement analysis: a narrative review, Arch. Comput. Methods Eng. 29 (2022) 4915–4958, https://doi.org/10.1007/s11831-022-09757-0. [14] K.A. Inkol, J. McPhee, Assessing control of fixed-support balance recovery in wearable lower-limb exoskeletons using multibody dynamic modelling, in: Proceedings of the IEEE RAS and EMBS International Conference on Biomedical Robotics and Biomechatronics 2020-November, 2020, pp. 54–60, https://doi. org/10.1109/BIOROB49111.2020.9224430. [15] C. Quental, M. Azevedo, J. Ambr´ osio, S.B. Gonçalves, J. Folgado, Influence of the musculotendon dynamics on the muscle force-sharing problem of the shoulder—A fully inverse dynamics approach, J. Biomech. Eng. 140 (2018), https://doi.org/10.1115/1.4039675. [16] S.B. Gonçalves, P. Flores, M.T. da Silva, On the use of mixed coordinates for simultaneous determination of joint angles and kinematic consistent positions, in: MMT Symposium, Guimar˜ aes, Portugal, 2024. [17] A.R.C. Oliveira, S.B. Gonçalves, M.A. De Carvalho, M.T.d. Silva, Development of a musculotendon model within the framework of multibody systems dynamics, Comput. Methods Appl. Sci. 42 (2016), https://doi.org/10.1007/978-3-319-30614-8_10. [18] T.M. Malaquias, S.B. Gonçalves, M.T. da Silva, A three-dimensional multibody model of the human ankle-foot complex, Mech. Mach. Sci. 24 (2015) 445–453, https://doi.org/10.1007/978-3-319-09411-3_47. [19] N. Monteiro, M.T. da Silva, J. Folgado, J. Melancia, Structural analysis of the intervertebral discs adjacent to an interbody fusion using multibody dynamics and finite element cosimulation, Multibody Syst. Dyn. 25 (2010) 245–270, https://doi.org/10.1007/S11044-010-9226-7, 2010 25:2. [20] M. Busch, B. Schweizer, Coupled simulation of multibody and finite element systems: an efficient and robust semi-implicit coupling approach, Arch. Appl. Mech. 82 (2012) 723–741, https://doi.org/10.1007/S00419-011-0586-0. [21] F. Guedes de Melo, S.B. Gonçalves, P. Areias, M.T. Silva, Analysis of the foot-ground contact using an MSD-FEM Co-simulation approach, in: Giulio Rosati, Alessandro Gasparetto, Marco Ceccarelli (Eds.), New Trends in Mechanism and Machine Science: Proceedings of EuCoMeS 2024, Springer Cham, 2024, pp. 54–62, https://doi.org/10.1007/978-3-031-67295-8_7. [22] C.D. Twigg, D.L. James, Many-worlds browsing for control of multibody dynamics, in: Proceedings of the ACM SIGGRAPH Conference on Computer Graphics, 2007, https://doi.org/10.1145/1275808.1276395. [23] A. Parra, A.J. Rodriguez, A. Zubizarreta, J. Perez, Validation of a real-time capable multibody vehicle dynamics formulation for automotive testing frameworks based on simulation, IEEE Access. 8 (2020) 213253–213265, https://doi.org/10.1109/ACCESS.2020.3040232. [24] D. Negrut, A. Tasora, M. Anitescu, H. Mazhar, T. Heyn, A. Pazouki, Solving large multibody dynamics problems on the GPU, GPU Comput. Gems Jade Ed. (2012) 269–280, https://doi.org/10.1016/B978-0-12-385963-1.00020-4. [25] J. Cuadrado, D. Dopico, M. Gonzalez, M.A. Naya, A combined penalty and recursive real-time formulation for multibody dynamics, J. Mech. Design, Trans. ASME 126 (2004) 602–608, https://doi.org/10.1115/1.1758257. [26] W. Schiehlen, Multibody system dynamics: roots and perspectives, Multibody Syst. Dyn. 1 (1997) 149–188, https://doi.org/10.1023/A:1009745432698. [27] J. Cuadrado, D. Dopico, M.A. Naya, M. Gonzalez, Real-time multibody dynamics and applications, CISM Int. Centre Mech. Sci., Courses Lectures 507 (2009) 247–311, https://doi.org/10.1007/978-3-211-89548-1_6. [28] P. Nikravesh, Computer-aided Analysis of Mechanical systems,, 1st ed, Prentice Hall, New Jersey, 1988, p. 07632. [29] J.G. de Jalon, E. Bayo, Kinematic and Dynamic Simulation of Multibody Systems: the Real-Time Challenge, Springer Verlag, New York, 1994. [30] I. Roupa, S.B. Gonçalves, M.T. da Silva, Kinematics and dynamics of planar multibody systems with fully Cartesian coordinates and a generic rigid body, Mech. Mach. Theory. 180 (2023) 105–134, https://doi.org/10.1016/j.mechmachtheory.2022.105134. [31] P.E. Nikravesh, An overview of several formulations for multibody dynamics, in: D. Talabua, T. Roche (Eds.), Product Engineering: Eco-Design, Technologies and Green Energy, Springer, Netherlands, Dordrecht, 2005, pp. 189–226, https://doi.org/10.1007/1-4020-2933-0_13. [32] F. Marques, I. Roupa, M.T. Silva, P. Flores, H.M. Lankarani, Examination and comparison of different methods to model closed loop kinematic chains using lagrangian formulation with cut joint, clearance joint constraint and elastic joint approaches, Mech. Mach. Theory. 160 (2021) 104294, https://doi.org/ 10.1016/j.mechmachtheory.2021.104294. [33] J.G. Jal´ on, Twenty-five years of natural coordinates, Multibody Syst. Dyn. 18 (2007) 15–33. [34] J.G. De Jal´ on, A. Callejo, A straight methodology to include multibody dynamics in graduate and undergraduate subjects, Mech. Mach. Theory. 46 (2011) 168–182, https://doi.org/10.1016/j.mechmachtheory.2010.09.008. [35] J. García de Jal´ on, M.D. Guti´ errez-L´ opez, Multibody dynamics with redundant constraints and singular mass matrix: existence, uniqueness, and determination of solutions for accelerations and constraint forces, Multibody Syst. Dyn. 30 (2013) 311–341, https://doi.org/10.1007/s11044-013-9358-7. [36] J. Cuadrado, J. Cardenal, E. Bayo, Modeling and solution methods for efficient real-time simulation of multibody dynamics, Multibody Syst. Dyn. 1 (1997) 259–280, https://doi.org/10.1023/A:1009754006096. [37] D.S. Bae, E.J. Haug, A recursive formulation for constrained mechanical system dynamics: part III, Parallel Processor Implementation, Mech. Struct. Mach. 15 (1987) 359–382, https://doi.org/10.1080/08905458808960263. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 33
[38] X. Yu, A. Mikkola, Y. Pan, J.L. Escalona, The explanation of two semi-recursive multibody methods for educational purpose, Mech. Mach. Theory. 175 (2022) 104935, https://doi.org/10.1016/j.mechmachtheory.2022.104935. [39] A.A. Shabana, Dynamics of Multibody Systems, Cambridge university press, 2020, https://doi.org/10.1017/CBO9781107337213. [40] E.J. Haug, Computer Aided Kinematics and Dynamics of Mechanical Systems, 1, Basic Methods, Allyn & Bacon, Inc., 1989. [41] J.G. De Jal´ on, J. Unda, A. Avello, Natural coordinates for the computer analysis of multibody systems, Comput. Methods Appl. Mech. Eng. 56 (1986) 309–327. [42] P.E. Nikravesh, H.A. Affifi, Construction of the equations of motion for multibody dynamics using point and joint coordinates, Comput. Aided Anal. Rigid Flexible Mech. Syst. (1994) 31–60, https://doi.org/10.1007/978-94-011-1166-9_2. [43] J.A. Ambr´ osio, M. Tavares da Silva, A biomechanical multibody model with a detailed locomotion muscle apparatus, Adv. Comput. Multibody Syst. (2005) 155–184. [44] M.T. Gameiro, P. Silva, Modelaç˜ ao e Simulaç˜ ao Sistem´ atica em Coordenadas Cartesianas Totais de Sistemas Multicorpo. Actas Do Congresso De M´ etodos Num´ ericos Em Engenharia, Junho, Porto, Portugal, 2007, p. 2007, 13-15. [45] S. Uhlar, P. Betsch, A rotationless formulation of multibody dynamics: modeling of screw joints and incorporation of control constraints, Multibody Syst. Dyn. 22 (2009) 69–95, https://doi.org/10.1007/s11044-009-9149-3. [46] C.M. Pappalardo, A natural absolute coordinate formulation for the kinematic and dynamic analysis of rigid multibody systems, Nonlinear. Dyn. 81 (2015) 1841–1869. [47] C.M. Pappalardo, D. Guida, On the Lagrange multipliers of the intrinsic constraint equations of rigid multibody mechanical systems, Arch. Appl. Mech. 88 (2018) 419–451, https://doi.org/10.1007/s00419-017-1317-y. [48] C.M. Pappalardo, D. Guida, Dynamic analysis of planar rigid multibody systems modeled using natural absolute coordinates, Appl. Comput. Mech. 12 (2018). [49] I. Roupa, S.B. Gonçalves, M. Tavares da Silva, EZMOTION – A computational tool to perform dynamic analysis of planar (Bio) mechanical systems, in: 5th Meeting of the Young Researchers of LAETA, May 5-6, Lisboa, Portugal, 2022. [50] I. Roupa, R. Peneque, S.B. Gonçalves, M. Tavares da Silva, Calculation of the reaction and muscle forces in squat and lunge exercises – comparison between a static optimization technique and a muscle reduction approach, in: DSM - 2nd Portuguese Conference on Multibody Systems Dynamics, Guimar˜ aes, Portugal, 2022. [51] I. Roupa, S.B. Gonçalves, M. Tavares da Silva, Kinematic analysis of planar biomechanical models using mixed coordinates, in: ECCOMAS Thematic Conference on Multibody Dynamics, Budapest, Hungary, 2021. [52] J.G. de Jal´ on, J. Cuadrado, A. Avello, J.M. Jimenez, Kinematic and dynamic simulation of rigid and flexible systems with fully cartesian coordinates, Comput.- Aided Anal. Rigid Flex. Mech. Syst. (1994) 285–323, https://doi.org/10.1007/978-94-011-1166-9_9. [53] M. Tavares da Silva, Human Motion Analysis Using Multibody Dynamics and Optimization Tools, Universidade T´ ecnica de Lisboa - Instituto Superior T´ ecnico, 2003. [54] M. Gonz´ alez, F. Gonz´ alez, A. Luaces, J. Cuadrado, A collaborative benchmarking framework for multibody system dynamics, Eng. Comput. 26 (2010) 1–9, https://doi.org/10.1007/s00366-009-0139-0. [55] IFToMM technical committee for multibody dynamics, library of computational benchmark problems, (2022). https://www.iftomm-multibody.org/benchmark/ browse/ (accessed January 10, 2025). [56] F. Amirouche, Fundamentals of Multibody dynamics: Theory and Applications, Springer Science & Business Media, 2007. [57] W.C. Rheinboldt, Methods for Solving Systems of Nonlinear Equations, SIAM, 1998. [58] R. Soram, S. Roy, S.R. Singh, M. Khomdram, S. Yaikhom, S. Takhellambam, On the rate of convergence of Newton-Raphson method, Int. J. Eng. Science (IJES) 2 (2013) 5–12. [59] D. Dopico, ´ A.L. Varela, A. Luaces, Augmented lagrangian index-3 semi-recursive formulations with projections: kinematics and dynamics, Multibody Syst. Dyn. 52 (2021) 377–405, https://doi.org/10.1007/s11044-020-09771-9. [60] K. Augustynek, A. Urba´ s, Numerical investigation on the constraint violation suppression methods efficiency and accuracy for dynamics of mechanisms with flexible links and friction in joints, Mech. Mach. Theory. 181 (2023) 105211, https://doi.org/10.1016/J.MECHMACHTHEORY.2022.105211. [61] P. Flores, M. Machado, E. Seabra, M. Tavares da Silva, A parametric study on the Baumgarte stabilization method for forward dynamics of constrained multibody systems, J. Comput. Nonlinear. Dyn. 6 (2011) 011019, https://doi.org/10.1115/1.4002338. [62] J. Baumgarte, Stabilization of constraints and integrals of motion in dynamical systems, Comput. Methods Appl. Mech. Eng. 1 (1972) 1–16, https://doi.org/ 10.1016/0045-7825(72)90018-7. [63] F. Marques, A.P. Souto, P. Flores, On the constraints violation in forward dynamics of multibody systems, Multibody Syst. Dyn. 39 (2017) 385–419, https://doi. org/10.1007/s11044-016-9530-y. [64] M. Khoshnazar, M. Dastranj, A. Azimi, M.M. Aghdam, P. Flores, Application of the Bezier integration technique with enhanced stability in forward dynamics of constrained multibody systems with Baumgarte stabilization method, Eng. Comput. 40 (2024) 1559–1573, https://doi.org/10.1007/S00366-023-01884-X. [65] E. Paraskevopoulos, N. Potosakis, S. Natsiavas, An augmented lagrangian formulation for the equations of motion of multibody systems subject to equality constraints, Procedia Eng. 199 (2017) 747–752, https://doi.org/10.1016/j.proeng.2017.09.037. [66] E. Bayo, A. Avello, Singularity-free augmented lagrangian algorithms for constrained multibody dynamics, Nonlinear. Dyn. 5 (1994) 209–231, https://doi.org/ 10.1007/BF00045677. [67] B. Ruzzeh, J. K¨ ovecses, A penalty formulation for dynamics analysis of redundant mechanical systems, J. Comput. Nonlinear. Dyn. 6 (2011) 1–12, https://doi. org/10.1115/1.4002510. [68] L. Yang, S. Xue, W. Yao, Application of Gauss principle of least constraint in multibody systems with redundant constraints, Proc. Inst. Mech. Eng., Part K: J. Multi-Body Dyn. 235 (2021) 150–163, https://doi.org/10.1177/1464419320975301. [69] J. Cuadrado, D. Dopico, M.A. Naya, M. Gonzalez, Penalty, semi-recursive and hybrid methods for MBS real-time dynamics in the context of structural integrators, Multibody Syst. Dyn. 12 (2004) 117–132, https://doi.org/10.1023/B:MUBO.0000044421.04658.de. [70] J.C. García Orden, S. Conde Martín, Controllable velocity projection for constraint stabilization in multibody dynamics, Nonlinear. Dyn. 68 (2012) 245–257, https://doi.org/10.1007/s11071-011-0224-y. [71] S.S. Kim, M.J. Vanderploeg, A general and efficient method for dynamic analysis of mechanical systems using velocity transformations, J. Mech. Design, Trans. ASME 108 (1986) 176–182, https://doi.org/10.1115/1.3260799. [72] A. Avello, J.M. Jim´ enez, E. Bayo, J.G. de Jal´ on, A simple and highly parallelizable method for real-time dynamic simulation based on velocity transformations, Comput. Methods Appl. Mech. Eng. 107 (1993) 313–339, https://doi.org/10.1016/0045-7825(93)90072-6. [73] D. Dopico, F. Gonz´ alez, J. Cuadrado, J. K¨ ovecses, Determination of holonomic and nonholonomic constraint reactions in an index-3 augmented lagrangian formulation with velocity and acceleration projections, J. Comput. Nonlinear. Dyn. 9 (2014), https://doi.org/10.1115/1.4027671. [74] C.W. Walton, E.C. Steeves, New matrix theorem and its application for establishing independent coordinates for complex dynamical systems with constraints, NASA-Tech Rep. R-326 (1969). [75] N.K. Mani, E.J. Haug, K.E. Atkinson, Application of singular value decomposition for analysis of mechanical system dynamics, J Mech. Trans. Automation Design 107 (1985) 82–87. [76] S.K. Ider, F.M.L. Amirouche, Coordinate reduction in the dynamics of constrained multibody systems, New Approach, Am. Soc. Mech. Eng. (Paper) (1989). [77] R.L. Wang, J.T. Huston, A comparison of analysis methods of redundant multibody systems, Mech. Res. Commun. 16 (1989) 175–182. [78] J.T. Wang, R.L. Huston, Computational methods in constrained multibody dynamics: matrix formalisms, Comput. Struct. 29 (1988) 331–338, https://doi.org/ 10.1016/0045-7949(88)90267-2. [79] E. Pennestrì, P.P. Valentini, Coordinate reduction strategies in multibody dynamics: a review, Conf. Multibody Syst. Dyn. (2004) 1–17. [80] R.A. Wehage, E.J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained, J. Mech. Design 104 (1982) 247–255. [81] S.S. Kim, M.J. Vanderploeg, QR decomposition for state space representation of constrained mechanical dynamic systems, ASME J. Mech., Trans., Automation Design 108 (1986) 183–188, https://doi.org/10.1115/1.3260800. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 34
[82] E.J. Haug, Computer Aided Kinematics and Dynamics of Mechanical Systems, Allyn and Bacon Boston, 2021. [83] M.A. Serna, R. Avil´ es, J. García de Jal´ on, Dynamic analysis of plane mechanisms with lower pairs in basic coordinates, Mech. Mach. Theory. 17 (1982) 397–403, https://doi.org/10.1016/0094-114X(82)90032-5. [84] D.S. Lopes, M.T. Silva, J.A. Ambr´ osio, Tangent vectors to a 3-D surface normal: a geometric tool to find orthogonal vectors based on the Householder transformation, Computer-Aided Design 45 (2013) 683–694, https://doi.org/10.1016/J.CAD.2012.11.003. [85] H. Tan, L. Li, Q. Huang, Z. Jiang, Q. Li, Y. Zhang, D. Yu, Influence of two kinds of clearance joints on the dynamics of planar mechanical system based on a modified contact force model, Sci. Rep. 13 (2023) 1–27, https://doi.org/10.1038/s41598-023-47315-1, 2023 13:1. [86] K.H. Chang, e-Design: Computer-Aided Engineering Design, Academic Press, 2016. [87] X. Iriarte, J. Bacaicoa, A. Plaza, J. Aginaga, A unified analytical disk cam profile generation methodology using the Instantaneous center of rotation for educational purpose, Mech. Mach. Theory. 196 (2024) 105625, https://doi.org/10.1016/J.MECHMACHTHEORY.2024.105625. [88] P. Masarati, M.J.U. Quro, A. Zanoni, Projection continuation for minimal coordinate set formulation and singularity detection of redundantly constrained system dynamics, Multibody Syst. Dyn. 61 (2023) 453–480, https://doi.org/10.1007/S11044-023-09930-8. [89] I. Roupa, S.B. Gonçalves, M.T. Silva, Dynamic analysis of planar multibody systems with fully cartesian coordinates, in: Proceedings of International Conference on Multibody System Dynamics, Lisbon, Portugal, 2018, p. 2018. June 24-28. S.B. Gonçalves et al. Mechanism and Machine Theory 209 (2025) 105955 35