Dynamic modeling and performance analysis of a lower-mobility parallel robot with a rotatable platform

1.   Introduction

Parallel robots possess features of high stiffness, high load-bearing capacity, low cumulative error, and high motion accuracy ^[1]. They have been widely used in processing and manufacturing, aerospace, and medical fields ^{[2,3,4,5,6,7]}. However, they also have some disadvantages, such as smaller working spaces, more singularities, and difficulty in dynamic solutions. These limitations have restricted the development of parallel robots to some extent.

Due to their structural features, parallel robots have several singularities. Near the singularity, the robot's performance sharply drops. In more critical cases, it may lead to control loss ^[1]. Therefore, when designing and applying parallel robots, singular dislocations should be avoided, especially type Ⅱ singularities ^[8]. The redundancy method is a widely-used method that permits the avoidance of such singularities ^{[9,10,11,12,13]}. This method is further divided into the kinematic ^[9,10] and the drive redundancy methods ^[11,12,13]. Compared to the prior, the latter is more operational because it does not introduce a new DOF in the elimination process.

With the industry's rapid development, demand for the dynamic performance of parallel robots has been increasing. To obtain superior dynamic characteristics, parallel robots are often designed to be light. However, their flexible components are prone to elastic deformations during high-speed motion, potentially causing motion errors in the moving platform in addition to elastic vibration or resonance in the mechanism. This consequently results in the parallel mechanism's kinematic and dynamic performance degradation. Therefore, it is of the utmost importance to consider the elastic deformation of components during the modeling process, in addition to establishing a flexible model. To simplify this study, we discretized the components and used the following main methods: 1) the Finite Element Method (FEM) ^[14,15]; 2) the Assumed Mode Method (AMM) ^[16,17]; 3) the finite section method ^[18]; and 4) the lumped parameter method ^[19]. The assumed mode method could be further combined with the modal truncation technique to truncate the dominant low-order modes. After that, we can use the LaGrangian formulation and Kane's equation to establish the system's dynamic equations with lower dimensions. Therefore, the assumed mode method is suitable for the construction of active control strategies for flexible robots.

Academic research on the treatment of rigid-flexible kinematic relationships of flexible bodies has led to the development of various modeling methods. The most widely used of these are the following: 1) The Kinematic Elastic Dynamics (KED); 2) the Floating Frame of Reference (FFRF); and 3) the Absolute Nodal Coordinate Formulation (ANCF). The floating frame of reference proved to be a reliable method of establishing the rigid-flexible coupled dynamics model.

Recently, there have been some studies on the dynamic modeling of flexible parallel robots. For example, Liang Dong et al ^[20] studied a multi-drive 2 DOF parallel robot and established a systematic rigid-flexible coupling dynamics model by combining the assumed mode method with the Augmented Lagrange method. They also conducted a study on the active control strategy. İDer, S. K et al ^[21] established a dynamic model for flexible joints and studied the inverse dynamic control of flexible joint parallel robots. Sheng et al ^[22] developed a dynamics model for a 3-RRR parallel robot with a flexible intermediate linkage and investigated the vibration features. Cammarata et al ^[23] proposed an MSA/CMS framework to describe the complete and simplified model of fully flexible parallel robots' dynamic elasticity.

We present in this paper a novel 4-RRR redundantly-driven parallel robot. First, we discretized the flexibility of the system's arbitrary rod, j, based on the assumed mode method. Then, we established the flexible rod's dynamic equations. Using the Lagrangian formulation, we established the end-effector's rigid platform's dynamic equations. We also categorized the eight rods into four branch chains to establish the open-chain system's dynamic equations. We established the complete system's dynamic equations on the system's closed-loop vector equations and transformed them into the closed-loop system dynamic equations in the form of independent generalized system coordinates.

We regarded the robot's eight rods as flexible bodies while the end-effector platform as a rigid body. They thus formed a rigid-flexible coupled system and were numerically simulated. In these simulations, the dynamic performance was discussed under 3 types of drive-force distribution: a) a non-redundant drive; b) a redundant drive with minimum input torque; and c) a redundant drive with minimum input energy. Based on the inverted rigid body dynamics model, we respectively calculated three sets of corresponding driving moments. Taking the previously mentioned modes' driving moments as feedforward inputs, numerical calculations were conducted in MATLAB to obtain the system's dynamic response. To provide a reference for control and strategy, we then compared and analyzed the dynamic performance differences under these three modes.

We developed in this paper a dynamic model of a rigid-flexible coupled parallel manipulator with a rigid end-effector. Moreover, we analyzed the model's dynamic performance under three driving modes. The modular modeling method was conducive to the generalization to other models; additionally, it was helpful for the spatial structure's modeling. The established rigid-flexible coupling system dynamics model has proved to be beneficial for the optimization of light design and dynamic performance of parallel robots. Moreover, it is conducive to the development of subsequent control strategies.

