Dynamic analysis of planar rigid multibody systems modeled using natural absolute coordinates

This paper deals with the dynamic simulation of rigid multibody systems described with the use of two-dimensional natural absolute coordinates. The computational methodology discussed in this investigation is referred to as planar Natural Absolute Coordinate Formulation (NACF). The kinematic representation used in the planar NACF is based on a vector of generalized coordinates that includes two translational coordinates and four rotational parameters. In particular, the set of natural absolute coordinates is employed for describing the global location and the geometric orientation relative to the general configuration of a planar rigid body. The kinematic description utilized in the planar NACF is based on the separation of variable principle. Therefore, a constant symmetric positive-definite mass matrix and a zero inertia quadratic velocity vector associated with the centrifugal and Coriolis inertia effects enter in the formulation of the equations of motion. However, since a redundant set of rotational parameters is used in the kinematic description of the planar NACF for defining the geometric orientation of a rigid body, the introduction of a set of intrinsic normalization conditions is necessary for the mathematical formulation of the algebraic constraint equations. Thus, the intrinsic constraint equations associated with the natural absolute coordinates must be properly taken into account in addition to the extrinsic constraint equations that model the kinematic pairs which form the mechanical joints. This investigation discusses in details the mathematical derivation and the numerical implementation of the multibody system differential-algebraic equations of motion elaborated in the context of the planar NACF. For this purpose, simple geometric considerations are employed in the paper to develop the algebraic equations associated with the intrinsic and extrinsic constraints, whereas the fundamental principles of classical mechanics are utilized for the formal deduction of the dynamic equations. By using the augmented formulation, the index-three form of the differential-algebraic equations of motion is reduced to the corresponding index-one counterpart in order to be able to apply the Udwadia-Kalaba approach for the analytical calculation of the multibody system generalized acceleration vector. Furthermore, in the numerical implementation of the equations of motion based on the planar NACF, the direct correction method is utilized for stabilizing the algebraic constraint equations at both the position and velocity levels. The direct correction approach represents a new methodology recently developed in the field of multibody system dynamics for treating the algebraic constraint equations leading to physically correct and numerically stable dynamic simulations. A standard numerical integration algorithm is employed for obtaining an approximate solution of the nonlinear dynamic equations derived by using the planar NACF. The numerical implementation of a general-purpose multibody computer program based on the planar NACF is demonstrated in the paper considering four simple benchmark examples of rigid multibody systems. c © 2018 University of West Bohemia. All rights reserved.


