A framework for establishing constraint Jacobian matrices of planar rigid-flexible-multibody systems

1.   Introduction

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

where $\boldsymbol 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. $\ddot{\boldsymbol{q}}$, the second time derivative of $\boldsymbol q$, denotes the generalized accelerations. $\boldsymbol{M}$ represents a generalized mass matrix. $f$ represents generalized forces on the system. $\boldsymbol \lambda$ consists of Lagrange multipliers. $\boldsymbol\varPhi$ describes the constraint equations caused by the connection relation of all joint elements in the system. $\boldsymbol \varPhi_{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

where $\boldsymbol q$, $\dot{\boldsymbol{q}}$ and $\ddot{\boldsymbol{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 $\boldsymbol{ \ddot{q}}$ for a mechanical multibody system by utilizing its state vectors $(q, \boldsymbol{\bar{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

where $\boldsymbol 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.

2.   Constraint violation correction approach

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

where $\boldsymbol c$ denotes constraints on position level and $\boldsymbol c_{\boldsymbol{q}}(\boldsymbol 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 $\boldsymbol c_{\boldsymbol{q}} = \boldsymbol{ \dot{c}}_{\boldsymbol{q}}$ exists. The correction formula can be expressed as

If the system has m constraint equations and n generalized coordinates, the constraint Jacobian matrix is given as

3.   Form of constraint equations of a secondary joint

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 ${j}$th closed loop labeled $L_{j}$ 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

where $\boldsymbol r_{\mathrm{OB}}$ denotes the position of the base point w.r.t. the inertial frame. $\boldsymbol r_{\mathrm{OM}}$ denotes the position of the moving point. $\boldsymbol {\dot{r}}_{\text {OB }}\left(\boldsymbol {\dot{r}}_{\text {OM }}\right)$ is the absolute velocity of the base (moving) point. Note that if the vector in the plane is $\boldsymbol{r} = \left[\begin{array}{ll} x & y \end{array}\right]^{\mathrm{T}}$, then its cross product matrix $\boldsymbol {\tilde{r}} = \left[\begin{array}{ll} -y & x \end{array}\right]$ is defined. Required variables $\boldsymbol{r}_{\mathrm{OM}}, \boldsymbol{r}_{\mathrm{OB}}, \dot{\boldsymbol{r}}_{\mathrm{OM}}, \dot{\boldsymbol{r}}_{\mathrm{OB}}$ 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.

4.   Intermediate variables of partial derivatives of constraints with respect to state vectors

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 $\partial \boldsymbol {c} / \partial q_{i}(i = 1 \cdots n)$ in Eq (6), there is

For convenience, two separate paths are discussed here. If element $i$ is associated with the moving point of the pin secondary joint, variables $\boldsymbol{r}_{\mathrm{OB}}$ are not related to generalized coordinate $q_i$. The partial derivative of $\boldsymbol{r}_{\mathrm{OM}}$ with respect to $q_i$ in this path is $\boldsymbol{r}_{\mathrm{OM}, 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

Similarly, on the other path, there is

Intermediate variables $\boldsymbol{r}_{\mathrm{OM}, q_i}, \boldsymbol{r}_{\mathrm{OB}, q_i}$ 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.

4.1.   Jacobian intermediate for a primary joint

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. $\boldsymbol{r}$ and $\boldsymbol{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

where subscript $\mathrm{O}$ denotes the inertial coordinate system. Subscripts $\mathrm{R}$ and $\mathrm{T}$ denote the body-fixed coordinate system with origin $\mathrm{R}$ and $\mathrm{T}$, respectively. $\boldsymbol{r}_{\mathrm{OP}}$ denotes the position of the P point w.r.t. the inertial frame. $\boldsymbol l_{\mathrm{RT}}$ describes position vector from the R point to the T point described in the body-fixed coordinate system with the R point. Cosine matrix $\boldsymbol{A}_{\mathrm{OP}}$ for the coordinate transformation from the body-fixed coordinate system with the P point to the inertial frame reads as

Generalized coordinates $q_i$ in the equation only impact the parameters of the element, namely, $\boldsymbol{l}_{\mathrm{RT}}$ and $\boldsymbol{A}_{\mathrm{RT}}$. Therefore, variables in Eq (12) varying with generalized coordinate $\boldsymbol{q}_{H_i}$ can be expressed as

where the partial derivatives with respect to the generalized coordinates are indicated as subscripts, denoted as $\partial \boldsymbol{r}_{\mathrm{OM}} / \partial q_i = \boldsymbol{r}_{\mathrm{OM}, q}$.

Variables $\boldsymbol{l}_{\mathrm{RT}, q}, \boldsymbol{A}_{\mathrm{RT}, 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

Variables in group Ⅱ can be described here. The pin joint is rotated around the z-axis by an angle $\theta_{\mathrm{RT}}$, and it has no relative slip along the z-axis, $\boldsymbol{l}_{\mathrm{RT}} = \mathit{\pmb{0}} \cdot \boldsymbol{A}_{\mathrm{RT}}$ is a direction cosine matrix corresponding to a simple rotation around the z-axis in the body-fixed frame. The corresponding $\boldsymbol{A}_{\mathrm{RT}}, \boldsymbol{l}_{\mathrm{RT}}$ can be described as

Variables $\boldsymbol{l}_{\mathrm{RT}}, \boldsymbol{A}_{\mathrm{RT}}$ with respect to generalized coordinates can be read as

The above equations are substituted into Eq (14). The kinematic quantities influenced by the generalized coordinates become

where $\boldsymbol{r}_{\mathrm{TB}} = \boldsymbol{A}_{\mathrm{OT}} \boldsymbol{l}_{\mathrm{TB}}$. Meanwhile, the relationship $\boldsymbol{A}_{\mathrm{OT}} \boldsymbol{D} \boldsymbol{A}_{\mathrm{OT}}^{\mathrm{T}} = \boldsymbol{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.

4.2.   Jacobian intermediate for a flexible body

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

where subscript f denotes a floating coordinate system. $\boldsymbol{\rho}_{\mathrm{R}}\left(\boldsymbol{\rho}_{\mathrm{T}}\right)$ denotes a deformed displacement vector from the origin of floating frame $O_{\mathrm{f}}$ to the $\mathrm{R}(\mathrm{T})$ end described in the floating frame. Variables $\boldsymbol{r}_{\mathrm{OR}}, \boldsymbol{A}_{\mathrm{OR}}, \theta_{\mathrm{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. $\boldsymbol{A}_{\mathrm{OT}}^{\mathrm{T}} \boldsymbol{r}_{\mathrm{TP}}$ described in the body-fixed frame are independent of generalized coordinates and velocities, respectively. Only quantities $\boldsymbol{A}_{\mathrm{RO}}, \boldsymbol{\rho}_{\mathrm{R}}, \boldsymbol{\rho}_{\mathrm{T}}$ and $\boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{T}}$ change with modal coordinates. Kinematic quantities for position level affected by the modal coordinates become

Observing the above equation, intermediate variables $\boldsymbol{A}_{\mathrm{RO}}, \boldsymbol{q}_{\mathrm{f}}, \boldsymbol{\rho}_{\mathrm{R}, \boldsymbol{q}_{\mathrm{f}}}, \boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{T}, \boldsymbol{q}_{\mathrm{f}}}, \boldsymbol{\rho}_{\mathrm{T}, \boldsymbol{q}_{\mathrm{f}}}$ are selected in group II. The quantities $\boldsymbol{A}_{\mathrm{RO} \mathrm{f}, \boldsymbol{q}_{\mathrm{f}}}, \boldsymbol{\rho}_{\mathrm{R}, \boldsymbol{q}_{\mathrm{f}}}, \boldsymbol{\rho}_{\mathrm{T}, \boldsymbol{q}_{\mathrm{f}}}, \boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{T}, \boldsymbol{q}_{\mathrm{f}}}$ 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

The quantities of the flexible body will be presented here, as shown in Figure 9. Radius vectors $\boldsymbol{\rho}_{\mathrm{R}}, \boldsymbol{\rho}_{\mathrm{T}}$ and angular deformation vector $\theta_{L, N_k}$ of the $N_k$ node could be described as

where $\theta_{L, N_k}$ describes the attitude of the flexible body moving in the plane. $\boldsymbol{u}_{\mathrm{O}_{\mathrm{f}} \mathrm{R}}\left(\boldsymbol{u}_{\mathrm{O}_{\mathrm{f}} \mathrm{T}}\right)$ represents the original coordinate of node $\mathrm{R}(\mathrm{T})$ in the undeformed configuration. The modal analysis and synthesis are employed to discretize deformation displacement $\boldsymbol{\varPhi}_{\boldsymbol{u}, \mathrm{R}} \boldsymbol{q}_{\mathrm{f}}\left(\boldsymbol{\varPhi}_{\boldsymbol{u}, \mathrm{~T}} \boldsymbol{q}_{\mathrm{f}}\right). \boldsymbol{\varPhi}$ is a shape function, i.e., modal function matrix in modal analysis. Subscript $u(\theta)$ denotes deformation in the displacement (angle). The directional cosine matrix of the flexible body from the body-fixed frame $o_{\mathrm{N}_k} x_{\mathrm{N}_k} y_{\mathrm{N}_k}$ to the floating frame is represented as $\boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{N}_{{k}}}$. Directional cosine matrix $\boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{N}_{{k}}}$ can be expressed as

Directional cosine matrix $A_{\mathrm{O}_{\mathrm{f}} \mathrm{N}_k}$ cannot be reduced using the same method as that in a primary joint. However, the appearance of the directional cosine matrix $\boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{N}_k}$ in the equation is inevitably multiplied by a variable independent of the modal coordinates. A mathematical relationship has been discovered, which is

The proof process is as follows. Since $\boldsymbol{u}$ is an unrelated variable, $\partial \boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{N}_k} \boldsymbol{u} / \partial t$ can be obtained by the direct derivative rule as follows:

The chain rule can be used for the same variable $\partial \boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{N}_k} \boldsymbol{u} / \partial t$, giving

Combining Eqs (25) and (26) expressed in different forms, there is

Substituting Eq (27) into the original formula, the relationships between relevant variables $\boldsymbol{A}_{\mathrm{Ro}_{\mathrm{f}}, q}$, $\boldsymbol{A}_{\mathrm{O}_{\mathrm{f}} \mathrm{T}, q}, \boldsymbol{u}$ and $\boldsymbol{q}_{\mathrm{f}}$ are

The intermediate variables in group II could be expressed as

Substituting Eq (30) into Eq (20), variables $\boldsymbol{r}_{\mathrm{OP}, q}, \boldsymbol{A}_{\mathrm{OP}, q}$ affected by state vectors of a flexible body can be read as

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.

5.   Semi-analytic solution for constraint Jacobian matrix

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:

5.1.   Composition of the system Jacobian matrix

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 $L_j$ 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 $\boldsymbol{c}_{L_j, q_{H_j}}$ for generalized coordinate $q_{H_i}$ is first demonstrated. Assume that pin joint $i$ is associated with the base point of the secondary joint. Intermediate variables $\boldsymbol{r}_{\mathrm{OB}, q}$ in group I can be expressed by variables in group Ⅱ as

Intermediate variables under the effect of the primary joint become

Constraint Jacobian matrix $\boldsymbol{c}_{L_j, q_{H_i}}$ could be expressed as

Submitting Eq (17) into the above equation, we have

In the other path, pin joint $i$ is associated with the moving point. The position of the moving point can be expressed as

Constraint Jacobian matrix $\boldsymbol{c}_{L_j, q_{H_i}}$ in this path could be expressed as

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 $L_j$ is obtained. If the closed-loop subsystem $L_j$ is without any flexible bodies, the corresponding Jacobian matrix can be written as

The assembly process can be applied for flexible bodies. Jacobian matrix $\boldsymbol c_{L_j, q_{F_i}}$ 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 $\boldsymbol{r}_{\mathrm{OB}, \boldsymbol{q}_{\mathrm{f}}}$. Intermediate quantities can be read as

According to Eq (30), the partial derivative of $\boldsymbol{r}_{\mathrm{OB}}$ with respect to the modal coordinates can be written as

According to the constraint equation of the secondary pin joint shown in Eq (7), constraint Jacobian matrix $\boldsymbol{c}_{L_j, q_{F_i}}$ reads as

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 $\boldsymbol{q}_{F_i}$ can be written as

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 $L_j$, the corresponding Jacobian matrix consists of that for each element, written as

The Jacobian matrix for all primary joints is described as $\boldsymbol{c}_{L_j, \boldsymbol q_H}$, and the Jacobian matrix for flexible bodies is described as $\boldsymbol{c}_{L_j, q_F}$. The Jacobian matrix $\boldsymbol{c}_{L_j, q}$ can be composed as

where Jacobian matrix $\boldsymbol{c}_{L_j, \boldsymbol q_H}$ is composed of the Jacobian matrix for each primary joint swept through the path, i.e., $\boldsymbol{c}_{L_j, \boldsymbol{q}_H} = \left[\begin{array}{llll}\boldsymbol{c}_{L_j, \boldsymbol{q}_{H_1}} & \boldsymbol{c}_{L_j, q_{H_2}} & \cdots & \boldsymbol{c}_{L_j, \boldsymbol{q}_{H_n}}\end{array}\right]$. Jacobian matrix $\boldsymbol{c}_{L_j, \boldsymbol{q}_F}$ is composed of the Jacobian matrix for each flexible body swept through the path, i.e., $\boldsymbol{c}_{L_j, \boldsymbol{q}_F} = \left[\begin{array}{llll}\boldsymbol{c}_{L_j, \boldsymbol{q}_{{1}}} & \boldsymbol{c}_{L_j, \boldsymbol{q}_{F_2}} & \cdots & \boldsymbol{c}_{L_j, \boldsymbol{q}_{F_n}}\end{array}\right]$.

According to the topological structure of a multibody system, Jacobian matrix $\boldsymbol{c}_q$ can be obtained by assembling the Jacobian matrices of all sub-closed loops in the system, denoted as

Then, the system Jacobian matrix can be expressed as

The specific programming procedure for the system Jacobian matrix can be described below.

5.2.   Assembling the system Jacobian matrix

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 $j$th secondary joint is solved by sweeping through the primary joint elements and flexible body elements along the path from the base point of the $j$th 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 $\boldsymbol{c}_{L_j, \boldsymbol{q}_{H_i}}$ and $c_{L_j, q_{\digamma_i}}$ are obtained. According to the location of the generalized coordinates, constraint Jacobian matrices $\boldsymbol{c}_{L_j, \boldsymbol{q}_{H_j}}$ and $\boldsymbol{c}_{L_j, \boldsymbol{q}_{f_i}}$ are composed to obtain constraint Jacobian matrix $\boldsymbol{c}_{L_j, q}$ of the $j$th secondary joint. System Jacobian matrix $\boldsymbol{c}_q$ is obtained ultimately by assembling constrained Jacobian matrices $\boldsymbol{c}_{L_{j, q}}(j = 1 \ldots 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.

6.   Direct differentiation for Jacobian matrix

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 $\boldsymbol{q}_j$ is put in the corresponding row. To obtain the solution for $\partial \boldsymbol{r}_{\mathrm{OM}} / \partial q_i$, the direct differentiation method is applied. The parameters $\partial \boldsymbol{r}_{\mathrm{OT}_G} / \partial q_i$ and $\partial A_{\mathrm{OT}_G} / \partial q_i$ 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 $\partial \boldsymbol{r}_{\mathrm{OR}_j} / \partial q_i$ and $\partial A_{\mathrm{OR}_j} / \partial q_i$ of the element $j$ can be determined. The kinematic relations of the element $j$ can be used to express $\partial \boldsymbol{r}_{\mathrm{OT}_j} / \partial q_i, \partial A_{\mathrm{OT}_j} / \partial q_i$ for the tip end, which can be written as

The process is iterated until the values of $\partial \boldsymbol{r}_{\text {OT }} / \partial q_i$ and $\partial A_{\text {OT }} / \partial q_i$ 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 $\partial \boldsymbol{r}_{\mathrm{OM}} / \partial q_i$. This process is then repeated cyclically to obtain $\partial r_{\mathrm{oM}} / \partial q_i(i = 1 \cdots a)$ for different generalized coordinates. The solution for $\partial \boldsymbol{r}_{\mathrm{OB}} / \partial q_i(i = a \cdots n)$ is obtained in a similar fashion as that of $\partial \boldsymbol{r}_{\mathrm{OM}} / \partial q_i(i = 1 \cdots 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.

7.   Constraint violation correction approach

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.

8.   Application example using the approach

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 $q_1 = q_3 = \pi / 3, q_2 = -\pi / 3, q_4 = q_4 = 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.

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.

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.

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.

9.   Conclusions

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.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

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).

Conflict of interest

The authors declare that they have no conflict of interest.