2.   Model

We present in this paper a planar 3 DOF parallel robot design as shown in Figure 1. It contains a total of four branch chains, each of which contains two rods and a square end-effector platform. The four branch chains are connected to the end-effector platform's four corners. Among all the branch chains, the rods labeled 1–4 are active while the 5–8 ones are driven. We placed the robot in the inertial coordinate system $O - x - y$. Let ${A_i}$ denote the solid joint coordinates between the i^th branch chain and the inertial coordinate system, its coordinates denoted as $({x_{{A_i}}}, {y_{{A_i}}})$, and its angle as ${\phi _i}(i = 1, \cdots, 4)$. Let ${B_i}$ denote the position of the joint between the two rods of the i^th branch chain, its coordinate denoted as $({x_{{B_i}}}, {y_{{B_i}}})$, and its angle as ${\phi _i}(i = 5, \cdots, 8)$. Let ${C_i}$ denote the position of the joint at the end of the follower rod of the i^th branch chain, its coordinate denoted as $({x_{{C_i}}}, {y_{{C_i}}})$. The operating platform's center point's position is denoted as P, having coordinates of $({x_P}, {y_P})$. We denote the angle between the line of the center point P, ${C_i}$, and the $O - x$ axis as ${\phi _i}$.

3.   The components' dynamics model

3.1.   The flexible rod j's dynamic modeling

3.1.1.   The discretization of an arbitrary flexible rod j

According to the floating coordinate frame method, the motion of an arbitrary flexible rod j can be decomposed into two parts. One is a large range of rigid bodies solidly connected to the floating coordinate system and the other is a small range of elastic deformation motions relative to it. Thus, we selected the arbitrary rod's near-frame endpoint as the origin of its floating coordinate system. The x-direction is the rod's axial direction while the y-direction is its tangential one. At this point, we simulated the deformation of the flexible rod j as the cantilever beam's deformation. Figure 2 shows this deformation while moving.

In Figure 2, $O - x - y$ represents the global coordinate system. The origin O is fixed in the plane of motion. ${O_j} - {x_j} - {y_j}$ represents the floating coordinate system fixed on the flexible rod j. The position vector of the origin ${O_j}$ in the overall coordinate system is ${{\boldsymbol{r}}_{{O_j}}} = {\left[{\begin{array}{*{20}{c}} {{x_{{O_j}}}} & {{y_{{O_j}}}} \end{array}} \right]^T}$. Before the flexible rod's motion deformation, ${P_j}$ is a point on it and its position vector in the floating coordinate system ${O_j} - {x_j} - {y_j}$ is ${{\boldsymbol{r}}_{O, {P_j}}} = {\left[{\begin{array}{*{20}{c}} {\bar x} & 0 \end{array}} \right]^T}$. After motion deformation, the deformation vector diameter of the point ${P_j}$ in the floating coordinate system becomes ${{\boldsymbol{\delta }}_{{P_j}}} = {\left[{\begin{array}{*{20}{c}} {{v_{{P_j}}}} & {{w_{{P_j}}}} \end{array}} \right]^T}$. The position vector of the point ${P_j}$ in the overall coordinate system is denoted as ${{\boldsymbol{r}}_{{P_j}}}$. According to Fig. 2, the following position-vector relationship can be obtained as follows:

In the equation, ${\boldsymbol{R}}({\phi _j}) = \left[{\begin{array}{*{20}{c}} {\cos {\phi _j}} & { - \sin {\phi _j}} \\ {\sin {\phi _j}} & {\cos {\phi _j}} \end{array}} \right]$ denotes the rotation transformation matrix of the floating coordinate system ${O_j} - {x_j} - {y_j}$, in which the flexible rod j is located with respect to the overall coordinate system $O - x - y$, where ${\phi _j}$ is the relative angle of rotation.

The following equation can be obtained using the assumed mode method to discretize the flexible rod j:

In the equation, ${v_{{P_j}}}(\bar x, t)$ and ${w_{{P_j}}}(\bar x, t)$ denote the axial and transverse elastic deformations on the flexible rod j at a distance $\bar x$ from the origin under the floating coordinate system ${O_j} - {x_j} - {y_j}$ at the moment of motion t, respectively.

${{\boldsymbol{\psi }}_{{v_j}, {k_v}}}(\bar x)$ denotes the ${k_v}$ order axial modal shape function.

${{\boldsymbol{\psi }}_{{w_j}, {k_w}}}(\bar x)$ denotes the transverse modal shape function of the ${k_w}$ order.

${q_{{v_j}, {k_v}}}(t)$ denotes the axial modal time-varying coordinates.

