© Y. Sheng and X. Chen Published by EDP Sciences 2024

1 Introduction

Spatial parallel mechanism [1] has some competitive advantages compared with other mechanisms, such as more compact mechanical structure, smaller kinematic inertia, higher motion accuracy and good isotropy [2], it has been successfully applied in the field of coordinate measuring mechanism. The use of spatial parallel mechanism improves the measurement accuracy and mechanism stability of coordinate measuring mechanism [3–5]. Spherical joints are extensively used linking devices with the advantages of higher transmission efficiency, but existence of clearance between joint elements is inevitable because of assembly requirements. Spherical clearance joints exert a great influence on the dynamic behavior of spatial parallel coordinate measuring mechanism, reducing measurement accuracy and causing unpredictable vibration [6]. Studying the influence of clearance joints on spatial parallel coordinate measuring mechanism can be demonstrated by the dynamic behavior under different research parameters [7]. Based on this, it is critical to establish a dynamics modeling and analysis method applied for spatial parallel coordinate measuring mechanism, and the method can also enable the research of the effect of spherical clearance joints on dynamic behavior of spatial parallel mechanisms.

In the past, research about dynamic of spatial parallel mechanism is imperfect, and research direction is mainly focused on the plane multi-link mechanism [8–11]. Li et al. [12] developed a comprehensive methodology for predicting wear characteristic in crank slider mechanism with multiple clearance joints, and found a result that changing boundary conditions have large effect on wear characteristics of the mechanism. Bai et al [13] researched the influence of joints with axial and radial clearances on the dynamic response of planar mechanical systems, and discussed dynamic response of hybrid clearance joints to crank-slider mechanism with different examples. Chen et al. [14] compared the experimental data with numerical simulations to investigate the effect of clearance joints on the dynamical characteristics of planar multibody systems, and verified the correctness of modeling and analysis methods. Wu et al. [15] analyzed the system dynamics of planar crank-slider mechanism with multiple clearance joints by studying the correlation dimension and bifurcation action, and found that the response of system has different sensitivities to different clearance joints. Wang et al. [16] explored dynamic behavior of planar four-bar mechanism considering different clearance joint mounting positions and number of clearance joints, and inferred that the larger the number of clearance joints, the worse the stability of mechanism.

The research on chaotic characteristics of spatial parallel mechanisms is immature, and the current research mainly focuses on the dynamic response of spatial parallel mechanisms [17–20]. Erkaya et al. [21] explored influence of vibration generated by clearance joint on chaos of mechanism, concluded that vibration reduces dynamic performance and making rod flexible can reduce the negative impact which clearance bring to mechanism. Hou et al. [22] explored dynamic response of four-bar mechanism with rigid rod, flexible rod and flexible socket, and inferred that flexible rod and flexible socket can weaken adverse impact on the system produced by clearance joint. Wang et al. [23] built dynamic model of parallel mechanism with spherical clearance and flexible input rod, found a result that making the rod flexible is very necessary to study the influence of clearance joint on mechanism. Zhang et al. [24] analyzed the dynamics of a new parallel mechanism using implicit integration algorithm and Baumgart stabilization technology, and verified the validity of numerical results of the dynamic response using ADAMS simulation software. Chen et al. [25] used the principle of virtual work to analyze dynamics of a new type of over constrained parallel mechanism, and verified the validity of dynamic simulation results through the physical model simulation in Simscape and rigid flexible coupling model simulation in Adams.

In this paper, because it is immature and imperfect in the field of researching dynamic response and chaotic characteristics of spatial parallel mechanisms at the same time, a universal dynamics modeling and analysis method applied for spatial parallel coordinate with multiple spherical clearance joints was proposed. In Section 2, motion collision process of spherical clearance joints is described in detail, kinematic model of spherical clearance joint is established, and different force models are introduced to analyze contact force at the contact collision point. In Section 3, kinematic model of spatial parallel coordinate measuring mechanism is established, and accompanying motion of mechanism is explained, and clearance dynamic model based on Lagrange multiplier technology is established. In Section 4, the correctness of clearance dynamics modeling method is verified by comparison with Adams virtual simulation, and taking clearance value and friction coefficient as single input variables, effect of spherical clearance joints on response and chaotic characteristics of mechanism is studied.

2 Motion description of spherical clearance joint

2.1 Kinematic model

Spatial parallel coordinate measuring mechanism (SPCMM) is made up of three identical kinematic chains connected to a mobile platform, SPCMM can realize data measurement of three degrees of freedom (one translation measurement and two rotation measurement). Every branch consists of a driving rod and a driven rod, driving rod is fixed connection with the frame through the translational joint, the driving rod is fixed connection with the driven rod through spherical joint, and driven rod is fixed connection with the terminal platform through rotating joint.

Figure 1 is the three-dimensional virtual model of SPCMM. Figure 2 is structural sketch of SPCMM, P_i represents translational joint, S_i represents spherical joint, R_i represents rotating joint, O_A-XYZ is world coordinate system fixed connection with the frame which is expressed as {A}, its vertical distance from point P is a₁, o_i is object coordinate system fixed connection with other moving components which is expressed as {o_i}, and {o_i}(i = 1 ~ 6) is fixed at the centroid position of each component, and o₇ is established at the center position of R₁ and R₂ on the terminal platform. Length h_i(i = 1 ~ 3) is the initial installation position of the translational joint from point P.

Spherical joints S1, S₂  and S₃ are clearance joints, with the ball socket of the ball joint fixedly connected to the driving rod, and the ball joint ball fixedly connected to the driven rod. Existing clearance makes the position constraint of spherical joint invalid, and the ball can get a freely motion in six degrees of freedom, therefore, the modeling of spherical clearance joint can be seen in Figure 3. Coordinate system {A} is world coordinate system, Point Z_i and Z_j are center points of the socket and ball, R and r is radius of socket and ball. Vector ^An_n and ^At_t are normal vector and tangential vector of plane where the contact point is located. Vector e_ij is eccentric vector of center points in spherical joint. Vector F_N and F_T are force at the point of impact.

The position coordinate about origin of object coordinate system fixed connection with driving rod in world coordinate system is

$${}^{\mathit{A}}\mathit{P}_{\mathit{o}\mathit{i}}={[{}^{\mathit{A}}\mathit{P}_{\mathit{o}\mathit{i}\mathit{x}}{}^{\mathit{A}}\mathit{P}_{\mathit{o}\mathit{i}\mathit{y}}{}^{\mathit{A}}\mathit{P}_{\mathit{o}\mathit{i}\mathit{z}}]}^{\text{T}}(\mathit{i}=\mathrm{1,2,3})$$(1)

The Euler matrix of object coordinate system in world coordinate system is written as

$$\begin{array}{l}{}^{\mathit{A}}\mathit{R}_{\mathit{i}}=\mathit{R}(\mathit{x},\alpha )\cdot \mathit{R}(\mathit{z},\beta )\cdot \mathit{R}(\mathit{y},\gamma )\\ =\left[\begin{array}{l}\mathit{c}\beta \cdot \mathit{c}\gamma -\mathit{s}\beta \mathit{c}\beta \cdot \mathit{s}\gamma \\ \mathit{s}\alpha \cdot \mathit{s}\gamma +\mathit{c}\alpha \cdot \mathit{c}\gamma \cdot \mathit{s}\beta \mathit{c}\alpha \cdot \mathit{c}\beta \mathit{c}\alpha \cdot \mathit{s}\beta \cdot \mathit{s}\gamma -\mathit{c}\gamma \cdot \mathit{s}\alpha \\ \mathit{s}\alpha \cdot \mathit{s}\gamma +\mathit{c}\alpha \cdot \mathit{c}\gamma \cdot \mathit{s}\beta \mathit{c}\alpha \cdot \mathit{c}\beta \mathit{c}\alpha \cdot \mathit{s}\beta \cdot \mathit{s}\gamma -\mathit{c}\gamma \cdot \mathit{s}\alpha \end{array}\right]=[\mathit{u}\mathit{v}\mathit{w}]=\left[\begin{array}{l}{\mathit{u}}_{\mathit{x}}{\mathit{v}}_{\mathit{x}}{\mathit{w}}_{\mathit{x}}\\ {\mathit{u}}_{\mathit{y}}{\mathit{v}}_{\mathit{y}}{\mathit{w}}_{\mathit{y}}\\ {\mathit{u}}_{\mathit{z}}{\mathit{v}}_{\mathit{z}}{\mathit{w}}_{\mathit{z}}\end{array}\right].\end{array}$$(2)

