Multibody Dynamics of Nonsymmetric Planar 3PRR Parallel Manipulator with Fully Flexible Links

Abstract

This paper presents the implementation of the floating frame of reference formulation to model the flexible multibody dynamics of a nonsymmetric planar 3PRR parallel manipulator. All of the links, including the moving platform, of the manipulator under study are assumed flexible whereas the joints are assumed rigid. Using the Euler-Bernoulli beam, the flexibility of the links is modeled by using the Rayleigh-Ritz and finite element approximations. In both approximations, fixed-free boundary conditions are applied to the elastic coordinates of the links. These boundary conditions enable the evaluation of the elastic displacement at a link tip coincident with the end-effector of the manipulator which is of interest in the high precision robotics application. Both the approximations were compared by applying two different types of loads to the manipulator. It is shown that the elastic displacements obtained by using both the approximations have an agreement with a slight difference in the magnitude. In addition, the sensitivity analysis shows that the rigidity of the manipulator is much affected by the in-plane depth of the manipulator links’ cross section.

1. Introduction

No machine components are perfectly rigid. Flexibility exists in the machine components regardless of the magnitude. Components with high rigidity are commonly assumed to be rigid. A flexible body assumption becomes significant for components with low rigidity, such as those with slender or thin components. Parallel manipulators (PMs) are typically built by using slender or thin components to maintain their low weight advantage while they provide an acceptable stiffness. When a parallel manipulator is used for a high-precision application such as machining, a small elastic displacement matters and should be taken care of in order to achieve the expected accuracy.

The dynamics of a flexible multibody system undergoing small elastic displacements can be represented by using either an analytical (exact) model, discretized model or finite segment (lumped parameter) model. The analytical model employs a continuous model which only applies to very simple geometry, that is, bar or beam. The continuous model results in an infinite number of degrees of freedom which is not feasible for control purpose. The discretized models include the Rayleigh-Ritz model, sometimes called the assumed modes model (AMM), and the finite element model (FEM). The drawback of the Rayleigh-Ritz model is its inconvenience to handle a complex geometry. On the other, the finite element model can handle a complex geometry by using an appropriate element. The finite segment model represents a link as several rigid body links interconnected by springs and dampers to incorporate the flexibility of the link. A list of works utilizing the aforementioned models for the flexible multibody dynamics of manipulators can be found in References [1,2].

The finite element model of a flexible multibody dynamics can be categorized broadly into two approaches—the incremental finite element modeling [3,4] and the floating frame of reference formulation (FFRF) [5,6,7,8]. In the former approach, the shape function can only describe a small (infinitesimal) rotation and a convected coordinate system is used in order to deal with the finite rigid body motion. Accordingly, this approach involves a linearization of the equations of motion. This approach cannot provide an exact modeling of the rigid body dynamics, that is, exact rigid body inertia, in the presence of rigid body rotation. Furthermore, this approach does not guarantee zero strains in an arbitrary rigid body motion [9]. These drawbacks are overcome by the latter approach in which the displacements are presented as reference (rigid) and elastic coordinates. In the latter approach, which is used in this work, the elastic displacements are defined in floating frames attached to the bodies. To avoid the linearization of the equations of motion, an intermediate element coordinate system is introduced in this approach.

When the FFRF is used in a constrained flexible multibody dynamics, the constraints are defined for the rigid and elastic components. While the constraints for the rigid components are formulated in an identical manner to those in the rigid multibody dynamics, the constraints for the elastic components can be defined by using a certain boundary condition. Using the assumed modes model, the pinned-free [10,11,12,13], pinned-pinned [14,15], fixed-pinned [16] and fixed-free [17] boundary conditions have been used in symmetric planar PMs. Using the finite element model based on the FFRF or FEM-FFRF for short, Zhang et al. [18] imposed the pinned-free boundary condition to the flexible links of a symmetric planar 3RRR PM while Cammarata and Pappalardo [19] applied pinned-rolled, pinned-pinned and fixed-pinned boundary conditions to a similar PM. Using FEM-FFRF, Wang and Mills [20] applied both pinned-free and fixed-free boundary conditions to the flexible links of a symmetric 3PRR PM. When the finite element model is used in an incremental approach, such boundary conditions are not used [21,22,23].

Table 1. Previous works studying the flexible body dynamics of parallel manipulators (PMs) based on floating frame of reference formulation (FFRF).

Reference	PM Topology	Symmetry	Moving Platform	Boundary Condition	Method
Kang and Mills [10]	planar 3PRR	Yes	rigid	Pinned-free RR	AMM
Zhang et al. [11]	planar 3PRR	Yes	rigid	Pinned-free RR	AMM
Zhang et al. [12]	planar 3PRR	Yes	rigid	Pinned-free RR	AMM
Sheng et al. [13]	planar 3RRR	Yes	rigid	Pinned-free RR	AMM
Zhang et al. [14]	planar 3PRR	Yes	rigid	Pinned-pinned RR	AMM
Zhang et al. [15]	planar 3PRR	Yes	rigid	Pinned-pinned RR	AMM
Firoozabadi et al. [16]	planar 3RPR	Yes	rigid	Fixed-pinned PR	AMM
Ebrahimi and Firoozabadi [17]	planar 3RPR	Yes	rigid	Fixed-free PR	AMM
Zhang et al. [18]	planar 3RRR	Yes	rigid	Fixed-free RR	FEM
Cammarata and Pappalardo [19]	planar 3RRR	Yes	rigid	Pinned-rolled RR, pinned-pinned RR, fixed-pinned RR	FEM
Wang and Mills [20]	planar 3PRR	Yes	rigid	Pinned-free RR, fixed-free RR	FEM

summarizes the aforementioned use of various boundary conditions. It can be noticed from the table that all the investigated PMs are symmetric, which means that all the legs are identical and configured in a symmetric manner and all the studies assumed that the moving platform is rigid. The flexible links under the studies are either RR links (which means that the links have revolute joints at both ends) or PR links (which means that the links have prismatic joint at one end and revolute joint at another end). While the flexible PR links were modeled by applying fixed-pinned and fixed-free boundary conditions, the flexible RR links were modeled by using various boundary conditions namely pinned-free, pinned-rolled, pinned-pinned, fixed-free and fixed-pinned boundary conditions.

This work investigates the flexible multibody dynamics of a nonsymmetric planar 3PRR PM by using FFRF [7,8] with the following distinguished aspects:

• While the previous works as summarized in Table 1 studied symmetric PMs, this work investigates a nonsymmetric planar 3PRR manipulator previously discussed in References [24,25]. In the previous works with symmetric topologies, only one branch, that is, one leg, can be used to evaluate the dynamics of the manipulator. As this work deals with a non-symmetric PM in which all the three legs are not identical and they are not configured in a symmetric manner, all the legs should be included in the dynamics analysis of the manipulator.
• While the previous works as summarized in Table 1 assumed that the moving platform is rigid whereas another work such as [26] assumed that only the moving platform is flexible while the legs are rigid, this work assumes that the moving platform and all the legs are flexible. Hence, the links are called fully flexible.
• While the previous works as summarized in Table 1 used only a single method, that is, either AMM or FEM, this work uses both the Rayleigh-Ritz and finite element methods in a comparative study. It should be noticed that the finite element model used here is based on the FFRF, not the incremental method.
• While many existing works evaluated the elastic displacements of an intermediate or middle points of the links, this work evaluates the elastic displacement of the end-effector of the manipulator due to its significance in robotics application.

In a nutshell, this work has two-fold goals. First, it aims at evaluating the pose of a link tip of a nonsymmetric planar 3PRR PM, which represents the end-effector, due to the elasticity of the links. By knowing the deviation of the end-effector pose from its rigid body pose due to the link elasticity, a compensation can be performed to increase the accuracy of the manipulator motion. This is crucial for a high precision application such as machining. Second, it aims at comparing the elastic displacements of the manipulator at hand obtained by Rayleigh-Ritz and finite element methods based on FFRF. A simple geometry of the manipulator which can be modeled by beams is used to enable the utilization of both the methods. After an acceptable agreement between the two methods is verified in this study, the use of only finite element method with 3D elements for a similar manipulator with a complex geometry, in which Rayleigh-Ritz method is no more applicable, can be performed in the future. For this future work, a nonsymmetric planar 3PRR manipulator with a complex geometry which will be described in Section 3 was built for machining application.

For convenience, the formulation of FFRF using the Rayleigh-Ritz and finite element methods will be presented here briefly. First, the kinematics of a general flexible multibody system is presented. Second, the components of the flexible body dynamics which include the inertia matrix, the quadratic velocity vector, the external force vector and the elastic force vector are described. These components should be defined for each body and subsequently an assembly should be performed to obtain the components at the system level. For simplicity, friction forces are not taken into consideration. Since the system is constrained, the constraint equations should also be defined. Finally, by incorporating the constraint equations into the system equations of motion, a constrained flexible body dynamics for the system is composed. In order to solve the obtained flexible body dynamics, a set of boundary conditions should be imposed to the system. In this work, the elastic displacements of all of the flexible RR links in the manipulator at hand are solved by applying a fixed-free boundary condition to the links. The links are modeled by using Euler-Bernoulli beam model. The limitation of this link model is its lack of applicability to complex geometry. In a case where links with complex geometry should be used, one should apply another model such as finite element model. However, the finite element model for complex geometry typically results in a huge number of degrees of freedom (DOFs). For a pure analysis, this is acceptable. But for control purpose, the computation involving the huge number of DOFs is intensive while a control should be performed in real time. In this case, either a model which provides a reasonably small number of DOFs should be used or a model order reduction should be performed to a model which provides the huge number of DOFs to obtain a smaller number of DOFs.

2. Floating Frame of Reference Formulation

2.1. Kinematics

Referring to Figure 1, the pose of the body i in the global frame is given by $q^i$ which consists of the rigid coordinates $q_r{}^i$ and the flexible coordinates $q_f{}^i$:

$$q^i=\left[\begin{array}{c}{q}_r{}^i\\ q_f{}^i\end{array}\right]=\left[\begin{array}{c}{r}_o{}^i\\ {\theta }^i\\ q_f{}^i\end{array}\right],$$