Introduction
In modern scientific literature, the study of the dynamic behavior of mechanical systems constrained by kinematic joints is referred to as multibody system dynamics [17].Multibody mechanical systems are formed by collections of rigid and/or flexible continuum bodies interconnected between each other by means of mechanical joints, passive force elements, and active control actuators [39,40].Common examples of mechanical systems widely employed in engineering applications that can be modeled using the multibody system approach are machines, mechanisms, vehicles, robots, space structures, and biomechanical systems [8,9,63,69,75].The mechanical components of a multibody system can be connected in open-loop or closed-loop configurations and may also come into contact between one another or with the surrounding environment [45,50,70,77].The resultant evolution in time of a multibody system originates from the action of external force fields and applied forces, is influenced by the possible presence of contact and friction forces, and is heavily determined by the motion restrictions imposed by the kinematic joints [25,26,83,84].
One distinguishing feature that characterizes the motion of a general multibody system is the presence of large displacements, large finite rotations, and possibly small and/or large deformations [4].Therefore, the correct description of the rotational motion using a coordinate parameterization that is free of kinematic singularities if of fundamental importance for effectively performing dynamic simulations of complex multibody systems [71,100].On the other hand, the mechanical deformations of the continuum bodies and of the kinematic pairs that form a general multibody system can significantly affect the resulting motion of the system itself.The dynamic behavior of such complex systems is governed by nonlinear differential-algebraic equations of motion resulting from the inherently nonlinear nature of the mutual relationships between the multibody system components [98].Thus, it becomes necessary to adopt general and robust analysis approaches capable of correctly capturing geometric nonlinearities in order to systematically formulate and numerically solve the equations of motion and the algebraic equations of a general multibody system [3,47,51,74,91].This issue is particularly relevant for the optimal control problem of rigid-flexible multibody systems in which the standard algorithms based on the linearization of the equations of motion are inadequate [43,53,66,67].To this end, appropriate numerical techniques and formulation procedures are needed to obtain realistic numerical solutions for describing the general motion of flexible multibody systems such as, for example, the Floating Frame of Reference Formulation (FFRF) and the Absolute Nodal Coordinate Formulation (ANCF) [85].The FFRF is a formulation approach that allows for modeling large translations and large finite rotations of flexible multibody systems which exhibit small deformations, whereas in the ANCF computational framework both small and large deformations can be accurately described [61,86].One of the prominent features of both the FFRF and the ANCF is represented by the possibility to achieve an accurate description of complex curved geometry and the development of effective computational algorithms for the numerical resolution of the differential-algebraic equations of motion based on a nonincremental solution approach [49,54,62].Therefore, both the FFRF and the ANCF computational procedures for the description of flexible multibody systems can be readily used in conjunction with the analytical methodologies typically employed for the dynamic analysis of rigid multibody systems [60,72,73].
In recent years, multibody system dynamics has emerged as an independent research field in which several formulation strategies and numerical procedures have been developed for the systematic derivation and the numerical resolution of the differential-algebraic equations of motion that characterize the dynamic behavior of a general multibody system [79,95].In particular, there is a considerable difference between the formulation strategies employed for modeling rigid multibody systems and the formulation approaches adopted for the description of flexible multibody systems [82,87].The formulation procedures used for the simulation of the dynamic behavior of rigid multibody systems can be subdivided into three general categories, namely the redundant coordinate formulation, the recursive coordinate approach, and the semirecursive method.The redundant coordinate formulation is one of the principal approaches used in multibody system dynamics for modeling mechanical engineering applications.This formulation employs a nonminimal set of generalized coordinates for the kinematic description of a multibody system and, therefore, the dimension of the coordinate parameterization is larger than the strict number of degrees of freedom of the system itself [55].On the contrary, the recursive coordinate approach is widely used in robotics and makes use of a minimal set of generalized coordinates that are specific for each topology of the kinematic joints of a given multibody system.The semi-recursive method, on the other hand, represents a hybrid strategy in which the set of generalized coordinates includes both redundant and recursive variables [12].As extensively discussed in the multibody literature, the three general categories of formulation techniques mentioned before have distinct advantages and drawbacks of a theoretical nature or in terms of computational convenience [5,11,44,48].
The redundant coordinate formulation is of main interest for this investigation and represents the computational strategy employed in this paper.Following the concept underlying the redundant formulation technique for modeling the kinematics of two-dimensional and threedimensional multibody systems composed only of rigid bodies, numerous computational methodologies and computer codes were developed for the dynamic analysis of rigid multibody systems [88].In particular, the two principal formulation strategies that adopt the redundant formulation approach are referred to as Reference Point Coordinate Formulation (RPCF) and Natural Coordinate Formulation (NCF) [24].The RPCF is a robust and reliable formulation approach that is widely used for performing dynamic simulations of complex multibody systems [42].The kinematic description employed in the RPCF is based on a nonminimal set of generalized coordinates and can be distinguished into two types of coordinate formulations according to the sets of rotational coordinates used to define the orientation of a general rigid body.In this respect, one can distinguish the RPCF with Euler angles, that is based on the Cartesian coordinates of a reference point together with set of three Euler angles in a three-dimensional space, and the RPCF with Euler parameters, which uses the reference point coordinates as translational variables and the four components of a unit quaternion for the identification the body orientation in the space [1,58,59].The main advantage of the RPCF is the fact that the algebraic constraint equations that mathematically model the joint constraints can be expressed in a simple and general form even in the case of complex kinematic pairs [56].On the other hand, the NCF is an effective formulation procedure which can be used in conjunction with specially tailored numerical algorithms for performing real-time dynamic simulations of rigid multibody systems.The kinematic representation used in the NCF deliberately avoids the use of angular coordinates as rotational coordinates and exploits the separation of spatial-dependent variables and time-dependent coordinates leading to a constant mass matrix and zero centrifugal and Coriolis generalized inertia forces [22,23].However, the main challenge encountered in the computer implementation of the NCF is the definition of the algebraic equations that model the joint constraints which are difficult to obtain in a systematic manner.This paper, on the other hand, is based a formulation approach called Natural Absolute Coordinate Formulation (NACF) in which the advantages of both the RPCF and the NCF are ideally combined [64].Indeed, the NACF represents a robust formulation approach for the dynamic analysis of rigid multibody systems that has a clear geometric interpretation, involves simple algebraic manipulations, and allows for a systematic derivation of the differential-algebraic multibody equations of motion that can be readily implemented in a general-purpose multibody computer program [65].This investigation is focused on the development of a methodology for performing nonlinear dynamic simulations of two-dimensional rigid multibody systems modeled by using the set of natural absolute coordinates as generalized coordinates.The formulation approach discussed in this paper is referred to as planar Natural Absolute Coordinate Formulation (NACF) and represents the two-dimensional counterpart of the spatial NACF recently developed for the kinematic and dynamic analysis of three-dimensional rigid multibody systems.In general, the NACF computational framework is suitable for modeling rigid-flexible multibody systems in combination with the ANCF method and this computational strategy can effectively lead to an unified geometry/analysis approach that facilitates the integration of computer-aided design and analysis (I-CAD-A) [46,81].In a two-dimensional space, the set of natural absolute coordinates is composed only of Cartesian coordinates of appropriate Euclidean vectors.To this end, for identifying the general configuration of a two-dimensional rigid body, the Cartesian coordinates of a specified reference point are used as translational coordinates, while the components of two unit vectors endowed with the direction cosines of a body-fixed reference frame are employed as rotational coordinates.Therefore, the natural absolute coordinates employed in this investigation are natural in the sense that they include only generalized Cartesian coordinates and absolute in the sense that they are relative to geometric vectors associated with the general configuration of a rigid body represented with respect to an inertial frame of reference.In the planar NACF, the underlying kinematic representation is based on the separation of variable principle and leads to a substantial simplification of the mathematical derivation of the dynamic equations.In fact, a matrix of shape functions can be defined for a general rigid body employing the set of natural absolute coordinates leading to a kinematic representation featuring the separation of space-dependent variables and time-dependent variables.It is shown in the paper that the position, velocity, and acceleration fields of a general rigid body are vector functions which can be respectively written as linear combinations of the natural absolute coordinates, velocities, and accelerations.
Employing the formulation strategy discussed in this work, the algebraic equations that describe the geometric restrictions imposed by the joint constraints can be readily obtained in a concise form.In the mathematical formulation of the constraint equations in the context of the planar NACF, the algebraic constraints are conceptually divided into two families, namely the intrinsic constraint equations and the extrinsic constraint equations.While the intrinsic constraint equations are a mathematical representation of the physical rigidity of a given body that forms the multibody system, the extrinsic constraint equations formulate analytically the relationships between the bodies of the multibody system that originate from the presence of the kinematic pairs forming the mechanical joints.Furthermore, by using the D'Alembert-Lagrange principle of virtual work together with the Lagrange multiplier method the differential-algebraic equations of motion for a general multibody system can be systematically obtained in the framework of the planar NACF employing a standard assembly procedure [21].The mass matrix of a general rigid body modeled employing the set of natural absolute coordinates is a positive-definite symmetric constant matrix and, therefore, in the planar NACF the inertia quadratic velocity vector that includes the centrifugal and Coriolis inertia effects is a null vector.In particular, the collocation of the body-fixed reference system in correspondence of the principal axes of inertia associated with the center of mass of the rigid body leads to a diagonal mass matrix which implies a minimal inertia coupling of the equations of motion.It is important to note that having a constant mass matrix as well as a zero inertia quadratic velocity vector is advantageous for performing effective and efficient dynamic simulations.When compared with other multibody computational procedures based on the redundant formulation approach, in the planar NACF the extrinsic algebraic equations which mathematically model the kinematic joints assume a simple analytical form which in some important cases is linear and, therefore, can be eliminated at the preprocessing stage leading to a set of differential-algebraic equations of motion featuring a significantly reduced dimension [94].
The formulation procedure employed in the planar NACF leads to a sparse structure of the differential-algebraic equations of motion which can be readily implemented in a generalpurpose multibody computer program.However, the set of generalized coordinates used in the planar NACF is based on a nonminimal set of rotational parameters and, consequently, in the numerical solution of the equations of motion, an effective method for enforcing the algebraic equations at the position, velocity, and acceleration levels is needed in order to take into account the intrinsic constraint equations associated with the normalization conditions of the direction cosines as well as the extrinsic constraint equations that model the kinematic pairs of the mechanical joints.For this purpose, in this investigation, the Udwadia-Kalaba equations are used for computing the generalized acceleration vector of a rigid multibody system whereas the direct correction approach is employed for the numerical stabilization of the algebraic equations [2,13].The augmented formulation is used to transform the index-three form of the differential-algebraic equations of motion into the corresponding index-one form in order to obtain a standard formulation of the dynamic equations amenable to be treated with the use of the Udwadia-Kalaba approach to the analytical dynamics [92].The specific form of Udwadia-Kalaba equations employed in this investigation comes from the original formulation of the fundamental equations of constrained motion developed by Udwadia and Kalaba, and are appropriately modified in order to exploit the simplified structure of the dynamic equations that result from the development of the planar NACF [35].On the other hand, the direct correction approach is a methodology recently developed in the field of rigid multibody system dynamics for satisfying the algebraic equations of a general rigid multibody system [52].The key idea of the direct correction methodology is to devise effective correction terms for the entire sets of generalized coordinates and velocities in order to ensure the fulfillment of the algebraic equations at both the position and velocity levels [99].In analogy to the well-known generalized coordinate partitioning procedure, the direct correction method allows for the selection of a constraint tolerance that is enforced in the algorithm for the algebraic equations at each time step [96].Since the direct correction approach does not require the identification of dependent and independent coordinates, it is generally more efficient than the generalized coordinate partitioning method while conserves the robustness of the partitioning algorithm [97].
In this paper, the direct correction method is revisited in the framework of the planar NACF and is combined with the Udwadia-Kalaba method in order to obtain a reliable multibody solver for the differential-algebraic dynamic equations.Four numerical examples are considered in the paper in order to evaluate by means of numerical experiments the effectiveness and the efficiency of the formulation procedure analyzed in this work.In the numerical resolution of the differential-algebraic equations of motion of the four numerical examples considered as illustrative examples, the sixth-order explicit linear multistep method called Adams-Bashforth algorithm is used to march forward the numerical solution on the time grid and, at the same time, the direct correction approach is used to minimize the drift phenomenon of the algebraic constraint equations at both the position and velocity levels.The numerical results obtained in the paper employing the planar NACF are compared with the numerical results computed with the use of the analytical procedure recently developed for modeling two-dimensional multibody systems called planar RPCF with Euler parameters and considering the conventional modeling approach referred to as planar RPCF with Euler angles.A very good matching between the numerical results computed employing the planar NACF, the planar RPCF with Euler parameters, and the planar RPCF with Euler angles is found confirming the effectiveness of the computational methodology developed in this work.
This paper is organized as follows.In Section 2, the kinematic representation employed in the planar NACF for modeling two-dimensional rigid multibody systems is described.For this purpose, the set of natural absolute coordinates is defined and the separation of variable principle is used for obtaining the matrix of shape functions of a general two-dimensional rigid body.Section 3 encompasses a detailed description of the algebraic constraint equations used for modeling the intrinsic normalization conditions of a generic rigid body and the extrinsic constraint equations relative to the kinematic pairs that form the joint constraints of a general multibody system.In particular, the algebraic equations corresponding to the extrinsic constraints are subdivided in this section into four subsets of constraint equations called basic constraints in order to facilitate the construction of the mathematical equations relative to complex kinematic pairs.In Section 4, the fundamental principles of classical mechanics are used to derive the equations of motion of a rigid multibody system in the framework of the planar NACF.To this end, the D'Alembert-Lagrange principle of virtual work is utilized in conjunction with the Lagrange multiplier technique for obtaining the differential-algebraic form of the dynamic equations.Section 5 discusses the computational algorithm used in this investigation for the enforcement the algebraic constraints and the numerical solutions of the nonlinear dynamic equations.For this purpose, the direct correction approach is used to stabilize the algebraic constraints at both the position and velocity levels and the Udwadia-Kalaba method is employed for enforcing the constraint equations at the acceleration level.In Section 6, the computer implementation of a general-purpose multibody computer code developed in MATLAB based on the planar NACF considering four simple numerical examples of two-dimensional rigid multibody systems is demonstrated.It is shown in this section that the numerical results calculated by using the planar NACF are consistent with the numerical results computed employing the planar RPCF with Euler parameters as well as the planar RPCF with Euler angles.Finally, Section 7 provides a summary of the paper and the conclusions drawn from this investigation.