The position coordinate about center point of ball and socket in world coordinate system is written as follow

$$\{\begin{array}{l}{}^{\mathit{A}}\mathit{P}_{\mathit{Z}\mathit{i}}={}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {}^{\mathit{o}}\mathit{p}_{\mathit{Z}\mathit{i}}+{}^{\mathit{A}}\mathit{P}_{\mathit{o}\mathit{i}}\\ {}^{\mathit{A}}\mathit{P}_{\mathit{Z}\mathit{j}}={}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {}^{\mathit{o}}\mathit{p}_{\mathit{Z}\mathit{j}}+{}^{\mathit{A}}\mathit{P}_{\mathit{o}\mathit{i}}\end{array}.$$(3)

Vector ^op_zi and ^op_Zj are position coordinate of center point of ball and socket in object coordinate system and they can be expressed as

$$\{\begin{array}{l}{}^{\mathit{o}}\mathit{p}_{\mathit{Z}\mathit{i}}={}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\left[}^{\frac{{\mathit{L}}_1}{2}}\\ {}^{\mathit{o}}\mathit{p}_{\mathit{Z}\mathit{j}}={}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {[}^{-}\end{array}.$$(4)

Symbol L₁ and L₂ are the length of driving rod and driven rod.

The eccentric vector ^Ae_ij of spherical joint is written as

$${\hspace{0.17em}}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}{=}^{\mathit{A}}{\mathit{P}}_{\mathit{Z}\mathit{j}}{-}^{\mathit{A}}{\mathit{P}}_{\mathit{Z}\mathit{i}}$$(5)

The amplitude of eccentric vector ^Ae_ij at clearance joint is expressed as

$${\mathit{e}}_{\mathit{i}\mathit{j}}=\sqrt{{\hspace{0.17em}}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}^{\text{T}\bullet \mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}}=\sqrt{{\hspace{0.17em}}^{\mathit{A}}\mathit{e}\mathit{i}\mathit{j}{\left(1,1\right)}^2{+}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}{\left(2,1\right)}^2{+}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}{\left(3,1\right)}^2}.$$(6)

The relative impact velocity of center points between ball and socket is

$${\dot{\mathit{e}}}_{\mathit{i}\mathit{j}}=\frac{{\hspace{0.17em}}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}\left(1,1\right){\bullet }^{\mathit{A}}{\dot{\mathit{e}}}_{\mathit{i}\mathit{j}}\left(1,1\right){+}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}\left(2,1\right){\bullet }^{\mathit{A}}{\dot{\mathit{e}}}_{\mathit{i}\mathit{j}}\left(2,1\right){+}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}\left(3,1\right){\bullet }^{\mathit{A}}{\dot{\mathit{e}}}_{\mathit{i}\mathit{j}}\left(3,1\right)}{{\mathit{e}}_{\mathit{i}\mathit{j}}}.$$(7)

The normal vector of the plane where the impact point C_i and D_j are located is

$${\hspace{0.17em}}^{\mathit{A}}{\mathit{n}}_{\mathit{n}}=\frac{{\hspace{0.17em}}^{\mathit{A}}{\mathit{e}}_{\mathit{i}\mathit{j}}}{{\mathit{e}}_{\mathit{i}\mathit{j}}}$$(8)

The clearance value of the spherical joint is

$${\mathit{c}}_{\mathit{i}\mathit{j}}=\mathit{R}-\mathit{r}.$$(9)

Pseudo-penetration depth between spherical joint elements is expressed as

$${\delta }_{\mathit{i}\mathit{j}}={\mathit{e}}_{\mathit{i}\mathit{j}}-{\mathit{c}}_{\mathit{i}\mathit{j}}.$$(10)

According to the value of pseudo-penetration depth, the specific motion state of spherical joint elements can be judged as follow

$$\{\begin{array}{l}{\delta }_{\mathit{i}\mathit{j}}>0,\text{in occurred}\\ {\delta }_{\mathit{i}\mathit{j}}=0,\text{Just contact or separation}\\ {\delta }_{\mathit{i}\mathit{j}}<0,\text{no contact and free motion}\end{array}$$(11)

when spherical joint elements are in contact and collision state, the position coordinate of potential collision points C_i and D_j in world coordinate system are

$$\{\begin{array}{l}{}^{\mathit{A}}\mathit{P}_{\mathit{C}\mathit{i}}={}^{\mathit{A}}\mathit{P}_{\mathit{Z}\mathit{i}}+\mathit{R}\cdot {}^{\mathit{A}}\mathit{n}_{\mathit{n}}\\ {}^{\mathit{A}}\mathit{P}_{\mathit{D}\mathit{j}}={}^{\mathit{A}}\mathit{P}_{\mathit{Z}\mathit{j}}+\mathit{r}\cdot {}^{\mathit{A}}\mathit{n}_{\mathit{n}}\end{array}.$$(12)

Velocity vector of collision point in world coordinate system is

$$\{\begin{array}{l}{}^{\mathit{A}}\mathit{V}_{\mathit{C}\mathit{i}}={}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{i}}+\mathit{R}\cdot {}^{\mathit{A}}{\displaystyle \stackrel{.}{\mathit{n}}}_{\mathit{n}}\\ {}^{\mathit{A}}\mathit{V}_{\mathit{D}\mathit{j}}={}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{j}}+\mathit{r}\cdot {}^{\mathit{A}}{\displaystyle \stackrel{.}{\mathit{n}}}_{\mathit{n}}\end{array},$$(13)where, ^AV_Zi and ^AV_Zj are velocity vector of center point in world coordinate system, they can be written as

$$\{\begin{array}{l}{}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{i}}={}^{\mathit{A}}\mathit{W}_{\mathit{o}\mathit{i}}\times {}^{\mathit{A}}\mathit{r}_{\mathit{Z}\mathit{i}}+{}^{\mathit{A}}\mathit{V}_{\mathit{o}\mathit{i}}\\ {}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{j}}={}^{\mathit{A}}\mathit{W}_{\mathit{o}\mathit{i}}\times {}^{\mathit{A}}\mathit{r}_{\mathit{Z}\mathit{j}}+{}^{\mathit{A}}\mathit{V}_{\mathit{o}\mathit{i}}\end{array},$$(14)

where, ^AW_oi is the Euler velocity vector, ^Ar_Zi is position vector of center point Z_i in object coordinate system which transformed from object coordinate system to world coordinate system, ^AV_oi is velocity vector of the origin of object coordinate system in world coordinate system.

Normal component of relative collision velocity is

$${}^{\mathit{A}}\mathit{V}_{\mathit{n}}=(({}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{j}}-{}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{i}}{)}^{\text{T}}\cdot {}^{\mathit{A}}\mathit{n}_{\mathit{n}})\cdot {}^{\mathit{A}}\mathit{n}_{\mathit{n}},$$(15)

Tangential component of relative collision velocity is

$${}^{\mathit{A}}\mathit{V}_{\mathit{t}}=({}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{j}}-{}^{\mathit{A}}\mathit{V}_{\mathit{Z}\mathit{i}}{)}^{\text{T}}-{}^{\mathit{A}}\mathit{V}_{\mathit{n}}.$$(16)

Amplitude of tangential component of relative collision velocity can be expressed as

$${\mathit{V}}_{\mathit{t}}=\sqrt{{}^{\mathit{A}}\mathit{V}_{\mathit{t}}^{\text{T}}\times {}^{\mathit{A}}\mathit{V}_{\mathit{t}}}$$(17)

So tangential vector of plane where impact point C_i and D_j are located is

$${}^{\mathit{A}}\mathit{t}_{\mathit{t}}=\frac{{}^{\mathit{A}}\mathit{V}_{\mathit{t}}}{{\mathit{V}}_{\mathit{t}}}.$$(18)

Fig. 1 Characteristic map of SPCMM.

Fig. 2 The kinematic scheme of SPCMM.

Fig. 3 Spherical clearance joint model.

2.2 Normal contact force and tangential friction force model

Lankarani-Nikravesh (L-N) contact force model [26] has been confirmed by many experiments, and this model has been widely used [27–29]. L-N contact force model performs well in dealing with contact problems with low energy dissipation and low collision velocities at high coefficient of restitution [30, 31]. Therefore, L-N contact force model was chosen in this paper to deal with normal contact force at collision point in spherical clearance joins of spatial parallel mechanism. Normal contact force model based on L-N contact force model can be written as