(1)

where r_oⁱ is the position vector of Oⁱ with respect to the global frame whereas θⁱ is the orientation vector of the body i with respect to the global frame.

Point P in Figure 1 denotes the undeformed point whereas point P′ denotes the deformed point. The position of the deformed point P′ in the body i frame, uⁱ, is the vectorial sum of the position of the undeformed point P in the body frame i, u_oⁱ and the deformation vector u_f ⁱ:

$$u^i=u_o{}^i+u_f{}^i\mathrm{where}{u}_f{}^i=S^i{q}_f{}^i.$$

(2)

The matrix S is called the shape matrix whereas q_f is the vector of elastic coordinates which combines the flexible coordinates of all bodies in a multibody system, that is,

$$q_f={\left[\begin{array}{cccc}{q}_f{}^1& q_f{}^2& \cdots & q_f{}^n\end{array}\right]}^T,$$

(3)

where n is the number of bodies in the multibody system. The matrix S is space dependent whereas the vector q_f is time dependent.

The position vector of the deformed point P′ in the global frame is given by:

$$r_{P^{\prime }}^i=r_o{}^i+R^i\left(u_o{}^i+S^i{q}_f{}^i\right),$$

(4)

where Rⁱ is the rotation matrix from the body frame i to the global frame.

The velocity vector of the deformed point P′ with respect to the global frame is obtained by differentiating (4) with respect to time, that is,

$${\dot{r}}_{P^{\prime }}^i={\dot{r}}_o{}^i+{\dot{R}}^i\left(u_o{}^i+S^i{q}_f{}^i\right)+R^i{S}^i{\dot{q}}_f{}^i,$$

(5)

where a dot above a letter denotes differentiation with respect to time. Using the transport theorem, the time derivative of the rotation matrix Rⁱ can be omitted in the velocity expression. Accordingly, the velocity vector of the deformed point P′ with respect to the global frame can be written in terms of the angular velocity vector ${\omega }^i$ as follows:

$${\dot{r}}_{P^{\prime }}^i={\dot{r}}_o{}^i+{\tilde{u}}^i{}^T{\omega }^i+R^i{S}^i{\dot{q}}_f{}^i,$$

(6)

where ${\tilde{u}}^i$ is the skew-symmetric matrix of the vector uⁱ.

Upon differentiating Equation (6) with respect to time, the acceleration vector of the deformed point P′ with respect to the global frame can be obtained by utilizing the transport theorem which avoids the time derivative of the rotation matrix Rⁱ:

$${\ddot{r}}_{P^{\prime }}^i={\ddot{r}}_o{}^i+\left({\omega }^i\times \left({\omega }^i\times u^i\right)\right)+\left({\alpha }^i\times u^i\right)+\left(2{\omega }^i\times \left(R^i{\dot{u}}^i\right)\right)+R^i{\ddot{u}}^i,$$

(7)

where a double dot above a letter denotes double differentiation with respect to time and ${\alpha }^i$ denotes the angular acceleration vector of body i.

2.2. Inertia

The inertia matrix of body i can be derived from its kinetic energy expression:

$$T^i=\frac{1}{2}{\displaystyle \underset{V^i}{\int }{\rho }^i{\dot{r}}^i{}^T{\dot{r}}^{i}d{V}^i},$$

(8)

where ${\rho }^i$ is the mass density of body i and ${\dot{r}}^i$ denotes ${\dot{r}}_{P^{\prime }}^i$ by omitting the subscript P′ for brevity.

The inertia matrix of the body i is then given by the following symmetric matrix:

$$M^i=\left[\begin{array}{cc}{m}_{rr}{}^i& m_{rf}{}^i\\ m_{fr}{}^i& m_{ff}{}^i\end{array}\right]=\left[\begin{array}{ccc}{m}_{r_o{r}_o}{}^{{}^i}& m_{r_o\theta }{}^i& m_{r_{o}f}{}^i\\ & m_{\theta \theta }{}^i& m_{\theta f}{}^i\\ sym& & m_{ff}{}^i\end{array}\right]$$

(9)

where the elements of the inertia matrix are defined by the following volume integrals:

$$m_{r_o{r}_o}{}^i={\displaystyle \underset{V^i}{\int }{\rho }^{i}Id{V}^i}$$

(10)

$$m_{r_o\theta }{}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{B}^{i}d{V}^i};\begin{array}{c}\end{array}where\begin{array}{c}\end{array}{B}^i=-R^i{\tilde{u}}^i{G}^i\begin{array}{c}\end{array}and\begin{array}{c}\end{array}{\omega }^i=G^i{\dot{\theta }}^i$$

(11)

$$m_{r_{o}f}{}^i=R^i{\displaystyle \underset{V^i}{\int }{\rho }^i{S}^{i}d{V}^i}$$

(12)

$$m_{\theta \theta }{}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{B}^i{}^T{B}^{i}d{V}^i}$$

(13)

$$m_{\theta f}{}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{B}^i{}^T{R}^i{S}^{i}d{V}^i}$$

(14)

$$m_{ff}{}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{S}^i{}^T{S}^{i}d{V}^i}.$$

(15)

For a planar case, which is the case in this work, the integrals Equations (10)–(15) are given by the following:

$$m_{r_o{r}_o}{}^i=\left[\begin{array}{cc}{m}^i& 0\\ 0& m^i\end{array}\right]\mathrm{where}{m}^i\mathrm{is}\mathrm{the}\mathrm{mass}\mathrm{of}\mathrm{the}\mathrm{body}i$$

(16)

$$m_{r_o\theta }{}^i=R_{\theta }{}^i\left(I_1{}^i+{\overline{S}}^i{q}_f{}^i\right),\mathrm{where}{R}_{\theta }{}^i=\frac{\partial R^i}{\partial {\theta }^i}$$

(17)

$$I_1{}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{u}_o{}^{i}d{V}^i}$$

(18)

$${\overline{S}}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{S}^{i}d{V}^i}$$

(19)

$$m_{r_{o}f}{}^i=R^i{\overline{S}}^i$$

(20)

$$m_{\theta \theta }{}^i={\left(m_{\theta \theta }{}^i\right)}_{r_o{r}_o}+{\left(m_{\theta \theta }{}^i\right)}_{r_{o}f}+{\left(m_{\theta \theta }{}^i\right)}_{ff},\mathrm{where}$$

(21)

$${\left(m_{\theta \theta }{}^i\right)}_{r_o{r}_o}={\displaystyle \underset{Vi}{\int }{\rho }^i{u}_o{{}^i}^T{u}_o{}^i}d{V}^i$$

(22)

$${\left(m_{\theta \theta }{}^i\right)}_{r_{o}f}=2\left({\displaystyle \underset{Vi}{\int }{\rho }^i{u}_o{{}^i}^T{S}^{i}d{V}^i}\right)q_f{}^i$$

(23)

$${\left(m_{\theta \theta }{}^i\right)}_{ff}=q_f{{}^i}^T{m}_{ff}{q}_f{}^i$$

(24)

$$m_{ff}{}^i={\displaystyle \underset{Vi}{\int }{\rho }^i{S}^i{}^T{S}^{i}d{V}^i}$$

(25)

$$m_{\theta f}{}^i={\displaystyle \underset{V^i}{\int }{\rho }^i{u}_o{{}^i}^T\tilde{I}{S}^{i}d{V}^i}+q_f{}^i\left({\overline{S}}_{12}-{\overline{S}}_{21}\right),\mathrm{where}$$

(26)

$$\tilde{I}=\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$$

(27)

$$m_{ff}{}^i={\overline{S}}_{11}{}^i+{\overline{S}}_{22}{}^i.$$

(28)

Notice that and ${\overline{S}}_{kl}$ in Equations (26) and (28) is given by:

$${\overline{S}}_{kl}={\displaystyle \underset{Vi}{\int }{\rho }^i{S}_k{{}^i}^T{S}_l{}^{i}d{V}^i};\begin{array}{c}\end{array}k,l=1,2.$$

(29)

2.3. Force Vectors

The non-inertial forces in the constrained flexible body dynamics are provided by the quadratic velocity vector, gravitational force vector, external force vector, elastic force vector and constraint force vector.

2.3.1. Quadratic Velocity Vector

The quadratic velocity vector represents the Coriolis and centrifugal force vectors in the dynamics. For a planar case, the quadratic velocity vector of body i is:

$$Q_v{}^i=\left[\begin{array}{c}{\left(Q_v{}^i\right)}_r\\ {\left(Q_v{}^i\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left(Q_v{}^i\right)}_{r_o}\\ {\left(Q_v{}^i\right)}_{\theta }\\ {\left(Q_v{}^i\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left({\dot{\theta }}^i\right)}^2{R}^i\left(\overline{S}{q}_f{}^i+I_1{}^i\right)-2{\dot{\theta }}^i{R}_{\theta }{}^i{\overline{S}}^i{\dot{q}}_f{}^i\\ -2{\dot{\theta }}^i{\dot{q}}_f{}^{iT}\left(m_{ff}{}^i{q}_f{}^i+I_o{}^i\right)\\ {\left({\dot{\theta }}^i\right)}^2\left(m_{ff}{}^i{q}_f{}^i+I_o{}^i\right)+2{\dot{\theta }}^i{\tilde{S}}^i{\dot{q}}_f{}^i\end{array}\right]$$

(30)

The matrix $I_1{}^i$ is given by Equation (18) whereas

$$I_o{}^i={\left({\overline{I}}_{11}{}^i+{\overline{I}}_{22}{}^i\right)}^T,$$

(31)

where:

$${\overline{I}}_{kl}{}^i={\displaystyle \underset{Vi}{\int }{\rho }^i{x}_k{}^i{S}_l{}^{i}d{V}^i};\begin{array}{c}\end{array}k,l=1,2.$$

(32)

2.3.2. Gravitational Force Vector

The gravitational force vector of body i can be written as follows:

$$Q_g{}^i=\left[\begin{array}{c}{\left(Q_g{}^i\right)}_r\\ {\left(Q_g{}^i\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left(Q_g{}^i\right)}_{r_o}\\ {\left(Q_g{}^i\right)}_{\theta }\\ {\left(Q_g{}^i\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left(F_g{{}^i}^{T}I\right)}^T\\ {\left(F_g{}^{iT}{R}_{\theta }{}^i{u}^i\right)}^T\\ {\left(F_g{}^T{R}^i{S}^i\left(\xi ={\xi }_{COM}\right)\right)}^T\end{array}\right],$$

(33)

where F_g is the gravitational force vector expressed in the global frame. For a planar case, F_g on the body i can be written as:

$$F_g{{}^i}^T=\left[\begin{array}{c}0\\ -m^{i}g\end{array}\right],$$

(34)

whereas $S^i\left(\xi ={\xi }_{COM}\right)$ represents $S^i$ at a value of $\xi$ at which the center of mass of the body is located. For example, if the center of mass is at the middle of the body i, then we should use $S^i\left(\xi =0.5\right)$ as $\xi$ represents the length of the body i in the scale from 0 to 1.

2.3.3. External Force Vector

The external force vector of body i corresponds to any external force applied to a coordinate. Since an external force is typically applied to a rigid coordinate, the elements of the external force vector corresponding to the elastic coordinates are zeros. This can be written as follows:

$$Q_e{}^i=\left[\begin{array}{c}{\left(Q_e{}^i\right)}_r\\ {\left(Q_e{}^i\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left(Q_e{}^i\right)}_{r_o}\\ {\left(Q_e{}^i\right)}_{\theta }\\ {\left(Q_e{}^i\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left(Q_e{}^i\right)}_{r_o}\\ {\left(Q_e{}^i\right)}_{\theta }\\ 0\end{array}\right].$$

(35)

2.3.4. Elastic Force Vector

The elastic force vector of body i can be derived from the elastic potential energy expression:

$$U^i=\frac{1}{2}{q}_f{}^{iT}{K}_{ff}{}^i{q}_f{}^i.$$

(36)

The stiffness matrix of the body i is given by:

$$K^i=\left[\begin{array}{cc}{K}_{rr}{}^i& K_{rf}{}^i\\ K_{fr}{}^i& K_{ff}{}^i\end{array}\right]=\left[\begin{array}{ccc}{K}_{r_o{r}_o}{}^i& K_{r_o\theta }{}^i& K_{r_{o}f}{}^i\\ K_{\theta r_o}{}^i& K_{\theta \theta }{}^i& K_{\theta f}{}^i\\ K_{fr}{}^i& K_{f\theta }{}^i& K_{ff}{}^i\end{array}\right]=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& K_{ff}{}^i\end{array}\right].$$

(37)

Using Euler-Bernoulli beam model, the elastic potential energy of the body i can be written as:

$$U^i=\frac{1}{2}{\displaystyle \underset{0}{\overset{L^i}{\int }}\left(EI{\left(\frac{{\partial }^2{f}_2}{\partial x^2}\right)}^2+EA{\left(\frac{\partial f_1}{\partial x}\right)}^2\right)dx},$$

(38)

where E, I, A and L denote the modulus of elasticity, the area moment of inertia, the cross section area and the length of the beam, whereas f₂ and f₁ are the transversal and axial elastic displacements of the beam, respectively.

2.3.5. Constraint Forces

From the position constraint equations $\mathsf{\Gamma }$, we can obtain the velocity constraint equations by differentiating the position constraint equations with respect to time:

$$\frac{d}{dt}\mathsf{\Gamma }=B(q,t)\dot{q}=0.$$

(39)

Alternatively, the velocity constraint matrix B can also be obtained by deriving the position constraint equations with respect to the generalized coordinates q. For a body i:

$$B^i=\frac{\partial {\mathsf{\Gamma }}^i}{\partial q^i}.$$

(40)

Notice that in the case of flexible body dynamics, the generalized coordinates consist of not only the rigid coordinates but also the elastic coordinates. The constraint forces are given by the product between the transpose of the matrix B and the Lagrange multipliers λ.

2.4. Equations of Motion

The equations of motion (EOMs) in the body i level are given by the following differential-algebraic equations (DAEs):

$$M^i{\ddot{q}}^i+K^i{q}^i+B^i{}^T\lambda =Q_v{}^i+Q_g{}^i+Q_e{}^i$$

(41)

$$B^i{\dot{q}}^i=0.$$

(42)

The EOMs in the system level are represented by the following DAEs obtained by assembling all the EOMs corresponding to all the bodies in the system:

$$M\ddot{q}+Kq+B^T\lambda =Q_v+Q_g+Q_e$$

(43)

$$B\dot{q}=0.$$

(44)

3. Kinematics and Frame Assignment

The manipulator discussed in this paper is shown in Figure 2. It consists of a fixed base, three legs, a moving platform and several joints. The three legs and the moving platform will be referred to as bodies or links. The joints used in an order starting from the fixed base are prismatic (P) joints, revolute (R) joints which connect the sliders with the legs and other revolute (R) joints which connect the legs with the moving platform. Therefore, the manipulator is called the 3PRR PM. The prismatic joints are actuated, have translation along X-axis and are implemented by using sliders moving along the base through guideway. As a result of this configuration, the manipulator provides three DOFs namely translation in X-axis, translation in Y-axis and rotation about the Z-axis. The placement of the three sliders in a parallel manner along the X-axis implies a non-symmetric configuration of the planar 3PRR PM. The translational motion of the sliders along the X-axis does not only enable the moving platform to move upward and downward as well as to tilt but also creates more workspace. This is suitable for an application requiring an elongated workspace which cannot be provided by a symmetric configuration of a planar 3PRR PM. Furthermore, two neighboring legs of the PM at hand are put in a crossing configuration to avoid singularity as well as maximize the useful workspace and the stiffness of the manipulator.

A hybrid manipulator composed of two non-symmetric planar 3PRR PMs, as shown in Figure 3, was developed. The two similar PMs are conjugated—one on the bottom and another one on the perpendicular side. This gives a spatial workspace with a redundancy in the translation along the sliders motion direction. This translational redundancy provides a larger workspace in the sliders motion direction. Both the lower and upper platforms use planar PMs as they are easier to design than spatial PMs and more geometrical constraints on the DOFs are expected to minimize the error due to joint errors. This hybrid manipulator can be used for five-axis machining by putting a worktable on the lower module and holding a tool spindle on the upper module. To avoid a redundancy between the rotation of the spindle and the rotation of the upper module, the tool spindle is oriented horizontally. As a result, a five-axis machine is composed. The five-axis machine has three orthogonal translational degrees of freedom as well as two rotational degrees of freedom.

The flexible body dynamics of the manipulator at hand is formulated based on the frame assignment shown in Figure 4. As the problem is planar, the coordinate systems are given in a 2D space. The global frame is indicated by O-XY frame whereas four local frames corresponding to the four bodies of the manipulator are indicated Oⁱ -x₁ⁱ x₂ⁱ frames. Since the local frames are attached exactly at the joint positions, no intermediate joint frames are required. This accordingly simplifies the problem. The position of the joints is indicated by the letters A, B, C, D, E and F. These letters are used when defining the position constraint equations.

4. Boundary Conditions for the Elastic Coordinates

What meant by the boundary conditions here are some reference conditions exclusively corresponding with the elastic coordinates only. They are also commonly called the reference conditions. These boundary conditions should be defined to eliminate the redundancy in the resulting linear system so that a unique solution, that is, unique displacement, can be obtained. Hence, these boundary conditions are not defined to constrain the manipulator in order to achieve the expected rigid body motions which are already constrained by the constraint equations representing the joints in the manipulator. In other words, the elastic coordinates should be reduced through the boundary conditions in order to convert the system from a singular system to a non-singular one and therefore make it solvable. This reduction is conducted by suppressing some of the elastic coordinates based on the boundary conditions. The determination of the boundary conditions depends on the expected behavior of the elastic displacement in a multibody system. In general, different boundary conditions give different solutions. These boundary conditions include fixed-pinned, fixed-free (cantilever), pinned-pinned, pinned-rolled (simply supported) and pinned-free boundary conditions.

The fixed-pinned boundary conditions do not allow any type of elastic displacement at the initial end of the body i but allow the slope at the final end. The fixed-rolled boundary conditions do not allow any type of elastic displacement at the initial end of the body i but allow the axial elastic displacement and the slope at the final end. The fixed-free boundary conditions do not allow any type of elastic displacement at the initial end of the body i but allow all types of elastic displacement at the final end. The pinned-pinned boundary conditions imply that only slopes are allowed at the two ends of the body i. The pinned-rolled boundary conditions imply that the axial elastic displacement is allowed at the final end of the body i beside the slopes at both the ends. The pinned-free boundary conditions indicate that the slope at the initial end of the body i is allowed while the axial elastic displacement, transversal elastic displacement and slope are all allowed at the final end.

Since the translational elastic displacement at the final end of each body is of interest, at least one translational DOF at the final end should be allowed. If fixed-pinned or pinned-pinned boundary conditions are used, no translational elastic displacement is allowed at the final end. As the result, there will be no elastic displacement at the end-effector of the manipulator at hand as it is defined at a link final end. This will be useless as our interest is the elastic displacement at the end-effector. However, in some other applications where one is interested in the elastic displacement at a certain intermediate point along the link, the fixed-pinned or pinned-pinned boundary conditions can be used.

Using either the fixed-rolled or the pinned-rolled boundary conditions, the axial displacement at the final point is allowed. However, the transversal elastic displacement is not allowed. For slender links such as the links of the manipulator at hand, the transversal elastic displacement is typically much more significant than the axial displacement. Using either the fixed-free or the pinned-free boundary conditions, both the axial and transversal displacements at the final point is allowed. If the pinned-free boundary condition is used, the system is still singular, which physically means underconstrained, unless an additional constraining force is added. Based on this consideration, the fixed-free boundary conditions are used in this work. This implies that all the DOFs at the initial end are suppressed while all the DOFs at the final end are retained.

5. Flexible Body Dynamics Based on Rayleigh-Ritz Shape Matrix

Using Rayleigh-Ritz shape matrix, the elastic displacement vector of body i in a two-dimensional space is given by:

$$u_f{}^i=\left[\begin{array}{c}{u}_{f_1}{}^i\\ u_{f_2}{}^i\end{array}\right]=S_b{q}_f{}^i,$$

(45)

where $u_{f_1}{}^i$ and $u_{f_2}{}^i$ represent the axial and transversal displacement of any point between the end points of the body i. Hence, it is clear that Equation (45) is a displacement interpolation function.

The matrix S_b for a planar case is the following Rayleigh-Ritz shape matrix:

$$S_b=\left[\begin{array}{cccccc}1-\xi & 0& 0& \xi & 0& 0\\ 0& 1-3{\xi }^2+2{\xi }^3& l\left(\xi -2{\xi }^2+{\xi }^3\right)& 0& 3{\xi }^2-2{\xi }^3& l\left({\xi }^3-{\xi }^2\right)\end{array}\right],$$

(46)

whereas q_f ⁱ is the vector of elastic coordinates of the two endpoints of the body i in the 2D space:

$$q_f{}^i={\left[\begin{array}{cccccc}{q}_{f1}{}^i& q_{f2}{}^i& q_{f3}{}^i& q_{f4}{}^i& q_{f5}{}^i& q_{f6}{}^i\end{array}\right]}^T.$$

(47)

The first three coordinates of q_f ⁱ belong to one end point of the body i whereas the remaining three coordinates belong to the other end point. For each end point, the three coordinates represent the axial elastic displacement, the transversal elastic displacement and the slope about an axis perpendicular to the plane. For the manipulator at hand, all the elastic coordinates are shown in Figure 5. Notice that a coordinate subscript represents the order of the coordinate which represents the degree of freedom, whereas the superscript represents the body to which the coordinate belongs. Since each link has 6 coordinates while the manipulator has 4 links, the total number of the coordinates is 24.

They can be written as a single vector $q_f$:

$$q_f={\left[\begin{array}{cccc}{q}_f{}^1& q_f{}^2& q_f{}^3& q_f{}^4\end{array}\right]}^T,$$

(48)

where:

$$q_f{}^1=\left[\begin{array}{c}{q}_{f1}{}^1\\ q_{f2}{}^1\\ q_{f3}{}^1\\ q_{f4}{}^1\\ q_{f5}{}^1\\ q_{f6}{}^1\end{array}\right];\begin{array}{cc}& \end{array}{q}_f{}^2=\left[\begin{array}{c}{q}_{f1}{}^2\\ q_{f2}{}^2\\ q_{f3}{}^2\\ q_{f4}{}^2\\ q_{f5}{}^2\\ q_{f6}{}^2\end{array}\right];\begin{array}{cc}& \end{array}{q}_f{}^3=\left[\begin{array}{c}{q}_{f1}{}^3\\ q_{f2}{}^3\\ q_{f3}{}^3\\ q_{f4}{}^3\\ q_{f5}{}^3\\ q_{f6}{}^3\end{array}\right];\begin{array}{cc}& \end{array}{q}_f{}^4=\left[\begin{array}{c}{q}_{f1}{}^4\\ q_{f2}{}^4\\ q_{f3}{}^4\\ q_{f4}{}^4\\ q_{f5}{}^4\\ q_{f6}{}^4\end{array}\right].$$

(49)

Due to the imposed boundary conditions, the elastic coordinates of the bodies of the manipulator are reduced to the following:

$$q_f{}^1=\left[\begin{array}{c}{q}_{f4}{}^1\\ q_{f5}{}^1\\ q_{f6}{}^1\end{array}\right];\begin{array}{cc}& \end{array}{q}_f{}^2=\left[\begin{array}{c}{q}_{f4}{}^2\\ q_{f5}{}^2\\ q_{f6}{}^2\end{array}\right];\begin{array}{cc}& \end{array}{q}_f{}^3=\left[\begin{array}{c}{q}_{f4}{}^3\\ q_{f5}{}^3\\ q_{f6}{}^3\end{array}\right];\begin{array}{cc}& \end{array}{q}_f{}^4=\left[\begin{array}{c}{q}_{f4}{}^4\\ q_{f5}{}^4\\ q_{f6}{}^4\end{array}\right].$$

(50)

Accordingly, the columns of the shape matrix corresponding to the eliminated elastic coordinates are also omitted. Hence, the shape matrix for each body is reduced to the following:

$$S=\left[\begin{array}{ccc}\xi & 0& 0\\ 0& 3{\xi }^2-2{\xi }^3& L\left({\xi }^3-{\xi }^2\right)\end{array}\right].$$

(51)

Subsequently, the inertia matrix of each body is obtained by using Equations (16)–(29) whereas the quadratic velocity vector for each body is obtained by using Equations (30)–(32). The gravitational force vector is defined by assuming that the center of mass (COM) of each body is at the middle of the body, therefore the matrix Sⁱ in Equation (35) is defined by using $\xi =0.5$. The external force vector is obtained by using Equation (35) where the elements corresponding to the rigid coordinates are given by the forces or torques applied to the rigid coordinates. This is written as follows:

$${\left(Q_e\right)}_r=\left[\begin{array}{cccccccccccc}{Q}_{x_1}& 0& Q_{{\alpha }_1}& Q_{x_2}& 0& Q_{{\alpha }_2}& Q_{x_3}& 0& Q_{{\alpha }_3}& Q_x& Q_y& Q_{\theta }\end{array}\right],$$

(52)

where Q_x1, Q_x2 and Q_x3 are the actuation forces given by the linear motors whereas Q_α1, Q_α2 and Q_α3 are the external torques applied at the lower revolute joints. Without a gravity compensation applied at the lower revolute joints, Q_α1, Q_α2 and Q_α3 are zero. However, if gravity compensation torques are applied at the lower revolute joints then Q_α1, Q_α2 and Q_α3 are not zero. The elements Q_x, Q_y and Q_θ represent the external wrench applied to the moving platform. The values of Q_x, Q_y and Q_θ are zero when there is no external wrench applied to the moving platform.

To obtain the elastic force vector, first we define the vector u_f ⁱ for each body i by using the reduced shape matrix S:

$$u_f{}^i=\left[\begin{array}{c}{u}_{f_1}{}^i\\ u_{f_2}{}^i\end{array}\right]=S{q}_f{}^i.$$

(53)

Subsequently, we find the first derivative of $u_{f_1}{}^i$ and the second derivative of $u_{f_2}{}^i$ with respect to x and substitute them into Equation (38) which can also be written as follows:

$$U^i=\frac{1}{2}{\displaystyle \underset{0}{\overset{L^i}{\int }}\left(\left[\begin{array}{cc}\frac{\partial f_1}{\partial x}& \frac{{\partial }^2{f}_2}{\partial x^2}\end{array}\right]\left[\begin{array}{cc}EA& 0\\ 0& EI\end{array}\right]{\left[\begin{array}{cc}\frac{\partial f_1}{\partial x}& \frac{{\partial }^2{f}_2}{\partial x^2}\end{array}\right]}^T\right)dx}.$$

(54)

Upon the derivation, Equation (54) yields the following:

$$U^i=\frac{1}{2}{\left[\begin{array}{c}{q}_{f4}{}^i\\ q_{f5}{}^i\\ q_{f6}{}^i\end{array}\right]}^T{K}_{ff}{}^i\left[\begin{array}{c}{q}_{f4}{}^i\\ q_{f5}{}^i\\ q_{f6}{}^i\end{array}\right]$$

(55)

from which we take the matrix K_ff ⁱ and substitute it into Equation (37) to get the stiffness matrix of each body.

After obtaining all the aforementioned inertia matrices, stiffness matrices and force vectors, we assemble them in a system level as the following where subscripts r and f refer to rigid and flexible coordinates, respectively:

$$M_{rr}=\left[\begin{array}{cccc}{m}_{rr}{}^1& & & 0\\ & m_{rr}{}^2& & \\ & & m_{rr}{}^3& \\ 0& & & m_{rr}{}^4\end{array}\right]$$

(56)

$$M_{rf}=\left[\begin{array}{cccc}{m}_{rf}{}^1& & & 0\\ & m_{rf}{}^2& & \\ & & m_{rf}{}^3& \\ 0& & & m_{rf}{}^4\end{array}\right]$$

(57)

$$M_{fr}=M_{rf}{}^T$$

(58)

$$M_{ff}=\left[\begin{array}{cccc}{m}_{ff}{}^1& & & 0\\ & m_{ff}{}^2& & \\ & & m_{ff}{}^3& \\ 0& & & m_{ff}{}^4\end{array}\right]$$

(59)

$$K_{ff}=\left[\begin{array}{cccc}{K}_{ff}{}^1& & & 0\\ & K_{ff}{}^2& & \\ & & K_{ff}{}^3& \\ 0& & & K_{ff}{}^4\end{array}\right]$$

(60)

$$Q_g=\left[\begin{array}{c}{\left(Q_g\right)}_r\\ {\left(Q_g\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left[\begin{array}{cccc}{\left(Q_g{}^1\right)}_r& {\left(Q_g{}^2\right)}_r& {\left(Q_g{}^3\right)}_r& {\left(Q_g{}^4\right)}_r\end{array}\right]}^T\\ {\left[\begin{array}{cccc}{\left(Q_g{}^1\right)}_f& {\left(Q_g{}^2\right)}_f& {\left(Q_g{}^3\right)}_f& {\left(Q_g{}^4\right)}_f\end{array}\right]}^T\end{array}\right]$$

(61)

$$Q_e=\left[\begin{array}{c}{\left(Q_e\right)}_r\\ {\left(Q_e\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left[\begin{array}{cccc}{\left(Q_e{}^1\right)}_r& {\left(Q_e{}^2\right)}_r& {\left(Q_e{}^3\right)}_r& {\left(Q_e{}^4\right)}_r\end{array}\right]}^T\\ {\left[\begin{array}{cccc}{\left(Q_e{}^1\right)}_f& {\left(Q_e{}^2\right)}_f& {\left(Q_e{}^3\right)}_f& {\left(Q_e{}^4\right)}_f\end{array}\right]}^T\end{array}\right]$$

(62)

$$Q_v=\left[\begin{array}{c}{\left(Q_v\right)}_r\\ {\left(Q_v\right)}_f\end{array}\right]=\left[\begin{array}{c}{\left[\begin{array}{cccc}{\left(Q_v{}^1\right)}_r& {\left(Q_v{}^2\right)}_r& {\left(Q_v{}^3\right)}_r& {\left(Q_v{}^4\right)}_r\end{array}\right]}^T\\ {\left[\begin{array}{cccc}{\left(Q_v{}^1\right)}_f& {\left(Q_v{}^2\right)}_f& {\left(Q_v{}^3\right)}_f& {\left(Q_v{}^4\right)}_f\end{array}\right]}^T\end{array}\right].$$

(63)

The position constraints of the flexible body model of the manipulator are the following:

$${\mathsf{\Gamma }}_1:r_{o_1}{}_{,y}=0$$

(64)

$${\mathsf{\Gamma }}_2:r_{o_2}{}_{,y}=0$$

(65)

$${\mathsf{\Gamma }}_3:r_{o_3}{}_{,y}=0.$$

(66)

The constraints Equations (64)–(66) implies that all the rigid sliders always slide on the X-axis without any displacement in the y direction. The other position constraints are corresponding to the flexible links:

$${\mathsf{\Gamma }}_4,{\mathsf{\Gamma }}_5:r_{o_4}=r_{o_2}+R^2{u}_D$$

(67)

$$u_D{}^2=u_{o,D}{}^2+u_{f,D}{}^2$$

(68)

$$u_{o,D}{}^2={\left[\begin{array}{cc}{L}^2& 0\end{array}\right]}^T$$

(69)

$$u_{f,D}{}^2=S^2\left(\xi =1\right)q_f{}^2$$

(70)

$${\mathsf{\Gamma }}_6,{\mathsf{\Gamma }}_7:r_{o_4}+R^4{u}_E{}^4=r_{o_1}+R^1{u}_E{}^1,\mathrm{where}:$$

(71)

$$u_E{}^4=u_{o,E}{}^4+u_{f,E}{}^4$$

(72)

$$u_{o,E}{}^4={\left[\begin{array}{cc}{x}_{p2}& 0\end{array}\right]}^T$$

(73)

$$u_{f,E}{}^4=S^4\left(\xi =\frac{x_{p2}}{x_{p3}}\right)q_f{}^4$$

(74)

$$u_E{}^1=u_{o,E}{}^1+u_{f,E}{}^1$$

(75)

$$u_{o,E}{}^1={\left[\begin{array}{cc}{L}^1& 0\end{array}\right]}^T$$

(76)

$$u_{f,E}{}^1=S^1\left(\xi =1\right)q_f{}^1$$

(77)

$${\mathsf{\Gamma }}_8,{\mathsf{\Gamma }}_9:r_{o_4}+R^4{u}_F{}^4=r_{o_3}+R^3{u}_F{}^3,\mathrm{where}:$$

(78)

$$u_F{}^4=u_{o,F}{}^4+u_{f,F}{}^4$$

(79)

$$u_{o,F}{}^4={\left[\begin{array}{cc}{x}_{p3}& 0\end{array}\right]}^T$$

(80)

$$u_{f,F}{}^4=S^4\left(\xi =1\right)q_f{}^4$$

(81)

$$u_F{}^3=u_{o,F}{}^3+u_{f,F}{}^3$$

(82)

$$u_{o,F}{}^3={\left[\begin{array}{cc}{L}^3& 0\end{array}\right]}^T$$

(83)

$$u_{f,F}{}^3=S^3\left(\xi =1\right)q_f{}^3.$$

(84)

The velocity constraint matrix B is obtained by using Equation (40), that is, by differentiating all the nine position constraints above with respect to all the coordinates including the rigid coordinates and the elastic coordinates:

$$B_r=\left[\begin{array}{ccccccc}\frac{\partial {\mathsf{\Gamma }}_1}{\partial r_{o_1,x}}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial r_{o_1,y}}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial {\alpha }_1}& \cdots & \frac{\partial {\mathsf{\Gamma }}_1}{\partial r_{o_4,x}}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial r_{o_4,y}}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial \theta }\\ \frac{\partial {\mathsf{\Gamma }}_2}{\partial r_{o_1,x}}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial r_{o_1,y}}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial {\alpha }_1}& \cdots & \frac{\partial {\mathsf{\Gamma }}_2}{\partial r_{o_4,x}}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial r_{o_4,y}}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial \theta }\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ \frac{\partial {\mathsf{\Gamma }}_9}{\partial r_{o_1,x}}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial r_{o_1,y}}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial {\alpha }_1}& \cdots & \frac{\partial {\mathsf{\Gamma }}_9}{\partial r_{o_4,x}}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial r_{o_4,y}}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial \theta }\end{array}\right]$$

