
Constraint violation correction is an important research topic in solving multibody system dynamics. For a multibody system dynamics method which derives acceleration equations in a recursive manner and avoids overall dynamics equations, a fast and accurate solution to the violation problem is paramount. The direct correction method is favored due to its simplicity, high accuracy and low computational cost. This method directly supplements the constraint equations and performs corrections, making it an effective solution for addressing violation problems. However, calculating the significant Jacobian matrices for this method using dynamics equations can be challenging, especially for complex multibody systems. This paper presents a programmatic framework for deriving Jacobian matrices of planar rigid-flexible-multibody systems in a simple semi-analytic form along two paths separated by a secondary joint. The approach is verified by comparing constraint violation errors with and without the constraint violation correction in numerical examples. Moreover, the proposed method's computational speed is compared with that of the direct differential solution, verifying its efficiency. The straightforward, highly programmable and universal approach provides a new idea for programming large-scale dynamics software and extends the application of dynamics methods focused on deriving acceleration equations.
Citation: Lina Zhang, Xiaoting Rui, Jianshu Zhang, Guoping Wang, Junjie Gu, Xizhe Zhang. A framework for establishing constraint Jacobian matrices of planar rigid-flexible-multibody systems[J]. AIMS Mathematics, 2023, 8(9): 21501-21530. doi: 10.3934/math.20231096
[1] | Liandi Fang, Li Ma, Shihong Ding . Finite-time fuzzy output-feedback control for $ p $-norm stochastic nonlinear systems with output constraints. AIMS Mathematics, 2021, 6(3): 2244-2267. doi: 10.3934/math.2021136 |
[2] | Tan Zhang, Pianpian Yan . Asymmetric integral barrier function-based tracking control of constrained robots. AIMS Mathematics, 2024, 9(1): 319-339. doi: 10.3934/math.2024019 |
[3] | Jan Nordström, Fredrik Laurén, Oskar Ålund . An explicit Jacobian for Newton's method applied to nonlinear initial boundary value problems in summation-by-parts form. AIMS Mathematics, 2024, 9(9): 23291-23312. doi: 10.3934/math.20241132 |
[4] | Mingliang Zheng . Study on the symmetries and conserved quantities of flexible mechanical multibody dynamics. AIMS Mathematics, 2023, 8(11): 27969-27982. doi: 10.3934/math.20231430 |
[5] | Sani Aji, Poom Kumam, Aliyu Muhammed Awwal, Mahmoud Muhammad Yahaya, Kanokwan Sitthithakerngkiet . An efficient DY-type spectral conjugate gradient method for system of nonlinear monotone equations with application in signal recovery. AIMS Mathematics, 2021, 6(8): 8078-8106. doi: 10.3934/math.2021469 |
[6] | Yimeng Xi, Zhihong Liu, Ying Li, Ruyu Tao, Tao Wang . On the mixed solution of reduced biquaternion matrix equation $ \sum\limits_{i = 1}^nA_iX_iB_i = E $ with sub-matrix constraints and its application. AIMS Mathematics, 2023, 8(11): 27901-27923. doi: 10.3934/math.20231427 |
[7] | Zhengqi Zhang, Huaiqin Wu . Cluster synchronization in finite/fixed time for semi-Markovian switching T-S fuzzy complex dynamical networks with discontinuous dynamic nodes. AIMS Mathematics, 2022, 7(7): 11942-11971. doi: 10.3934/math.2022666 |
[8] | Ke Su, Wei Lu, Shaohua Liu . An area-type nonmonotone filter method for nonlinear constrained optimization. AIMS Mathematics, 2022, 7(12): 20441-20460. doi: 10.3934/math.20221120 |
[9] | Qiangqiang Zhang, Yiyan Han, Chuandong Li, Le You . Constraint impulsive consensus of nonlinear multi-agent systems with impulsive time windows. AIMS Mathematics, 2020, 5(4): 3682-3701. doi: 10.3934/math.2020238 |
[10] | Chunyan Liu . The traveling wave solution and dynamics analysis of the fractional order generalized Pochhammer–Chree equation. AIMS Mathematics, 2024, 9(12): 33956-33972. doi: 10.3934/math.20241619 |
Constraint violation correction is an important research topic in solving multibody system dynamics. For a multibody system dynamics method which derives acceleration equations in a recursive manner and avoids overall dynamics equations, a fast and accurate solution to the violation problem is paramount. The direct correction method is favored due to its simplicity, high accuracy and low computational cost. This method directly supplements the constraint equations and performs corrections, making it an effective solution for addressing violation problems. However, calculating the significant Jacobian matrices for this method using dynamics equations can be challenging, especially for complex multibody systems. This paper presents a programmatic framework for deriving Jacobian matrices of planar rigid-flexible-multibody systems in a simple semi-analytic form along two paths separated by a secondary joint. The approach is verified by comparing constraint violation errors with and without the constraint violation correction in numerical examples. Moreover, the proposed method's computational speed is compared with that of the direct differential solution, verifying its efficiency. The straightforward, highly programmable and universal approach provides a new idea for programming large-scale dynamics software and extends the application of dynamics methods focused on deriving acceleration equations.
Multibody system dynamics (MSD) guides practical engineering applications by studying the movement of systems with multiple linked elements. Several main approaches for MSD have been developed since the 1970s [1,2,3,4,5,6,7,8]. The absolute coordinate methods [9], known for their programmatic nature and ability to handle both tree and non-tree systems, are commonly utilized. The approaches employ global coordinates to describe the motion of the system and obtain a general form for the dynamical equations, written as
{M¨q+ΦTqλ=fΦ(q,t)=00, | (1) |
where q denotes the generalized coordinates, which consist of the position and orientation of each rigid body-fixed coordinate system in the global inertial coordinate system, as well as the position and orientation of the flexible body floating coordinate system in the global inertial coordinate system and the flexible deformation generalized coordinates. ¨q, the second time derivative of q, denotes the generalized accelerations. M represents a generalized mass matrix. f represents generalized forces on the system. λ consists of Lagrange multipliers. Φ describes the constraint equations caused by the connection relation of all joint elements in the system. Φq denotes the partial derivative of the constraint equations with respect to generalized coordinates.
However, when dealing with a complex multi-degree-of-freedom system, the solution of the dynamics can be quite challenging to derive. To overcome this issue, one technique utilized in structural dynamics is the static condensation method, which was extended by Alessandro Cammarata and colleagues to mechanisms and structures with internal joints [10]. Alternatively, one can use relative coordinates to establish the dynamics equations. The relative coordinate methods, developed by Robertson and Wittenburg [3,6], reduce the number of unknown variables and are used in tree structures with elegant types. Among the relative coordinate methods, highly programmatic and efficient approaches are coming to the fore, including Featherstone's algorithm [11,12], the reduced multibody system transfer matrix method (RMSTMM) [13,14], the Brandl method [15,16], etc. The approaches utilize hinge coordinates as generalized coordinates and avoid global dynamics equations with a system inertia matrix. They can fully describe the configuration of tree systems [17,18] in the form of ordinary differential equations, written as
ˉq=f(q,˙q,t), | (2) |
where q, ˙q and ¨q are the generalized coordinates, velocities and accelerations, which have analogous meanings to those given in Eq (1) but with much lower dimensions. The focus of these methods is to recursively derive generalized accelerations ¨q for a mechanical multibody system by utilizing its state vectors (q,ˉq). However, in the context of closed-loop systems, the hinge coordinates are no longer independent [19]. To overcome this limitation, a set of differential algebraic equations (DAEs) is established by adding a set of algebraic equations that represent the constraint of cut-joints [4,7], as follows
{¨q=f(q,˙q,t)c(q)=00, | (3) |
where c is assembled by the constraint equations of all cut joints. Despite the approaches' effectiveness, numerical solutions are prone to divergence, as indicated by the red line in Figure 1.
Ensuring the accuracy of solutions for DAEs is known to be numerically challenging. To address this issue, various methods have been developed and broadly classified into three main categories: (i) constraint stabilization approaches, (ii) coordinate partitioning methods and (iii) direct correction methods [20,21]. Constraint stabilization approaches involve developing and modifying the Baumgarte method, which introduces the control feedback theory. However, one problem with the Baumgarte stabilization method is the ambiguity of the stabilization parameters [22]. To address this issue, Lin and Huang [23] proposed a method that used the Runge-Kutta algorithm to determine the stabilization parameters. Other approaches, such as those by Zhang [24] and Flores [25], assisted with stabilization coefficient determination by combining it with Taylor expansion, but the methods lack catholicity.
Another approach is the coordinate partitioning method, which involves dividing the generalized coordinates into independent and dependent sets. Numerical integration is then performed for the independent generalized coordinates, while the constraint equations are solved for the dependent generalized coordinates [26]. However, this method is numerically inefficient, as it necessitates frequent changes to the set of independent coordinates, limiting its effectiveness [27]. Decomposing the Jacobi matrix, currently available as triangular decomposition [28], singular value decomposition [29] or orthogonal trigonometric decomposition [30], is the technique used to accomplish this method.
The direct correction method, which is popular in correcting violation errors, introduces position and velocity corrections based on direct integration [31]. Over the years, a variety of direct correction techniques have been developed [20,32,33,34,35]. Yoon et al. [36] established a direct correction method that directly corrects the values of state variables, leading to a better fit with the constraint equations. However, this method is limited to the position level. Yu and Chen [37] corrected violation errors by implementing constraints on both position and velocity levels. The effectiveness of this approach was demonstrated through a simple case, which was compared to the standard formulation and the Baumgarte method. The direct correction method offers a more accurate solution compared to stabilization methods and is less computationally costly than coordinate partitioning methods [38]. In dealing with constraint violations arising from the DAEs, a fast and accurate solution to violation errors is crucial. The direct correction method is favored due to its simplicity, high accuracy and low computational cost.
For dynamics methods involving global dynamics equations, the analytic form of the constraint Jacobian matrix can be obtained explicitly [39]. However, for methods avoiding global dynamics equations and using hinge coordinates, the value of the constraint Jacobian matrix becomes uncertain dealing with constraint violation problems in closed-loop systems. To address this issue, the direct differentiation method can be used in computational programming by taking partial derivatives of each quality encountered along the path rather than solving the Jacobian matrix in its specific form. As the complexity of the closed-loop system increases, the constraint Jacobian matrix of the system in the direct correction method becomes more complex. Nevertheless, this approach is computationally intensive and time-consuming when dealing with large and complex systems, which limits the advantages of the direct correction method.
In this paper, semi-analytic Jacobian matrices of the direct correction method are derived to solve the divergence of algorithms. The correction method is first introduced in Section 2. A further attempt is to build a general framework for constraint Jacobian matrices. The form of constraint equations in the framework is shown in Section 3. The correlated variables affected by generalized coordinates in the constraint equations are encapsulated in the class of a secondary joint (cut-joint in a closed-loop system [40]). In Section 4, the variables of a secondary joint can be described by primary joints (other joints except for the secondary joint in the closed-loop system [40]) or flexible bodies. Among the correlation variables described, the ones affected by generalized coordinates are encapsulated in the corresponding class of the primary joint or the flexible body. The Jacobian matrices between different classes are programmatically derived in Section 5. The quantities of the secondary joint are cycled in the closed loop, sweeping through the classes of the primary joint or the flexible body. The solved Jacobian matrix is then substituted back into the constraint equations. Section 6 gives the pseudo-code for solving Jacobian matrices by direct differentiation. The corrected generalized coordinates and velocities are obtained using the Newton-Raphson formula in Section 7. Section 8 offers numerical examples validating the program. Finally, the concluding remarks are given in Section 9.
The generalized coordinates of the closed-loop system are not independent, causing constraint violations. Corrected generalized coordinates and velocities can be obtained by the Newton-Raphson method. The Newton-Raphson formula of the constraint equation on position level can be written as
00=c(q+Δq)≈c(q)+cq(q)Δq, | (4) |
where c denotes constraints on position level and cq(q) is the constraint Jacobian matrix of the constraint equations with respect to the generalized coordinates. Here, the generalized coordinates include the modal coordinates and the coordinates of the primary joint elements.
Upon differentiating the aforementioned equation, the correction equation for the generalized velocity is obtained, wherein the well-known relation cq=˙cq exists. The correction formula can be expressed as
00=˙c(q,˙q+Δ˙q)≈˙c(q,˙q)+˙c˙qΔ˙q=˙c(q,˙q)+cqΔ˙q. | (5) |
If the system has m constraint equations and n generalized coordinates, the constraint Jacobian matrix is given as
∂c∂q=[∂c1∂q1∂c1∂q2⋯∂c1∂qn∂c2∂q1∂c2∂q2⋯∂c2∂qn⋮⋮⋱⋮∂cm∂q1∂cm∂q2⋯∂cm∂qn]. | (6) |
The constraint equations of a secondary joint on both position and velocity levels can be established. The description of secondary joints, as shown in Figure 2, is based on the works of Roberson [40]. A span system is generated by cutting the joint which is called the secondary joint. The joints still in the span system are called primary joints. A joint connects two bodies of the system. One of the two bodies is called the base body and the other one is the moving body of the joint. The points B and M are fixed on the secondary joint, connecting the base body and the moving body, respectively.
For example, a planar pin joint in the jth closed loop labeled Lj is cut. The relative displacement of the pin joint moving in the plane is restricted. The constraint equation on position level could be expressed as
cLj=[rOM−rOB]Lj=00, | (7) |
˙cLj=[˙rOM−˙rOB]Lj=00, | (8) |
where rOB denotes the position of the base point w.r.t. the inertial frame. rOM denotes the position of the moving point. ˙rOB (˙rOM ) is the absolute velocity of the base (moving) point. Note that if the vector in the plane is r=[xy]T, then its cross product matrix ˜r=[−yx] is defined. Required variables rOM,rOB,˙rOM,˙rOB in Eqs (7) and (8) will be encapsulated in the class of the secondary joint, as shown in Figure 3. Function CalculateD0ConstraintEq() denotes constraints on position level. Analogically, the function CalculateD1ConstraintEq() gathers the process of constraints on velocity level. The Jacobian matrix, also known as the matrix of partial derivatives, is an essential tool for making accurate corrections in the system. The process of solving the constraint Jacobian matrix is put into the function CalculateConstraintJacobian. The process of solving the constraint Jacobian matrix is shown in the following.
The constraint equations established are ultimately affected since the position of the secondary joint is affected by state vectors of a system. The effect of each state vector in the system on the constraint equations will be considered here. For the pin secondary joint mentioned above, focusing on a specific column ∂c/∂qi(i=1⋯n) in Eq (6), there is
∂c∂qi=∂rOM∂qi−∂rOB∂qi. | (9) |
For convenience, two separate paths are discussed here. If element i is associated with the moving point of the pin secondary joint, variables rOB are not related to generalized coordinate qi. The partial derivative of rOM with respect to qi in this path is rOM,q. According to Eq (9), the partial derivatives of the constraint equation on position level with respect to the generalized coordinates in this path can be written as
cLj,qi=[rOM,q]. | (10) |
Similarly, on the other path, there is
cLj,qj=[rOB,q]. | (11) |
Intermediate variables rOM,qi,rOB,qi are collected in group I in a unified form. The generalized coordinates of a multi-rigid-flexible body system include the joint coordinate and the modal coordinates of flexible bodies. These are discussed below.
For elucidating the effect of generalized coordinates on constraint equations, subscript P can be replaced with base B or moving M point of the secondary joint. A primary joint moving in a plane connects to bodies with a root end described as R and a tip end described as T. r and A are chosen to represent the position level quantities considering the different types of primary joints. As shown in Figures 4 and 5, the kinematic quantities of the P node can be expressed as
{rOP=rOR+rRT+rTP=rOR+AORlRT+AORARTlTPAOP=AORARTATP, | (12) |
where subscript O denotes the inertial coordinate system. Subscripts R and T denote the body-fixed coordinate system with origin R and T, respectively. rOP denotes the position of the P point w.r.t. the inertial frame. lRT describes position vector from the R point to the T point described in the body-fixed coordinate system with the R point. Cosine matrix AOP for the coordinate transformation from the body-fixed coordinate system with the P point to the inertial frame reads as
AOP=[cosθOP−sinθOPsinθOPcosθOP]. | (13) |
Generalized coordinates qi in the equation only impact the parameters of the element, namely, lRT and ART. Therefore, variables in Eq (12) varying with generalized coordinate qHi can be expressed as
{rOP,q=AORlRT,q+AORART,qlTP,AOP,q=AORART,qATP, | (14) |
where the partial derivatives with respect to the generalized coordinates are indicated as subscripts, denoted as ∂rOM/∂qi=rOM,q.
Variables lRT,q,ART,q changing with the generalized coordinate of a primary joint are encapsulated in group II of the primary joint class. The unified modeling language (UML) diagram [41] of a primary joint is provided in Figure 6. For convenience, the above derivation is loaded in the function PreprocessConstraintJacobian of the primary joint.
A smooth pin joint moving in a plane is taken as an example to illustrate the main process in detail, as shown in Figure 6. Pin joint i rotates around the z-axis of the body-fixed coordinate system with the T point, which has only one generalized coordinate. The generalized coordinate can be described as
qHi=[θz]. | (15) |
Variables in group Ⅱ can be described here. The pin joint is rotated around the z-axis by an angle θRT, and it has no relative slip along the z-axis, lRT=00⋅ART is a direction cosine matrix corresponding to a simple rotation around the z-axis in the body-fixed frame. The corresponding ART,lRT can be described as
{lRT=00,ART=Az(θz), | (16) |
Variables lRT,ART with respect to generalized coordinates can be read as
{∂lRT∂θz=00,∂ART∂θz=Az(θz)D=ARTD, | (17) |
The above equations are substituted into Eq (14). The kinematic quantities influenced by the generalized coordinates become
{∂rOP∂θz=˜rTP,∂AOP∂θz=AOTDATOTAOP=DAOP, | (18) |
where rTB=AOTlTB. Meanwhile, the relationship AOTDATOT=D can be proven.
For further insight into partial derivatives of constraints with respect to state vectors of a primary joint, Eq (18) is elucidated. Only intermediate variables in group Ⅱ are used for subsequent assembly in the framework for the system Jacobian matrix.
The tip and root ends, denoted by R and T, describe the points on a flexible body connected to joints, as shown in Figure 7. Under the assumption of small deformation, dynamics equations of flexible bodies are described by the widely used floating coordinate system. The motion of the body is regarded as a superposition of the large range motion of the floating coordinate system and the small elastic deformation [42,43,44,45].
Using the kinematic quantities of a flexible body i, the kinematic quantities of node P for position level can be described as
{rOP=rOR+rRO+rOfT+rTP=rOR−AORAROfρR+AORAROfρT+AOTATOTrTP,AOP=AORAROfAOfTATP, | (19) |
where subscript f denotes a floating coordinate system. ρR(ρT) denotes a deformed displacement vector from the origin of floating frame Of to the R(T) end described in the floating frame. Variables rOR,AOR,θOR do not vary with the modal coordinates because they are transferred from the element associated to the flexible body, similar to the primary joint. ATOTrTP described in the body-fixed frame are independent of generalized coordinates and velocities, respectively. Only quantities ARO,ρR,ρT and AOfT change with modal coordinates. Kinematic quantities for position level affected by the modal coordinates become
{rOP,qf=−AORAROf,qfρR−AORAROfρR,qf+AORARO,qfρT+AORAROfρT,qf+AORAROf,qfAOTATOTrTP+AORAROAOfT,qfATOTrTPAOP,qf=AORARO,qfAOfTATP+AORAROAOfT,qfATP. | (20) |
Observing the above equation, intermediate variables ARO,qf,ρR,qf,AOfT,qf,ρT,qf are selected in group II. The quantities AROf,qf,ρR,qf,ρT,qf,AOfT,qf vary with the modal coordinates. The others vary with the generalized velocities. The UML diagram of a flexible body is given in Figure 8. The above derivation is loaded in the function PreprocessConstraintJacobian.
For flexible body i, the dynamics equations are derived based on the floating frame of reference. Component modal synthesis methods are applied to reduce the order of the dynamics equations of the flexible bodies [42,43,44,45]. The modal coordinates of the flexible body are denoted as
qFi=[qf]. | (21) |
The quantities of the flexible body will be presented here, as shown in Figure 9. Radius vectors ρR,ρT and angular deformation vector θL,Nk of the Nk node could be described as
{ρR=uOfR+Φu,RqfρT=uOfT+Φu,TqfθL, Nk=Φθ,Nkqf | (22) |
where θL,Nk describes the attitude of the flexible body moving in the plane. uOfR(uOfT) represents the original coordinate of node R(T) in the undeformed configuration. The modal analysis and synthesis are employed to discretize deformation displacement Φu,Rqf(Φu, Tqf).Φ is a shape function, i.e., modal function matrix in modal analysis. Subscript u(θ) denotes deformation in the displacement (angle). The directional cosine matrix of the flexible body from the body-fixed frame oNkxNkyNk to the floating frame is represented as AOfNk. Directional cosine matrix AOfNk can be expressed as
AOfNk=[cosθL, Nk−sinθL, NksinθL, NkcosθL, Nk]. | (23) |
Directional cosine matrix AOfNk cannot be reduced using the same method as that in a primary joint. However, the appearance of the directional cosine matrix AOfNk in the equation is inevitably multiplied by a variable independent of the modal coordinates. A mathematical relationship has been discovered, which is
∂AOfNku∂qf=∂AOfNku∂θL, Nk∂θL, Nk∂qf=DAOfNkuΦθ,Nk. | (24) |
The proof process is as follows. Since u is an unrelated variable, ∂AOfNku/∂t can be obtained by the direct derivative rule as follows:
∂AOfNku∂t=DAOfNk˙θL, Nku | (25) |
The chain rule can be used for the same variable ∂AOfNku/∂t, giving
∂AOfNku∂t=∂AOfNku∂θL, Nk∂θL, Nk∂t=∂AOfNku∂θL, Nk˙θL, Nk. | (26) |
Combining Eqs (25) and (26) expressed in different forms, there is
∂AOfNku∂θL, Nk=DAOfNku. | (27) |
Substituting Eq (27) into the original formula, the relationships between relevant variables ARof,q, AOfT,q,u and qf are
∂AOfTu∂qf=∂AOfTu∂θL,OfT∂θL,OfT∂qf=DAOfTuΦθ,T, | (28) |
∂AROfu∂qf=∂AROfu∂θL.OfR∂θL.OfR∂qf=−DAROfuΦθ,R. | (29) |
The intermediate variables in group II could be expressed as
{∂ρR∂qf=Φu,R,∂ρT∂qf=Φu, T,∂AROfu∂qf=−DAROuΦθ,R,∂AOfTu∂qf=DAOfTuΦθ,T. | (30) |
Substituting Eq (30) into Eq (20), variables rOP,q,AOP,q affected by state vectors of a flexible body can be read as
{∂rOP∂qf=D(AOOfρR−AOOfρT−rTP)Φθ,R−AOOfΦu,R+AOOfΦu,T+DrTPΦθ,T,∂AOP∂qf=−AORDARPΦθ,R+AOOfDAOfPΦθ,T. | (31) |
The existence of Eq (31) makes it possible to see the partial derivatives of the kinematic quantities clearly and also facilitates the composition of the constraint Jacobian matrix in the following. Nonetheless, this step does not appear in the actual programmatic implementation but only facilitates understanding of the solution concepts.
For a system with arbitrary topologies, the Jacobian matrix is a given temporary matrix changed with time. The calculation of the Jacobian matrix becomes very complicated, so a framework is established to simplify computation:
A programmatic derivation for the constraint Jacobian matrices is described using intermediate variables in groups I and Ⅱ. It is known that the constraint Jacobian matrix for correcting the generalized coordinates is identical to correcting generalized velocities. A closed-loop subsystem Lj in which the secondary pin joint j is associated with a primary pin joint i is taken as an example to explain in detail. The process of assembling Jacobian matrix cLj,qHj for generalized coordinate qHi is first demonstrated. Assume that pin joint i is associated with the base point of the secondary joint. Intermediate variables rOB,q in group I can be expressed by variables in group Ⅱ as
rOB,q=AORlRT,q+AORART,qlTB. | (32) |
Intermediate variables under the effect of the primary joint become
∂rOB∂qHi=AOR∂ART∂θzlTB, | (33) |
Constraint Jacobian matrix cLj,qHi could be expressed as
cLj,qHi=[−AOR(∂ART/∂θz)lTB]Lj, | (34) |
Submitting Eq (17) into the above equation, we have
cLj,qHj=−˜rTB. | (35) |
In the other path, pin joint i is associated with the moving point. The position of the moving point can be expressed as
rOM,q=AORlRT,q+AORART,qlTM. | (36) |
Constraint Jacobian matrix cLj,qHi in this path could be expressed as
cLj,qHi=[AOR(∂ART,q/∂θz)lTM]Lj=[˜rTM]Lj. | (37) |
Observing Eqs (35) and (37), the Jacobian matrix of the hinge i can be solved by calling the intermediate variables of the corresponding hinge element. Putting the calculated Jacobian matrix into the corresponding position, the Jacobian matrix for the closed-loop subsystem Lj is obtained. If the closed-loop subsystem Lj is without any flexible bodies, the corresponding Jacobian matrix can be written as
cLj,q=[cLj,qH1⋯cLj,qHn]. | (38) |
The assembly process can be applied for flexible bodies. Jacobian matrix cLj,qFi is derived as follows, if flexible body element i is associated with the base point of the secondary joint. According to Eq (20), kinematic quantities of the flexible body element are used to express intermediate quantities rOB,qf. Intermediate quantities can be read as
rOB,qf=−AORAROf,qfρR−AORAROρR,qf+AORAROOf,qfρT+AORAROρT,qf+AORARO,qfAOTATOTrTB+AORAROAOT,qfATOTrTB. | (39) |
According to Eq (30), the partial derivative of rOB with respect to the modal coordinates can be written as
∂rOB∂qf=−AOR(∂AROfρR/∂qf)−AOOf(∂pR/∂qf)+AOR(∂AROfρT/∂qf)+AOOf(∂ρT/∂qf)+AOR(∂AROrTB/∂qf)+AORAROf(∂ATOOfrTB/∂qf) | (40) |
According to the constraint equation of the secondary pin joint shown in Eq (7), constraint Jacobian matrix cLj,qFi reads as
cLj,qFi=[AOR(∂AROfρR/∂qf)+AOO(∂ρR/∂qf)−AOR(∂AROfρT/∂qf)−AOOf(∂ρT/∂qf)−AOR(∂AROrTB/∂qf)−AORARO(∂ATOOfrTB/∂qf)]=[D(AOOfρR−AOOfρT−rTB)Φθ,R−AOOfΦu,R+AOOfΦu,T+DrTBΦθ,T]. | (41) |
Similarly, if flexible body i is associated with the moving point of the secondary joint, point P will be replaced with M in Eq (31). The partial derivative of the constraint equation on position level with respect to modal coordinates qFi can be written as
cLj,qFi=[AORf(∂AROρR/∂qf)+AOOf(∂ρR/∂qf)−AOR(∂AROfρT/∂qf)−AOOf(∂ρT/∂qf)−AOR(∂AROfrTB/∂qf)−AORARO(∂ATfOrTM/∂qf)]=[D(AOOfρR−AOOfρT−rTB)Φθ,R−AOOfΦu,R+AOOfΦu, T+DrTMΦθ,T]. | (42) |
The Jacobian matrix of the flexible body i can be solved by calling the intermediate variables of a flexible body. For the closed-loop subsystem Lj, the corresponding Jacobian matrix consists of that for each element, written as
cLj,q=[cLj,qH1⋯cLj,qHncLj,qF1⋯cLj,qFn]. | (43) |
The Jacobian matrix for all primary joints is described as cLj,qH, and the Jacobian matrix for flexible bodies is described as cLj,qF. The Jacobian matrix cLj,q can be composed as
cLj,q=[cLj,qHcLj,qF], | (44) |
where Jacobian matrix cLj,qH is composed of the Jacobian matrix for each primary joint swept through the path, i.e., cLj,qH=[cLj,qH1cLj,qH2⋯cLj,qHn]. Jacobian matrix cLj,qF is composed of the Jacobian matrix for each flexible body swept through the path, i.e., cLj,qF=[cLj,q1cLj,qF2⋯cLj,qFn].
According to the topological structure of a multibody system, Jacobian matrix cq can be obtained by assembling the Jacobian matrices of all sub-closed loops in the system, denoted as
cq=[cL1,qcL2,q⋮cLn,q]. | (45) |
Then, the system Jacobian matrix can be expressed as
cq=[cL1,qcL2,q⋮cLn,q]=[cL1,qHcL1,qFcL2,qHcL2,qF⋮⋮cLn,qHcLn,qF]. | (46) |
The specific programming procedure for the system Jacobian matrix can be described below.
The universalized composition process of the system Jacobian matrix based on the intermediate variables in groups I and II can be implemented by computer, as shown in Algorithm 1. The Jacobian matrix for the jth secondary joint is solved by sweeping through the primary joint elements and flexible body elements along the path from the base point of the jth secondary joint element to the root element. In this process, the PreprocessConstraintJacobian() function in the classes of the primary joint element and the flexible body element is called to solve the intermediate variables in group I. The same process takes place on the other path. Corresponding constraint Jacobian matrices cLj,qHi and cLj,qϝi are obtained. According to the location of the generalized coordinates, constraint Jacobian matrices cLj,qHj and cLj,qfi are composed to obtain constraint Jacobian matrix cLj,q of the jth secondary joint. System Jacobian matrix cq is obtained ultimately by assembling constrained Jacobian matrices cLj,q(j=1…n) according to Eq (46). The pseudo-code of the computational process is illustrated in Algorithm 1.
In the programmatic derivation, multiple calculations of the intermediate variables in group II are avoided. The required variables can be obtained in the first calculation and called continuously in the subsequent calculations.
Algorithm 1. Computation of system Jacobian matrix |
Function assemblingsystemJacobianmatrix (i:int, ω:double, s:int) |
{Return Jacobian matrix cq of the system} |
{ cL,q is a Jacobian matrix of the secondary joint element.} |
{ SecondaryJointElement is an indicator of the secondary joint element.} |
{ PrimaryJointElement is an indicator of the primary joint element.} |
{ FlexibleBodyElement is an indicator of the flexible body element.} |
for k = 1 to the number of secondary joint elements do |
for i from element connected to the BasePoint of secondary joint k to the root element in closed-loop k do |
if element i is a primary joint element then |
SecondaryJointElement[k]. cL,q [i] ← SecondaryJointElement[k].Calculate |
ConstraintJacobian() ← PrimaryJointElement[i].PreprocessConstraintJacobian() |
end |
if element i is a flexible body element then |
SecondaryJointElement[k]. cL,q [i] ← SecondaryJointElement[k]. Calculate |
ConstraintJacobian() ← FlexibleBodyElement[i].PreprocessConstraintJacobian() |
end |
end |
for j from element connected to the MovingPoint of secondary joint k to the root element in closed-loop k do |
if element j is a primary joint element then |
SecondaryJointElement[k]. cL,q [j] ← SecondaryJointElement[k].Calculate |
ConstraintJacobian() ← PrimaryJointElement[j].PreprocessConstraintJacobian() |
end |
if element j is a flexible body element then |
SecondaryJointElement[k]. cL,q [j] ← SecondaryJointElement[k]. Calculate |
ConstraintJacobian() ← FlexibleBodyElement[j].PreprocessConstraintJacobian() |
end |
end |
cq [k]←SecondaryJointElement[k]. cL,q |
return cq |
The direct differentiation method obtains the Jacobian matrix by differentiating the constraint equations with respect to state variables in accordance with the chain rule of differentiation, which obtains an accurate calculation result.
The Jacobian matrix can be calculated row-wise. The same row is the partial derivatives of the constraints with respect to the same generalized coordinate. The partial derivatives in the same row are passed along the transfer path from the system root end to the position of the secondary joint recursively.
The pin joint is still used as an example with the constraint equation shown in Eq (7). In this case, the partial derivative of r with respect to the generalized coordinates qj is put in the corresponding row. To obtain the solution for ∂rOM/∂qi, the direct differentiation method is applied. The parameters ∂rOTG/∂qi and ∂AOTG/∂qi of the root element of the system can be solved using initial conditions. There exists a relationship where the parameter of the point connected to the root element is equal to the parameter of the root point of the element j to which it is connected. With this information, the quantities ∂rORj/∂qi and ∂AORj/∂qi of the element j can be determined. The kinematic relations of the element j can be used to express ∂rOTj/∂qi,∂AOTj/∂qi for the tip end, which can be written as
rOT=rOR+AORlRT | (47) |
AOT=AORART | (48) |
The process is iterated until the values of ∂rOT /∂qi and ∂AOT /∂qi for the tip end connected element to point M are obtained. The geometry of the transfer is such that the parameter of the M point is equivalent to the parameter of the tip end, leading to the determination of ∂rOM/∂qi. This process is then repeated cyclically to obtain ∂roM/∂qi(i=1⋯a) for different generalized coordinates. The solution for ∂rOB/∂qi(i=a⋯n) is obtained in a similar fashion as that of ∂rOM/∂qi(i=1⋯a).
By analogy, the partial derivatives are calculated separately for all the generalized coordinates in the system. The system Jacobian matrix is obtained, where the partial derivatives are placed in the appropriate rows. The direct differentiation method for Jacobian matrices uses a for loop to achieve recursion, and the algorithm is shown below.
Algorithm 2. Computation of system Jacobian matrix. |
Function assemblingsystemJacobianmatrix (i:int, ω:double, s:int) |
{Return Jacobian matrix cq of the system} |
{ cL,q is a Jacobian matrix of the secondary joint element.} |
{ SecondaryJointElement is an indicator of the secondary joint element.} |
{ PrimaryJointElement is an indicator of the primary joint element.} |
{ FlexibleBodyElement is an indicator of the flexible body element.} |
{ BodyElement is an indicator of the body element.} |
{ Index1 (Index2) is an indicator of the generalized coordinate.} |
for i from 0 to the number of generalized coordinates do |
Index1 + = the number of generalized coordinate of the element[i] |
Ground. ∂r/∂qi.setzero(). |
Ground. ∂A/∂qi.setzero(). |
for j from 0 to the number of body elements do |
Index2 + = the number of generalized coordinate of the element[j] |
if Index2 = Index1 then |
PrimaryJointElement[j].PreprocessConstraintJacobian() |
FlexibleBodyElement[j].PreprocessConstraintJacobian() |
end |
PrimaryJointElement[j] ∂r/∂qi(∂A/∂qi) ←Ground. ∂r/∂qi(∂A/∂qi). |
BodyElement[j] ∂r/∂qi(∂A/∂qi) ←PrimaryJointElement[j]. ∂r/∂qi(∂A/∂qi). |
end |
for k = 1 to the number of secondary joint elements do |
SecondaryJointElement[k]. cL,q ←SecondaryJointElement[k]. ∂r/∂qi(∂A/∂qi) |
end |
end |
for k = 1 to the number of secondary joint elements do |
cq [k]←SecondaryJointElement[k]. cL,q |
end |
return cq |
As illustrated in Figure 10, the complete correction procedure is described as follows:
(i) Uncorrected generalized coordinates and velocities are obtained during the numerical solution. Then, constraint equations on both position and velocity levels are established in the class of secondary joint elements. The established constraint equations of each secondary joint are combined into the system constraint equations.
(ii) If any values in the system constraint equations are greater than the tolerance error 10-10, the system will be penalized to correct the generalized coordinates and velocities. Algorithm 1 or Algorithm 2 is used to assemble the system Jacobian matrix. Using Eq (4), the corrected generalized coordinates can be solved. The value of the constraint equation is calculated again. If any value of the constraint violation errors is still greater than the tolerance error, continue with the correction. Otherwise, the correction is stopped.
(iii) The system constraint equations on velocity level are obtained by assembling the constraint equations on velocity level of each secondary joint. If the constraint violation errors on velocity level are greater than the tolerance error, the generalized velocity of the system is corrected using Eq (5). Then, the value of the constraint equation is calculated again. If the constraint violation errors are less than the tolerance error, the corrected generalized coordinates and velocity are obtained.
The simplicity of engineering models, such as the four-link mechanism, can have a wide range of practical applications, such as frames and rotating dampers [46,47]. Therefore, the ability to simulate these models quickly and accurately is crucial for engineering. A four-link mechanism is used to demonstrate the validity of the method. The system is depicted in Figure 11 and is subject only to gravity. All rigid links are assumed to be identical in length (1 meter) and mass (1 kg), and the secondary joint C is removed to create a spanning system. Element 3 is a flexible beam with two modes. Dynamic calculations were performed using the RMSTMM method, and initial conditions were set at q1=q3=π/3,q2=−π/3,q4=q4=0. A fourth-order Runge-Kutta algorithm with a very small step size (0.01 s) was employed for the numerical solution.
Figure 12 compares the constraint violations on x-direction of secondary joint C with and without the correction. The comparison confirms the correctness of the correction and the derived Jacobian matrix. Figure 13 shows violations of secondary joint C on x-direction by different methods. The results demonstrate that both the direct differential method and the proposed approach are successful in eliminating constraint violations. However, as anticipated, the proposed method is faster in terms of computational speed, as demonstrated by the comparison in Table 1 for a simple four-link mechanism.
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 19.036 | 19.766 | 24.875 |
100 | 38.017 | 38.531 | 44.547 |
200 | 79.593 | 80.985 | 87.516 |
In contrast to the previous example which had only a single closed loop, the upcoming example features an increased number of elements, closed loops and flexible body elements as depicted in Figures 14 and 15. The body elements 5, 11 and 21 are flexible and have two order modes, and the other elements are still rigid. The system is moved by gravity only. The dynamics of the system is simulated using the RMSTMM.
The result shows a divergent trend without constraints, as shown in Figure 16. The violations obtained by the Direct differentiation methods and proposed approach are shown in Figure 17, and both effectively keep the violations within the tolerance error. As shown in table 2, as the system calculation time increases, so does the CPU time used by the system. What is clear is that the longer the system computation time is, the longer the CPU time used by the direct differentiation method compared to the proposed method.
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 73.125 | 73.469 | 88.593 |
100 | 146.969 | 148.625 | 177.89 |
200 | 296.516 | 311.297 | 371.703 |
The third numerical example is shown in Figures 18 and 19 to verify the effectiveness of the programmatic stabilization solution. The example with rigid and flexible bodies is labeled with thick solid arrows and circles. Joint 19 is a secondary joint affiliated with the corresponding spanning tree. For primary joint, 1, 3, 5, 7, 9 are pin joints, and 11, 13, 15, 17 are prismatic joints. The others are rigid body elements labeled 2, 4, 8, 10, 12, 14, 16 and 18. The quantities of the body elements are the same as in the previous example but in the shape of rods.
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 34.922 | 35.641 | 42.968 |
100 | 75 | 78.235 | 96.391 |
200 | 144.75 | 151.297 | 181.25 |
The example uses the RMSTMM for dynamics analysis, too. If the correction is not added, the violation error on the x-direction of secondary joint 19 shows a clear divergence trend in Figure 20. The direct correction method is imposed while solving the Jacobian matrix using the direct differentiation method and the proposed method, respectively. The variation of the violation error is shown in Figure 21. Both methods are workable for correcting violation errors. An attempt is made below to compare the computation speed.
It can be found that the proposed approach is faster. The longer the system simulation time is, the longer the CPU time used by the direct differentiation method compared to the proposed method.
For computational speed assessment, the proposed approach and the direct differentiation are compared for a different total number of bodies in Figure 22. The abscissa represents the count of connected bodies within a single closed-loop system. The specific configuration of this system is illustrated in Figure 19. It is worth noting that the number of rods in the system shown in Figure 19 is no longer limited to 7. Instead, it includes various possibilities, such as 100, 200, 300 and 500 rods, respectively.
Upon comparing the two methods in Figure 22, it becomes evident that the proposed strategy outperforms the direct differentiation method in terms of computational time. The results indicate that the proposed strategy is far more efficient than the direct differentiation method.
This work proposes a framework for programmatical deriving of a system Jacobin matrix, which has significant advantages for correcting violations in the efficient nonlinear dynamics method focused on acceleration level. Simultaneously, Jacobian matrices in the direct correction method on both position and velocity levels are derived simply semi-analytically. The class for the system Jacobian matrix is designed, leading to a better-structured framework. Intermediate variables in the class of a secondary joint described by variables in the class of primary joints or flexible body elements can be used to assemble the Jacobian matrix. Next, the Jacobian matrix is assembled according to the topological structure of the system. The Newton-Raphson iteration method is used to obtain the corrected generalized coordinates and velocities. The proposed method is sufficiently accurate and faster than the widely used direct differentiation method. Numerical examples are used to verify the proposed approach's effectiveness by comparing the constraint violation errors of multi-rigid-flexible body systems with different methods. Due to their simplicity, the proposed method is well suited for rapid, draft-style modeling of multibody systems.
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
This work was supported by the National Natural Science Foundation of China (Grant No: 11902158), the Natural Science Foundation of Jiangsu Province (Grant No: BK20190438) and the National Natural Science Foundation of China (Grant No: 92266201).
The authors declare that they have no conflict of interest.
[1] | F. M. Amirouche, Computational methods in multibody dynamics, Englewood Cliffs, NJ: Prentice-Hall, 1992. |
[2] | J. G. De Jalon, E. Bayo, Kinematic and dynamic simulation of multibody systems: The real-time challenge, Springer, 2012. https://doi.org/10.1007/978-1-4612-2600-0 |
[3] | W. Jens, Dynamics of systems of rigid bodies, Berlin: Springer, 2013. https://doi.org/10.1007/978-3-322-90942-8 |
[4] | Y. Liu, Z. Pan, X. Ge, Dynamics of multibody systems, Beijing: Higher Education Press, 2014. |
[5] | P. E. Nikravesh, Computer-aided analysis of mechanical systems, Upper Saddle River: Prentice-Hall, 1988. |
[6] | R. E. Roberson, R. Schwertassek, Dynamics of multibody systems, Berlin: Springer, 2012. https://doi.org/10.1007/978-3-642-86464-3 |
[7] | A. A. Shabana, Dynamics of multibody systems, New York: Cambridge University Press, 2020. https://doi.org/10.1017/9781108757553 |
[8] | S. Werner, Multibody systems handbook, Berlin: Springer, 1990. https://doi.org/10.1007/978-3-642-50995-7 |
[9] |
W. M. Silver, On the equivalence of Lagrangian and Newton-Euler dynamics for manipulators, Int. J. Rob. Res., 1 (1982), 60–70. https://doi.org/10.1177/027836498200100204 doi: 10.1177/027836498200100204
![]() |
[10] |
A. Cammarata, R. Sinatra, P. D. Maddìo, Static condensation method for the reduced dynamic modeling of mechanisms and structures, Arch. Appl. Mech., 89 (2019), 2033–2051. https://doi.org/10.1007/s00419-019-01560-x doi: 10.1007/s00419-019-01560-x
![]() |
[11] | R. Featherstone, Rigid body dynamics algorithms, Springer, 2014. https://doi.org/10.1007/978-1-4899-7560-7 |
[12] |
F. I. T. Petrescu, Advanced dynamics processes applied to an articulated robot, Processes, 10 (2022), 640. https://doi.org/10.3390/pr10040640 doi: 10.3390/pr10040640
![]() |
[13] |
X. Rui, J. Zhang, X. Wang, B. Rong, B. He, Z. Jin, Multibody system transfer matrix method: The past, the present, and the future, Int. J. Mech. Syst. Dyn., 2 (2022), 3–26. https://doi.org/10.1002/msd2.12037 doi: 10.1002/msd2.12037
![]() |
[14] |
R. Xue, D. Bestle, Reduced multibody system transfer matrix method using decoupled hinge equations, Int. J. Mech. Syst. Dyn., 1 (2021), 12. https://doi.org/10.1002/msd2.12026 doi: 10.1002/msd2.12026
![]() |
[15] |
H. Brandl, R. Johanni, M. Otter, A very efficient algorithm for the ssimulation of robots and similar multibody systems without inversion of the mass matrix, IFAC Proc., 19 (1986), 95–100. https://doi.org/10.1016/S1474-6670(17)59460-4 doi: 10.1016/S1474-6670(17)59460-4
![]() |
[16] |
A. Cammarata, R. Sinatra, P. D. Maddio, Interface reduction in flexible multibody systems using the floating frame of reference formulation, J. Sound Vib., 523 (2022). https://doi.org/10.1016/j.jsv.2021.116720 doi: 10.1016/j.jsv.2021.116720
![]() |
[17] | K. S. Anderson, Recursive derivation of explicit equations of motion for efficient dynamic/control simulation of large multibody systems, Stanford University, 1990. |
[18] |
A. Jain, G. Rodriguez, Recursive flexible multibody system dynamics using spatial operators, J. Guid. Control Dyn., 15 (1992), 1453–1466. https://doi.org/10.2514/3.11409 doi: 10.2514/3.11409
![]() |
[19] | H. Brandl, An algorithm for the simulation of multibody systems with kinematic loops, Proc.7th World Congr. Theory Mach. Mech., 1987. |
[20] |
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 doi: 10.1007/s11044-016-9530-y
![]() |
[21] | P. E. Nikravesh, Some methods for dynamic analysis of constrained mechanical systems: a survey, Berlin: Springer, 1984. https://doi.org/10.1007/978-3-642-52465-3_14 |
[22] |
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 doi: 10.1016/0045-7825(72)90018-7
![]() |
[23] |
S. T. Lin, J. N. Huang, Stabilization of baumgarte's method using the Runge-Kutta approach, J. Mech. Design, 124 (2002). https://doi.org/10.1115/1.1519277 doi: 10.1115/1.1519277
![]() |
[24] |
P. Zhang, A Stabilization of constraints in the numerical solution of Euler-Lagrange equation, Chin. J. Eng. Math., 20 (2003), 13–18. https://doi.org/10.3969/j.issn.1005-3085.2003.04.003 doi: 10.3969/j.issn.1005-3085.2003.04.003
![]() |
[25] |
P. Flores, M. Machado, E. Seabra, T. Miguel, A parametric study on the baumgarte stabilization method for forward dynamics of constrained multibody systems, J. Comput. Nonlinear Dyn., 6 (2009), 011019. https://doi.org/10.1115/1.4002338 doi: 10.1115/1.4002338
![]() |
[26] |
K. T. Wehage, R. A. Wehage, B. Ravani, Generalized coordinate partitioning for complex mechanisms based on kinematic substructuring, Mech. Mach. Theory, 92 (2015), 464–483. https://doi.org/10.1016/j.mechmachtheory.2015.06.006 doi: 10.1016/j.mechmachtheory.2015.06.006
![]() |
[27] |
P. Fisette, B. Vaneghem, Numerical integration of multibody system dynamic equations using the coordinate partitioning method in an implicit Newmark scheme, Comput. Methods Appl. Mech. Eng., 135 (1996), 85–105. https://doi.org/10.1016/0045-7825(95)00926-4 doi: 10.1016/0045-7825(95)00926-4
![]() |
[28] | E. J. Haug, J. Yen, Generalized coordinate partitioning methods for numerical integration of differential-algebraic equations of dynamics, In: Real-time integration methods for mechanical system simulation, Springer, 1991. https://doi.org/10.1007/978-3-642-76159-1_5 |
[29] |
R. Singh, P. Likins, Singular value decomposition for constrained dynamical systems, J. Appl. Mech., 52 (1985), 943–948. https://doi.org/10.1115/1.3169173 doi: 10.1115/1.3169173
![]() |
[30] |
S. Kim, M. Vanderploeg, QR decomposition for state space representation of constrained mechanical dynamic systems, J. Mech. Design, 108 (1986), 183–188. https://doi.org/10.1115/1.3260800 doi: 10.1115/1.3260800
![]() |
[31] |
Q. Yu, J. Hong, A new violation correction method for constraint multibody systems, Chin. J. Theor. Appl., 30 (1998), 300–306. https://doi.org/10.6052/0459-1879-1998-3-1995-130 doi: 10.6052/0459-1879-1998-3-1995-130
![]() |
[32] |
G. Lyu, R. Liu, Errors control of constraint violation in dynamical simulation for constrained mechanical systems, J. Comput. Nonlinear Dyn., 14 (2019). https://doi.org/10.1115/1.4042493 doi: 10.1115/1.4042493
![]() |
[33] |
X. Xu, J. Luo, Z. Wu, Extending the modified inertia representation to constrained rigid multibody systems, J. Appl. Mech., 87 (2020), 011010. https://doi.org/10.1115/1.4045001 doi: 10.1115/1.4045001
![]() |
[34] |
J. Zhang, D. Liu, Y. Liu, A constraint violation suppressing formulation for spatial multibody dynamics with singular mass matrix, Multibody Syst Dyn., 36 (2016), 87–110. https://doi.org/10.1007/s11044-015-9458-7 doi: 10.1007/s11044-015-9458-7
![]() |
[35] |
L. Zhang, X. Rui, J. Zhang, J. Gu, H. Zheng, T. Li, Study on transfer matrix method for the planar multibody system with closed-loops, J. Comput. Nonlinear Dyn., 16 (2021). https://doi.org/10.1115/1.4052433 doi: 10.1115/1.4052433
![]() |
[36] |
S. Yoon, R. M. Howe, D. T. Greenwood, Geometric elimination of constraint violations in numerical simulation of lagrangian equations, J. Mech. Design, 116 (1994), 1058–1064. https://doi.org/10.1115/1.2919487 doi: 10.1115/1.2919487
![]() |
[37] |
Y. Q, I. M. Chen, A direct violation correction method in numerical simulation of constrained multibody systems, Comput. Mech., 26 (2000), 52–57. https://doi.org/10.1007/s004660000149 doi: 10.1007/s004660000149
![]() |
[38] | J. Hong, Computational multibody system dynamics, Beijing: Higher Education Press, 1999. |
[39] | D. Negrut, A. Dyer, Adams/solver primer, MSC Software Ann Arbor, 2004. |
[40] |
J. Wittenburg, Dynamics of mulitibody systems-a brief review, Space Humanity, 1989, 89–92. https://doi.org/10.1016/B978-0-08-037877-0.50015-6 doi: 10.1016/B978-0-08-037877-0.50015-6
![]() |
[41] | R. James, The unified modeling language reference manual, Addison-Wesley Professional, 2006. |
[42] |
F. Liu, J. Zhang, Q. Hu, A modified constraint force algorithm for flexible multibody dynamics with loop constraints, Nonlinear Dyn., 90 (2017), 1885–1906. https://doi.org/10.1007/s11071-017-3770-0 doi: 10.1007/s11071-017-3770-0
![]() |
[43] |
H. Lu, X. Rui, Y. Ding, Y. Chang, Y. Chen, J. Ding, X. Zhang, A hybrid numerical method for vibration analysis of linear multibody systems with flexible components, Appl. Math. Model., 101 (2022), 748–771. https://doi.org/10.1016/j.apm.2021.09.015 doi: 10.1016/j.apm.2021.09.015
![]() |
[44] |
Y. Lu, Z. Chang, Y. Lu, Y. Wang, Development and kinematics/statics analysis of rigid-flexible-soft hybrid finger mechanism with standard force sensor, Robot. Comput. Integr. Manuf., 67 (2021), 101978. https://doi.org/10.1016/j.rcim.2020.101978 doi: 10.1016/j.rcim.2020.101978
![]() |
[45] |
J. Zhang, X. Rui, F. Liu, Q. Zhou, L. Gu, Substructuring technique for dynamics analysis of flexible beams with large deformation, J. Shanghai Jiaotong Univ., 22 (2017), 562–569. https://doi.org/10.1007/s12204-017-1875-8 doi: 10.1007/s12204-017-1875-8
![]() |
[46] |
A. E. Nabawy, A. A. Abdelrahman, W. S. Abdalla, A. M. Abdelhaleem, S. S. Alieldin, Analysis of the dynamic behavior of the double wishbone suspension system, Int. J. Appl. Mech., 11 (2019), 1950044. https://doi.org/10.1142/S1758825119500443 doi: 10.1142/S1758825119500443
![]() |
[47] |
B. Zhang, Z. Li, Mathematical modeling and nonlinear analysis of stiffness of double wishbone independent suspension, J. Mech. Sci. Technol., 35 (2021), 5351–5357. https://doi.org/10.1007/s12206-021-1107-x doi: 10.1007/s12206-021-1107-x
![]() |
1. | Xizhe Zhang, Xiaoting Rui, Jianshu Zhang, Lina Zhang, Junjie Gu, Use of a path-following method for finding static equilibria of multibody systems modeled by the reduced transfer matrix method, 2024, 1384-5640, 10.1007/s11044-024-09996-y |
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 19.036 | 19.766 | 24.875 |
100 | 38.017 | 38.531 | 44.547 |
200 | 79.593 | 80.985 | 87.516 |
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 73.125 | 73.469 | 88.593 |
100 | 146.969 | 148.625 | 177.89 |
200 | 296.516 | 311.297 | 371.703 |
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 34.922 | 35.641 | 42.968 |
100 | 75 | 78.235 | 96.391 |
200 | 144.75 | 151.297 | 181.25 |
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 19.036 | 19.766 | 24.875 |
100 | 38.017 | 38.531 | 44.547 |
200 | 79.593 | 80.985 | 87.516 |
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 73.125 | 73.469 | 88.593 |
100 | 146.969 | 148.625 | 177.89 |
200 | 296.516 | 311.297 | 371.703 |
Simulation time (s) | CPU time (s) | ||
Without stabilization | Proposed approach | Direct differentiation | |
50 | 34.922 | 35.641 | 42.968 |
100 | 75 | 78.235 | 96.391 |
200 | 144.75 | 151.297 | 181.25 |