$${\mathit{F}}_{\mathit{N}}={\mathit{K}}_{\mathit{i}\mathit{j}}{\delta }_{\mathit{i}\mathit{j}}^{\mathit{n}}\left[1+\frac{3(1-{{\mathit{c}}_{\mathit{e}}}^2){{\displaystyle \dot{\delta }}}_{\mathit{i}\mathit{j}}}{4{{\displaystyle \dot{\delta }}}_{\mathit{i}\mathit{j}0}}\right],$$(19)

Where, K_ijK_ij is the stiffness coefficient, n is the power exponent of material and δ_ij is pseudo-penetration depth, ${{\displaystyle \dot{\delta }}}_{\mathit{i}\mathit{j}}$ is pseudo-penetration velocity, it comes from the derivative of δ_ij with respect to time, ${{\displaystyle \dot{\delta }}}_{\mathit{i}\mathit{j}0}$ is initial value of pseudo-penetration velocity, c_e is coefficient of restitution, υ_i and υ_j are Poisson's ratio of socket and ball, E_i and E_j are elastic modulus, R is radius of socket and r is radius of ball.

In equation (19),${\mathit{K}}_{\mathit{i}\mathit{j}}=\frac{4}{3({\delta }_{\mathit{i}}+{\delta }_{\mathit{j}})}\sqrt{\frac{\mathit{R}\mathit{r}}{\mathit{R}-\mathit{r}}}$ and in extended formula of K_ij, δ_i and δ_j are fixed values related to material properties, and${\delta }_{\mathit{i}}=\frac{1-{\upsilon }_{\mathit{i}}^2}{{\mathit{E}}_{\mathit{i}}}$, ${\delta }_{\mathit{j}}=\frac{1-{\upsilon }_{\mathit{j}}^2}{{\mathit{E}}_{\mathit{j}}}$.

In order to describe tangential force at collision point, the tangential friction model need to be established. Among the tangential force models for multi-body system dynamics, the Coulomb friction model is more common [32]. However, when the velocity approaches zero, this model does not describe the tangential friction at collision point well. Ambrosio [33] proposed an improved Coulomb friction model by introducing dynamic correction coefficients, improved Coulomb friction model compensates for the reduced stability of numerical integration due to the direction changed of friction force when collision velocity is zero and also the improved Coulomb friction model can well describe the dry friction and the viscous phenomena in spherical clearance joints at relatively low speeds. The improved Coulomb friction model can be expressed as

$${\mathit{F}}_{\mathit{T}}=-{\mathit{c}}_{\mathit{f}}{\mathit{c}}_{\mathit{g}}{\mathit{F}}_{\mathit{N}}\frac{{}^{\mathit{A}}\mathit{t}_{\mathit{t}}}{\left|{}^{\mathit{A}}\mathit{t}_{\mathit{t}}\right|},$$(20)

Where c_f is sliding friction coefficient, c_g is dynamic correction coefficient, ^At_t is scalar of tangential velocity.

$${\mathit{c}}_{\mathit{g}}=\{\begin{array}{l}0\mathit{i}\mathit{f}\left|{}^{\mathit{A}}\mathit{t}_{\mathit{t}}\right|<{\mathit{v}}_{\alpha }\\ \frac{\left|{}^{\mathit{A}}\mathit{t}_{\mathit{t}}\right|-{\mathit{v}}_{\alpha }}{{\mathit{v}}_{\theta }-{\mathit{v}}_{\alpha }}\mathit{i}\mathit{f}{\mathit{v}}_{\alpha }\le \left|{}^{\mathit{A}}\mathit{t}_{\mathit{t}}\right|\le {\mathit{v}}_{\theta }\\ 1\mathit{i}\mathit{f}\left|{}^{\mathit{A}}\mathit{t}_{\mathit{t}}\right|>{\mathit{v}}_{\theta }\end{array}.$$(21)

Symbols υ_a and υ_θ are threshold of given limiting velocity.

According to equations (19) and (20), the contact force ball exerts on socket can be expressed as

$${\mathit{F}}_{\mathit{i}}={\mathit{F}}_{\mathit{N}}\cdot {}^{\mathit{A}}\mathit{n}_{\mathit{n}}+{\mathit{F}}_{\mathit{T}}\cdot {}^{\mathit{A}}\mathit{t}_{\mathit{t}}.$$(22)

Reaction force of socket to ball is

$${\mathit{F}}_{\mathit{j}}=-{\mathit{F}}_{\mathit{i}}.$$(23)

Under the action of contact force, the moment at the centroid of components i and j can be written as

$$\{\begin{array}{l}{\mathit{M}}_{\mathit{i}}={(}^{\mathit{A}}{\mathit{P}}_{\mathit{C}\mathit{i}}{-}^{\mathit{A}}{\mathit{P}}_{\mathit{O}\mathit{i}})\times {\mathit{F}}_{\mathit{i}}\\ {\mathit{M}}_{\mathit{j}}={(}^{\mathit{A}}{\mathit{P}}_{\mathit{D}\mathit{j}}{-}^{\mathit{A}}{\mathit{P}}_{\mathit{O}\mathit{j}})\times {\mathit{F}}_{\mathit{j}}\end{array}.$$(24)

3 Dynamic model of SPCMM with spherical clearance joints

3.1 Kinematic model

Kinematic modeling of SPCMM is the basis for solving clearance dynamics, therefore, it is critical to build an appropriate kinematic model.

Position coordinate of origin of object coordinate system {0₇} in world coordinate system is O_Ao₇ = [x y z]^T.

Position coordinates of the spherical joints S_i in world coordinate system are

$${}^{\mathit{A}}\mathit{P}_{\mathit{S}\mathit{i}}={[{}^{\mathit{A}}\mathit{P}_{\mathit{S}\mathit{i}\mathit{x}}{}^{\mathit{A}}\mathit{P}_{\mathit{S}\mathit{i}\mathit{y}}{}^{\mathit{A}}\mathit{P}_{\mathit{S}\mathit{i}\mathit{z}}]}^{\text{T}}.$$(25)

Position coordinates of rotating joints R_i in world coordinate system are

$${}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}}={[{}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}\mathit{x}}{}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}\mathit{y}}{}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}\mathit{z}}]}^{\text{T}}.$$(26)

Vectors of driven rods in world coordinate system are

$${\mathit{S}}_{\mathit{i}}{\mathit{R}}_{\mathit{i}}={}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}}-{}^{\mathit{A}}\mathit{P}_{\mathit{S}\mathit{i}}.$$(27)

Position coordinates of the rotating joints R_i in object coordinate system{o₇} are

$${}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}}={[{}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}\mathit{x}}{}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}\mathit{y}}{}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}\mathit{z}}]}^{\text{T}}..$$(28)

The transformational vector of the above vector changed from object coordinate system {o₇} to world coordinate system {A} is

$${}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}}={[{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}\mathit{x}}{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}\mathit{y}}{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}\mathit{z}}]}^{\text{T}}.$$(29)

Vector ε_i is axis vector of rotation joints, it can be written as

$${\epsilon }_1={\left[0\hspace{0.17em}1\hspace{0.17em}0\right]}^{\text{T}},{\epsilon }_2={\left[0\hspace{0.17em}1\hspace{0.17em}0\right]}^{\text{T}},{\epsilon }_3={\left[0\hspace{0.17em}0\hspace{0.17em}1\right]}^{\text{T}}.$$(30)

According to the constraint relationship at the rotating joints R_i in the parallel mechanism, S_iR_i is always perpendicular to ε_i, constraint equation can be written as

$${\mathit{S}}_{\mathit{i}}{\mathit{R}}_{\mathit{i}}\cdot {ϵ}_{\mathit{i}}=0.$$(31)

By solving equation (31), the accompanying motion of the mechanism can be obtained. The Euler matrix of object coordinate system under world coordinate system can be written as

$${\mathit{R}}_{\mathit{i}}=\mathit{T}(\mathit{z},{\mathit{h}}_{\mathit{i}})\cdot \mathit{T}(\mathit{x}ʹ,-{\mathit{L}}_1)(\mathit{i}=\mathrm{1,2,3}).$$(32)

The Euler transformation process from point O_A to point R_i can be written as

$${\mathit{R}}_{\mathit{i}}=\mathit{R}(\mathit{x},{\alpha }_{\mathit{i}})\mathit{T}(\mathit{z}ʹ,{\mathit{h}}_{\mathit{i}})\cdot \mathit{R}(\mathit{z}ʹʹ,\beta )\cdot \mathit{R}(\mathit{y}ʹʹʹ,\gamma )\mathit{T}(\mathit{x}ʹʹʹʹ,{\mathit{L}}_2)(\mathit{i}=\mathrm{4,5,6}).$$(33)