${q_{{w_j}, {k_v}}}(t)$ denotes the transverse modal time-varying coordinates.

${k_v}$ and ${k_w}$ denote axial and transverse modal orders, respectively.

${{\boldsymbol{\psi }}_{{v_j}}}(\bar x)$ denotes the axial mode shape function matrices.

${{\boldsymbol{\psi }}_{{w_j}}}(\bar x)$ denotes the transverse mode shape function matrices.

${{\boldsymbol{q}}_{{f_{{v_j}}}}}(t)$ denotes the column matrix of the axial modal time-varying coordinates.

${{\boldsymbol{q}}_{{f_{{w_j}}}}}(t)$ denotes the column matrix of the transverse modal time-varying coordinates.

For this paper's 4-RRR parallel robot, we can model the rods as a cantilever beam model; i.e.: with fixed-free boundary conditions. Based on structural dynamics ^[24], the modal shape functions ${{\boldsymbol{\psi }}_{{v_j}, {k_v}}}(\bar x)$ and ${{\boldsymbol{\psi }}_{{w_j}, {k_w}}}(\bar x)$ can be obtained with specific expressions as follows:

In the equation,

After discretizing the flexible rod j in the assumed mode method, we can express the position vector of the arbitrary point P_j after the rod's kinematic deformation in the overall coordinate system as:

In the equation,

By taking the first order derivative of equation (4) with respect to time, we can obtain the velocity vector at the global coordinates of point P_j as:

In the equation, ${{\boldsymbol{S}}_j} = \left[{\begin{array}{*{20}{c}} {{{\boldsymbol{I}}_{2 \times 2}}} & { - {\boldsymbol{R}}({\phi _j}){\boldsymbol{\tilde K}}\left({{{\boldsymbol{r}}_{O, {P_j}}} + {{\boldsymbol{\psi }}_j}(\bar x){{\boldsymbol{q}}_{{f_j}}}} \right)} & {{\boldsymbol{R}}({\phi _j}){{\boldsymbol{\psi }}_j}(\bar x)} \end{array}} \right] \in {\mathbb{R}^{2 \times (3 + {n_v} + {n_w})}}$, ${{\boldsymbol{q}}_j} = {\left[{\begin{array}{*{20}{c}} {{\boldsymbol{r}}_{{O_i}}^T} & {{\phi _j}} & {{\boldsymbol{q}}_{{f_j}}^T} \end{array}} \right]^T} \in {\mathbb{R}^{(3 + {n_v} + {n_w}) \times 1}}$ are defined as a column matrix of the flexible rod's generalized coordinates. It contains both rigid generalized coordinates, ${{\boldsymbol{r}}_{{O_j}}}$ and ${\phi _j}$, for a wide range of variation and flexible generalized coordinates, ${{\boldsymbol{q}}_{{f_j}}}$, for a small one.

3.1.2.   The arbitrary flexible rod's kinetic energy

By integrating the flexible rod's microelements' kinetic energy, we obtained its kinetic energy at moment t as:

In the equation, $\rho$ denotes the flexible rod j's density; A denotes its cross-sectional area; ${L}_{j}$ denotes its length; ${{\boldsymbol{M}}_{t, j}}({q_j}) = \int_0^{{L_j}} {\rho A{\boldsymbol{S}}_j^T} {{\boldsymbol{S}}_J}$ denotes its mass matrix, which is a function of ${{\boldsymbol{q}}_j}$.

3.1.3.   The arbitrarily flexible rod's potential energy

In this paper, the 4-RRR parallel robot always moves in the horizontal plane, and the gravitational potential energy is constant; thus, we can nullify the latter. This implies that we only need to calculate the mechanism's elastic potential energy. Neglecting the rod's shear deformation and considering only the deformation and bending ones, we can obtain its potential energy as:

In the equation, E denotes the flexible rod's material's modulus of elasticity. I denotes the cross-sectional's area moment of inertia, and ${{\boldsymbol{K}}_{t, j}}$ denotes its stiffness matrix.

3.1.4.   The flexible rod's kinetic equations

We can express the Lagrangian formulation in the following form:

Where L is the difference between the kinetic and potential energies, $L = {T_{t, j}} - {U_{t, j}}$, ${{\boldsymbol{q}}_i}$ is the generalized coordinate, ${{\boldsymbol{\dot q}}_i}$ is the generalized velocity, and ${{\boldsymbol{\tau }}_i}$ is the generalized force corresponding to the generalized coordinate.

Upon collation, the flexible rod's kinetic equation can be obtained as:

In the equation, ${{\boldsymbol{M}}_{t, j}}({{\boldsymbol{q}}_j})$ and ${{\boldsymbol{K}}_{t, j}}$ denote the flexible rod's mass and stiffness matrices, respectively.${{\boldsymbol{F}}_{t, j}}$ denotes its generalized external moment column matrix.${{\boldsymbol{C}}_{t, j}}({{\boldsymbol{q}}_j}, {{\boldsymbol{\dot q}}_j})$ denotes its centrifugal Force/Coriolis force matrix with the following expressions:

where I denotes the unit matrix and $⨂$ the Kronecker Product symbol.

3.2.   The flexible rod's kinetic equations in the system's generalized coordinate form

For this paper's system, we can express the near-frame end joint coordinates ${{\boldsymbol{r}}_{{B_i}}}$ of the driven rod ${B_i}{C_i}$ of the arbitrary branch chain i in terms of the drive rod's angle of rotation and deformation modal coordinates, ${A_i}{B_i}$. Therefore, we chose all the rods' angles of rotation and modal coordinates as the generalized coordinates. Since the end effector is a rigid body, its center point coordinate, ${{\boldsymbol{r}}_P}$, and its rotation angle, ${\varphi _1}$, are chosen as its generalized coordinates. By grouping all flexible rods with the end effector, a total of eight flexible rods and one rigid operator are included. The whole system's generalized coordinates are denoted as ${{\boldsymbol{q}}^s}$ (non-independent coordinates) and are specified below:

To represent the flexible bar's generalized coordinates ${{\boldsymbol{q}}_j}$ in the system's generalized coordinates, it is necessary to express each quantity in ${{\boldsymbol{q}}_j}$ as ${{\boldsymbol{q}}^s}$.

The arbitrary flexible rod's angular coordinates of rotation, ${\phi _j}$ can be expressed as:

In the equation, ${{\boldsymbol{J}}_{{\phi _j}}}$ denotes the transfer matrix between the flexible rod's angular coordinates of rotation and the system's generalized coordinates, which is a Boolean matrix; i.e.: its elements are either 0 or 1. We specify it as:

Similarly, we can express the arbitrary flexible rod's array of modal coordinates, ${{\boldsymbol{q}}_{{f_j}}}$, as:

In the equation, ${{\boldsymbol{J}}_{{q_{{f_j}}}}}$ denotes the transfer matrix between the array of the arbitrary flexible rod's modal coordinates and the system's generalized coordinates. It is also a Boolean matrix specified as:

Due to the mechanism's configuration, we cannot express the flexible rod's near-frame end joint coordinates, ${{\boldsymbol{r}}_{{O_j}}}$, as a general formula. Thus, it is necessary to divide each rod's near-frame end joint coordinates into two categories expressed as a function of the system's generalized coordinates, ${{\boldsymbol{q}}^s}$, respectively:

Taking each of these two equations first order derivative with respect to time t, we can obtain:

In the equation,

According to equations (12), (14), and (17), we can express the rod's generalized coordinate column matrix's first-order derivative, ${{\boldsymbol{\dot q}}_j}$, in terms of the system's generalized coordinates' first-order derivative, ${{\boldsymbol{\dot q}}^s}$, as:

Inputting (18) into (9), we can express it in the form of the system's generalized coordinates as follows:

where ${{\boldsymbol{D}}_j}$ is a Boolean matrix specified as:

To organize equation (19) symmetrically, we multiply both ends by ${\boldsymbol{J}}_j^T$, leading to:

In the equation,

where ${\boldsymbol{M}}_{t, j}^s({\boldsymbol{q}}_j^s)$, ${\boldsymbol{C}}_{t, j}^s({\boldsymbol{q}}_j^s, {\boldsymbol{\dot q}}_j^s)$, ${\boldsymbol{K}}_{t, j}^s$, and ${\boldsymbol{F}}_{t, j}^s$ respectively denote the flexible rod's mass, centrifugal force/Coriolis force, stiffness, and generalized external moment column matrices expressed in the form of the system's generalized coordinates.

3.3.   The end effector's kinetic equations

The end effector has a square structure and high stiffness. We can consider its motion to be purely rigid. We denoted its position and posture vector by ${{\boldsymbol{X}}_P} = {\left[{\begin{array}{*{20}{c}} {{x_P}} & {{y_P}} & {{\varphi _1}} \end{array}} \right]^T}$, its position vector as ${{\boldsymbol{r}}_P} = {\left[{\begin{array}{*{20}{c}} {{x_P}} & {{y_P}} \end{array}} \right]^T} = {{\boldsymbol{S}}_{P1}}{{\boldsymbol{X}}_P}$, and its rotation angle as ${\varphi _1} = {{\boldsymbol{S}}_{P2}}{{\boldsymbol{X}}_P}$, where ${{\boldsymbol{S}}_{P1}}$ and ${{\boldsymbol{S}}_{P2}}$ are the transformation matrices expressed as: ${{\boldsymbol{S}}_{P1}} = \left[{\begin{array}{*{20}{l}} 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}} \right]$ and ${{\boldsymbol{S}}_{P2}} = \left[{\begin{array}{*{20}{c}} 0 & 0 & 1 \end{array}} \right]$.