Reference kinematics
In this section, the kinematic representation employed in the framework of the planar NACF is described.To this end, the position, velocity, and acceleration fields of a general body i of a planar multibody system composed of rigid bodies are derived in terms of natural absolute coordinates.For a generic rigid body i, the set of natural absolute coordinates is formed by the Cartesian coordinates of a preassigned reference point on the rigid body i together with the components of the unit vectors associated with the direction cosines of a body-fixed reference system.In fact, in the kinematic description based on the planar NACF, two types of reference frames are employed.The first type of reference system is an inertial frame of reference that is called global reference system and serves as a common standard for the kinematic representation used in the planar NACF.The second type of reference system is a local frame of reference that is attached to a given point O i of the rigid body i.The origin O i of the body-fixed reference system represents the reference point associated with the rigid body i.The local reference system is also called floating reference frame since it translates and rotates with the generic rigid body i of the multibody system.The global position vector of the reference point O i of the rigid body i measured with respect to the inertial reference system is denoted with R i and is defined as: where R i 1 and R i 2 are the Cartesian coordinates of the reference point O i of the rigid body i defined with respect to the inertial reference frame.Consider an arbitrary point P i on a general rigid body i.The local position vector of a generic point P i on the rigid body i defined with respect to the floating reference system is denoted with ūi and is given by: where xi and ȳi are the Cartesian coordinates of an arbitrary point P i defined with respect to the local reference frame pertaining to the rigid body i. Employing fundamental geometric concepts of the mechanics of rigid bodies, the global position vector of a generic point P i on the rigid body i measured with respect to the inertial reference system is denoted with r i and can be written as: where R i is the global position vector of the reference point O i measured with respect to the inertial coordinate frame, A i is the rotation matrix which defines the orientation of the floating reference system with respect to the global coordinate system, and ūi is the local position vector of an arbitrary point P i on the rigid body i defined in the floating reference system associated with the body i.The rotation matrix A i transforms local vector quantities expressed in the floating reference system into their global counterparts represented in the inertial coordinate system [19].In a two-dimensional space, the rotation matrix A i can be formally obtained by using the set of direction cosines as follows: where α i and β i are the unit vectors formed by the direction cosines associated with the floating reference frame of the rigid body i and are given by: where α i 1 and α i 2 are the direction cosines referred to the unit vector directed along the horizontal axis of the floating reference frame of the body i, while β i 1 and β i 2 are the direction cosines associated with the unit vector oriented along the vertical axis of the floating reference system of the body i.The direction cosines that identify the unit vectors α i and β i can also be grouped to form an orientation vector denoted with δ i which is defined as: Considering an absolute angular displacement of the rigid body i denoted with θ i and evaluated counterclockwise from the horizontal axis of the inertial frame of reference, the two-dimensional set of direction cosines can be computed as follows: Since the two unit vectors α i and β i formed by the direction cosines represent a nonminimal set of orientation parameters, they must satisfy an appropriate set of normalization conditions.Indeed, the normalization conditions of the orientation parameters based on the direction cosines can be readily written as: where ϕ i is an intrinsic constraint vector which identifies the rigidity constraint equations of the body i which ensure that the rotation matrix A i is an orthogonal matrix.In fact, it can be proved that the normalization conditions of the orientation vectors α i and β i mathematically formulate the physical rigidity of a general body i that belongs to the rigid multibody system.Therefore, the normalization conditions of the orientation vectors formed by the direction cosines represent a set of intrinsic constraint equations.In the planar NACF, the general configuration of a rigid body i pertaining to a multibody system can be univocally identified by using a vector of natural absolute coordinates denoted with e i which contains the translational coordinate vector R i and the orientation parameter vector δ i of the rigid body i and is defined as: The vector of natural absolute coordinates e i represents the generalized coordinate vector employed in the planar NACF.By using the set of natural absolute coordinates as generalized coordinates, the position vector of an arbitrary point P i on the rigid body i can be expressed as follows: where e i is the vector of natural absolute coordinates and S i is the matrix of linear shape functions associated with the rigid body i defined as: where I is the identity matrix and H i is a constant matrix that depends on the Cartesian coordinates of a generic point P i defined with respect to the floating frame of reference on the rigid body i as follows: where xi and ȳi are the Cartesian coordinates of an arbitrary point P i represented with respect to the local reference frame associated with the rigid body i.The mathematical structure of the constant matrix of shape functions S i associated with the rigid body i demonstrates that the global position vector r i of an arbitrary point P i can be written as a linear combination of the global position vector of the body reference point R i with the unit vectors associated with the direction cosines α i and β i that identify the orientation of the floating reference frame with respect to the inertial reference frame.In particular, the constant coefficients of the linear combination are the Cartesian coordinates xi and ȳi of an arbitrary point P i defined with respect to the local reference frame and, therefore, the position field r i of a rigid body can be explicitly written as: Consequently, in the kinematic representation of the general geometric configuration of a generic rigid body i based on the planar NACF, it is apparent that the use of the set of natural absolute coordinates leads to the separation of the spatial-dependent variables xi and ȳi from the time-dependent variables R i , α i , and β i .Therefore, the planar NACF features the principle of separation of variables.The mathematical property of separation of variables considerably facilitates the kinematic description of a general rigid body i and leads to a substantial simplification of the derivation of the multibody system equations of motion in the context of the planar NACF.It is important to note that the position field r i of a general rigid body i is a linear function of the vector of natural absolute coordinates e i that serve as generalized coordinates in the planar NACF.Furthermore, employing the separation of variable technique, the virtual displacement field δr i , the velocity field ṙi , and the acceleration field ri of a general rigid body i can be respectively expressed by using the set of natural absolute coordinates as follows: Thus, in the planar NACF, it is clear that the virtual displacement field δr i of the rigid body i is a linear vector function of the natural absolute virtual displacement vector δe i , the velocity field ṙi of the rigid body i is a linear vector function of the natural absolute velocity vector ėi , and the acceleration field ri of the rigid body i is a linear vector function of the natural absolute acceleration vector ëi .On the other hand, the magnitude of the angular velocity vector defined in the floating reference system ωi can be written as a linear combination of the time derivative of the vector of orientation coordinates δi as: where the angular velocity transformation matrix Ḡi is defined as follows: The mathematical derivation in terms of natural absolute coordinates of the position, velocity, and acceleration fields of a general rigid body i completes the analytical description of the reference kinematics developed in the framework of the planar NACF.