(85)

$$B_f=\left[\begin{array}{cccccccccc}\frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f4}{}^1}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f5}{}^1}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f6}{}^1}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f4}{}^2}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f5}{}^2}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f6}{}^2}& \cdots & \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f4}{}^4}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f5}{}^4}& \frac{\partial {\mathsf{\Gamma }}_1}{\partial q_{f6}{}^4}\\ \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f4}{}^1}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f5}{}^1}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f6}{}^1}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f4}{}^2}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f5}{}^2}& \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f6}{}^2}& \cdots & \frac{\partial {\mathsf{\Gamma }}_2}{\partial q_{f4}{}^4}& \frac{\partial {\mathsf{\Gamma }}_5}{\partial q_{f5}{}^4}& \frac{\partial {\mathsf{\Gamma }}_5}{\partial q_{f6}{}^4}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f4}{}^1}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f5}{}^1}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f6}{}^1}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f4}{}^2}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f5}{}^2}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f6}{}^2}& \cdots & \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f4}{}^4}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f5}{}^4}& \frac{\partial {\mathsf{\Gamma }}_9}{\partial q_{f6}{}^4}\end{array}\right]$$

(86)

$$B=\left[\begin{array}{cc}{B}_r& B_f\end{array}\right]$$

(87)

Finally, the assembled equations of motion are given by:

$$\left[\begin{array}{cc}{M}_{rr}& M_{rf}\\ M_{fr}& M_{ff}\end{array}\right]\left[\begin{array}{c}{q}_r\\ q_f\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& K_{ff}\end{array}\right]\left[\begin{array}{c}{q}_r\\ q_f\end{array}\right]+\left[\begin{array}{c}{B}_r{}^T\\ B_f{}^T\end{array}\right]\lambda =\left[\begin{array}{c}{\left(Q_v\right)}_r\\ {\left(Q_v\right)}_f\end{array}\right]+\left[\begin{array}{c}{\left(Q_g\right)}_r\\ {\left(Q_g\right)}_f\end{array}\right]+\left[\begin{array}{c}{\left(Q_e\right)}_r\\ {\left(Q_e\right)}_f\end{array}\right]$$

(88)

To solve the system DAEs, the coordinate partitioning method was applied as in the rigid body dynamics presented earlier. In this case, the coordinates are partitioned into the following dependent and independent coordinates which are indicated by subscripts d and i, respectively:

$$q=\left[\begin{array}{c}{q}_{d(9x1)}\\ q_{i(15x1)}\end{array}\right]=\left[\begin{array}{c}{\left[\begin{array}{ccccccccc}{r}_{o_{1,x}}{}^1& r_{o_{1,y}}{}^1& {\alpha }_1& r_{o_{2,x}}{}^2& r_{o_{2,y}}{}^2& {\alpha }_3& r_{o_{3,x}}{}^3& r_{o_{3,y}}{}^3& {\alpha }_3\end{array}\right]}^T\\ {\left[\begin{array}{ccccccc}{r}_{o_{4,x}}{}^4& r_{o_{4,y}}{}^4& \theta & q_{f1}& q_{f2}& q_{f3}& q_{f4}\end{array}\right]}^T\end{array}\right]$$

(89)

Using again the native names of the rigid coordinates and identifying zero coordinates, we can write Equation (89) as:

$$q=\left[\begin{array}{c}{q}_{d(9x1)}\\ q_{i(15x1)}\end{array}\right]=\left[\begin{array}{c}{\left[\begin{array}{ccccccccc}{x}_1& 0& {\alpha }_1& x_2& 0& {\alpha }_3& x_3& 0& {\alpha }_3\end{array}\right]}^T\\ {\left[\begin{array}{ccccccc}x& y& \theta & q_{f1}& q_{f2}& q_{f3}& q_{f4}\end{array}\right]}^T\end{array}\right]$$

(90)

The descriptor form of the DAEs is:

$$\left[\begin{array}{cc}M& B^T\\ B& 0\end{array}\right]\left[\begin{array}{c}\ddot{q}\\ \lambda \end{array}\right]=\left[\begin{array}{c}{Q}_v+Q_g+Q_e-Kq\\ -\dot{B}\dot{q}\end{array}\right]$$

(91)

which can be partitioned into the following:

$$\left[\begin{array}{ccc}{M}_{dd}& M_{di}& B_d{}^T\\ M_{id}& M_{ii}& B_i{}^T\\ B_d& B_i& 0\end{array}\right]\left[\begin{array}{c}{\ddot{q}}_d\\ {\ddot{q}}_i\\ \lambda \end{array}\right]=\left[\begin{array}{c}{\left(Q_v\right)}_d+{\left(Q_g\right)}_d+{\left(Q_e\right)}_d-K_{d}q\\ {\left(Q_v\right)}_i+{\left(Q_g\right)}_i+{\left(Q_e\right)}_i-K_{i}q\\ -\dot{B}\dot{q}\end{array}\right]$$

(92)

Upon manipulating Equation (92), the dynamics can be written as follows:

$$\overline{M}{\ddot{q}}_i=\overline{F}$$