The end effector's kinetic energy, T_p , consists of both translational and rotational kinetic energies expressed as follows:

Where ${m_P}$ denotes the end effector's mass, ${{\boldsymbol{I}}_P}$ its rotational inertia, and ${{\boldsymbol{M}}_P}$ its mass matrix, which is expressed as ${{\boldsymbol{M}}_P} = {m_P}{\boldsymbol{S}}_{P1}^T{{\boldsymbol{S}}_{P1}} + {I_P}{\boldsymbol{S}}_{P2}^T{{\boldsymbol{S}}_{P2}}$.

As shown in Figure 3, we theoretically divided the end manipulator's four jointed sections and analyzed them separately. Since the end effector constantly is in a pure rigid motion, we can set both its gravitational and elastic potential energies as zero. Moreover, its kinetic energy can be plugged into the Lagrangian formulation to obtain its kinetic equation:

In the equation,

${{\boldsymbol{F}}_P}$ denotes the column matrix of the equivalent external force on the center of mass of the end-effector and its external torque.

To represent its kinetic equations as a system in a generalized coordinate form, it is necessary to map ${{\boldsymbol{\ddot X}}_P}$ onto ${{\boldsymbol{\ddot q}}^s}$. ${{\boldsymbol{\ddot X}}_P}$ can represented as:

In the equation,

Substituting (23) into equation (22) allows us to write it as:

Multiplying both sides of (24) by ${\boldsymbol{J}}_{{{\boldsymbol{X}}_P}}^T$ gives the end effector system's dynamic equations in the generalized coordinate form:

In the equation,

Here, ${\boldsymbol{M}}_P^S$ and ${\boldsymbol{F}}_P^S$ denote the end effector's mass and generalized external moment column matrices in the form of the system's generalized coordinates, respectively.

4.   The system's kinetic equations

4.1.   The kinetic equations of open chain systems

For this paper's 4-RRR parallel robot system, we expressed all of the components' dynamic equations using the system's generalized coordinates. We can obtain the open chain system's dynamic equations by grouping all of the components:

In the equation,

${{\boldsymbol{M}}^s}({{\boldsymbol{q}}^s})$, ${{\boldsymbol{C}}^s}({{\boldsymbol{q}}^s}, {{\boldsymbol{\dot q}}^s})$, ${{\boldsymbol{K}}^s}$, and ${{\boldsymbol{F}}^s}$ respectively denote the open chain system's mass, centrifugal force/Coriolis force, stiffness, and generalized external moment column matrices.

4.2.   The kinetic equations of closed chain systems

The system dynamic equations obtained above are in the form of an open chain; i.e.: the four branch chains are hypothetically severed from the end table. To obtain the complete closed chain system dynamic equations, we must add back the closed chain constraint to the equation.

We can write four sets of closed chain constraint equations (i = 1, 2, 3, 4) based on the closed chain vector relationship:

In the equation, ${\varphi _i}$ denotes the angle between the line from P to ${C_i}$ and the $O - x$ axis; ${L_{P{C_i}}}$ denotes the length of the line from P to ${C_i}$.