Algebraic constraints
In this section, a comprehensive description of the algebraic equations necessary for the definition of the constraint equations in the context of the planar NACF is analyzed.In particular, the algebraic equations that mathematically represent the kinematic constraints are classified according to their physical nature into two general groups, namely the set of intrinsic constraint equations and the set of extrinsic constraint equations [90].The intrinsic constraint equations are peculiar for each set of rotational parameters employed as rotational coordinates in a general multibody formulation and are a mathematical expression of the body physical rigidity.In fact, in the planar NACF, the intrinsic constraint equations arise from the normalization conditions of the unit vectors based on the set of direction cosines that are employed as orientation parameters in the formulation of the equations of motion of multibody systems composed of two-dimensional rigid bodies.On the other hand, while enforcing the intrinsic constraint equations is of fundamental importance in order to ensure the rigidity of each body i pertaining to the multibody system, the extrinsic constraint equations originate from the presence of the kinematic joints and taking into account the extrinsic constraints is necessary for modeling the mechanical restrictions on the motion of the rigid bodies that form the multibody system [89].It is important to note that the correct formulation of the extrinsic constraint equations that mathematically represent the kinematic pairs of the mechanical joints of a general multibody system composed of rigid bodies is a fundamental step for performing reliable kinematic and dynamic analysis in the framework of the planar NACF.For this purpose, the algebraic equations for modeling complex kinematic joints can be derived from a simple set of extrinsic constraint equations which are called in this investigation basic constraints.The basic constraint equations associated with the extrinsic constraints serve as a set of elementary building blocks and allows for the formulation of the constraint equations of more complex kinematic pairs in a two-dimensional space such as, for example, the generalized joint, the rigid joint, the revolute joint, the prismatic joint, and others.In this work, the basic equations for the extrinsic constraints are classified into the following four types: generalized constraints, position constraints, orthogonality constraints of the first kind, and orthogonality constraints of the second kind.The accurate description of the algebraic equations used in the planar NACF for modeling the intrinsic constraints and the four types of basic equations corresponding to the extrinsic constraints is discussed in details in this section.

Intrinsic constraint equations for the orientation vectors
In this section, the algebraic equations relative to the intrinsic constraints of a general rigid body i modeled in a two-dimensional space with the use of natural absolute coordinates are developed.In the planar NACF, the intrinsic constraint equations allow for enforcing the rigidity of a generic body i and originate from the normalization conditions of the unit vectors composed of direction cosines which form a redundant set of orientation parameters.For a generic rigid body i, the intrinsic constraint equations can be written as follows: where ϕ i is the intrinsic constraint vector which mathematically models the rigidity property of the body i, whereas α i and β i are the unit vectors that contain the direction cosines representing the orientation of the floating reference system associated with the rigid body i.A virtual change of the intrinsic constraint vector ϕ i can be written as: where ϕ i e i is the Jacobian matrix of the intrinsic constraint vector ϕ i computed with respect to the vector of natural absolute coordinates e i of the rigid body i which is defined as follows: In addition, the first and the second time derivatives of the intrinsic constraint vector ϕ i can be readily calculated as: where Q i d,ϕ is the intrinsic constraint quadratic velocity vector that absorbs the terms which are quadratic in the generalized velocities arising from the second time derivative of the intrinsic constraint equations associated with the normalization conditions of the orientation vectors α i and β i based on the direction cosines.The intrinsic constraint quadratic velocity vector Q i d,ϕ is defined as: When an index-one form of the multibody system equations of motion is employed for performing dynamic simulations using the planar NACF, the analytical calculation of the intrinsic constraint quadratic velocity vector Q i d,ϕ associated with each rigid body i of the multibody system is required for the imposition of the intrinsic constraint equations of the body orientation parameters at the acceleration level.

Basic algebraic equations for the extrinsic generalized constraints
In this section, the algebraic constraint equations associated with the basic extrinsic constraints that mathematically describe the generalized constraints for a kinematic pair which interconnects two planar rigid bodies are developed in the context of the planar NACF.Consider a general set of generalized constraint equations labeled with the index k and connecting two generic rigid bodies of the multibody system identified by the indices i and j.The generalized constraints represent a set of extrinsic kinematic constraint equations between the natural absolute coordinates of two bodies i and j involved in the kinematic pair k and remove a preassigned number of relative degrees of freedom between the two rigid bodies.The basic algebraic equations that describe the kinematic restrictions of the generalized constraints are given by: where ψ k is the extrinsic constraint vector which identifies the basic algebraic equations for the generalized constraints, B i and B j are appropriate Boolean matrices, e i and e j are the natural absolute coordinate vectors of the rigid bodies i and j involved in the kinematic pair k, and d i,j is a constant vector.A virtual change of the extrinsic constraint vector ψ k yields the following equations: where ψ k e k is the Jacobian matrix of the extrinsic constraint vector ψ k calculated with respect to the vector of natural absolute coordinates e k involved in the kinematic pair k of the basic algebraic equations associated with the generalized constraints.The the Jacobian matrix ψ k e k is defined as: Furthermore, the first and the second time derivatives of the extrinsic constraint vector ψ k of the basic algebraic equations relative to the generalized constraints can be easily calculated as: where Q k d,ψ is the extrinsic constraint quadratic velocity vector that absorbs the terms which are quadratic in the generalized velocities written in terms of natural absolute coordinates.For the basic constraint equations of the generalized constraints, the constraint quadratic velocity vector Q k d,ψ is a zero vector: In the planar NACF, the analytical calculation of the extrinsic constraint quadratic velocity vector Q k d,ψ relative to the basic algebraic equations of the generalized constraints is necessary for the numerical implementation of the index-one form of the multibody system equations of motion.

Basic algebraic equations for the extrinsic position constraints
In this section, the algebraic constraint equations relative to the basic extrinsic constraints that mathematically describe the position constraints for a kinematic pair which interconnects two planar rigid bodies are developed in the framework of the planar NACF.Consider a general set of position constraints labeled with the index k that connect two points P i and P j belonging to two generic rigid bodies of the multibody system identified by the indices i and j.The position constraints remove two relative degrees of freedom between the bodies i and j involved in the kinematic pair k.The basic algebraic equations that describe the kinematic conditions of the position constraints can be written as: where ψ k is the extrinsic constraint vector which identifies the basic algebraic equations for the position constraints, r i and r j are the position vectors of the two points P i and P j on the two rigid bodies i and j involved in the kinematic pair k, and p i,j is a constant vector.A virtual change of the extrinsic constraint vector ψ k leads to the following equations: where ψ k e k is the Jacobian matrix of the extrinsic constraint vector ψ k computed with respect to the vector of natural absolute coordinates e k involved in the kinematic pair k of the basic algebraic equations relative to the position constraints.The Jacobian matrix ψ k e k is given by: where S i and S j are the matrices of shape functions evaluated in correspondence of the two points P i and P j on the two rigid bodies i and j involved in the kinematic pair k.Moreover, the first and the second time derivatives of the extrinsic constraint vector ψ k of the basic algebraic equations associated to the position constraints can be readily computed as: where Q k d,ψ is the extrinsic constraint quadratic velocity vector which absorbs the terms that are quadratic in the generalized velocities expressed in terms of natural absolute coordinates.For the basic constraint equations of the position constraints, the extrinsic constraint quadratic velocity vector Q k d,ψ is a zero vector: In the planar NACF, the analytical computation of the extrinsic constraint quadratic velocity vector Q k d,ψ associated with the basic algebraic equations of the position constraints is required for the numerical implementation of the index-one form of the multibody system equations of motion.