(93)

$$\overline{M}=M_{ii}-M_{id}{B}_d{}^{-1}{B}_i+B_i{}^T{\left(B_d{}^T\right)}^{-1}{M}_{dd}{B}_d{}^{-1}{B}_i-B_i{}^T{\left(B_d{}^T\right)}^{-1}{M}_{di}$$

(94)

$$\overline{F}={\overline{Q}}_v+{\overline{Q}}_{\dot{q}}+{\overline{Q}}_g+{\overline{Q}}_e+{\overline{F}}_s$$

(95)

$${\overline{Q}}_v={\left(Q_v\right)}_i-B_i{}^T{\left(B_d{}^T\right)}^{-1}{\left(Q_v\right)}_d$$

(96)

$${\overline{Q}}_{\dot{q}}=M_{id}{B}_d{}^{-1}\dot{B}\dot{q}-B_i{}^T{\left(B_d{}^T\right)}^{-1}{M}_{dd}{B}_d{}^{-1}\dot{B}\dot{q}$$

(97)

$${\overline{Q}}_g={\left(Q_g\right)}_i-B_i{}^T{\left(B_d{}^T\right)}^{-1}{\left(Q_g\right)}_d$$

(98)

$${\overline{Q}}_e={\left(Q_e\right)}_i-B_i{}^T{\left(B_d{}^T\right)}^{-1}{\left(Q_e\right)}_d$$

(99)

$${\overline{F}}_s=-K_{i}q+B_i{}^T{\left(B_d{}^T\right)}^{-1}{K}_{d}q$$

(100)

The forward dynamics is given by:

$${\ddot{q}}_i={\overline{M}}^{-1}\overline{F}$$

(101)

For the purpose of numerical integration, Equation (101) can be presented in the following state space representation:

$$\left[\begin{array}{c}{\dot{Y}}_1\\ {\dot{Y}}_2\\ ⋮\\ {\dot{Y}}_{15}\end{array}\right]=\left[\begin{array}{c}{Y}_{16}\\ Y_{17}\\ ⋮\\ Y_{30}\end{array}\right]$$

(102)

$${\ddot{q}}_i=\left[\begin{array}{c}{\dot{Y}}_{16}\\ {\dot{Y}}_{17}\\ ⋮\\ {\dot{Y}}_{30}\end{array}\right]=\left[\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{\theta }\\ {\ddot{q}}_f\end{array}\right]$$

(103)

$$\left[\begin{array}{c}{\dot{Y}}_1\\ {\dot{Y}}_2\\ ⋮\\ {\dot{Y}}_{30}\end{array}\right]={\left[\begin{array}{cc}{I}_{15x15}& 0_{15x15}\\ 0_{15x15}& {\overline{M}}_{15x15}\end{array}\right]}^{-1}\left[\begin{array}{c}{Y}_{16}\\ ⋮\\ Y_{30}\\ {\overline{F}}_{15x1}\end{array}\right]$$

(104)

Numerically integrating Equation (104) gives the pose of the moving platform as well as the elastic displacements q_f and their corresponding velocities. The acceleration of the moving platform can be computed by imposing the obtained pose and velocities of the moving platform to Equation (101). The velocities and accelerations of the sliders can be obtained based on the coordinate partitioning method as the following:

$${\dot{q}}_d=B_d{}^{-1}\left(-B_i{\dot{q}}_i\right)$$

(105)

$${\ddot{q}}_d=B_d{}^{-1}\left(\dot{B}\dot{q}-B_i{\ddot{q}}_i\right)$$

(106)

6. Finite Element Formulation

The finite element formulation of the floating frame of reference formulation is different from the common finite element modeling as the former separates the elastic coordinates from the rigid (undeformed) coordinates. In other words, the elastic displacement is presented exclusively by using each body frame as reference. Using the finite elements, each body is discretized into a number of elements. The displacement field w is given by:

$$w=Se,$$

(107)

where S and e are the element shape function and the nodal coordinates, respectively. The element shape function as well as the number and types of the nodal coordinates depends on the element type. The nodal coordinates e can be written as:

$$e=e_o+e_f,$$

(108)

where e_o is the undeformed nodal coordinates whereas e_f is the elastic nodal coordinates.

In the finite element formulation of the floating frame of reference formulation, four coordinate systems are used as shown in Figure 6—(1) the global coordinate system X₁X₂X₃, (2) the body i coordinate system $x_1{}^i{x}_2{}^i{x}_3{}^i$, (3) the element coordinate system $x_1{}^{ij}{x}_2{}^{ij}{x}_3{}^{ij}$ and (4) the intermediate element coordinate system $x_{i1}{}^{ij}{x}_{i2}{}^{ij}{x}_{i3}{}^{ij}$. The global coordinate system is inertially fixed. The body i coordinate system, one in each body i, is not rigidly attached to a point in the body. The element coordinate system, one in each element, has an origin rigidly attached to a point in the element. As a result, it translates and rotates with the element. The intermediate element coordinate system, one in each element, has an origin rigidly attached to the origin of the body i coordinate system Oⁱ. This coordinate system has a fixed orientation with respect to the body i coordinate system. Its orientation is selected such that it is initially parallel to the axes of the element coordinate system.

The displacement field in the element j expressed in the body i coordinate system is:

$$u^{ij}=N^{ij}{q}_n{}^{ij},$$

(109)

where:

$$N^{ij}=C^{ij}{S}^{ij}{\overline{C}}^{ij}{B}_1{}^{ij}$$

(110)

and $q_n{}^{ij}$ is the vector of nodal coordinates in the body i, that is,

$$q_n{}^i={\left[\begin{array}{ccccc}{e}_1{}^i& e_2{}^i& \cdots & e_{n_n-1}{}^i& e_{n_n}{}^i\end{array}\right]}^T$$

(111)

For a 2D beam element, that is, a beam element with only longitudinal and transversal translations and a rotation about its longitudinal axis,

$$C^{ij}=\left[\begin{array}{cc}\cos {\beta }^{ij}& -\sin {\beta }^{ij}\\ \sin {\beta }^{ij}& \cos {\beta }^{ij}\end{array}\right]$$

(112)

$${\overline{C}}^{ij}=\left[\begin{array}{cc}{\overline{C}}_1{}^{ij}& 0\\ 0& {\overline{C}}_1{}^{ij}\end{array}\right]$$

(113)

$${\overline{C}}_1{}^{ij}=\left[\begin{array}{ccc}\cos {\beta }^{ij}& \sin {\beta }^{ij}& 0\\ -\sin {\beta }^{ij}& \cos {\beta }^{ij}& 0\\ 0& 0& 1\end{array}\right]$$

(114)

where β^ij is the fixed angle between the intermediate element coordinate system and the body coordinate system. In other words, the matrix C^ij is the constant rotation matrix between the intermediate element coordinate frame and the body coordinate frame. In a planar case, the matrix C^ij is a 2 × 2 elementary rotation matrix such as that given in Equation (112). In a spatial case, the matrix C^ij is a 3 × 3 rotation matrix.

The matrix $B_1{}^{ij}$ is a Boolean matrix defined for each element with respect to the body i nodal coordinates $q_n{}^i$ to establish a connectivity between elements in that body:

$$q_n{}^{ij}=e^{ij}=B_1{}^{ij}{q}_n{}^i.$$

(115)

Similar to Equation (2), the body i nodal coordinates $q_n{}^i$ can be written as:

$$q_n{}^i=q_o{}^i+q_f{}^i,$$

(116)

where $q_o{}^i$ is the vector of undeformed nodal coordinates in the body i whereas $q_f{}^i$ is the vector of elastic nodal deformations in the body i before any boundary condition is applied. Due to a set of applied boundary conditions, the vector $q_f{}^i$ is reduced to the vector ${\widehat{q}}_f{}^i$. This reduction can be conveniently written as:

$$q_f{}^i=B_2{}^i{\widehat{q}}_f{}^i,$$

(117)

where the matrix $B_2{}^i$ is a Boolean matrix.

6.1. Formulation at the Element Level

The inertia matrix of an element is given by the following symmetric matrix:

$$M^{ij}=\left[\begin{array}{cc}{m}_{rr}{}^{ij}& m_{rf}{}^{ij}\\ sym& m_{ff}{}^{ij}\end{array}\right]=\left[\begin{array}{ccc}{m}_{r_o{r}_o}{}^{ij}& m_{r_o\theta }{}^{ij}& m_{r_{o}f}{}^{ij}\\ & m_{\theta \theta }{}^{ij}& m_{\theta f}{}^{ij}\\ sym& & m_{ff}{}^{ij}\end{array}\right]$$

(118)

where:

$$m_{r_o{r}_o}{}^{ij}=\left[\begin{array}{cc}{m}^{ij}& 0\\ 0& m^{ij}\end{array}\right]\mathrm{for}\mathrm{planar}\mathrm{case}.$$

(119)

$$m_{r_o\theta }{}^{ij}=R_{\theta }{}^i{C}^{ij}{\overline{S}}^{ij}{\overline{C}}^{ij}{e}^{ij}$$

(120)

$$m_{r_{o}f}{}^{ij}=R^i{C}^{ij}{\overline{S}}^{ij}{\overline{C}}^{ij}$$

(121)

$$m_{ff}{}^{ij}={\overline{C}}^{ij}{}^T{S}_{ff}{}^{ij}{\overline{C}}^{ij}$$

(122)

$$m_{\theta \theta }{}^{ij}=e^{ij}{}^T{m}_{ff}{}^{ij}{e}^{ij}$$

(123)

$$m_{\theta f}{}^{ij}=e^{ij}{}^T{\overline{C}}^{ij}{}^T{\tilde{S}}^{ij}{\overline{C}}^{ij}$$

(124)

For a 2D beam element with length l^ij, the matrices $S_{ff}$, ${\overline{S}}^{ij}$ and ${\tilde{S}}^{ij}$ are given by:

$$S_{ff}{}^{ij}=m^{ij}\left[\begin{array}{cccccc}1/3& & & & & sym\\ 0& 13/35& & & & \\ 0& \frac{11}{210}{l}^{ij}& \frac{{(l^{ij})}^2}{105}& & & \\ 1/6& 0& 0& 1/3& & \\ 0& 9/70& \frac{13}{420}{l}^{ij}& 0& 13/35& \\ 0& -\frac{13}{420}{l}^{ij}& -\frac{{(l^{ij})}^2}{140}& 0& -\frac{11}{210}{l}^{ij}& \frac{{(l^{ij})}^2}{105}\end{array}\right]$$

(125)

$${\overline{S}}^{ij}=\frac{m^{ij}}{12}\left[\begin{array}{cccccc}6& 0& 0& 6& 0& 0\\ 0& 6& l^{ij}& 0& 6& -l^{ij}\end{array}\right]$$

(126)

$${\tilde{S}}^{ij}=\frac{m^{ij}}{60}\left[\begin{array}{cccccc}0& 21& 3l^{ij}& 0& 9& -2l^{ij}\\ -21& 0& 0& -9& 0& 0\\ -3l^{ij}& 0& 0& -2l^{ij}& 0& 0\\ 0& 9& 2l^{ij}& 0& 21& -3l^{ij}\\ -9& 0& 0& -21& 0& 0\\ 2l^{ij}& 0& 0& 3l^{ij}& 0& 0\end{array}\right]$$

(127)

The stiffness matrix of an element is:

$$K^{ij}=\left[\begin{array}{ccc}0& 0& 0\\ & 0& 0\\ sym& & K_{ff}{}^{ij}\end{array}\right]$$

(128)

where:

$$K_{ff}{}^{ij}={\overline{C}}^{ij}{\overline{K}}_{ff}{}^{ij}{\overline{C}}^{ij}$$

(129)

$${\overline{K}}_{ff}{}^{ij}=\frac{E^{ij}{I}^{ij}}{l^{ij}}\left[\begin{array}{cccccc}\frac{A}{I^{ij}}& & & & & sym\\ 0& \frac{12}{{(l^{ij})}^2}& & & & \\ 0& 6/l^{ij}& 4& & & \\ -\frac{A}{I^{ij}}& 0& 0& \frac{A}{I^{ij}}& & \\ 0& -\frac{12}{{(l^{ij})}^2}& -6/l^{ij}& 0& \frac{12}{{(l^{ij})}^2}& \\ 0& 6/l^{ij}& 2& 0& -6/l^{ij}& 4\end{array}\right]$$

(130)

6.2. Formulation at the Body Level

The inertia matrix of a body i, before the application of any boundary condition, is the following symmetric matrix:

$$M^i=\left[\begin{array}{cc}{m}_{rr}{}^i& m_{rf}{}^i\\ sym& m_{ff}{}^i\end{array}\right]={\displaystyle \sum _{j=1}^{n_e}{M}^{ij}}$$

(131)

The stiffness matrix of a body i, before the application of any boundary condition, is:

$$K^i=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& K_{ff}{}^i\end{array}\right]$$

(132)

where:

$$K_{ff}{}^i={\displaystyle \sum _{j=1}^{n_e}{K}_{ff}{}^{ij}}.$$

(133)

The process of assembly indicated in Equations (131) and (133) is conducted by overlapping the elements of the matrices corresponding to the nodes of connectivity between elements, as well-known in the finite element method. The notation n_e in Equations (131) and (133) denotes the number of elements in the body i. Upon the application of a set of boundary conditions, the inertia and stiffness matrices are reduced to ${\widehat{M}}^i$ and ${\widehat{K}}^i$. This reduction can be easily conducted by suppressing the rows and columns corresponding to the constrained elastic coordinates. It should be noticed that the assembly of the elements and the suppression of the rows and columns corresponding to the constrained elastic coordinates are conducted at the body level.

The quadratic velocity vector, the gravitational force vector and the external force vector in the body i consist of their rigid and elastic components and are defined as in Equations (30), (33) and (35), respectively. The matrix B for the body i can be obtained by using Equation (40).

6.3. Formulation at the System Level

The inertia and stiffness matrices of the system after the application of the boundary conditions are:

$$\widehat{M}=\left[\begin{array}{cc}{M}_{rr}& {\widehat{M}}_{rf}\\ {\widehat{M}}_{fr}& {\widehat{M}}_{ff}\end{array}\right]=\left[\begin{array}{cccccccc}{m}_{rr}{}^1& 0& \cdots & 0& {\widehat{m}}_{rf}{}^1& 0& \cdots & 0\\ 0& m_{rr}{}^2& \cdots & 0& 0& {\widehat{m}}_{rf}{}^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & m_{rr}{}^{n_b}& 0& 0& \cdots & {\widehat{m}}_{rf}{}^{n_b}\\ & & & & {\widehat{m}}_{ff}{}^1& 0& \cdots & 0\\ & & & & 0& {\widehat{m}}_{ff}{}^1& \cdots & 0\\ & & & & ⋮& ⋮& \ddots & ⋮\\ sym& & & & 0& 0& \cdots & {\widehat{m}}_{ff}{}^{n_b}\end{array}\right]$$

(134)

$$\widehat{K}=\left[\begin{array}{cc}0& 0\\ 0& {\widehat{K}}_{ff}\end{array}\right]=\left[\begin{array}{ccccc}0& 0& 0& \cdots & 0\\ 0& {\widehat{K}}_{ff}{}^1& 0& \cdots & 0\\ 0& 0& {\widehat{K}}_{ff}{}^2& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & {\widehat{K}}_{ff}{}^{n_b}\end{array}\right]$$

(135)

The matrices $\widehat{M}$ and $\widehat{K}$ have an identical dimension. The notation n_b in Equations (134) and (135) denotes the number of bodies in the system.

The EOMs of the system are given by DAEs as in Equations (43) and (44). The expanded expression of the EOMs showing its rigid and elastic components is provided in Equation (88).

7. Flexible Body Dynamics Based on Finite Element Approximation

Using a 2D Euler-Bernoulli beam model as in the Rayleigh-Ritz approximation previously discussed, let’s discretize each of the four bodies in the manipulator into four elements with equal lengths, that is,

$$l^{i1}=l^{i2}=l^{i3}=l^{i4}=\frac{L^i}{4}.$$

(136)

As a result, there are five nodes in each body. All the nodes along with the nodal coordinates in each body i, before the application of any boundary condition, is depicted in Figure 7. The nodes 2–4 are the nodes of connectivity at which two adjacent elements are connected.

As the manipulator at hand has four bodies (links), the nodal coordinates in the system are given by:

$$q_n={\left[\begin{array}{cccc}{q}_n{}^1& q_n{}^2& q_n{}^3& q_n{}^4\end{array}\right]}^T.$$

(137)

The sub-vector $q_n{}^i$ in Equation (142) has a dimension of 5 × 3 = 15 as there are 5 nodes in each body i with 3 DOFs in each node, that is,

$$q_n{}^1=\left[\begin{array}{c}{e}_1{}^1\\ e_2{}^1\\ e_3{}^1\\ ⋮\\ e_{14}{}^1\\ e_{15}{}^1\end{array}\right];\begin{array}{cc}& \end{array}{q}_n{}^2=\left[\begin{array}{c}{e}_1{}^2\\ e_2{}^2\\ e_3{}^2\\ ⋮\\ e_{14}{}^2\\ e_{15}{}^2\end{array}\right];\begin{array}{cc}& \end{array}{q}_n{}^3=\left[\begin{array}{c}{e}_1{}^3\\ e_2{}^3\\ e_3{}^3\\ ⋮\\ e_{14}{}^3\\ e_{15}{}^3\end{array}\right];\begin{array}{cc}& \end{array}{q}_n{}^4=\left[\begin{array}{c}{e}_1{}^4\\ e_2{}^4\\ e_3{}^4\\ ⋮\\ e_{14}{}^4\\ e_{15}{}^4\end{array}\right]$$

(138)

The dynamics of the manipulator at hand is constrained by nine position constraint equations. Three of them are position constraints corresponding to the three rigid sliders as given in Equations (64)–(66). The rest of the position constraints are as follow:

$${\mathsf{\Gamma }}_4,{\mathsf{\Gamma }}_5:r_{o_4}=r_{o_2}+R^2{u}_D{}^2$$

(139)

$${\mathsf{\Gamma }}_6,{\mathsf{\Gamma }}_7:r_{o_4}+R^4{u}_E{}^4=r_{o_1}+R^1{u}_E{}^1$$

(140)

$${\mathsf{\Gamma }}_8,{\mathsf{\Gamma }}_9:r_{o_4}+R^4{u}_F{}^4=r_{o_3}+R^3{u}_F{}^3,$$

(141)

where:

$$u_E{}^1=N^1\left(q_o{}^1+B_2{}^1{q}_f{}^1\right)$$

(142)

$$u_D{}^2=N^2\left(q_o{}^2+B_2{}^2{q}_f{}^2\right)$$

(143)

$$u_F{}^3=N^3\left(q_o{}^3+B_2{}^3{q}_f{}^3\right)$$

(144)

$$u_E{}^4=N_E{}^4\left(q_o{}^4+B_2{}^4{q}_f{}^4\right)$$

(145)

$$u_F{}^4=N_F{}^4\left(q_o{}^4+B_2{}^4{q}_f{}^4\right)$$

(146)

$$N^1=C^1{S}^1\left(\xi =1\right){\overline{C}}^1{B}_1{}^1$$

(147)

$$N^2=C^2{S}^2\left(\xi =1\right){\overline{C}}^2{B}_1{}^2$$

(148)

$$N^3=C^3{S}^3\left(\xi =1\right){\overline{C}}^3{B}_1{}^3$$

(149)

$$N_E{}^4=C^4{S}^4\left(\xi =x_{p2}/x_{p3}\right){\overline{C}}^4{B}_1{}^4$$

(150)

$$N_F{}^4=C^4{S}^4\left(\xi =1\right){\overline{C}}^4{B}_1{}^4$$

(151)

$$q_o{}^1={\left[\begin{array}{ccccccccccccccc}0& 0& 0& l^{11}& 0& 0& l^{12}& 0& 0& l^{13}& 0& 0& l^{14}& 0& 0\end{array}\right]}^T$$