Simplified equation (33), and simplified results can be written as

$$\left[\begin{array}{c}{}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}\mathit{x}}\\ {}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}\mathit{x}}\\ {}^{\mathit{A}}\mathit{P}_{\mathit{R}\mathit{i}\mathit{x}}\end{array}\right]=\left[\begin{array}{c}{\mathit{L}}_2\cdot \cos {\beta }_{\mathit{i}}\cdot \cos {\gamma }_{\mathit{i}}\\ {\mathit{L}}_2\cdot (\sin {\alpha }_{\mathit{i}}\cdot \sin {\gamma }_{\mathit{i}}+\cos {\alpha }_{\mathit{i}}\cdot \cos {\gamma }_{\mathit{i}}\cdot \sin {\beta }_{\mathit{i}})\\ {\mathit{h}}_1-{\mathit{L}}_2\cdot (\cos {\alpha }_{\mathit{i}}\cdot \sin {\gamma }_{\mathit{i}}-\cos {\gamma }_{\mathit{i}}\cdot \sin {\alpha }_{\mathit{i}}\cdot \sin {\beta }_{\mathit{i}})\end{array}\right].$$(34)

Solved equation (34), the Euler angle γ_i and β_i can be calculated as

$${\gamma }_{\mathit{i}}=\mathit{a}\sin (({\mathit{F}}_{\mathit{i}}\cdot \sin {\alpha }_{\mathit{i}}-{\mathit{G}}_{\mathit{i}}\cdot \cos {\alpha }_{\mathit{i}})/{\mathit{L}}_1),$$(35)

$${\beta }_{\mathit{i}}=\mathit{a}\tan (({\mathit{F}}_{\mathit{i}}\cdot \cos {\alpha }_1+{\mathit{G}}_{\mathit{i}}\cdot \sin {\alpha }_1)/({\mathit{N}}_{\mathit{i}})).$$(36)

Symbol F_i, G_i and N_i are simplification parameters.

The kinematics closed-loop vector equation of mechanical system is established as

$${\mathit{O}}_{\mathit{A}}{\mathit{o}}_7={}^{\mathit{A}}\mathit{P}_{\mathit{S}\mathit{i}}+{\mathit{S}}_{\mathit{i}}{\mathit{R}}_{\mathit{i}}-{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}}.$$(37)

According to the coordinate correspondence, the following equations can be obtained

$$\{\begin{array}{l}\mathit{x}=\mathit{f}({\mathit{P}}_1)+{\mathit{S}}_1{\mathit{R}}_1(\mathrm{1,1})+{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}\mathit{x}}\\ \mathit{y}=\mathit{f}({\mathit{P}}_2)+{\mathit{S}}_2{\mathit{R}}_2(\mathrm{2,1})+{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}\mathit{y}}\\ \mathit{z}=\mathit{f}({\mathit{P}}_3)+{\mathit{S}}_3{\mathit{R}}_3(\mathrm{3,1})+{}^{\mathit{A}}\mathit{P}_{\mathit{o}7\mathit{R}\mathit{i}\mathit{z}}\end{array},$$(38)

where, symbol f(P_i) is driving function, by simplified equation (38), f(P_i) can be expressed as

$$\left\{\begin{array}{l}\mathit{f}\left({\mathit{P}}_3\right)=\frac{-\left(2\bullet \mathit{C}\bullet \mathit{D}-2\bullet \mathit{A}\right)+\sqrt{{\left(2\bullet \mathit{C}\bullet \mathit{D}-2\bullet \mathit{A}\right)}^{\text{2}}\text{-}4\bullet \left({\mathit{C}}^2+1\right)\bullet \left({\mathit{A}}^2-\mathit{B}+{\mathit{D}}^2\right)}}{2\bullet ({\mathit{C}}^2+1)}\hfill \\ \mathit{f}\left({\mathit{P}}_1\right)={\mathit{b}}_1\bullet {\mathit{w}}_{\mathit{z}}+\left(\left(\mathit{f}\left({\mathit{P}}_3\right)-\mathit{x}\bullet \mathit{t}\mathit{a}\mathit{n}\beta \right)\bullet \mathit{s}\mathit{i}\mathit{n}\beta \bullet \mathit{t}\mathit{a}\mathit{n}\gamma -\mathit{x}\bullet \mathit{c}\mathit{o}\mathit{s}\beta \bullet \mathit{t}\mathit{a}\mathit{n}\gamma \right)+\sqrt{{\mathit{L}}_1^{\text{2}}-{\left({\mathit{b}}_1\bullet {\mathit{w}}_{\mathit{x}}+\mathit{x}\right)}^2-{\left({\mathit{b}}_1\bullet {\mathit{w}}_{\mathit{y}}+\mathit{x}\bullet \mathit{t}\mathit{a}\mathit{n}\beta \right)}^2}\hfill \\ \mathit{f}\left({\mathit{P}}_2\right)={\mathit{b}}_1\bullet {\mathit{w}}_{\mathit{z}}-\left(\left(\mathit{f}\left({\mathit{P}}_3\right)-\mathit{x}\bullet \mathit{t}\mathit{a}\mathit{n}\beta \right)\bullet \mathit{s}\mathit{i}\mathit{n}\beta \bullet \mathit{t}\mathit{a}\mathit{n}\gamma -\mathit{x}\bullet \mathit{c}\mathit{o}\mathit{s}\beta \bullet \mathit{t}\mathit{a}\mathit{n}\gamma \right)+\sqrt{{\mathit{L}}_2^{\text{2}}-{\left(\mathit{x}-{\mathit{b}}_1\bullet {\mathit{w}}_{\mathit{x}}\right)}^2-{\left(\mathit{x}\bullet \mathit{t}\mathit{a}\mathit{n}\beta -{\mathit{b}}_1\bullet {\mathit{w}}_{\mathit{y}}\right)}^2}\hfill \end{array}\right).$$(39)

The symbol A, B, C and D is simplification parameters.

The position coordinate O_Ao₇ = [x y z]^T of terminal platform is a function of variable x, β and γ. So linear velocity v_p and angular velocity w_p of terminal platform can be expressed by independent variable parameter${\displaystyle \dot{\mathit{H}}}=[{\displaystyle \dot{\mathit{x}}}{\displaystyle \dot{\beta }}{\displaystyle \dot{\gamma }}]$, it can be expressed as

$${\mathit{v}}_{\mathit{p}}=\left[\begin{array}{ccc}\frac{\partial \mathit{x}}{\partial \mathit{x}}& \frac{\partial \mathit{x}}{\partial \beta }& \frac{\partial \mathit{x}}{\partial \gamma }\\ \frac{\mathit{d}\gamma }{\partial \mathit{x}}& \frac{\partial \gamma }{\partial \beta }& \frac{\partial \gamma }{\partial \gamma }\\ \frac{\partial \mathit{z}}{\partial \mathit{x}}& \frac{\partial \mathit{z}}{\partial \beta }& \frac{\partial \mathit{z}}{\partial \gamma }\end{array}\right]\bullet \left[\begin{array}{c}\dot{\mathit{x}}\\ \dot{\beta }\\ \dot{\gamma }\end{array}\right]={\mathit{J}}_{1\mathit{A}}\left[\begin{array}{c}\dot{\mathit{x}}\\ \dot{\beta }\\ \dot{\gamma }\end{array}\right]$$(40)

Matrix J_1A is linear velocity transformation matrix.

$$\begin{array}{l}{\mathit{w}}_{\mathit{p}}=\mathit{R}(\mathit{x},\alpha )\cdot \left[\begin{array}{l}1\\ 0\\ 0\end{array}\right]\cdot {\displaystyle \dot{\alpha }}+\mathit{R}(\mathit{z},\beta )\cdot \left[\begin{array}{l}0\\ 0\\ 1\end{array}\right]\cdot {\displaystyle \dot{\beta }}+\mathit{R}(\mathit{y},\gamma )\cdot \left[\begin{array}{l}0\\ 1\\ 0\end{array}\right]\cdot {\displaystyle \dot{\gamma }}\\ =\left[\begin{array}{ccc}1& 0& -\sin \beta \\ 0& -\sin \alpha & \cos \alpha \cdot \cos \beta \\ 0& \cos \alpha & \cos \beta \cdot \sin \alpha \end{array}\right]\cdot \left[\begin{array}{l}{\displaystyle \dot{\alpha }}\\ {\displaystyle \dot{\beta }}\\ {\displaystyle \dot{\gamma }}\end{array}\right]\\ ={\mathit{J}}_{2\mathit{A}}\cdot \left[\begin{array}{l}{\displaystyle \dot{\alpha }}\\ {\displaystyle \dot{\beta }}\\ {\displaystyle \dot{\gamma }}\end{array}\right]\end{array}$$(41)