Basic algebraic equations for the extrinsic orthogonality constraints of the first kind
In this section, the algebraic constraint equations associated with the basic extrinsic constraints that mathematically describe the orthogonality constraints of the first kind for a kinematic pair which interconnects two planar rigid bodies are developed in the context of the planar NACF.Consider a general set of orthogonality constraints of the first kind labeled with the index k that involve two geometric directions P i 1 O i = P i 1 − O i and P j 2 O j = P j 2 − O j belonging to two generic rigid bodies of the multibody system identified by the indices i and j.The orthogonality constraints of the first kind remove one relative degree of freedom between the bodies i and j involved in the kinematic pair k.The basic algebraic equations that describe the kinematic restrictions of the orthogonality constraints of the first kind are given by: where ψ k is the extrinsic constraint vector which identifies the basic algebraic equations for the orthogonality constraints of the first kind, v i 1 and v j 2 are the unit vectors associated with the two geometric directions P i 1 O i and P j 2 O j on the two rigid bodies i and j involved in the kinematic pair k and identify two axes of the kinematic joint, and a i,j is a constant scalar quantity.A virtual change of the extrinsic constraint vector ψ k yields the following equations: where ψ k e k is the Jacobian matrix of the extrinsic constraint vector ψ k calculated with respect to the vector of natural absolute coordinates e k involved in the kinematic pair k of the basic algebraic equations associated with the orthogonality constraints of the first kind.The Jacobian matrix ψ k e k is defined as: where N i 1 and N j 2 are constant matrices associated with the unit vectors v i 1 and v j 2 corresponding to the two geometric directions P i 1 O i and P j 2 O j on the two rigid bodies i and j involved in the kinematic pair k which are respectively defined as follows: where H i 1 and H j 2 are constant matrices functions of the local positions of the two generic points P i 1 and P j 2 on the corresponding rigid bodies i and j involved in the kinematic pair k.The two unit vectors v i 1 and v j 2 associated with the two geometric directions P i 1 O i and P j 2 O j on the two rigid bodies i and j can be respectively written by using the constant matrices N i 1 and N j 2 as: Furthermore, the first and the second time derivatives of the extrinsic constraint vector ψ k of the basic algebraic equations relative to the orthogonality constraints of the first kind can be easily calculated as: where Q k d,ψ is the extrinsic constraint quadratic velocity vector which absorbs the terms that are quadratic in the generalized velocities written in terms of natural absolute coordinates.For the basic constraint equations of the orthogonality constraints of the first kind, the extrinsic constraint quadratic velocity vector Q k d,ψ is given by: In the planar NACF, the analytical calculation of the extrinsic constraint quadratic velocity vector Q k d,ψ relative to the basic algebraic equations of the orthogonality constraints of the first kind is necessary for the numerical implementation of the index-one form of the multibody system equations of motion.

Basic algebraic equations for the extrinsic orthogonality constraints of the second kind
In this section, the algebraic constraint equations associated with the basic extrinsic constraints that mathematically describe the orthogonality constraints of the second kind for a kinematic pair which interconnects two planar rigid bodies are developed in the framework of the planar NACF.Consider a general set of orthogonality constraints of the second kind labeled with the index k that involve the geometric direction P i 1 O i = P i 1 − O i belonging to the generic rigid body i and two points P i and P j belonging to two generic rigid bodies of the multibody system identified by the indices i and j.The orthogonality constraints of the second kind remove one relative degree of freedom between the bodies i and j involved in the kinematic pair k.The basic algebraic equations that describe the kinematic conditions of the orthogonality constraints of the second kind can be written as: where ψ k is the extrinsic constraint vector which identifies the basic algebraic equations for the orthogonality constraints of the second kind, v i 1 is the unit vector associated with the geometric direction P i 1 O i on the rigid body i in the kinematic pair k which identify the axis of the kinematic joint, r i and r j are the position vectors of the two points P i and P j on the two rigid bodies i and j involved in the kinematic pair k, and b i,j is a constant scalar quantity.A virtual change of the extrinsic constraint vector ψ k leads to the following equations: where ψ k e k is the Jacobian matrix of the extrinsic constraint vector ψ k computed with respect to the vector of natural absolute coordinates e k involved in the kinematic pair k of the basic algebraic equations relative to the orthogonality constraints of the second kind.The Jacobian matrix ψ k e k is given by: where N i 1 is a constant matrix relative to the unit vector v i 1 corresponding to the geometric direction P i 1 O i on the rigid body i in the kinematic pair k, whereas S i and S j are the matrices of shape functions evaluated in correspondence of the two points P i and P j on the two rigid bodies i and j involved in the kinematic pair k.Moreover, the first and the second time derivatives of the extrinsic constraint vector ψ k of the basic algebraic equations associated with the orthogonality constraints of the second kind can be readily computed as: where Q k d,ψ is the extrinsic constraint quadratic velocity vector which absorbs the terms that are quadratic in the generalized velocities expressed in terms of natural absolute coordinates.For the basic constraint equations of the orthogonality constraints of the second kind, the extrinsic constraint quadratic velocity vector Q k d,ψ is defined as: In the planar NACF, the analytical computation of the extrinsic constraint quadratic velocity vector Q k d,ψ relative to the basic algebraic equations of the orthogonality constraints of the second kind is required for the numerical implementation of the index-one form of the multibody system equations of motion.