Set ${{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_1}$, ${{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_2}$, ${{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_3}$, and ${{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_4}$ in a matrix to obtain the closed chain constraint matrix equation for the system:

The above equation does not explicitly include the time t for which the first-order derivative with respect to time is obtained.

In the equation, ${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_q^s({{\boldsymbol{q}}^s}, t) = \frac{{\partial {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}^s}({{\boldsymbol{q}}^s}, t)}}{{\partial {{({{\boldsymbol{q}}^s})}^T}}}$ denotes the system's closed chain-constrained Jacobi matrix, which is a sparse matrix. ${{\boldsymbol{\dot q}}^s}$ denotes the derivative of the system's generalized coordinate array.

We reintroduce the generalized binding forces, ${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_q^{sT}\lambda$, at the joint cut-offs (joints C1, C2, C3, C4) into the open chain system's dynamic equations to obtain the complete closed chain system dynamic equations:

In the equation, $\lambda \in {\mathbb{R}^{8 \times 1}}$, is a Lagrangian multiplier indicating the generalized binding force's magnitude at the joint cut-off (joints C1, C2, C3, C4). ${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_q^{sT}$ is the transpose of the system constraint Jacobi matrix, ${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_q^s({q^s}, t)$, whose column vectors indicate the generalized binding force's direction at the joint cut-off (joints C1, C2, C3, C4).

4.3.   The kinetic equations of closed chain systems written in the form of system independent generalized coordinates

The system dynamic equations obtained above contain the system generalized coordinate array, ${{\boldsymbol{q}}^s}$, with a total of $\left({11 + 8({n_v} + {n_w})} \right)$ coordinate variables that are mutually coupled and non-independent. To solve the system's dynamic equations, we must obtain them in the form of the system's independent generalized coordinates. Considering that this paper's parallel robot is in the horizontal plane and contains two translational and one rotational DOFs, we can utilize the position and posture array ${{\boldsymbol{X}}_P} = {\left[{\begin{array}{*{20}{c}} {{x_P}} & {{y_P}} & {{\varphi _1}} \end{array}} \right]^T}$ at point P as the rigid system's independent coordinates. The mechanism consists of 8 flexible rods.

All the flexible rods' modal coordinates can be used as the flexible system's independent coordinates. Therefore, we can define the system's independent generalized coordinates as:

The first 3 terms in the equation are the end effector's rigid independent coordinates and the last 8 terms are the independent modal coordinates of the 8 flexible rods.

To represent the system's non-independent generalized coordinate array, ${{\boldsymbol{q}}^s}$, as a system-independent one, ${{\boldsymbol{\bar q}}^s}$, we divide all terms in ${{\boldsymbol{q}}^s}$ into 3 classes and rewrite them as:

In this equation, ${{\boldsymbol{q}}_\phi } = {\left[{\begin{array}{*{20}{c}} {{\phi _1}} & {{\phi _2}} & \cdots & {{\phi _8}} \end{array}} \right]^T}$, ${{\boldsymbol{q}}_f} = {\left[{\begin{array}{*{20}{c}} {{\boldsymbol{q}}_{{f_1}}^T} & {{\boldsymbol{q}}_{{f_2}}^T} & \cdots & {{\boldsymbol{q}}_{{f_8}}^T} \end{array}} \right]^T}$, ${{\boldsymbol{q}}_f} = {\left[{\begin{array}{*{20}{c}} {{\boldsymbol{q}}_{{f_1}}^T} & {{\boldsymbol{q}}_{{f_2}}^T} & \cdots & {{\boldsymbol{q}}_{{f_8}}^T} \end{array}} \right]^T}$, and ${{\boldsymbol{q}}_P} = {\left[{\begin{array}{*{20}{c}} {{x_P}} & {{y_P}} & {{\varphi _1}} \end{array}} \right]^T}$ denote all the rods' angular array, modal coordinate array, and the end effector position and posture array, respectively.

Following the above description, we can write the system-independent generalized coordinates as:

Taking the first order derivative with respect to time for both ends of the system's closed chain constraint matrix equation (28) gives:

In the equation,

Reforming equation (34) gives:

Substituting the above equation into the first-order derivative ${{\boldsymbol{\dot q}}^s} \triangleq {\left[{\begin{array}{*{20}{c}} {{\boldsymbol{\dot q}}_\phi ^T} & {{\boldsymbol{\dot q}}_f^T} & {{\boldsymbol{\dot q}}_P^T} \end{array}} \right]^T}$ of the system's non-independent generalized coordinate array, we get:

In the equation,

At this point, the system's non-independent generalized coordinate array's first-order derivative ${{\boldsymbol{\dot q}}^s}$ has been expressed by the first-order derivative ${{\boldsymbol{\dot {\bar q}}}^s}$ of the system's independent generalized coordinate array. Plugging (36) in (29) gives:

Since ${{\boldsymbol{\dot {\bar q}}}^s}$ is a system independent generalized coordinate, the following consistently holds for the arbitrary ${{\boldsymbol{\dot {\bar q}}}^s}$ in the system workspace:

Hereby, we can define G as the system's closed-chain constrained Jacobi's zero-space matrix.

To eliminate the generalized constraint, ${\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_q^{sT}\lambda$, which is difficult for the closed-chain system's dynamic equations, we multiply ${{\boldsymbol{G}}^T}$ with both sides of equation (30), resulting in the following:

Finally, plugging (36) in (39) gives the closed chain system's dynamic equations in the independent generalized coordinate form:

In the equation,

5.   The dynamic response's numerical simulation

Based on the closed chain system's dynamics model in the form of system-independent generalized coordinates, we conducted a dynamics simulation for a 4-RRR parallel robot containing two parts when moving at high speed. One part is a large range of rigid motions while the other is a small range of flexible ones. Since the two influence each other, we utilized the ode23tb solver to solve the differential equations and simplify the model for efficiency. We took all flexible rods as the first order of modalities ^[25]. To compare the robot's dynamic performance under different drive torque assignments, we first proposed different optimization schemes for the distribution of the rigid body model's active joint drive torques. We then inputted them into the rigid-flexible coupled 4-RRR parallel robot. The positive kinematic solution was then obtained to determine the system's dynamic response.

The 4-RRR parallel robot has 3 DOFs and 4 active joints. If a non-redundant drive scheme is used, we can select any 3 active joints as inputs, where the input torque solution is unique and deterministic. Otherwise, if a redundant drive scheme is used, then there is an infinite number of possibilities for the allocation of input moments. The latter can be broadly divided into three categories: the first is the minimum input moment, i.e.: the minimum Euclidean norm ${F_{torque}} = {\text{min}}\sqrt {\sum\limits_{i = 1}^4 {\tau _i^2} }$ for the four active moments. The second is the minimum input energy, i.e.: the minimum Euclidean norm ${F_{power}} = {\text{min}}\sqrt {\sum\limits_{i = 1}^4 {{{\left({\tau {\phi _i}} \right)}^2}} }$ for the four active joints' input power. The third is the minimum power of some of the robot's internal forces. However, since the dynamics model's Lagrangian derivation only involves energy and not internal forces, we will not discuss this optimization scheme.

In this simulation, the planned path was the following: the end effector's center point P moves from point (0, 0) to the point (1.2mm, -17.5mm) in a spiral trajectory within 0-0.4s, while ${\varphi _1} = \frac{2}{3}\pi$ remains constant. We chose the eccentric thread because the spiral trajectory has a gradually increasing velocity and acceleration on the one hand, which facilitates dynamic analysis. On the other hand, it permits us to study the system's response under such an eccentric trajectory as an eccentric thread. The trajectory, velocity, and acceleration of point P are shown in figures 4, 5, and 6, respectively.

Table 1 demonstrates the parameters of a set of 4-RRR parallel robots. We assumed that all bars and end manipulators were of uniformly-textured aluminum and that all joints were lightweight ideal bearings; i.e.: they are massless and frictionless, etc.

Based on the previously-presented planning trajectory, we discuss the dynamic performance under 3 types of drive force distributions, namely mode a: a non-redundant drive (drive joints A1, A2, A3), mode b: a redundant drive with a minimum input moment, and mode c: a redundant drive with minimum input energy. According to the inverse rigid body dynamics model, their drive moments are calculated as shown in Figure 7.

In MATLAB, the drive moments of these drive modes were used as prior inputs for the numerical calculations of the system's dynamic response.

Figure 8 shows the end elastic deformation of rods 1–8 in the redundant drive modes. Since the axial deformation is much smaller than the transverse one, we only calculated the prior one here. To improve the calculation efficiency, the first-order mode was intercepted for all flexible rods for calculation. (1) For the active rod, the elastic deformation of the flexible rod's end in the non-redundant driving mode had the largest vibration amplitude in each small range, which is about 2–3 times bigger compared to that of the other two redundant driving modes. (2) For the active rod, the flexible rod's end's elastic deformation was in the non-redundant drive mode. The vibration amplitude in the period between (0–0.35s) was slightly larger than that of the other two redundant drive modes; the vibration amplitude in the period between (0.35–0.4s) was much larger than that of the other two redundant drive modes. (3) For all three drive modes, the driven rod's deformation was much smaller than that of the drive rod. (4) For the middle 4 driven rods, the deformation magnitudes under the 3 drive force modes did not greatly differ. (5) The redundant drive with the smallest input torque had a slightly smaller maximum deformation for rods 2, 3, and 4 compared to the redundant drive with the smallest input energy. (6) Since the robot's 1–4 rod was the active rod, its elastic deformation was affected by the driving force, demonstrating a trend fluctuation; the robot's 5–8 rod was the driven rod and its elastic deformation fluctuation was irregular. The deformation of the remaining flexible rods was equivalent. Overall, the deformation of the flexible rods in the two redundant drive modes was much smaller than that in the non-redundant drive mode. Compared to the two redundant drive modes, the redundant drive with the smallest input torque had a slightly smaller deformation of flexible rods compared to the redundant drive with the smallest input energy.

Figure 9 shows point P's actual trajectory in different drive modes. Figure 10 shows its displacement in each direction, the desired displacement error, the end table's angle of rotation, and the desired angle error in different drive modes. The redundant drive with the lowest input torque and the one with the lowest energy input have relatively small deviations compared to the non-redundant drive.

Figure 11 shows the error between the actual motion and desired angles of rods 1–8 in different drive modes. We can observe that the deformation of each flexible rod increases sharply from 0.2s in all three drive modes. Additionally, the degree of deformation is greater in the non-redundant drive mode compared to that in the remaining two redundant modes.

Figure 12 shows the actual kinematic angular velocities of rods 1–8 in different drive modes compared to the desired angular velocities. It is shown that for the active rod, the flexible rod in the non-redundant drive mode's angular velocity has the largest vibration amplitude in all small ranges, which is approximately 2–3 times greater than that of the other two redundant drive modes. Moreover, for the driven rod, the flexible rod in the non-redundant drive mode's angular velocities in some periods are slightly greater than the other two modes (for example, rod 5 at 0.35–0.4s).

Figure 13 shows the actual kinematic angular acceleration of rods 1–8 in different drive modes compared to the desired angular acceleration. The amplitude of vibration is slightly greater in some periods than in the other two modes (for example, rod 5 at 0.35–0.4s). In general, the flexible rod's angular acceleration vibration tendency is identical to that of the angular velocity, albeit at a higher frequency and with a greater amplitude.

6.   Multi-Body software simulation validation

We based the simulation validation on a parametric model of a rigid-flexible coupled system for a 4-RRR parallel robot built on Adams. It was modeled using an approach similar to finite segment discretization, in which we discretized the eight rods by dividing each one into nine rigid body subsections, each linked by a massless torsion spring. As shown in Figure 14, the end-effector platform was treated as a rigid body, i.e.: it was not divided ^[18].

We conducted the simulation in Adams using the numerical calculation's same trajectory and the simulation duration set to 0.4s. Since the workspace is in the horizontal plane, we set the gravity to 0.

Without loss of generality, the numerical calculation contained three torque distribution modes. The simulation was conducted with the lowest torque redundant drive. Figure 15 shows the trajectory displacement curve at the end console's center point P, which was obtained from the numerical calculation and the Adams simulation. Figure 16 shows the velocity profile at P. We can observe that the point's displacements and velocities obtained using the two methods are in good correspondence and have only some minor deviations. Figure 17 shows each rod's angular velocity profile obtained from the numerical calculation and the Adams simulation (taking rods 1–4 as an example). As observed, the angular velocities of each rod obtained by the two methods are in general agreement; their deviations are bound by 15%. These differences are due to the different means of discretization used by the two methods; the numerical calculation uses the assumed mode method, while the Adams simulation uses the finite section method, resulting in a slight variation in the number of cells.

7.   Conclusion

In this paper, we designed and fabricated a planar 4-RRR parallel robot by combining the assumed mode and floating coordinate frame methods to discretize an arbitrary, flexible rod j. The modal truncation method was used to intercept the rod's first-order modal state. Then, assembling the two flexible rods of an arbitrary branch chain i, one can further obtain the dynamics model of a single chain.

Additionally, we introduced the closed chain constraint to obtain its complete dynamic equations. We utilized the zero-space method that eliminates the generalized binding terms in the equations. There are seven functions in MATLAB that could be used to solve for the numerical solution of differential equations: ode45, ode23, ode113, ode15s, ode23s, ode23t, and ode23tb. Among them, ode23tb is realized by TR -- BDF2 because it is based on the implicit Runge -- Kutta formula, which is more adequate for solving rigid equations. Accordingly, we used the ode23td solver to calculate the established model and obtain the dynamic response. After that, we used the Adams software to verify the results. From that, the following could be concluded:

(1) The mechanism has planar rotational degrees of freedom, allowing it to move in the desired position and posture without being limited to a translational movement like the conventional planar parallel robots. This allows for more complex movements.

(2) We proposed an effective modular rigid-flexible coupling modeling method based on the assumed mode one. We simplified it to enhance the accuracy and computational efficiency, which is of value for models with more flexible components and larger computational volumes.

(3) Different torque distribution schemes for the drive joints were used as feedforward inputs to obtain different dynamic performances. The results show the following: firstly, the redundant drive mode can suppress the rod's elastic deformation to a certain extent and effectively improve the robot end-effector's trajectory tracking accuracy, which is an important guideline for the practical control of future robots. Secondly, the minimum input torque redundant drive has slightly less flexible rod deformation compared to the minimum input energy redundant one. However, they are almost equivalent.

Furthermore, the dynamic model of the fully flexible link system proposed in this study is expected to aid in the optimization of lightweight design and dynamic performance of parallel robots, which are expected to be used in special processing scenarios and packaging industries.

Work outlook

(1) We did not consider the effect of joint gaps in this paper's modeling, which affects the actual model's accuracy. Thus, we need to consider this in our future research.

(2) In our future work, we should further investigate the expansion of the concept of parallel robots and their theoretical analysis method's applications to spatially parallel mechanisms.

Acknowledgments

The authors would like to appreciate the editors and the reviewers for their significant contributions to the improvement of the paper. Thanks to Xi 'an Jiaotong University for their help in this study. In this paper, the MATLAB calculations and ADAMS dynamics simulations were completed at Xi 'an Jiaotong University.

Conflict of interest

The authors declare no conflict of interest.