The matrix J_2A is angular velocity transformation matrix.

Velocity vector of rotating joints in world coordinate system is written as

$${}^{\mathit{A}}\mathit{R}_{\mathit{i}\mathit{v}}={\mathit{v}}_{\mathit{p}}+{\mathit{w}}_{\mathit{p}}\times ({}^{\mathit{A}}\mathit{R}_7\cdot {}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}}).$$(42)

Also, the acceleration vector of rotating joints in world coordinate system can be written as

$${}^{\mathit{A}}\mathit{R}_{\mathit{i}\mathit{a}}={{\displaystyle \dot{\mathit{v}}}}_{\mathit{p}}+{{\displaystyle \dot{\mathit{w}}}}_{\mathit{p}}\times ({}^{\mathit{A}}\mathit{R}_7\cdot {}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}})+{\mathit{w}}_{\mathit{p}}\times ({\mathit{w}}_{\mathit{p}}\times ({}^{\mathit{A}}\mathit{R}_7\cdot {}^{\mathit{o}7}\mathit{P}_{\mathit{R}\mathit{i}})).$$(43)

3.2 Dynamic model with spherical clearance joints

The generalized coordinates of each component are selected to accurately describe the pose changes of the mechanism, it is expressed as

$$\mathit{Q}=({{\mathit{Q}}_1}^{\mathit{T}}{{\mathit{Q}}_2}^{\mathit{T}}{{\mathit{Q}}_3}^{\mathit{T}}{{\mathit{Q}}_4}^{\mathit{T}}{{\mathit{Q}}_5}^{\mathit{T}}{{\mathit{Q}}_6}^{\mathit{T}}{{\mathit{Q}}_7}^{\mathit{T}}),$$(44)

where ${\mathit{Q}}_{\mathit{i}}=({\mathit{O}}_{\mathit{i}}^{\mathit{T}}{\mathit{\phi }}_{\mathit{i}}^{\mathit{T}})$,O_i = (x_i y_i z_i)^T, φ_i = (α_i β_i γ_i)^T(i = 1, 2, ..., 7), O_i is position coordinate of the origin of object coordinate system in world coordinate system, ϕ_i is Euler angle vector of object coordinate system in world coordinate system.

Spherical joints restrict three-dimensional movement and provide three-dimensional rotation, so its constraint equation is written as

$${\Phi }_{\mathit{S}}={\mathit{O}}_{\mathit{A}}{\mathit{o}}_{\mathit{i}}+{}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{o}}_{\mathit{i}}{\mathit{c}}_{\mathit{i}}-{\mathit{O}}_{\mathit{A}}{\mathit{o}}_{\mathit{j}}-{}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{o}}_{\mathit{j}}{\mathit{c}}_{\mathit{j}}=0_{3\times 1},$$(45)

where ^OA_oi and ^OA_oj are position coordinate of origin of object coordinate system {o_i} and {o_j} in world coordinate system{A}, o_ic_i is position coordinate about center point of socket in object coordinate system {o_i}, o_jc_j is position coordinate about center point of ball in object coordinate system {o_j},^AR_i  is Euler transformation matrix.

Sliding joints only provide translation in one direction and limits the pose in the other five directions, its constraint equation can be written as

$${\Phi }_{\mathit{P}}=\left[\begin{array}{c}{{\mathit{d}}_{\mathit{z}\mathit{j}}}^{\text{T}}\cdot ({\mathit{O}}_{\mathit{A}}{\mathit{p}}_{\mathit{j}}-{\mathit{O}}_{\mathit{A}}{\mathit{o}}_{\mathit{i}}-{}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{o}}_{\mathit{i}}{\mathit{p}}_{\mathit{i}})\\ {{\mathit{d}}_{\mathit{y}\mathit{j}}}^{\text{T}}\cdot ({\mathit{O}}_{\mathit{A}}{\mathit{p}}_{\mathit{j}}-{\mathit{O}}_{\mathit{A}}{\mathit{o}}_{\mathit{i}}-{}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{o}}_{\mathit{i}}{\mathit{p}}_{\mathit{i}})\\ {{\mathit{d}}_{\mathit{y}\mathit{j}}}^{\text{T}}\cdot ({}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{d}}_{\mathit{x}\mathit{i}})\\ {{\mathit{d}}_{\mathit{z}\mathit{j}}}^{\text{T}}\cdot ({}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{d}}_{\mathit{x}\mathit{i}})\\ {{\mathit{d}}_{\mathit{z}\mathit{j}}}^{\text{T}}\cdot ({}^{\mathit{A}}\mathit{R}_{\mathit{i}}\cdot {\mathit{d}}_{\mathit{y}\mathit{i}})\end{array}\right]=0_{5\times 1},$$(46)

where ^AR_i is Euler matrix, P_i and P_j are center of two components that make up sliding joint, d_xi, d_yi, d_zi are mutually perpendicular unit vectors fixed at d_xj, d_yj, d_zj are mutually perpendicular unit vectors fixed at P_j, d_xi and d_xj are parallel to the translational direction of sliding joint. The first two equations restrict the movement of the sliding joint in the other two directions, and the last three equations restrict the rotation in three directions.

The rotating pair provides a rotational pose around the axis, so the movement in three directions and rotation in the other two directions should be restricted, and the constraint equation can be written as

$${\Phi }_R\left[\begin{array}{l}{\mathit{O}}_{\mathit{A}}{\mathit{o}}_{\mathit{i}}{+}^{\mathit{A}}{\mathit{R}}_{\mathit{i}}\bullet {\mathit{o}}_{\mathit{i}}{\mathit{r}}_{\mathit{i}}-{\mathit{O}}_{\mathit{A}}{\mathit{o}}_{\mathit{j}}{-}^{\mathit{A}}{\mathit{R}}_{\mathit{j}}\bullet {\mathit{o}}_{\mathit{j}}{\mathit{r}}_{\mathit{j}}\\ {(}^{\mathit{A}}{\mathit{R}}_{\mathit{i}}\bullet {\mathit{u}}_{\mathit{i}})\times {(}^{\mathit{A}}{\mathit{R}}_{\mathit{i}}\bullet {\mathit{u}}_{\mathit{j}})\end{array}\right]$$(47)

where ^AR_i is Euler matrices, u_i and u_j are rotational axis vectors of shaft and bearing, ^OA^o_i and ^OA^o_j are position coordinate of origin of object coordinate system in world coordinate system, r_i and r_j are center points of shaft and bearing

The rigid constraint equation of SPCMM including all joints can be written as

$$\Phi (\mathit{Q})=\left[\begin{array}{l}{{\Phi }_{\mathit{P}1}}^{\text{T}}{{\Phi }_{\mathit{P}2}}^{\text{T}}{{\Phi }_{\mathit{P}3}}^{\text{T}}{{\Phi }_{\mathit{S}1}}^{\text{T}}{{\Phi }_{\mathit{S}2}}^{\text{T}}\\ {{\Phi }_{\mathit{S}3}}^{\text{T}}{{\Phi }_{\mathit{R}1}}^{\text{T}}{{\Phi }_{\mathit{R}2}}^{\text{T}}{{\Phi }_{\mathit{R}3}}^{\text{T}}{{\Phi }_{\mathit{D}}}^{\text{T}}\end{array}\right]=0_{42\times 1},$$(48)

where ${\Phi }_{\mathit{D}}=\left[\begin{array}{l}\left|\mathit{P}{\mathit{P}}_1\right|-\mathit{f}({\mathit{P}}_1)\\ \left|\mathit{P}{\mathit{P}}_2\right|-\mathit{f}({\mathit{P}}_2)\\ \left|\mathit{P}{\mathit{P}}_3\right|-\mathit{f}({\mathit{P}}_3)\end{array}\right]=0_{3\times 1}$, f(P_i) is driving function of translation joints.

Formula (48) derivation of time t, the velocity constraint equation is written as

$${\Phi }_{\mathit{Q}}{\displaystyle \dot{\mathit{Q}}}=-{\Phi }_{\mathit{t}}\equiv \sigma$$(49)

where Φ_Q is Jacobean matrix, ${\displaystyle \dot{\mathit{Q}}}$ is generalized velocity vector, Φ_t means that constraint equation Φ(Q) derivate of time t, σ  is the other side of the equation.