Analytical dynamics
In this section, the analytical formulation of the dynamic equations in the framework of the planar NACF is discussed.By using the formulation approach developed in this section, the multibody equations of motion for a general rigid multibody system are systematically obtained employing a standard assembly procedure implemented in terms of natural absolute coordinates.For this purpose, the physical quantities that enter into the mathematical formulation of the differentialalgebraic equations of motion are derived using the set of natural absolute coordinates as generalized coordinates and exploiting the separation of variable property.In particular, the multibody system mass matrix, the generalized external force vector, the Jacobian matrix of the kinematic constraints, and the constraint quadratic velocity vector are analytically calculated in this section by using the D'Alembert-Lagrange principle of virtual work in conjunction with the Lagrange multiplier method.Considering a general rigid body i modeled in the planar NACF, the virtual work of the inertia forces is denoted with δW i i and can be expressed employing the kinematic representation based on the set of natural absolute coordinates as follows: where A i is the area of the two-dimensional body i, ρ i represents the body mass density, and M i denotes the mass matrix of the rigid body i.In the planar NACF, the mass matrix M i of a general two-dimensional rigid body i is a constant positive-definite symmetric matrix given by: where I represents the identity matrix, m i denotes the mass of the rigid body i, Ji O i ,x i and Ji are the first moments of mass of the two-dimensional body i computed with respect to the axes of the floating frame of reference having the point O i as origin, whereas Ji O i ,x i xi , Ji O i ,ȳ i ȳi , and Ji O i ,x i ȳi are the second moments of mass of the rigid body i calculated with respect to the axes of the body-fixed reference system.The first and second moments of mass of the rigid body i can be written in terms of the local Cartesian coordinates of the body center of mass G i and using the body mass moments of inertia as follows: where m i is the body mass, xi G i and ȳi G i represent the Cartesian coordinates of the center of mass G i of the rigid body i referred to the local coordinate system, while Īi are the mass moments of inertia associated with the body i and calculated with respect to the axes of the floating frame of reference having the point O i as origin, and Īi O i ,x i ȳi is the product of inertia of the rigid body i computed with respect to the axes of the body-fixed reference system.In the planar NACF, when the reference point O i is assumed coincident with the center of mass G i of the rigid body i and the axes of the floating frame of references coincide with the body principal axes of inertia, the constant mass matrix M i of the rigid body i has a special form characterized by a diagonal structure and can be written as: where m i is the mass of rigid body i, whereas Ji G i ,x i xi and Ji G i ,ȳ i ȳi are the second moments of mass of the rigid body i calculated with respect to the body principal axes of inertia referred to the center of mass G i of the rigid body i that are given by: where are the mass moments of inertia relative to the body i referred to the body principal axes of inertia having the body center of mass G i as origin of the coordinate frame.It is noteworthy to emphasize the fact that, in the planar NACF, the mass matrix M i of a general rigid body i is a constant positive-definite symmetric matrix featuring a full rank.Consequently, by using the formulation approach based on the set of natural absolute coordinates, the kinetic energy T i of a general rigid body i can be mathematically written as a quadratic form having constant coefficients that depends only on the body generalized velocities as follows: Therefore, since in the planar NACF the kinetic energy T i of a rigid body i is a pure quadratic form and the mass matrix M i is a constant square matrix, the inertia quadratic velocity vector Q i v that absorbs the centrifugal and Coriolis inertia terms which are quadratic in the generalized velocities vanishes.In fact, by using the Euler-Lagrange equations in the framework of the planar NACF, one can readily write: It is important to note that having a dynamic formulation of the equations of motion for a general rigid body i based on a constant mass matrix and zero Coriolis as well as centrifugal generalized inertia terms is a distinguishing feature of the planar NACF.This peculiar feature of the planar NACF is directly inherited from the separation of variable property employed in the kinematic description based on natural absolute coordinates.This important property of the planar NACF is advantageous for performing efficient and effective dynamic simulations for complex multibody mechanical systems and leads to a considerable simplification of the mathematical definition of the algebraic equations that represents the kinematic constraints.For instance, the body mass matrix can be transformed into an identity matrix by using a change of variables based on a set of Cholesky generalized coordinates leading to an optimal sparse structure of the multibody equations of motion.On the other hand, the virtual work of the generalized external forces δW i e that are applied on a general rigid body i can be easily written by using the set of natural absolute coordinates as follows: where A i is the area of the rigid body i, f i e is a general vector of external forces distributed on the two-dimensional body i, such as for instance the gravity force vector, and Q i e is the generalized external force vector acting on the rigid body i.In the planar NACF, the generalized external force vector Q i e is defined as: Furthermore, in the planar NACF, the virtual work of the generalized constraint forces δW i c,ϕ relative to the intrinsic constraint equations of the rigid body i pertaining to a general twodimensional multibody system can be written using the Lagrange multiplier method as: where λ i is the Lagrange multiplier vector associated with the intrinsic constraint equations of the two-dimensional body i of the rigid multibody system, ϕ i is the intrinsic constraint vector associated with the rigid body i, ϕ i e i is the Jacobian matrix of the intrinsic constraint vector ϕ i calculated with respect to the vector of natural absolute coordinates e i of the rigid body i, and Q i c,ϕ is the generalized constraint force vector resulting from the intrinsic constraints of the rigid body i.The generalized intrinsic constraint force vector Q i c,ϕ relative to the intrinsic constraint equations of the rigid body i is defined as follows: In the planar NACF, the virtual work of the generalized constraint forces δW k c,ψ associated with the the extrinsic constraint equations of the kinematic joint pertaining to a general kinematic pair k of the multibody system can be expressed using the Lagrange multiplier method as follows: where λ k is the Lagrange multiplier vector relative to the extrinsic constraint equations of the kinematic joint pertaining to a general kinematic pair k of the multibody system, ψ k is the extrinsic constraint vector relative to the kinematic joint k, ψ k e k is the Jacobian matrix of the extrinsic constraint vector ψ k calculated with respect to the vector of natural absolute coordinates e k of the rigid bodies involved in the kinematic pair k, and Q k c,ψ is the generalized constraint force vector resulting from the extrinsic constraints associated with the kinematic joint k.The generalized extrinsic constraint force vector Q k c,ψ associated with the the extrinsic constraint equations of the kinematic pair k is defined as: In the planar NACF, obtaining analytically the formal expressions of the rigid body mass matrix M i , the generalized external force vector Q i e , the intrinsic constraint generalized force vector Q i c,ϕ associated with a general rigid body i, and the extrinsic constraint generalized force vector Q k c,ψ relative to a generic kinematic pair k of the multibody system is necessary for assembling the complete set of differential-algebraic equations of motion.To this end, one can make use of the fundamental principles of analytical mechanics, such as the D'Alembert-Lagrange principle of virtual work, combined with the Lagrange multiplier method.Thus, the virtual work of all the forces acting on the multibody system can be expressed as follows: where δW i represents the total virtual work of the rigid body inertia forces, δW e denotes the total virtual work of the external forces applied to the multibody system, and δW c is the total virtual work of the intrinsic and extrinsic constraint forces that characterize the kinematic constraints of the mechanical system.Assuming that the number of rigid bodies that form the multibody system is N b and denoting with e the total vector of natural absolute coordinates that identify the general configuration of the multibody system, on can write the total virtual work of the rigid body inertia forces δW i in the planar NACF as follows: where M is the total mass matrix of the rigid multibody system which can be obtained by using a standard assembly procedure.The total virtual work of the external forces δW e can be readily obtained in the planar NACF as follows: where Q e is the total generalized external force vector associated with the complete set of external forces applied on the rigid multibody system which can be formulated employing a standard assembly procedure.Assuming a set of N b intrinsic algebraic constraints relative to the total set of rigid bodies that form the multibody system and a set of N c extrinsic algebraic constraints which model the joint constraints of all the kinematic pairs restricting the motion of the multibody system, on can write the total virtual work of the constraint forces δW c in the planar NACF by using the Lagrange multiplier technique as follows: where C is the total vector of algebraic constraint equations encompassing the complete set of intrinsic and extrinsic constraints, C e is the Jacobian matrix of the total vector of constraint equations C calculated with respect to the total vector of natural absolute coordinates of the multibody system e, and λ is the total vector of Lagrange multipliers associated with the entire set of algebraic constraints.In the planar NACF, the total vector of algebraic constraint equations C, the total Jacobian matrix of the kinematic constraints C e , and the total constraint quadratic velocity vector that absorbs the terms which are quadratic in the generalized velocities Q d are respectively defined as: where ϕ is the total vector of intrinsic constraint equations, ψ is the total vector of extrinsic constraint equations, ϕ e is the complete Jacobian matrix of the total intrinsic constraint vector ϕ computed with respect to the total vector of natural absolute coordinates of the multibody system e, ψ e is the complete Jacobian matrix of the total extrinsic constraint vector ψ computed with respect to the total vector of natural absolute coordinates of the multibody system e, Q d,ϕ is the total constraint quadratic velocity vector associated with the entire set of intrinsic constraints, and Q d,ψ is the total constraint quadratic velocity vector associated with the entire set of extrinsic constraints which can be obtained employing a standard assembly procedure.Assuming that the Lagrange multipliers satisfy the dependent subset of the system equations of motion and substituting the mathematical expressions of the total virtual works for the inertia, external, and constraint forces in the analytical formulation of the D'Alembert-Lagrange principle of virtual work, the total set of differential-algebraic equations of motion for the rigid multibody system can be developed in the context of the planar NACF to yield: As expected, in the planar NACF, the multibody system equations of motion form a set of index-three differential-algebraic dynamic equations in which the mass matrix is constant and the Coriolis and centrifugal generalized inertia terms are absent.Employing the augmented formulation approach, the index-three set of equations of motion obtained in the planar NACF can be readily transformed into a dynamically equivalent set of index-one equations of motion.
The augmented formulation represents an effective index reduction strategy which replaces the total vector of algebraic constraints with its second time derivative and leads to the following set of dynamic equations: The index-one form of dynamic equations of a general multibody system derived in terms of natural absolute coordinates and transformed using the augmented formulation can be rewritten in a compact matrix form as follows: The index-one form of the differential-algebraic equation of motion of a general rigid multibody system constrained by kinematic joints features an appropriate mathematical structure amenable to be treated with the use of the fundamental equations of constrained motion.