(152)

$$q_o{}^2={\left[\begin{array}{ccccccccccccccc}0& 0& 0& l^{21}& 0& 0& l^{22}& 0& 0& l^{23}& 0& 0& l^{24}& 0& 0\end{array}\right]}^T$$

(153)

$$q_o{}^3={\left[\begin{array}{ccccccccccccccc}0& 0& 0& l^{31}& 0& 0& l^{32}& 0& 0& l^{33}& 0& 0& l^{34}& 0& 0\end{array}\right]}^T$$

(154)

$$q_o{}^4={\left[\begin{array}{ccccccccccccccc}0& 0& 0& l^{41}& 0& 0& l^{42}& 0& 0& l^{43}& 0& 0& l^{44}& 0& 0\end{array}\right]}^T$$

(155)

The fixed-free elastic boundary conditions, identical to that in the Rayleigh-Ritz approximation, are used. These boundary conditions constrain all the elastic DOFs in the first nodes of each body and let the remaining elastic DOFs free. As a consequence, all the matrices and vectors for each body should be reduced by suppressing the rows (for the case of vectors) and rows and columns (for the case of matrices) corresponding to the constrained elastic DOFs. This results in a dimensional reduction of the elastic matrices and vectors for each body from 15 to 12. This suppression is expressed in Equation (117). For the fixed-free boundary conditions, this can be written for each body i as:

$$q_f{}^i=B_2{}^i{\widehat{q}}_f{}^i=\left[\begin{array}{cccc}{0}_{3x3}& 0_{3x3}& 0_{3x3}& 0_{3x3}\\ I_{3x3}& 0_{3x3}& 0_{3x3}& 0_{3x3}\\ 0_{3x3}& I_{3x3}& 0_{3x3}& 0_{3x3}\\ 0_{3x3}& 0_{3x3}& I_{3x3}& 0_{3x3}\\ 0_{3x3}& 0_{3x3}& 0_{3x3}& I_{3x3}\end{array}\right]\left[\begin{array}{c}{\left(e_f{}^i\right)}_4\\ {\left(e_f{}^i\right)}_5\\ {\left(e_f{}^i\right)}_6\\ {\left(e_f{}^i\right)}_7\\ ⋮\\ {\left(e_f{}^i\right)}_{15}\end{array}\right]$$

(156)

Accordingly, the dimension of the system is also reduced. The EOMs of the system are written as Equation (88) with the generalized coordinates q defined in Equations (89) and (90) but now the elastic coordinates for each body $q_f{}^i$ are given by Equation (156).

8. Numerical Results

8.1. Manipulator Parameters

Given the manipulator at hand with legs having length of L₁ = 486 mm, L₂ = 500 mm and L₃ = 507 mm. The moving platform has dimensions of x_p2 = x_p3 = x_p = 315 mm. Hence, two upper revolute joints are coincident. The COM of the legs are L_1G = 199 mm, L_2G = 244.8685 mm and L_3G = 238.9118 mm whereas the COM of the moving platform is x_PG = 165.9629 mm. The masses of all of the sliders, all of the three legs and the moving platform are m_s1 = m_s2 = m_s3 = 7 kg, m_L1 = 8.5 kg, m_L2 = 9.8470 kg, m_L3 = 10.8420 kg and m_p = 15.2425 kg, respectively. The mass moments of inertia of the legs with respect to their COM are given by I_G1 = 0.1673 kg m², I_G2 = 0.2051 kg m² and I_G3 = 0.2322 kg m², whereas the mass moment of inertia of the moving platform with respect to its COM is given by I_Gp = 0.1260 kg m². All the links have an identical modulus of elasticity E = 7.0 × 10¹⁰ Pa and an identical rectangular cross section of b × h = 150 × 50 mm. Assuming that the manipulator is on XY plane (and accordingly the gravity applies), the depth b of 150 mm is measured in z-direction (perpendicular to the paper) whereas the depth h of 50 mm is measured on XY plane (on the paper plane). It is assumed that there is no friction at all the joints.

8.2. Input Forces and Elastic Coordinate Boundary Condition

Two different types of loads namely sinusoidal and triangular pulse loads are applied to the active joints of the manipulator. The sinusoidal loads are given by:

$$Q_{a1}=10\sin 20t,\begin{array}{c}\end{array}{Q}_{a2}=-2\sin 20t,\begin{array}{c}\end{array}{Q}_{a3}=-2\sin 20t,$$

(157)

whereas the triangular pulse loads ramp up from zero at t = 0 to respectively 10 N, −2 N and −2 N for the three sliders at t = 0.1 s and ramp down to zero at t = 0.2 s. It should be noticed that the manipulator will keep moving after the triangular pulse loads ramp to zero due to the gravity. The elastic displacement is evaluated by applying fixed-free boundary conditions to each link.

8.3. Solutions

Since the DAEs representing the dynamics of the manipulator are transformed to ordinary differential equations (ODEs) through partitioning, the solution can be easily obtained by integrating the ODEs. The second order ODEs were transformed to a state space representation as described earlier and subsequently integrated numerically by using MATLAB ODE45. The integration was performed within 0 to 4 s. Figure 8 and Figure 9 show the end-effector position and orientation of the manipulator subject to the sinusoidal loads. Both figures depict the end-effector position and orientation of the manipulator modeled as rigid body and flexible body. The flexible body model used in Figure 8 is based on the Rayleigh-Ritz method whereas that used in Figure 9 is based on the finite element method.

For more clarity, Figure 10 shows the elastic displacements as well as the orientation of the end effector evaluated by using the two methods in the case of sinusoidal loads. Figure 11 shows the elastic displacements as well as the orientation of the end effector evaluated by using the two methods in the case of the triangular pulse loads. It should be noticed that the orientation plots shown in Figure 10 and Figure 11 indicate the absolute orientation values of the manipulator with flexibility of all the links taken into account, not the elastic displacements alone. It is shown that the flexibility of the manipulator links results in infinitesimal displacement in both the translational DOFs and the rotational DOF in addition to the rigid body displacement of the end-effector. This infinitesimal displacement should be taken into account when dealing with a high precision application such as machining. Furthermore, it is shown that the elastic displacement obtained by both the methods have some degree of agreement, with some difference in the magnitude.

It can also be seen that the elastic displacements in the y-direction is larger than that in the x-direction. The former is roughly twice the latter. This can be related to the manipulator stiffness. The larger elastic displacement in the y-direction might be caused by smaller manipulator stiffness in the y-direction. This can be due to less depth of the manipulator in the y-direction compared to its depth in the x-direction during the performed motion. In fact, the stiffness of the manipulator in the y-direction is at its best when the manipulator is standing tall and therefore has a large depth in the y-direction. In such a posture, the stiffness in the x-direction will be reduced. The varying stiffness corresponding with the changing posture is a typical characteristics of parallel manipulators.

Moreover, the magnitude of the elastic displacements which represents the flexibility of the manipulator is shown to be unrelated with the methods as it is shown that the finite element method gives fewer elastic displacements in the case of sinusoidal loading but it gives larger elastic translational displacements in the case of triangular pulse loading. This implies an inconsistent relation between the flexibility of the manipulator and both the methods given a variety of loadings.

8.4. Sensitivity Analysis

The sensitivity of the displacements x, y and θ with respect to the design parameters was analyzed by evaluating the gradients of those displacements with respect to the design parameters. As all the links of the manipulator are designed to have the same material, the inertial and elasticity properties of the links are dependent on the mass density $\rho$ and the modulus of elasticity E of the links, respectively. Other design parameters affecting both the mass and the elasticity of the links are the cross section geometry. As the links in this case have a rectangular cross section, these parameters are the width b and height h of the rectangular cross section. Thus, the design parameters are P = {E, $\rho$, b, h}. The gradients of the displacements with respect to the aforementioned design parameters were obtained by using finite difference approach. The difference between the displacements corresponding to perturbed design parameters and the displacements corresponding to unperturbed design parameters is divided by the amount of perturbation. Using SI units for all the design parameters, a perturbation of 1 × 10⁻⁵ was used for all the design parameters.

Figure 12, Figure 13 and Figure 14 show respectively the sensitivity of the displacements x, y and θ with respect to the design parameters. All the figures show that the displacements in all the DOFs are more sensitive to the dimension of the cross section than to both the modulus of elasticity and the mass density of the links material. Furthermore, the displacements in all the DOFs are more sensitive to the height h of the cross section than to the width b of the cross section. This is expected as the manipulator is working on the XY plane which is the plane where the height h is defined. In this case, the area moments of inertia of the links, which affect the elasticity of the links, have a cubic relation with the height h and have a linear relation with the width b of the cross section. Hence, a change in the height h will affect more the elasticity of the links than a change of the width b. Consequently, the rigidity of the links can be significantly increased by increasing the height h of the links’ cross section.

9. Conclusions

A more realistic model of a flexible parallel manipulator was presented in this paper by assuming that all the links, including the moving platform, are flexible. Since the moving platform connects the end-effector to the legs, in general its flexibility should be taken into account in order to obtain a more accurate end-effector position. Furthermore, the non-symmetry of the manipulator represents a more general topology of planar parallel manipulators. The simulation using both the Rayleigh-Ritz and finite element approximations, as expected, showed that the manipulator undergoes a larger displacement when the link flexibility is taken into account. The comparison between both approximations using Euler-Bernoulli beam to model the links, which is the main contribution of this paper, showed an agreement in the direction and rough magnitude of the elastic displacement. The agreement between the two methods verifies each other. This can serve as a base to use only the finite element method using FFRF framework in the future to evaluate the elastic displacements of the manipulator at hand with complex geometry of the links. In such a case, appropriate 3D finite elements should be used to adequately represent the complex geometry. The huge number of DOFs resulting from the use of 3D finite elements can be reduced, if required, by employing a model order reduction technique. Finally, based on the sensitivity analysis, the rigidity of the links of the manipulator at hand can be significantly increased by increasing the height h of the links’ cross section.