The equation (49) derivate of time t, acceleration constraint equation is written as

$${\Phi }_{\mathit{Q}}{\displaystyle \stackrel{..}{\mathit{Q}}}=-{({\Phi }_{\mathit{Q}}{\displaystyle \dot{\mathit{Q}}})}_{\mathit{Q}}{\displaystyle \dot{\mathit{Q}}}-2{\Phi }_{\mathit{Q}\mathit{t}}{\displaystyle \dot{\mathit{Q}}}-{\Phi }_{\mathit{t}\mathit{t}}\equiv \gamma ,$$(50)

where ${\displaystyle \stackrel{..}{\mathit{Q}}}$ represents generalized acceleration vector of mechanism which means ${\displaystyle \dot{\mathit{Q}}}$ derivate of time t, ${({\Phi }_{\mathit{Q}}{\displaystyle \dot{\mathit{Q}}})}_{\mathit{Q}}$ means ${\Phi }_{\mathit{Q}}{\displaystyle \dot{\mathit{Q}}}$ derivate of the generalized coordinate Q, Φ_Qt means Φ_Q derivate of time t, Φ_tt means Φ_t derivate of time t, γ is other side of the equation.

Mass matrix of SPCMM is followed

$${\mathit{M}}_{\mathit{Q}}=\mathit{d}\mathit{i}\mathit{a}\mathit{g}[{\mathit{N}}_1{\mathit{J}}_1{\mathit{N}}_2{\mathit{J}}_2\cdot \cdot \cdot {\mathit{N}}_7{\mathit{J}}_7],$$(51)

where N_k = diag[m_k m_k m_k] (k = 1, 2...7), the rotational inertia J_k = diag[I_xk I_yk I_zk] (k = 1, 2...7).

The above results are substituted into dynamic equation with Lagrange multiplier, equation may write as follow

$${\mathit{M}}_{\mathit{Q}}{\displaystyle \stackrel{..}{\mathit{Q}}}+{{\Phi }_{\mathit{Q}}}^{\text{T}}\lambda ={\mathit{F}}_{\mathit{Q}}$$(52)

where λ is Lagrange multiplier, it only related to the internal force and internal moment of connected rods, F_Q is generalized force vector and it included external force and external torque.

When SPCMM carried out with no load, the generalized force vector F_Q is expressed as

$${\mathit{F}}_{\mathit{Q}}={[{\mathit{M}}_1{\mathit{M}}_2\mathrm{...}{\mathit{M}}_7]}^{\text{T}},$$(53)

where M_k = [m_kg 0 0 0 0 0]^T (k = 1, 2...7), g is the acceleration of gravity.

Expression of formula (50) and (52) in form of differential algebra, the union matrix is shown as

$$\left[\begin{array}{l}{\mathit{M}}_{\mathit{Q}}{{\Phi }_{\mathit{Q}}}^{\text{T}}\\ {\Phi }_{\mathit{Q}}0\end{array}\right]\left[\begin{array}{l}{\displaystyle \stackrel{..}{\mathit{Q}}}\\ \lambda \end{array}\right]=\left[\begin{array}{l}{\mathit{F}}_{\mathit{Q}}\\ \gamma \end{array}\right].$$(54)

Because displacement constraint and velocity constraint of multibody system are not used in equation (54), so there may be a violated constraint. Baumgarte [34] proposed a method to solve this problem using feedback control theory to enhance the stability of numerical calculation, the new union matrix is written as

$$\left[\begin{array}{l}{\mathit{M}}_{\mathit{Q}}{{\Phi }_{\mathit{Q}}}^{\text{T}}\\ {\Phi }_{\mathit{Q}}0\end{array}\right]\left[\begin{array}{l}{\displaystyle \stackrel{..}{\mathit{Q}}}\\ \lambda \end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_{\mathit{Q}}\\ \gamma -2{\alpha }_{\mathit{c}}{\displaystyle \dot{\Phi }}-{{\beta }_{\mathit{c}}}^2\Phi \end{array}\right],$$(55)

where a_c and β_c are correction parameter bigger than zero, ${\displaystyle \dot{\Phi }}$ is that the rigid constraint equation Φ(Q) derivate of t.

In clearance dynamics model, the ball moves freely in the socket, the force constraint in spherical joint will replace the displacement constraint in rigid dynamics model, and the constraint equation and generalized force vector considering spherical clearance joints S₁, S₂ and S₃ are as follows:

$${}^{*}\Phi (\mathit{Q})=\left[{{\Phi }_{\mathit{P}1}}^{\text{T}}{{\Phi }_{\mathit{P}2}}^{\text{T}}{{\Phi }_{\mathit{P}3}}^{\text{T}}{{\Phi }_{\mathit{R}1}}^{\text{T}}{{\Phi }_{\mathit{R}2}}^{\text{T}}{{\Phi }_{\mathit{R}3}}^{\text{T}}{{\Phi }_{\mathit{D}}}^{\text{T}}\right]=0_{33\times 1},$$(56)

$${}^{\ast }\mathit{F}_{\mathit{Q}}={[{}^{\ast }\mathit{M}_1{}^{\ast }\mathit{M}_2\mathrm{...}{}^{\ast }\mathit{M}_7]}^{\text{T}}.$$(57)

In formula (57),

$${}^{\ast }\mathit{M}_{\mathit{k}}=[{\mathit{m}}_{\mathit{k}}\mathit{g}+{\mathit{F}}_{\mathit{i}\mathit{d}}^{\mathit{S}\mathit{w}}0\text{0}{\mathit{M}}_{\mathit{i}}^{\mathit{S}\mathit{w}}](\mathit{k}=\mathrm{1,3,5}\mathit{d}=\mathit{x},\mathit{y},\mathit{z}),$$

$${}^{\ast }\mathit{M}_7=[{\mathit{m}}_7\mathit{g}00\text{0 0 0}]$$

where the force ${\mathit{F}}_{\mathit{i}\mathit{d}}^{\mathit{S}\mathit{w}}$ is component of contact force of spherical clearance joint in world coordinate system, ${\mathit{M}}_{\mathit{i}}^{\mathit{S}\mathit{w}}$ is the moment of ball and socket from collision point to component centroid.

The union matrix with spherical clearance joints is written as

$$\left[\begin{array}{l}{\mathit{M}}_{\mathit{Q}}{}^{\ast }{{\Phi }_{\mathit{Q}}}^{\text{T}}\\ {}^{\ast }\Phi _{\mathit{Q}}0\end{array}\right]\left[\begin{array}{l}{\displaystyle \stackrel{..}{\mathit{Q}}}\\ \lambda \end{array}\right]=\left[\begin{array}{c}{}^{\ast }\mathit{F}_{\mathit{Q}}\\ \gamma -2\alpha {\displaystyle \dot{\Phi }}-{\beta }^2\Phi \end{array}\right].$$(58)

4 Dynamics behavior analysis of SPCMM

4.1 Solving process and applied parameters

Numerical technique is used to solve dynamic model of SPCMM with spherical clearance joints, and detailed structural parameters and simulation parameters is written in Tables 1 and 2. Simulation trajectory is given as following