Computational algorithm
In this section, the computational algorithm used in this paper for the analytical treatment and the numerical solution the equations of motion of a general rigid multibody system obtained in the framework of the planar NACF is illustrated.In particular, the generalized acceleration vector of the multibody system constrained by a general set of kinematic joints is explicitly obtained in terms of natural absolute coordinates employing the fundamental equations of constrained motion referred to as Udwadia-Kalaba equations.Furthermore, the direct correction approach for eliminating the violations of the algebraic constraint equations at both the position and velocity levels is employed in the computational algorithm for the numerical solution of the multibody system equations of motion.A schematic representation of the complete Fig. 1.Flowchart of the computational algorithm for the numerical solution of the multibody system equations of motion computational algorithm used in the solution procedure of the equations of motion is shown in Fig. 1.The fundamental equations of constrained motion represent an effective mathematical technique recently developed in the field of analytical mechanics and are referred to as Udwadia-Kalaba equations by the names of their discoverers.Udwadia and Kalaba developed this powerful and general analytical method by using the Gauss principle of least constraint as a fundamental principle of classical mechanics in conjunction with modern linear algebra techniques based on the Moore-Penrose pesudoinverse matrix [93].In the Udwadia-Kalaba equations, the system mass matrix M, the generalized external force vector Q e , the Jacobian matrix of the constraint equations C e , and the constraint quadratic velocity vector Q d are known matrix and vector quantities, whereas the generalized acceleration vector ë and the vector of Lagrange multipliers λ are unknown vector quantities that can be readily computed employing the fundamental equations of constrained motion.In the planar NACF, the Udwadia-Kalaba equations can be rewritten in a special form by using the index-one form of the multibody equations of motion as follows: where a is the generalized acceleration vector expressed in terms of natural absolute coordinates for the multibody system released from the algebraic constraints, ε represents the constraint generalized error vector that mathematically quantifies how much the generalized acceleration vector a relative to the unconstrained multibody system violates the second time derivative of the algebraic constraint equations, K denotes the kinetic matrix that characterizes the constrained dynamics of the rigid multibody mechanical system represented using natural absolute coordinates, F is the constraint feedback matrix that originates from the application of the kinematic constraints on the rigid bodies that from the multibody system modeled in the planar NACF, λ is the total vector of Lagrange multipliers associated with the algebraic constraints, Q c denotes the total generalized constraint force vector relative to the algebraic constraint equations, a c represents the additional generalized acceleration vector induced on the multibody system by the action of the algebraic constraints, and ë is the resultant generalized acceleration vector of the multibody system expressed in terms of natural absolute coordinates.In the Udwadia-Kalaba equations, the matrix K + denotes the Moore-Penrose pseudoinverse matrix of the kinetic matrix K of the multibody system which can be readily computed using a set of standard mathematical procedures of numerical linear algebra such as for example the SVD numerical technique.However, if the Jacobian matrix of the kinematic constraints C e has a full row rank, which means that there are no redundant constraints in the vector of the algebraic equations C, the generalized inverse matrix which defines the Moore-Penrose pseudoinverse matrix is identical to the regular inverse matrix and, therefore, the constraint feedback matrix F can be computed from kinetic matrix of the multibody system K using a regular numerical inversion procedure based on the LU factorization [14,76].Thus, by using the Udwadia-Kalaba equations, the total generalized acceleration vector ë and the total vector of Lagrange multipliers λ can be effectively calculated for a general rigid multibody system modeled with the use of natural absolute coordinates and, therefore, the regular mathematical structure of the state function associated with the multibody system can be recovered allowing for the use of a standard numerical integration scheme for the numerical resolution of the equations of motion.However, since the numerical solution of the index-one form of the equations of motion is prone to the drift phenomenon of the algebraic constraints, a constraint stabilization technique, such as for instance the generalized coordinate partitioning algorithm, the direct correction technique, the penalty method, or the Baumgarte stabilization method, must be utilized for controlling the violations of the constraint equations at the position and velocity levels in order to obtain consistent numerical results [6].In particular, in this investigation, the direct correction technique is employed to contrast the violations of the algebraic constraint equations.The direct correction method is an effective and efficient constraint stabilization technique which can be readily implemented in a general-purpose multibody code.Unlike the conventional constraint stabilization techniques that are based on a trial-and-error selection of the stabilization parameters, the direct correction methodology is able to reduce within a specified error tolerance the violations of the algebraic constraints at both the position and velocity levels without modifying the structure of the equations of motion [20].Furthermore, unlike the generalized coordinate partitioning method, the direct correction algorithm does not require for each time step the analysis of the system degrees of freedom in which the identification of the system dependent and independent coordinates is performed [41].This is a crucial aspect of the generalized coordinate partitioning algorithm that contributes significantly to the robustness of the method but, at the same time, involves a significant additional overhead in the algorithm leading to a computational intensive numerical procedure [7,57].The direct correction approach, on the other hand, produces effective correction terms that are computed for the entire vectors of generalized coordinates and velocities without distinction between independent and dependent variables and leads to a more efficient computational procedure when compared to the generalized coordinate partitioning method [18].In analogy with the generalized coordinate partitioning method, the algorithmic structure of the direct correction method is composed of two numerical procedures, namely an attenuation of the constraint violations at the position level and a subsequent stabilization of the constraint drift at the velocity level.In the first step of the direct correction algorithm, the constraint violations at the position level are eliminated adding a correction term to the generalized coordinate vector as follows: where e is the corrected configuration vector containing the multibody system generalized coordinates updated by the direct correction procedure, e 0 denotes the initial uncorrected vector of generalized coordinates associated with a general time step of the numerical simulation, and Δe represents the correction term for the system natural absolute coordinates.In the of the direct correction terms associated with the total vector of generalized coordinates, the rectangular matrix C + e denotes the Moore-Penrose pseudoinverse matrix of the Jacobian matrix of the algebraic constraints C e that is a rectangular matrix.This first step of the direction correction algorithm is based on an approximate solution of the constraint equations at the position level based on the Moore-Penrose pseudoinverse matrix.Thus, the direct correction term for the natural absolute coordinates Δe has a minimum norm and makes the least influence possible on the dynamic equations.It is important to note that, in the direct correction algorithm, the direct correction step on the generalized coordinates must be repeated iteratively until a prescribed tolerance is met in order to effectively control the constraint violations at the position level.On the other hand, in the second step of the direct correction algorithm, the corrected vector of generalized coordinates obtained from the first step of the algorithm is employed and the constraint violations at the velocity level are eliminated considering a correction term for the generalized velocity vector which can be readily computed as follows: where ė denotes the corrected vector of generalized velocities of the multibody systems expressed in terms of natural absolute coordinates updated by the direct correction procedure, ė0 is the initial uncorrected velocity vector of generalized velocities relative to a general step of the dynamic simulation, and Δė identifies the correction term for the generalized velocities of the rigid multibody system.In the expression of the direct correction terms associated with the total vector of generalized velocities, the vector C t represents the partial derivative of the constraint vector C computed with respect to the time.Even in this case, this second step of the direction correction algorithm is based on an approximate solution of the constraint equations at the velocity level based on the Moore-Penrose pseudoinverse matrix.As a result, the direct correction term for the natural absolute velocities Δė has a minimum norm and, therefore, makes the least influence possible on the equations of motion.In general, in the direct correction algorithm, it is not necessary to repeat iteratively the direct correction step on the generalized velocities because typically the magnitude of the constraint violation at the velocity level is of a minor entity.The analytical formulation of the differential-algebraic set of equations of motion in terms of natural absolute coordinates, the description of the formulation approach for the calculation of the multibody system generalized acceleration vector based on the Udwadia-Kalaba equations, and the illustration of the numerical procedure for the enforcement of the algebraic constraints at the position and velocity levels by using the direct correction method complete the discussion on the dynamic equations in the framework of the planar NACF.

Numerical results and discussion
In this section, four illustrative numerical examples are examined in order to demonstrate the use of the planar NACF for the kinematic and dynamic analysis of rigid multibody systems confined in a two-dimensional space.As discussed in details in the previous section of the paper, the flowchart of the computational algorithm developed in this investigation for the numerical solutions of the equations of motion is shown in models of the four mechanical systems considered as numerical examples, the local position of the origin of the reference system of each body is assumed to be coincident with the body center of mass, while the coordinate system of each body is considered parallel to the body principal axes of inertia.Tables 1-4 report the numerical values of the masses, of the mass moments of inertia, and of the gravity accelerations of the rigid bodies that form the multibody systems considered as illustrative examples.The four-bar linkage is composed of three rigid bodies and four revolute joints that are respectively located in the points A, B, C, and D of Fig. 2. In Fig. 2, the centers of mass of the three rigid bodies that form the four-bar linkage are respectively indicated with G 1 , G 2 , and G 3 .As shown in Fig. 2, a spring-damper element is collocated between the points B and C of the four-bar linkage.The slider-crank mechanism is formed by three rigid bodies, three revolute joints, and one prismatic joint which are respectively located in the points A, B, C, and G 3 of Fig. 3.In Fig. 3, the centers of mass of the three rigid bodies which Four-bar linkage     It is important to note that, for all the multibody models considered as numerical examples, the initial conditions ensure the fulfillment of the algebraic constraint equations at both position and velocity levels.In fact, the four rigid multibody systems begin their time evolutions having zero initial generalized velocities and the resulting motions are the effects of the combination of the gravity force field with the force fields exerted by the spring-damper elements.Employing the D'Alembert-Lagrange principle of virtual work in conjunction with the Lagrange multiplier technique, the differential-algebraic equations of motion of the four multibody systems assumed as numerical examples are formally derived in the framework of the planar NACF.In order to enforce the rigidity of all the bodies of the four multibody systems under examination, an appropriate set of intrinsic constraint equations formed by the normalization conditions of the orientation parameters is taken into account in the derivation of the equations of motion.Subsequently, the index-three form of the differential-algebraic equations of motion based on the planar NACF is transformed into the index-one counterpart by using the augmented formulation, while the generalized acceleration vectors of the constrained multibody systems are calculated employing the Udwadia-Kalaba method.To this end, the multibody system equations of motion based on the planar NACF are implemented and solved in a general-purpose computer program for the dynamic analysis of rigid multibody systems developed in the MATLAB simulation environment.The numerical integration algorithm implemented for the numerical solution of the resulting set of ordinary differential equations of motion is the Adams-Bashforth method, which is a six-order explicit linear multistep algorithm with a constant time step.Since the Adams-Bashforth algorithm is not a self-starting numerical integration procedure, a sixth-order Runge-Kutta scheme is used for the first six time steps of the dynamic simulations.However, after the first six time steps, the Adams-Bashforth algorithm is an effective numerical integration scheme Four-bar linkage which requires only one evaluation of the system state function at each time step for marching forward the numerical solution of the dynamic equations on the time grid and, therefore, efficient dynamic simulations can be performed.The numerical parameters used to perform the numerical integration of the system equations of motion and to carry out the iterative numerical procedure for enforcing the constraint equations are reported in Table 10.In particular, in order to obtain a numerical solution of the multibody system equations of motion which is physically correct and numerically stable, the direct correction method is used to satisfy the algebraic equations arising from the intrinsic normalization conditions of the rigid bodies as well as the extrinsic constraints associated with the kinematic joints.For this purpose, in analogy with the wellknown generalized coordinate partitioning technique, the direct correction method represents an effective and efficient constraint stabilization algorithm for enforcing the algebraic equations of the kinematic constraints at the position and velocity levels without altering the structure of the dynamic equations.The maximum norms of the constraint violations at the position and velocity levels for the complete time evolution of the four multibody systems considered as illustrative examples are reported in Table 11.As expected, for all the multibody systems examined in this section, the maximum norms of the constraint violations at the position and velocity levels are below the tolerance prescribed in the iterative procedure that is employed in the computer implementation of the direct correction algorithm used for the numerical stabilization of the algebraic constraints.Thus, the dynamic behaviors of the four illustrative examples predicted by means of numerical simulations are consistent with the geometric topologies of the rigid multibody systems and comply with the physics of the constrained mechanical systems.In order to carry out a comparative study for the multibody formulation developed in this paper, the numerical results obtained by using the planar NACF are compared with the numerical results obtained employing the recently developed planar RPCF with Euler parameters as well as the conventional methodology based on the planar RPCF with Euler angles [68].In Fig. 6, the vertical displacement of the centroid G 2 of the second body of the four-bar linkage is shown.Fig. 7 represents the horizontal displacement of the centroid G 3 of the third body of the slider- crank mechanism.In Fig. 8, the vertical displacement of the centroid G 3 of the third body of the double inverted pendulum is shown.Fig. 9 represents the horizontal displacement of the centroid G 5 of the fifth body of the quick return mechanism.In Figs.6-9, the circles denote the numerical solutions computed by using the planar NACF, the squares indicate the numerical solutions obtained employing planar RPCF with Euler parameters, and the diamonds represent the numerical solutions calculated with the use of the planar RPCF with Euler angles.A complete overlapping of the numerical results obtained using these three different formulation procedures is found and, therefore, there is a very good agreement between the dynamic behaviors predicted by using the proposed planar NACF, the planar RPCF with Euler parameters, and the planar RPCF with Euler angles.Therefore, the numerical results presented in this section demonstrate the analytical accuracy, the physical consistency, and the correct computer implementation of the planar multibody formulation based on natural absolute coordinates developed in this paper.

Conclusion
The principal research objective of the authors is the development of new, efficient, and effective analytical methods and/or computational procedures in order to obtain accurate dynamic models [15,[27][28][29]80], perform reliable parametric identifications using experimental data [16,[30][31][32]36], and devise optimal control strategies for rigid-flexible multibody systems leveraging on the deep connections between multibody dynamics, system identification, and control theory [10,33,34,37,38,78].In particular, an analytical method encapsulated in a computational framework for modeling rigid multibody systems in a two-dimensional space was elaborated in this investigation.The computational framework described in this paper is based on the use of a set of two-dimensional natural absolute coordinates and, therefore, is called pla-nar Natural Absolute Coordinate Formulation (NACF).As demonstrated in this investigation, in the context of the redundant formulation approaches for performing a nonlinear dynamic analysis of planar rigid multibody systems, the planar NACF represents a viable alternative to the well-known planar Reference Point Coordinate Formulation (RPCF).In the kinematic description of the planar NACF, the Cartesian coordinates of a reference point and the Cartesian components of two unit vectors associated with a body-fixed reference system are used for identifying the general configuration of a rigid body in a two-dimensional space.In virtue of the separation of variable principle featured by the kinematic description that characterizes the planar NACF, the algebraic equations that model the kinematic joints can be systematically obtained in terms of simple mathematical expressions derived by means of elementary geometric considerations.More importantly, in the planar NACF, the mass matrix of a general rigid body is a constant positive-definite symmetric matrix and, consequently, the inertia quadratic velocity vector that absorbs the centrifugal and Coriolis inertia terms vanishes.Since in the planar NACF the orientation of a rigid body is represented by using redundant rotational coordinates, an additional set of constraint equations enters in the analytical formulation of the equations of motion.This additional class of algebraic equations arises from the normalization conditions of the orientation parameters which are formed by a set of direction cosines and are referred to as intrinsic constraints in order to conceptually distinguish them from the extrinsic constraints deriving from the mathematical representation of the kinematic pairs that identify the mechanical joints.As discussed thoroughly in the paper, the planar NACF allows for obtaining a systematic definition of the analytical quantities necessary for the deduction of the differential-algebraic equations of motion which can be readily obtained for a general multibody system using a standard assembly procedure.In the paper, the index-three form of the dynamic equations is directly transformed into the corresponding index-one form by using the augmented formulation, while the generalized acceleration vector of a general multibody system is obtained employing an analytical formulation based on the Udwadia-Kalaba methodology.Moreover, the recently developed direct correction method is formulated in terms of natural absolute coordinates and is used in this work as a constraint stabilization technique for automatically counteracting the drift phenomenon of the algebraic constraint equations.Therefore, the combination of the Uwadia-Kalaba equations with the direct correction approach devised in this investigation ensures that the constraint equations are enforced at the position, velocity, and acceleration levels, leading to accurate numerical solutions of the differential-algebraic equations of motion that are consistent with the multibody system geometry and comply with the fundamental laws of classical mechanics.In order to demonstrate the effectiveness of the proposed computational methodology by using numerical experiments, a general-purpose multibody code was developed in the MATLAB simulation environment and four simple numerical examples of open-loop and closed-loop multibody systems were analyzed in the framework of the planar NACF.For this purpose, the numerical resolution of the differential-algebraic multibody system equations of motion was carried out with the use of the standard Adams-Bashforth numerical integration method and, at the same time, the direct correction approach was employed to adjust the resulting state vector at each time step.In order to perform a comparative analysis for the four numerical examples examined in the paper, the multibody system equations of motion were also formulated using the recently developed planar RPCF with Euler parameters as well as the conventional planar RPCF with Euler angles.A very good agreement was found between the numerical results obtained by using the planar NACF, the planar RPCF with Euler parameters, and the planar RPCF with Euler angles.

Table 1 .
Four-bar linkage inertia properties

Table 3 .
Double inverted pendulum inertia properties

Table 4 .
Quick return mechanism inertia properties Table 5 the spring undeformed length, the spring constant stiffness coefficient, and the damper viscous damping coefficient for each multibody system considered as illustrative examples are reported.Furthermore, Figs.2-5 represent the initial configurations for the four multibody systems considered as numerical examples.The initial configurations of all the rigid bodies that constitute the four multibody systems considered as numerical examples are respectively described in Tables 6-9.

Table 6 .
Initial configuration of the four-bar linkage

Table 7 .
Initial configuration of the slider-crank mechanism

Table 8 .
Initial configuration of the double inverted pendulum

Table 9 .
Initial configuration of the quick return mechanism

Table 10 .
Parameters of the numerical simulations Time step Time span Constraint tolerance

Table 11 .
Violations of the algebraic constraints in the planar NACF