$$\{\begin{array}{l}\mathit{x}=\sin (\pi \cdot \mathit{t})/40+0.203\\ \beta =\sin (\pi \cdot \mathit{t})/40\\ \gamma =\pi /12\end{array}$$(59)

According to Figure 4, the detailed numerical simulation process can be summarized as

• Solve initial displacement and velocity of each component at time=0 by inverse kinematics, these data will be loaded into dynamics program with multiple spherical clearance joints.

• Solve displacement constraint equations including all prismatic joints, spherical joints and rotating joints. The velocity and acceleration constraint equations can be obtained by taking time derivatives of the displacement constraint equations.

• Establish kinematics model of spherical clearance joint, and solve normal contact force and tangential contact force at the impact point.

• The force constraint in spherical clearance joint is introduced into Baumgarte stability equation instead of position constraint of component to enhance the stability of numerical simulation.

• Judge the contact state at the contact point and calculate collision depth, if collision depth is greater than 0, calculate normal contact force and tangential contact force, and generalized displacement, velocity and acceleration can be got at the same time.

• Perform iterative calculation until the end of motion time.

Fig. 4 Flow chart of numerical simulation.

4.2 Dynamic response of SPCMM with multiple spherical clearance joints

4.2.1 Dynamic response considering clearance values

It is inevitable that clearance will cause an unpredictable impact and reduce the performance of mechanism. Concreting dynamic response and chaotic characteristics of parallel mechanism under different clearances can provide theoretical guidance for selecting the appropriate clearance in production and assembly, and set the clearance value as the single variable, set other variables to fixed values and the friction coefficient is fixed at 0.2.

From Figure 5a, the displacement curves of the terminal platform in X direction has no obvious deviation under different clearance values, and three curves coincided together. In the Figure 5b, curves have obvious fluctuations within the 0–0.3 s, and the amplitude of fluctuation increases while clearance value becomes bigger. Maximum deviation is 0.012 rad when clearance value is 0.7 mm, maximum deviation is 0.008 rad when clearance value is 0.5 mm, and maximum deviation is 0.006 rad when clearance value is 0.3 mm. After 0.3 s, the curves basically coincide and mechanism start a stable operation state. It is a conclusion change of clearance values have a great influence on the displacement response in Beta in the initial stage of movement, curve deviation becomes bigger when matched a bigger clearance value.

From Figure 6a, It is obvious clearance values caused a big different fluctuations within 0–0.4 s. Maximum fluctuation peak is 0.186 m/s while matched a fixed clearance value 0.7 mm, maximum fluctuation peak is 0.167 m/s while matched a fixed clearance value 0.5 mm, maximum fluctuation peak is 0.147 m/s while matched a fixed clearance value 0.3 mm. In Figure 6b, and the fluctuation is greater when matched a bigger clearance value, maximum fluctuation peak is 1.511 rad/s while matched a fixed clearance value 0.7 mm, maximum fluctuation peak is 1.104 rad /s while matched a fixed clearance value 0.5 mm, maximum fluctuation peak is 0.906 rad /s while matched a fixed clearance value 0.3 mm, but the fluctuations weaken very quickly, the mechanical system becomes stable and the curves basically coincide when it is time 0.4 s. It can be concluded that clearance has a great effect on the velocity in X and Beta of mechanism, when the impact phenomenon in joints disappear and the spherical joints elements are in contact, the mechanism can work in a stable motion state.

From Figure 7a, acceleration fluctuation becomes larger while matched a greater clearance value, the peak acceleration in X direction is 667.5 m/s² when clearance value set a fixed value 0.7 mm, the peak acceleration in X direction is 446.2 m/s² when clearance value set a fixed value 0.5 mm, the peak acceleration is 430.7 m/s² when clearance value set a fixed value 0.3 mm. In Figure 7b, the acceleration in Beta is bigger than X, this means the clearance value does a great effect on acceleration in Beta, the initial fluctuation range of acceleration is very large, and then it weakened rapidly, until the time comes to 0.4 s the motion state becomes stable. The reason for this fluctuation rule of acceleration in X and Beta is that coupling motion and mutual interference occurs between the three spherical joints.

In Figure 8, contact forces fluctuate greatly within 0–0.21 s in different joints, violent collision occurred between joints elements at start of motion, and fluctuation peaks is different. Fluctuation peak is 809.7 N in S1 joint when clearance value matched a fixed value 0.3 mm, fluctuation peak is 463.3 N in S2 joint and fluctuation peak is 755.8 N in S3 joint, after 0.21 s, collision disappeared and mechanism are working steadily, here is a conclusion impact phenomenon in S1 is more obvious than impact phenomenon in S2 and S3. In Figure 9, the contact forces fluctuate greatly within 0–0.33 s in different joints S1, S2 and S3, after 0.33 s, collision disappeared and mechanism are working smoothly. When clearance value matched a fixed value 0.5 mm, fluctuation peak is 840 N in S1 joint, fluctuation peak is 1061 N in S2 joint and fluctuation peak is 982.7 N in S3 joint, it is concluded impact phenomenon in S2 is more obvious than impact phenomenon in S1 and S3. In Figure 10, contact forces fluctuate greatly within 0–0.39 s in different joints S1, S2 and S3, after 0.39 s, collision disappeared and mechanism are working smoothly. When clearance value matched a fixed value 0.7 mm, fluctuation peak is 1094 N in S1 joint, fluctuation peak is 917.4 N in S2 joint and fluctuation peak is 1176 N in S3 joint, it is concluded impact phenomenon in S3 is more obvious than impact phenomenon in S1 and S2. Compare contact forces considering different clearance values in joints S1, S2 and S3 and find that contact force increased while clearance values increased, it takes longer time for mechanism to reach stable motion state when matched a bigger clearance value.

Fig. 5 Displacement response in X and Beta under different clearance values.

Fig. 6 Velocity response in X and Beta under different clearance values.

Fig. 7 Acceleration response in X and Beta under different clearance values.

Fig. 8 Contact forces of S1, S2 and S3 with c is 0.3 mm.

Fig. 9 Contact forces of S1, S2 and S3 with c is 0.5 mm.

Fig. 10 Contact forces of S1, S2 and S3 with c is 0.7 mm.

4.2.2 Dynamic response considering friction coefficients

Oversized contact force existing in spherical clearance joint caused by violent collision may lead to an unbalanced motion state of the mechanism, and the main influence factor of the contact force is the friction coefficient, so it is particularly critical to research the dynamic response and chaotic characteristics of mechanism under different friction coefficients, here the clearance value takes a fixed 0.3 mm.

In Figure 11a, the displacement curves of different friction coefficients coincide, it can be concluded that different friction coefficients have little effect on the displacement curves in X direction. In Figure 11b, the curves fluctuate significantly before 0.3 s, and the fluctuation amplitude increases with increase of friction coefficient. Maximum deviation is 6.61 × 10⁻³ rad when friction coefficient is 0.15, maximum deviation is 6.77 × 10⁻³ rad when friction coefficient is 0.20, and maximum deviation is 6.90 × 10⁻³ rad when friction coefficient is 0.25. Fluctuation amplitude decreases gradually with movement of the mechanism, after 0.3 s, the curve does not fluctuate and the curves coincide. It can be concluded that friction coefficient has much effect on displacement curve at beginning of movement, the greater the friction coefficient, the greater the fluctuation of the curve.

From Figure 12a, the velocity curves in X fluctuate obviously within 0–0.3 s, and the peak of fluctuation reached 0.161 m/s, 0.156 m/s and 0.149 m/s separately while cf is 0.25, 0.20 and 0.15, and the curves fluctuate significantly before 0.3 s, the fluctuation amplitude decreases gradually within 0–0.3 s, and after 0.3 s, the curve has no obvious fluctuation. Here is a concluded that friction coefficient does a great effect on the velocity curve in X within 0–0.3 s, after 0.3 s, the friction coefficients have little effect on the velocity curve in X. In Figure 12b, the fluctuation of velocity curve in Beta is very obvious within 0–0.3 s, different friction coefficients lead to different fluctuation peaks, the fluctuation peak corresponding to friction coefficient 0.25 is 1.279 rad/s, the fluctuation peak corresponding to friction coefficient 0.20 is 1.277 rad/s, the fluctuation peak corresponding to friction coefficient 0.15 is 1.265 rad/s, after 0.3 s the fluctuation of curve is not obvious and curves basically coincide. It can be concluded that different friction coefficients have much effect on the velocity curve in Beta within 0–0.3 s, after 0.3 s, the friction coefficient has little effect on the curve, the mechanism reaches a stable motion state.

In Figure 13, the acceleration curves in X and Beta fluctuate violently within 0–0.3 s, after 0.3 s, the curves do not fluctuate and the curves present a stable periodic working state. In Figure 13a, biggest peak value is 533.2 m/s² while friction coefficient matched a fixed value 0.15, peak value is 430.7 m/s² while friction coefficient matched a fixed value 0.20, the peak value is 426 m/s² while friction coefficient matched a fixed value 0.25, Here is a concluded that the greater the friction coefficient, the greater the fluctuation peak. In Figure 13b, the biggest peak value is 2193 m/s² while friction coefficient is 0.20, the peak value is 1582 m/s² while friction coefficient matched a fixed value 0.25 and peak value is 993.1 m/s² while friction coefficient matched a fixed value 0.15. Here is a concluded that due to existence of coupled motion, influence of friction coefficient on acceleration in Beta is not following the rule the greater the friction coefficient, the greater the peak value of the curve, the acceleration curve in Beta has higher sensitivity than the acceleration curve in X when the friction coefficient changes.

From Figure 14, the contact force curves fluctuate greatly within 0–0.22 s, after 0.22 s, mechanism are working steadily. Peak value is 1057 N in joint S2, peak value is 966.4 N in joint S1 and the peak value is 866.8 N in joint S3, it is a conclusion collisions in S2 is more severe. From Figure 15, the contact force curves fluctuate greatly within 0–0.23 s, after 0.23 s, mechanism are working steadily. Peak value is 922.3 N in joint S2, peak value is 552.6 N in joint S1 and peak value is 444.7 N in joint S3, it is a conclusion collisions in S2 is more severe. In Figure 16, the contact force curves fluctuate greatly within 0–0.22 s, after 0.22 s, mechanism are working steadily. Peak value is 719.4 N in joint S1, peak value is 560.8 N in joint S2 and peak value is 500.2 N in joint S3, it is a conclusion collisions in S1 is more severe.

Fig. 11 Displacement response in X and Beta under friction coefficients.

Fig. 12 Velocity response in X and Beta under different friction coefficients.

Fig. 13 Acceleration response in X and Beta under different friction coefficients.

Fig. 14 Contact forces of S1 S2 and S3 with cf is 0.15

Fig. 15 Contact forces of S1 S2 and S3 with cf is 0.20.

Fig. 16 Contact forces of S1 S2 and S3 with cf is 0.25.

4.2.3 Reliability verification of dynamic model with multiple spherical clearance joints

To test the reliability of dynamic model for SPCMM with multiple spherical clearance joints, Matlab calculation results, Adams simulation results and ideal motion results without clearance joints were compared, the comparison results are shown in Figure 17.

Figure 17 shows simulation comparison curves between Adams and Matlab. Figure 17a and 17b shows the displacement curves in X and Beta, and Adams results, Matlab result and the no clearance results basically consistent. Figure 17c and 17d shows that velocity curves fluctuate greatly at the beginning of operation, the results solved by Matlab is bigger than the results solved by Adams, the curves are basically consistent when the mechanism reaches the stable state. Figure 17(e) and (f) shows that the acceleration curves fluctuate greatly at the beginning of operation, the peak fluctuation of Adams curve in X is 291 m/s², the peak fluctuation of Matlab curve in X is 319.2 m/s², the peak fluctuation of Adams curve in Beta is 1571 rad/s², the peak fluctuation of Matlab curve in Beta is 1754 rad/s², when the mechanism reaches the stable state, curves are basically consistent. It can be concluded that reliability of dynamics model of SPCMM with multiple spherical clearance joints is verified through comparison results.

Fig. 17 Simulation comparison curves of Beta between Adams and Matlab.

4.3 Chaotic characteristics of SPCMM with multiple spherical clearance joints

Through Phase diagram, Poincare map and Bifurcation diagram of terminal platform in X and Beta direction, effects of different clearance values and friction coefficients on chaotic characteristics of mechanism can be researched and analyzed intuitively.

From Figure 18a, the Phase diagrams in X coincide basically and the images are smooth ellipses. It can be concluded different clearance values have little effect on Poincare map in X, the mechanism is in periodic motion in X analyzed by Poincare map. In Figure 18b, the picture is presented as ellipses of different widths, and the width of the ellipse becomes bigger when increases clearance value, it can be concluded that with increase of clearance value, influence of clearance on Poincare map in Beta becomes more obvious.

In Figure 19a, the Poincare maps in Beta corresponding to different clearance values are presented as three points at different positions, and the larger the clearance value, the more right the point is. It can be concluded that mechanism is in quasi periodic motion at this moment. From Figure 19b, the image is represented as three circular areas, and circular area becomes bigger when increases clearance value, it can be concluded that clearance value has much influence on the Poincare map in Beta, the mechanism is in quasi periodic motion state in Beta.

In Figure 20a, the curve has no divergence and presents as a stable straight line, it is resulted that mechanism is in periodic motion state all time while clearance value increases from 0.01 mm to 1 mm. In Figure 20b, the image gradually diverges while increases clearance value, and it can be concluded the clearance value does a great effect on the chaotic characteristics in Beta direction, and mechanism varied from periodic motion state to chaotic motion state with increase of clearance value.

In Figure 21a, the image appears as three smooth ellipses which coincide together, it can be concluded that friction coefficient has little effect on the phase diagram in X, and the mechanism is in periodic motion under different friction coefficients. From Figure 21b, there are three ellipses with wider contours which coincide together, the contours of the ellipse do not become wider when the friction coefficient becomes bigger, here is a concluded that friction coefficient has less effect on the phase diagram in Beta, and mechanism is in quasi periodic motion.

From Figures 22a and (b), the Poincare maps corresponding to different friction coefficients are located in the same position, the closer the mapping points are to the center, the denser the distribution is. It can be concluded that changing friction coefficient has little effect on Poincare mapping in X and Beta, and the mechanism is in quasi periodic motion.

From Figure 23, the Bifurcation diagrams of terminal platform in X and Beta present two convergent straight lines when friction coefficient increases from 0.01 to 0.4, it can be concluded that the mechanism is always in periodic working state, friction coefficient has little effect on the Bifurcation diagram of terminal platform in X and Beta.

Fig. 18 Phase diagram of terminal platform.

Fig. 19 Poincare map of terminal platform.

Fig. 20 Bifurcation diagram of terminal platform.

Fig. 21 Phase diagram of terminal platform.

Fig. 22 Poincare map of terminal platform.

Fig. 23 Bifurcation diagram of terminal platform.

5 Conclusions

As an indispensable resource for studying dynamic behavior of parallel mechanism, dynamic response and chaos characteristics of spatial parallel coordinate measuring mechanism with multi spherical clearance joints are studied.

• A dynamic modeling method with multiple spherical clearance joints for spatial parallel coordinate measuring mechanism based on Lagrange multiplier technology is developed. The reliability of dynamics modeling method is verified by comparing the results from Adams and Matlab, also a detailed contact collision model of spherical clearance joints under different collision states are established.

• Dynamic response of spatial parallel coordinate measuring mechanism is studied by comparing displacement, velocity, acceleration curves and contact force curves. The clearance values are respectively taken as fixed values of 0.3 mm, 0.5 mm and 0.7 mm and friction coefficient are respectively taken as fixed values of 0.15, 0.20 and 0.25. Comparing the analysis results, it is found that dynamic response of mechanism becomes more violent with increase of clearance value and friction coefficient.

• Chaos characteristics of spatial parallel coordinate measuring mechanism are researched by comparing Phase diagram, Poincare map and bifurcation diagram. Comparing the analysis results, get a conclusion that clearance value has a great influence on chaotic characteristics of Beta direction, chaotic characteristics became more obvious when matched a bigger clearance value, and different friction coefficient has little effect on Chaos characteristics.

Funding

This research is supported by Shandong Provincial Natural Science Foundation (Grant no. ZR2022ME040).

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. The authors have no conflict to disclose.

Data availability statement

All data has already been reported in the manuscript.

Author contribution statement

Xiulong Chen came up with the idea and Yongchao Sheng wrote the article.

References

All Tables

All Figures

	Fig. 1 Characteristic map of SPCMM.
In the text

	Fig. 2 The kinematic scheme of SPCMM.
In the text

	Fig. 3 Spherical clearance joint model.
In the text

	Fig. 4 Flow chart of numerical simulation.
In the text

	Fig. 5 Displacement response in X and Beta under different clearance values.
In the text

	Fig. 6 Velocity response in X and Beta under different clearance values.
In the text

	Fig. 7 Acceleration response in X and Beta under different clearance values.
In the text

	Fig. 8 Contact forces of S1, S2 and S3 with c is 0.3 mm.
In the text

	Fig. 9 Contact forces of S1, S2 and S3 with c is 0.5 mm.
In the text

	Fig. 10 Contact forces of S1, S2 and S3 with c is 0.7 mm.
In the text

	Fig. 11 Displacement response in X and Beta under friction coefficients.
In the text

	Fig. 12 Velocity response in X and Beta under different friction coefficients.
In the text

	Fig. 13 Acceleration response in X and Beta under different friction coefficients.
In the text

	Fig. 14 Contact forces of S1 S2 and S3 with cf is 0.15
In the text

	Fig. 15 Contact forces of S1 S2 and S3 with cf is 0.20.
In the text

	Fig. 16 Contact forces of S1 S2 and S3 with cf is 0.25.
In the text

	Fig. 17 Simulation comparison curves of Beta between Adams and Matlab.
In the text

	Fig. 18 Phase diagram of terminal platform.
In the text

	Fig. 19 Poincare map of terminal platform.
In the text

	Fig. 20 Bifurcation diagram of terminal platform.
In the text

	Fig. 21 Phase diagram of terminal platform.
In the text

	Fig. 22 Poincare map of terminal platform.
In the text

	Fig. 23 Bifurcation diagram of terminal platform.
In the text