<figure><img src="image_1.png" /></figure>

HAL
open science

# Elasto-Dynamic Modeling of an Over-Constrained Parallel Kinematic Machine Using a Beam Model

Hélène Chanal, Aurélie Guichard, Benoît Blaysat, Stéphane Caro

► To cite this version:

Hélène Chanal, Aurélie Guichard, Benoît Blaysat, Stéphane Caro. Elasto-Dynamic Modeling of an Over-Constrained Parallel Kinematic Machine Using a Beam Model. Machines, 2022, 10 (3), pp.200. 10.3390/machines10030200. . hal-03604187

HAL Id: hal-03604187

https://hal.science/hal-03604187

Submitted on 18 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Elasto-Dynamic Modeling of an Over-Constrained Parallel Kinematic Machine Using a Beam Model

Hélène Chanal <fnref n="1" />,*, Aurélie Guichard <fnref n="2" />, Benoît Blaysat <fnref n="1" /> and Stéphane Caro <fnref n="3" />

**Abstract:** This article deals with the development of a simple model to evaluate the first natural frequencies of over-constrained parallel kinematic machines (PKMs). The simplest elasto-dynamic models are based on multi-body approaches. However, these approaches require an expression of the Jacobian matrices that may be difficult to obtain for complex PKMs. Therefore, this paper focuses on the determination of the global mass and stiffness matrices of an over-constrained PKM in stationary configurations without the use of Jacobian matrices. The PKM legs are modeled by beams. Because the legs are connected to a moving platform and the mechanism is over-constrained, constraint equations between the parameters that model the deformation of each leg are determined according to screw theory. The first natural frequencies and associated modes can then be determined. It should be noted that the proposed method can be easily used at the conceptual design stage of PKMs. The Tripteor X7 machine is used as an illustrative example and is characterized experimentally.

**Keywords:** parallel kinematic machine tool; over-constrained system; elasto-dynamic model; screw theory

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Chanal, H.; Guichard, A.; Blaysat, B.; Caro, S. Elasto-Dynamic Modeling of an Over-Constrained Parallel Kinematic Machine Using a Beam Model. *Machines* **2022**, *10*, 200. https://doi.org/10.3390/ machines10030200

Academic Editor: Antonio J. Marques Cardoso

Received: 4 February 2022

Accepted: 8 March 2022

Published: 10 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Parallel Kinematic Robots are used in many fields such as medicinal, aerospace, rehabilitation, and astronomy [1]. In industry, parallel kinematic robots are mainly used for pick-and-place operations [2]. However, a few Parallel Kinematic Machine Tools (PKMs) are used in the automotive or aeronautical industries for High-Speed Machining (HSM) operations [3]. Their dynamic performances, in terms of acceleration, are better for an equivalent motorization than Serial Kinematic Machine Tools (SKMs) due their closed-loop architecture [4]. However, the design of these closed-loop architectures imposes to control leg size that reduces mobile masses and leads to a loss of stiffness compared to SKM [5,6].

One way to improve PKM stiffness behavior is to use an over-constrained architecture. To this aim, an over-constrained PKM, named Exechon, was designed by Neumann [7], and the PCI-SCEMM company chose to integrate an Exechon robot into the design of Tripteor X7 machine tool (Figure 1). An over-constrained system is considered to be a hyperstatic system which is defined by the IFToMM terminology as a system in which the distribution of internal forces depends on the material properties of the members of the system [8]. Thus, an over-constrained PKM is defined as a PKM with common or redundant constraints that can be removed without changing the kinematic properties of the mechanism [9]. Thus, the study of an over-constrained mechanism behavior can be complex due to the generation of static indeterminacy and geometric constraints [10].

Previous work highlights the impact on part quality of PKM vibrations during machining [11]. Depending on the machined part position, marks appear on the machined surfaces. To minimize these vibrations, it is of prime importance to express a parameterized elasto-dynamic model for over-constrained PKMs. Such a model enables the prediction

<footnote n="1">¹ Clermont Auvergne INP, Institut Pascal, CNRS, Université Clermont Auvergne, F-63000 Clermont-Ferrand, France; benoit.blaysat@uca.fr</footnote>

<footnote n="2">² IUT Angers, Université Angers, 2 boulevard du Ronceray, F-49035 Angers, France; aurelie.guichard@univ-angers.fr</footnote>

<footnote n="3">³ CNRS, LS2N UMR CNRS 6004, 1 rue de la Noë, CEDEX 03, F-44321 Nantes, France; stephane.caro@ls2n.fr</footnote>

<footnote>* Correspondence: helene.chanal@sigma-clermont.fr</footnote>
of static and vibration behavior during the conceptual design stage. Moreover, as PKM behavior is anisotropic and changes regarding the tool pose, the elasto-dynamic model of the PKM should be fast to compute to quickly obtain its stiffness, its first natural frequencies and associated modes throughout the manipulator workspace.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Tripteor X7 Parallel Kinematic Machine Tools.</figcaption></figure>

Two main methods are classically used in the literature to express the elasto-dynamic model of an over-constrained PKM:

1. Finite Element Analysis (FEA): this method is the most accurate one, but it requires the exact definition of the mechanism components. It is usually used for validation purposes at the final design stage since it is time-consuming [12–15].

2. Multi-body approaches: PKM legs can be modeled, for example, using beam theory and their joint behavior relies on Virtual Joint Method (VJM) [16,17]. In the literature, the simplest elasto-dynamic models are based on multi-body approaches. This method can only be applied for robot architectures whose hypotheses are valid. Other formulations are introduced to consider large deformations of flexible manipulators [18] or to decrease simulation time for complex mechanisms such as Matrix Structural Analysis (MSA) [19]. The decreasing of the time simulation is generally based on the development of a methodology to merge stiffness and mass matrices of all elements of the mechanism [19,20].

In this article, multi-body approach, a robust but simple approach to obtain a param- eterized elasto-dynamic model, is preferred to the use of complex and time-consuming FEA. The prediction of vibration phenomena requires the determination of the first natural frequencies and their associated modes. In [16,17], this is performed by computing the global mass matrix and the global stiffness matrix. However, the computation of mass and stiffness matrices is complex in the case of over-constrained mechanisms such as Exechon PKM. Indeed, the complex coupling between displacement vector components due to over-constrained properties must be considered during the sub-systems matrix assembly stage. Thus, the application of this approach to over-constrained PKMs requires the deter- mination of Jacobian matrices to characterize kinematic dependencies due to closed-loop over-constrained mechanisms [16,21] or the definition of a complex methodology to merge stiffness matrices of each PKM element [19]. Those Jacobian matrices describe the kine- matic behavior of PKMs by derivation of the geometric model as in Germain's [16] and Zhang's [17] works. Geometric model or closed-loop constraints of complex PKM, such as Exechon robot, can sometimes not be computed in a systematic and straightforward way [16,22]. For example, the expression of the geometric model of an Exechon robot is based on a system of nonlinear equations that are not easy to establish [23]. Moreover, its solving requires a numerical optimization with the Newton-Raphson method.

This paper aims to introduce a reliable and simple method to compute global mass and stiffness matrices for over-constrained PKMs, under the assumption of small displacements and using screw theory. This methodology is relevant for high-stiffness architectures such as machine tools and industrial robots where flexibility generates only small displacements.
These mathematical tools allow defining geometric dependencies of leg movements without
the computation of Jacobian matrices [24]. Screw theory is used both:

1. To express the kinematic behavior of each PKM leg using the theory of Timoshenko beams [25] and

2. To determine the displacement constraints applied to each leg extremity due to over-constrained mechanism closed-loops [26].

There are different ways to parameterize the orientation of a rigid body such as quaternion [27]. The choice of screw theory ensures use of the same representation for the actuation and constraint wrenches applied to the moving platform by the actuators and the geometry of the mechanism.

This second point is the main contribution of our work and allows us to take into
account stress due to over-constrained with geometric constraint equations extracted only
from the screw theory application to the passive joints. With this method, the computation
of global mass and stiffness matrices is easy to implement. It assumes only geometrical
constraints between displacement vector components without the introduction of the
Jacobian matrices. Global mass and stiffness matrices are then computed to estimate PKM
Cartesian stiffness, the first natural frequencies and their associated modes. This work is
relevant at the embodiment design stage of the robot when only a parameterized CAD
modeling of the PKM under study is available. Indeed, at this stage, a fast computation
time is required to assess the capability of the PKM in terms of vibration behavior and to
determine the optimal architecture and dimensions of the robot under design. During the
next design stages, when accurate simulation results are needed, FEA can be preferred
although it is more computation intensive.

It should be noted that the revolute joints of the Tripteor legs are stiff; their stiffness is larger than 1.04 × 10⁹ N·m⁻¹ [28]. Moreover, the stiffness of the spherical wrist of the Tripteor is even more important. As a consequence, and for the sake of clarity, only the deflection of Tripteor X7 PKM legs is considered in what remains. The joint deflection can be taken into account by adding local stiffness such as with the Virtual Joint Method [29].

The paper is organized as follows: Section 2 presents the Tripteor X7 PKM and its parametrization. Section 3 describes the proposed method to express the elasto-dynamic model of PKMs and to compute their natural frequencies. As an illustrative example, the stiffness matrix, the natural frequencies and associated modes of the Tripteor X7 PKM are computed in Section 4. Finally, conclusions and future work are drawn in Section 5.

## 2. Tripteor X7 PKM

Tripteor X7 PKM is manufactured by the PCI-SCEMM company in France. It is a
hybrid PKM. Its architecture is a combination of a parallel mechanism and a serial wrist
mounted in series. Such a hybrid manipulator ensures the decoupling of translational and
rotational motion of the end-effector. Moreover, hybrid manipulators usually have a better
dexterity than fully parallel manipulators as explained in [30].

Tripteor X7 PKM is a six-axis PKM with the following actuated joint variables (q₁, q₂, q₃, q₄, q₅) shown in Figure 2 and one-DoF rotary table represented in Figure 1. The revolute joint variable of this rotary table is denoted as q₆. A parallel mechanism provides three Degrees of Freedom (DoFs) from the actuation of three ball screw systems (q₁, q₂ and q₃) and a serial wrist with two DoFs from the actuation of two revolute joints with direct and belt drives (q₄ and q₅). The elongations q₁, q₂ and q₃ of Leg 1, Leg 2 and Leg 3 provide translational motions of the moving platform along axis x<sub>b</sub>, y<sub>b</sub> and axis z<sub>b</sub> with induced rotational motions about axis x<sub>b</sub> and axis y<sub>m</sub>. The serial wrist ensures rotational motions of the end-effector.

The paper aims to introduce a new methodology to compute natural frequencies of
over-constrained closed-loop mechanism in the case of PKM. Thus, only the vibration
behavior of the parallel mechanism composed of Leg 1, Leg 2, Leg 3 and the mobile
platform (4) is studied in this paper (Figure 2). If a study of the complete PKM should be
carried, a simple assembly of global mass and stiffness matrices of the parallel mechanism
and the serial wrist can then be carried out.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Kinematic model of Tripteor X7 PKM.</figcaption></figure>

Each leg is connected to the mobile platform (4) by a revolute joint and to the fixed platform (0) by two revolute joints which form a universal joint (Figure 2). It is worth noting that Leg 2 has a supplementary revolute joint around its principal axis. Consequently, bending and traction-compression loads are prescribed to Leg 2 whereas Legs 1 and 3 are loaded in bending, traction-compression, and torsion [28]. The mobile platform is supposed to be rigid. Because stiffness and low frequencies are of primal importance in machining [11], an efficient and fast way to compute those frequencies is described in this paper. For this purpose, a local model of a PKM leg under bending, traction-compression, and torsion loads is introduced before the definition of constraint equations due to mechanism assembly.

**3. Elasto-Dynamic Modeling of PKMs**

This section introduces the novel elasto-dynamic modeling approach and its applica-
tion to PKMs. For the sake of pedagogy, we introduce our method for PKMs with telescopic
legs such as the Tripteor X7 (Figure 3). Aᵢ and Bᵢ are the leg extremities. Aᵢ is attached to
the PKM base and Bᵢ is fixed to the mobile platform. The coordinate system associated
with the fixed base is denoted as $\mathcal{R}_b = (O, \mathbf{x}_b, \mathbf{y}_b, \mathbf{z}_b)$ and with the mobile platform as
$\mathcal{R}_m = (P, \mathbf{x}_m, \mathbf{y}_m, \mathbf{z}_m)$. zᵢ is the unit vector along leg i direction.

This method can be applied to any PKMs whose legs can be modeled by beams.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Definition of used parameters for PKM with telescopic legs.</figcaption></figure>
### 3.1. Local Modeling of PKM Leg

First, each leg is modeled assuming the usual Euler-Bernoulli beam theory. It is worth noting that this model considers the usual mechanical loads, which are bending, traction-compression, and torsion conditions. Each leg is considered to be being connected at one of its ends via a given joint to a fixed part (Figure 4). This boundary condition depends on the joint type used for connecting the leg to the fixed base. In this sub-section, we first study the modeling and the parameter definition of a leg assuming generic boundary conditions to elaborate a leg stiffness matrix $K_e$ and a leg mass matrix $M_e$ in the leg coordinate system. With a generic beam model, leg kinematics are defined through six displacement functions, which correspond to the displacement of the neutral axis. $\mathbf{u}(x)$ is the vector that collects those displacement functions:

$$ \mathbf{u}(x)^T = [u(x) \ v(x) \ w(x) \ \theta_{\alpha}(x) \ \theta_{\beta}(x) \ \theta_{\gamma}(x)]^T \quad (1) $$

$u(x)$ is the displacement along the neutral axis, $v(x)$ and $w(x)$ are the displacements along **fi**-direction and **fl**-direction. $\theta_{\alpha}(x)$, $\theta_{\beta}(x)$, and $\theta_{\gamma}(x)$ correspond to the small rotations about **ff**-axis, **fi**-axis and **fl**-axis, respectively. Using a Bernoulli beam theory, cross-sectional displacements and rotations satisfy the following equations:

$$ \theta_{\gamma}(x) = \frac{dv}{dx}(x), \quad \theta_{\beta}(x) = -\frac{dw}{dx}(x) \quad (2) $$

<figure><img src="image_8.png" /><figcaption>Figure 4. Local leg model and displacement parameters.</figcaption></figure>

Each kinematic function corresponds to a force function. In the following, $\mathbf{f}_{\text{beam}}$ denotes the force vector and $\mathbf{m}_{\text{beam}}$ is the moment vector such that:

$$ \left\{ \begin{array}{l} \mathbf{f}_{\text{beam}}(x) = N(x)\mathbf{ff} + T_{\beta}(x)\mathbf{fi} + T_{\gamma}(x)\mathbf{fl} \\ \mathbf{m}_{\text{beam}}(x) = M_{t}(x)\mathbf{ff} + M_{\beta}(x)\mathbf{fi} + M_{\gamma}(x)\mathbf{fl} \end{array} \right. \quad (3) $$

where $N(x)$ is the axial force, $T_{\beta}(x)$ and $T_{\gamma}(x)$ are the shear forces, $M_{t}(x)$ is the torsion moment, and $M_{\beta}(x)$ and $M_{\gamma}(x)$ are the bending moments. As there is no distributed load, the governing equilibrium equations give:

$$ \left\{ \begin{array}{l} \frac{dN}{dx}(x) = 0, \quad \frac{dT_{\beta}}{dx}(x) = 0, \quad \frac{dT_{\gamma}}{dx}(x) = 0 \\ \frac{dM_{t}}{dx}(x) = 0, \quad \frac{dM_{\beta}}{dx}(x) - T_{\gamma}(x) = 0, \quad \frac{dM_{\gamma}}{dx}(x) + T_{\beta}(x) = 0 \end{array} \right. \quad (4) $$

In addition, the kinematic and force functions are related to each other as follows:

$$ N(x) = ES\frac{du}{dx}(x), \quad M_{t}(x) = GI_{G}\frac{d\theta_{\alpha}}{dx}(x), \quad M_{\beta}(x) = EI_{\beta}\frac{d\theta_{\beta}}{dx}(x), \quad M_{\gamma}(x) = EI_{\gamma}\frac{d\theta_{\gamma}}{dx}(x) \quad (5) $$
where E is the material Young's modulus, S is the cross-sectional area, and $I_G$, $I_β$ and $I_γ$ are moments of area about **ff-axis**, **fi-axis** and **fl-axis**.

Combining Equations (2), (4) and (5), we can conclude that $u(x)$ and $\theta_α(x)$ are linear functions, whereas $v(x)$ and $w(x)$ are polynomial functions of order 3. In a generic case, all kinematic functions can thus be defined through their nodal values expressed in Equation (6).

$$ \mathbf{u}(x) = \begin{bmatrix} u(0) + \frac{x}{q}(u(q) - u(0)) \\ v(0) + \theta_{\gamma}(0)x - \left(\frac{2\theta_{\gamma}(0) + \theta_{\gamma}(q)}{q} - \frac{3(v(q) - v(0))}{q^2}\right)x^2 + \left(\frac{\theta_{\gamma}(0) + \theta_{\gamma}(q)}{q^2} - \frac{2(v(q) - v(0))}{q^3}\right)x^3 \\ w(0) - \theta_{\beta}(0)x + \left(\frac{2\theta_{\beta}(0) + \theta_{\beta}(q)}{q} + \frac{3(w(q) - w(0))}{q^2}\right)x^2 - \left(\frac{\theta_{\beta}(0) + \theta_{\beta}(q)}{q^2} + \frac{2(w(q) - w(0))}{q^3}\right)x^3 \\ \theta_{\alpha}(0) + \frac{x}{q}(\theta_{\alpha}(q) - \theta_{\alpha}(0)) \\ \theta_{\beta}(0) - 2\left(\frac{2\theta_{\beta}(0) + \theta_{\beta}(q)}{q} + \frac{3(w(q) - w(0))}{q^2}\right)x + 3\left(\frac{\theta_{\beta}(0) + \theta_{\beta}(q)}{q^2} + \frac{2(w(q) - w(0))}{q^3}\right)x^2 \\ \theta_{\gamma}(0) - 2\left(\frac{2\theta_{\gamma}(0) + \theta_{\gamma}(q)}{q} - \frac{3(v(q) - v(0))}{q^2}\right)x + 3\left(\frac{\theta_{\gamma}(0) + \theta_{\gamma}(q)}{q^2} - \frac{2(v(q) - v(0))}{q^3}\right)x^2 \end{bmatrix} \quad (6) $$

$\mathbf{u}_{nod}$ denotes the vector collecting the nodal values of the kinematic functions (Figure 4) and is expressed as a $1 \times 12$ column vector:

$$ \mathbf{u}_{nod} = \begin{bmatrix} \mathbf{u}(0) \\ \mathbf{u}(q) \end{bmatrix} \quad (7) $$

Vectors $\mathbf{N}_u(x)$, $\mathbf{N}_v(x)$, $\mathbf{N}_w(x)$, $\mathbf{N}_{\theta_α}(x)$, $\mathbf{N}_{\theta_β}(x)$ and $\mathbf{N}_{\theta_γ}(x)$ express the nodal values as a function of $x$ coordinate bounded between 0 and $q$, i.e., $x \in [0, q]$:

$$ \left\{ \begin{aligned} u(x) &= \mathbf{N}_u(x)^T \mathbf{u}_{nod} \\ v(x) &= \mathbf{N}_v(x)^T \mathbf{u}_{nod} \\ w(x) &= \mathbf{N}_w(x)^T \mathbf{u}_{nod} \\ \theta_\alpha(x) &= \mathbf{N}_{\theta_\alpha}(x)^T \mathbf{u}_{nod} \\ \theta_\beta(x) &= \mathbf{N}_{\theta_\beta}(x)^T \mathbf{u}_{nod} \\ \theta_\gamma(x) &= \mathbf{N}_{\theta_\gamma}(x)^T \mathbf{u}_{nod} \end{aligned} \right. \quad (8) $$

From [31], the beam potential energy $E_p$ is written as:

$$ E_p = \frac{1}{2} \int_0^q \left[ E_S \left( \frac{du(x)}{dx} \right)^2 + G I_G \left( \frac{d\theta_\alpha(x)}{dx} \right)^2 + E I_\beta \left( \frac{d\theta_\beta(x)}{dx} \right)^2 + E I_\gamma \left( \frac{d\theta_\gamma(x)}{dx} \right)^2 \right] dx \quad (9) $$

From Equations (8) and (9), $E_p$ can be expressed in a compact form as a function of the beam stiffness matrix $\mathbf{K}_e$ and $\mathbf{u}_{nod}$ as follows:

$$ E_p = \frac{1}{2} \mathbf{u}_{nod}^T \mathbf{K}_e \mathbf{u}_{nod} \quad (10) $$

With matrix $\mathbf{K}_e$ taking the form of Equation (11).
$$ K_e = \begin{bmatrix} \frac{ES}{q} & 0 & 0 & 0 & 0 & -\frac{ES}{q} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{12EI_\gamma}{q^3} & 0 & 0 & 0 & \frac{6EI_\gamma}{q^2} & 0 & -\frac{12EI_\gamma}{q^3} & 0 & 0 & 0 & \frac{6EI_\gamma}{q^2} \\ 0 & 0 & \frac{12EI_\beta}{q^3} & 0 & -\frac{6EI_\beta}{q^2} & 0 & 0 & 0 & -\frac{12EI_\beta}{q^3} & 0 & -\frac{6EI_\beta}{q^2} & 0 \\ 0 & 0 & 0 & \frac{GI_G}{q} & 0 & 0 & 0 & 0 & -\frac{GI_G}{q} & 0 & 0 & 0 \\ 0 & 0 & -\frac{6EI_\beta}{q^2} & 0 & \frac{4EI_\beta}{q} & 0 & 0 & 0 & \frac{6EI_\beta}{q^2} & 0 & \frac{2EI_\beta}{q} & 0 \\ 0 & \frac{6EI_\gamma}{q^2} & 0 & 0 & 0 & \frac{4EI_\gamma}{q} & 0 & -\frac{6EI_\gamma}{q^2} & 0 & 0 & 0 & \frac{2EI_\gamma}{q} \\ -\frac{ES}{q} & 0 & 0 & 0 & 0 & \frac{ES}{q} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{12EI_\gamma}{q^3} & 0 & 0 & 0 & -\frac{6EI_\gamma}{q^2} & 0 & \frac{12EI_\gamma}{q^3} & 0 & 0 & 0 & -\frac{6EI_\gamma}{q^2} \\ 0 & 0 & -\frac{12EI_\beta}{q^3} & 0 & \frac{6EI_\beta}{q^2} & 0 & 0 & 0 & \frac{12EI_\beta}{q^3} & 0 & \frac{6EI_\beta}{q^2} & 0 \\ 0 & 0 & 0 & -\frac{GI_G}{q} & 0 & 0 & 0 & 0 & \frac{GI_G}{q} & 0 & 0 & 0 \\ 0 & 0 & -\frac{6EI_\beta}{q^2} & 0 & \frac{2EI_\beta}{q} & 0 & 0 & 0 & \frac{6EI_\beta}{q^2} & 0 & \frac{4EI_\beta}{q} & 0 \\ 0 & \frac{6EI_\gamma}{q^2} & 0 & 0 & 0 & \frac{2EI_\gamma}{q} & 0 & -\frac{6EI_\gamma}{q^2} & 0 & 0 & 0 & \frac{4EI_\gamma}{q} \end{bmatrix} \quad (11) $$

From [31], the beam kinetic energy $E_k$ is defined as:

$$ E_k = \frac{1}{2} \rho \int_0^q \iint_S \mathbf{v}_M^T \mathbf{v}_M dxdydz \quad (12) $$

where $\rho$ is the beam material density, and $\mathbf{v}_M$ is the point velocity along the beam expressed as:

$$ \mathbf{v}_M = (\dot{u}(x) + z\dot{\theta}_\beta(x) - y\dot{\theta}_\gamma(x))\mathbf{ff} + (\dot{v}(x) - z\dot{\theta}_\alpha(x))\mathbf{fi} + (\dot{w}(x) + y\dot{\theta}_\alpha(x))\mathbf{fl} \quad (13) $$

Combining Equations (12) and (13) with the hypothesis that **fi-axis** and **fl-axis** are principal axes of inertia, the kinematic energy $E_k$ writes:

$$ E_k = \frac{1}{2}\rho S \int_0^q (\dot{u}(x)^2 + \dot{v}(x)^2 + \dot{w}(x)^2) dx + \frac{1}{2}\rho I_{\beta\beta} \int_0^q (\dot{\theta}_\alpha(x)^2 + \dot{\theta}_\beta(x)^2) dx \\ + \frac{1}{2}\rho I_{\gamma\gamma} \int_0^q (\dot{\theta}_\alpha(x)^2 + \dot{\theta}_\gamma(x)^2) dx \quad (14) $$

with $I_{\beta\beta}$ and $I_{\gamma\gamma}$ the beam inertia about **fi-axis** and **fl-axis**. Equation (14) takes the following compact form:

$$ E_k = \frac{1}{2} \dot{\mathbf{u}}_{\mathrm{nod}}^T \mathbf{M}_e \dot{\mathbf{u}}_{\mathrm{nod}} \quad (15) $$

with $\dot{\mathbf{u}}_{\mathrm{nod}}$ the time derivative of vector $\mathbf{u}_{\mathrm{nod}}$ and $\mathbf{M}_e$ the mass matrix of beam expressed in Equation (16).

$$ \mathbf{M}_e = \begin{bmatrix} \frac{b}{3} & 0 & 0 & 0 & 0 & \frac{b}{6} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{13}{35}b + \frac{6\rho I_{\gamma\gamma}}{5q} & 0 & 0 & 0 & \frac{\rho I_{\gamma\gamma}}{10} + \frac{11c}{210} & 0 & \frac{9b}{70} - \frac{6\rho I_{\gamma\gamma}}{5q} & 0 & 0 & \frac{\rho I_{\gamma\gamma}}{10} - \frac{13c}{420} \\ 0 & 0 & \frac{6\rho I_{\beta\beta}}{-5q} + \frac{13}{35}b & 0 & -\frac{\rho I_{\beta\beta}}{10} - \frac{11c}{210} & 0 & 0 & 0 & \frac{9b}{70} - \frac{6\rho I_{\beta\beta}}{-5q} & 0 & -\frac{\rho I_{\beta\beta}}{10} + \frac{13c}{420} & 0 \\ 0 & 0 & 0 & a & 0 & 0 & 0 & 0 & \frac{a}{2} & 0 & 0 & 0 \\ 0 & 0 & -\frac{\rho I_{\beta\beta}}{-10} - \frac{11c}{210} & 0 & \frac{2\rho I_{\beta\beta}q}{15} + \frac{d}{105} & 0 & 0 & 0 & \frac{\rho I_{\beta\beta}}{10} - \frac{13c}{420} & 0 & -\frac{\rho I_{\beta\beta}q}{30} - \frac{d}{140} & 0 \\ 0 & \frac{\rho I_{\gamma\gamma}}{10} + \frac{11c}{210} & 0 & 0 & 0 & \frac{2\rho I_{\gamma\gamma}q}{15} + \frac{d}{105} & 0 & -\frac{\rho I_{\gamma\gamma}}{10} + \frac{13c}{420} & 0 & 0 & 0 & -\frac{\rho I_{\gamma\gamma}q}{30} - \frac{d}{140} \\ \frac{b}{6} & 0 & 0 & 0 & 0 & 0 & \frac{b}{3} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{9b}{70} - \frac{6\rho I_{\gamma\gamma}}{5q} & 0 & 0 & 0 & -\frac{\rho I_{\gamma\gamma}}{10} + \frac{13c}{420} & 0 & \frac{13}{35}b + \frac{6\rho I_{\gamma\gamma}}{5q} & 0 & 0 & 0 & -\frac{\rho I_{\gamma\gamma}}{10} - \frac{11c}{210} \\ 0 & 0 & \frac{9b}{70} - \frac{6\rho I_{\beta\beta}}{-5q} & 0 & \frac{\rho I_{\beta\beta}}{10} - \frac{13c}{420} & 0 & 0 & 0 & \frac{6\rho I_{\beta\beta}}{-5q} + \frac{13}{35}b & 0 & \frac{\rho I_{\beta\beta}}{10} + \frac{11c}{210} & 0 \\ 0 & 0 & 0 & \frac{a}{2} & 0 & 0 & 0 & 0 & a & 0 & 0 & 0 \\ 0 & 0 & -\frac{\rho I_{\beta\beta}}{-10} + \frac{13c}{420} & 0 & -\frac{\rho I_{\beta\beta}q}{30} - \frac{d}{140} & 0 & 0 & 0 & \frac{\rho I_{\beta\beta}}{10} + \frac{11c}{210} & 0 & \frac{2\rho I_{\beta\beta}q}{15} + \frac{d}{105} & 0 \\ 0 & \frac{\rho I_{\gamma\gamma}}{10} - \frac{13c}{420} & 0 & 0 & 0 & -\frac{\rho I_{\gamma\gamma}q}{30} - \frac{d}{140} & 0 & -\frac{\rho I_{\gamma\gamma}q}{10} - \frac{11c}{210} & 0 & 0 & 0 & \frac{2\rho I_{\gamma\gamma}q}{15} + \frac{d}{105} \end{bmatrix} \quad (16) $$

$$ \text{with } a = \frac{\rho(I_{\beta\beta} + I_{\gamma\gamma})q}{3}, \quad b = \rho Sq, \quad c = \rho Sq^2 \text{ and } d = \rho Sq^3. $$
To compute the global mass and stiffness matrices, the constraints between the dis- placements of the leg extremities are determined in the next section.

## 3.2. Constraint Equations

The mobile platform is considered to be rigid for the elaboration of the constraint equations between the leg extremities. We therefore express the movement of the mobile platform due to leg flexibility in stationary configurations of the PKM architecture using screw theory [26] (Figure 5).

<figure><img src="image_9.png" /><figcaption>Figure 5. Parametrization of mobile platform small displacement.</figcaption></figure>

The small displacement screw $\mathcal{T}_m$ of the mobile platform is influenced by displacement and rotation $\mathbf{u}_i(q_i)$ of the leg $i$ at point $B_i$ regarding the movements of passive joints between leg $i$ and mobile platform. A normalized screw is used to model passive joint movement [32]. Each movement between leg $i$ and mobile platform is model with a one-DoF motion-screw:

$$ S_{ik} = \begin{bmatrix} S_{Fik} \\ S_{Sik} \end{bmatrix} = \begin{cases} \begin{bmatrix} \mathbf{s}_{ik} \\ \mathbf{s}_{ik} \times \mathbf{r}_{ik} \end{bmatrix} & \text{if it is a rotational movement} \\ \begin{bmatrix} \mathbf{0} \\ \mathbf{s}_{ik} \end{bmatrix} & \text{if it is a prismatic movement} \end{cases} \qquad (17) $$

where $k = 1$ to $n_i$ ($n_i$ is the number of equivalent one-DoF joints between leg $i$ and mobile platform, $n_i = 1$ for the case study), $\mathbf{s}_{ik}$ is a unit vector along the axis of the screw $S_{ik}$, $\mathbf{r}_{ik}$ is a vector directed from any point on the screw axis to point $B_i$.

The computation of the movement transmitted by the leg to the mobile platform or vice versa can be done using the reciprocal screws of the motion-screws $S_{ik}$. There is a unique normalized reciprocal screw system $\mathcal{T}_{ik}^{\perp}$ of order 5 of $S_{ik}$ [32]. For a rotational joint, the reciprocal system is a 5-system which includes all rotational screws whose axes intersect the joint axis and all prismatic screws whose axes are perpendicular to joint axes and all the combinations of the above reciprocal screw [32]. The relation between screw $S_{ik}$ and the five one-DoF reciprocal screws $S_{ijk}^{\perp}$ is:

$$ \left( \begin{bmatrix} \mathbf{0} & \mathbf{I}_3 \\ \mathbf{I}_3 & \mathbf{0} \end{bmatrix} S_{ik} \right)^T S_{ijk}^{\perp} = 0 \qquad (18) $$

where $\mathbf{I}_3$ is the 3 × 3 identity matrix, $\mathbf{0}$ is the 3 × 3 zero matrix, $j = 1$ to 5. Then, $\mathcal{T}_{ik}^{\perp} = (S_{ijk}^{\perp}$ with $j = 1$ to 5).

The small displacement screw is defined as $\mathcal{T}_m = \begin{bmatrix} \mathbf{!} \\ \mathbf{d}_{B_i} \end{bmatrix}$ with $\mathbf{!}$ the rotation vector of the mobile platform and $\mathbf{d}_{B_i}$ the small displacement vector of point $B_i$ ($\mathbf{d}_{B_j} = \mathbf{d}_{B_i} + \mathbf{b}_{ij} \times \mathbf{!}$ with $\mathbf{b}_{ji}$ the vector from $B_j$ to $B_i$). For leg $i$, the relation between $\mathcal{T}_m$ and $\mathbf{u}_i(q_i)$ is written:

$$ \prod_{k=1}^{n_i} \left( \sum_{j=1}^{5} S_{ijk}^{\perp} \right) \left( \mathbf{u}_i(q_i) - \begin{bmatrix} \mathbf{R}_{\text{leg}_i \to \text{beam}} \mathbf{R}_{\text{R}_m \to \text{leg}_i} \mathbf{d}_{B_i} \\ \mathbf{R}_{\text{leg}_i \to \text{beam}} \mathbf{R}_{\text{R}_m \to \text{leg}_i} \mathbf{!} \end{bmatrix} \right) = \mathbf{0} \qquad (19) $$
where $\mathbf{R}_{leg_i \to beam}$ is the rotational matrix between the beam coordinate system and the leg i coordinate system, and $\mathbf{R}_{R_m \to leg_i}$ is the rotational matrix between the leg i coordinate system and the mobile platform coordinate system.

Constraint equations are defined by expressing Equation (19) for each leg and using small displacement screw properties.

In the same way, the nodal displacement $\mathbf{u}_i(0)$ is limited due to passive joints between the leg i and the fixed base. We note $S'_{ijk}$ the movement-screws between leg i and fixed base. Thus, the null nodal parameter is computed from:

$$ \left( \prod_{k=1}^{n'_i} \left( \sum_{j=1}^{5} S'_{ijk} \perp \right) \right) \mathbf{u}_i(0) = \mathbf{0} \qquad (20) $$

where $n'_i$ is the number of equivalent one DoF joints between leg i and fixed base. For the case study, $n'_1 = n'_3 = 2$, $n'_2 = 3$.

From Equations (19) and (20), constraint equations are grouped together as:

$$ \mathbf{L}_D \mathbf{u}_D + \mathbf{L}_I \mathbf{u}_I = \mathbf{0} \Rightarrow \mathbf{u}_D = -\mathbf{L}_D^{-1} \mathbf{L}_I \mathbf{u}_I \qquad (21) $$

where $\mathbf{u}_D$ is a vector composed of the set of beam nodal dependent parameters extracted from all $\mathbf{u}_i(q_i)$. Its size is equal to the rank of the equation system (19) for all legs. $\mathbf{u}_I$ is a vector composed of the set of beam nodal independent parameters, which is not included in $\mathbf{u}_D$. The choice of $\mathbf{u}_D$ and $\mathbf{u}_I$ is not unique. The components of matrices $\mathbf{L}_D$ and $\mathbf{L}_I$ are the factor of vectors $\mathbf{u}_D$ and $\mathbf{u}_I$ components in Equations (19) and (20).

Equation (21) corresponds to a simple formulation of the closed-loop geometric constraints for over-constrained PKM and is the major contribution of the paper. The determination of stiffness and vibratory behavior of a PKM is based on the definition of its global mass and stiffness matrices. These matrices are defined in the following section.

### 3.3. Mass Matrix Computation of Mobile Platform

The mobile platform and the serial wrist are considered to be a point-mass m located at its center of mass G. Thus, kinematic energy $E_{km}$ writes:

$$ E_{km} = \frac{1}{2} m \mathbf{v}_G^T \mathbf{v}_G \qquad (22) $$

$\mathbf{v}_G$ being the linear velocity of point G. Under the assumption of small displacement $\mathbf{v}_G$ is expressed as:

$$ \mathbf{v}_G = \dot{\mathbf{d}}_G = \dot{\mathbf{d}}_{B_1} + \mathbf{g}_1 \times \mathbf{\dot{i}} \qquad (23) $$

where $\mathbf{d}_G$ is the small displacement of point G, $\dot{\mathbf{d}}_G$ is its time derivative and $\mathbf{g}_1$ is the vector from G to $B_1$.

Equation (23) allows defining $E_{km}$ according to beam nodal parameters:

$$ E_{km} = \frac{1}{2} \left( \begin{pmatrix} \dot{\mathbf{u}}_D \\ \dot{\mathbf{u}}_I \end{pmatrix} \right)^T \mathbf{M}_{em} \left( \begin{pmatrix} \dot{\mathbf{u}}_D \\ \dot{\mathbf{u}}_I \end{pmatrix} \right) \qquad (24) $$

where $\mathbf{M}_{em}$ is the mobile platform mass matrix. To compute natural frequency of the system, this mass matrix is added to the legs mass matrix.

### 3.4. Cartesian Stiffness Computation

$\mathbf{K}$ is the global stiffness matrix of the PKM. It is obtained from the assembly of stiffness matrices $\mathbf{K}_e$ of the three legs such that:

$$ \mathbf{K} = \begin{bmatrix} \mathbf{K}_D & \mathbf{K}_{DI} \\ \mathbf{K}_{DI}^T & \mathbf{K}_I \end{bmatrix} \qquad (25) $$
where $\mathbf{K}_D$ is the stiffness matrix corresponding to dependent beam nodal parameters $\mathbf{u}_D$, and $\mathbf{K}_I$ corresponds to the independent ones $\mathbf{u}_I$. $\mathbf{K}_{DI}$ is the coupling matrix.

From Equation (21), a decomposition between dependent and independent beam nodal parameters is used to express the potential energy $E_p$:

$$E_p = \frac{1}{2} \begin{bmatrix} \mathbf{u}_D \\ \mathbf{u}_I \end{bmatrix}^T \begin{bmatrix} \mathbf{K}_D & \mathbf{K}_{DI} \\ \mathbf{K}_{DI}^T & \mathbf{K}_I \end{bmatrix} \begin{bmatrix} \mathbf{u}_D \\ \mathbf{u}_I \end{bmatrix} = \frac{1}{2} \mathbf{u}_I^T \tilde{\mathbf{K}} \mathbf{u}_I \quad (26)$$

Consequently, this allows the introduction of the condensed matrix form $\tilde{\mathbf{K}}$ of the global stiffness matrix of the studied PKM, which now only depends on the independent beam nodal parameters, such that:

$$\tilde{\mathbf{K}} = \mathbf{L}_I^T (\mathbf{L}_D^{-1})^T \mathbf{K}_D \mathbf{L}_D^{-1} \mathbf{L}_I - (\mathbf{K}_{DI}^T \mathbf{L}_D^{-1} \mathbf{L}_I + (\mathbf{K}_{DI}^T \mathbf{L}_D^{-1} \mathbf{L}_I)^T) + \mathbf{K}_I \quad (27)$$

The Cartesian stiffness matrix $\mathbf{K}_c$ can be introduced from the expression of the potential energy $E_p$ according to:

$$E_p = \frac{1}{2} \mathbf{u}_I^T \tilde{\mathbf{K}} \mathbf{u}_I = \frac{1}{2} \mathbf{ffix}^T \mathbf{K}_c \mathbf{ffix} \quad (28)$$

where $\mathbf{ffix}$ is the small displacement vector of point $O_m$ expressed in the base frame. Thus, $\mathbf{ffix} = \mathbf{R}_m^{-1}(\mathbf{d}_{B_i} + \mathbf{b}_{mi} \times !) = \mathbf{L}_{XI}\mathbf{u}_I$ where $\mathbf{R}_m$ is the rotation matrix from the base frame to the mobile platform frame and $\mathbf{b}_{mi}$ is the vector from $O_m$ to $B_i$. Finally, the Cartesian stiffness matrix of the PKM is computed from $\tilde{\mathbf{K}}$:

$$\mathbf{K}_c = (\mathbf{L}_{XI} \mathbf{L}_{XI}^T)^{-1} \mathbf{L}_{XI} \tilde{\mathbf{K}} \mathbf{L}_{XI}^T (\mathbf{L}_{XI} \mathbf{L}_{XI}^T)^{-1} \quad (29)$$

Note that $\mathbf{K}_c$ is a 3 × 3 matrix and its components are the stiffness along the PKM base frame vectors.

### 3.5. Natural Frequency Computation

The global stiffness matrix $\mathbf{M}$ of the PKM is obtained from the assembly of the mass matrices $\mathbf{M}_e$ and $\mathbf{M}_{em}$ of the three legs and mobile platform, respectively. It is expressed as:

$$\mathbf{M} = \begin{bmatrix} \mathbf{M}_D & \mathbf{M}_{DI} \\ \mathbf{M}_{DI}^T & \mathbf{M}_I \end{bmatrix} \quad (30)$$

where $\mathbf{M}_D$ is the mass matrix corresponding to dependent beam nodal parameters $\mathbf{u}_D$, whereas $\mathbf{M}_I$ corresponds to the independent ones $\mathbf{u}_I$. $\mathbf{M}_{DI}$ is the coupling matrices.

The decomposition between dependent and independent parameters is used to express the kinetic energy $E_k$ as follows:

$$E_k = \frac{1}{2} \begin{bmatrix} \mathbf{u}_D \\ \mathbf{u}_I \end{bmatrix}^T \begin{bmatrix} \mathbf{M}_D & \mathbf{M}_{DI} \\ \mathbf{M}_{DI}^T & \mathbf{M}_I \end{bmatrix} \begin{bmatrix} \mathbf{u}_D \\ \mathbf{u}_I \end{bmatrix} = \frac{1}{2} \mathbf{u}_I^T \tilde{\mathbf{M}} \mathbf{u}_I \quad (31)$$

Consequently, this allows the introduction of the condensed matrix form of the global mass matrix $\tilde{\mathbf{M}}$, which now only depends on the independent beam nodal parameters:

$$\tilde{\mathbf{M}} = \mathbf{L}_I^T (\mathbf{L}_D^{-1})^T \mathbf{M}_D \mathbf{L}_D^{-1} \mathbf{L}_I - (\mathbf{M}_{DI}^T \mathbf{L}_D^{-1} \mathbf{L}_I + (\mathbf{M}_{DI}^T \mathbf{L}_D^{-1} \mathbf{L}_I)^T) + \mathbf{M}_I \quad (32)$$

The eigenvalues $\lambda_{ev}$ and the associated eigenvector $\mathbf{u}_{ev}$ are determined from the spectral decomposition of matrix $\tilde{\mathbf{M}}^{-1}\tilde{\mathbf{K}}$. The eigenvalues are the solutions of the polynomial $\det(\tilde{\mathbf{M}}^{-1}\tilde{\mathbf{K}} - \lambda_{ev}\mathbf{I}) = 0$. The ith natural frequency $f_{0i}$ is expressed as a function of the ith
eignevalue $\lambda_{evi}$ as follows: $f_{0i} = \frac{\sqrt{\lambda_{evi}}}{2\pi}$. The natural modes are the eigenvectors associated with those natural eigenvalues.

## 4. Natural Frequencies and Modes of Tripteor X7

The methodology described in Section 3 is applied in this section to calculate the natural frequencies and associated modes of the Tripteor X7 PKM.

### 4.1. Identification of Null Legs Nodal Parameters

The null nodal parameters of $u_i(0)$ depend on the joint types between the legs and the base as shown in Figure 6. The first joint of Legs 1 and 3 with the base is a revolute joint about fl-axis, and the second joint is about fi-axis. Thus, for Legs 1 and 3, infinitesimal and reciprocal screw motion are the following:

$$S'_{i1} = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{pmatrix} \quad \text{and} \quad S'_{i2} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \quad \Rightarrow \quad \sum_{j=1}^{5} S'_{ij1} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 0 \end{pmatrix} \quad \text{and} \quad \sum_{j=1}^{5} S'_{ij2} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 0 \\ 1 \end{pmatrix} \quad (33)$$

<figure><img src="image_10.png" /><figcaption>Figure 6. Local parameters for Tripteor X7 legs.</figcaption></figure>

Thus, from Equation (20), only $\theta_{\beta i}(0)$ and $\theta_{\gamma i}(0)$ are not null.

For Leg 2, with the same methodology, we obtain that only $\theta_{\alpha2}(0)$, $\theta_{\beta2}(0)$ and $\theta_{\gamma2}(0)$ are not null (Figure 6b).

To compute the global mass and stiffness matrices, the motion constraints between the extremities of the three legs are explained hereafter.

### 4.2. Constraint Equations

The joints between Legs 1 and 3 and the mobile platform are revolute joints about ym-axis (fi-axis) and the joint between Legs 2 and the mobile platform is a revolute joint about xm-axis (fl-axis) (Figure 7). Thus, the reciprocal motion-screws in the leg coordinate system are:

$$\sum_{j=1}^{5} S_{1j1}^{\perp} = \sum_{j=1}^{5} S_{3j1}^{\perp} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 0 \\ 1 \end{pmatrix} \quad \text{and} \quad \sum_{j=1}^{5} S_{2j1}^{\perp} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 0 \end{pmatrix} \quad (34)$$
<figure><img src="image_11.png" /><figcaption>Figure 7. Parametrization of Tripteor mobile platform movement.</figcaption></figure>

From Equation (19), the equation associated with the motion constraints between the mobile platform and the three legs of Tripteor X7 are the following:

$$ \sum_{j=1}^{5} S_{1j1}^{\perp} \left[ \mathbf{u}_1(q_1) - \begin{bmatrix} \mathbf{R}_{\text{leg1}} \rightarrow \text{beam} \mathbf{R}_{\text{Rm}} \rightarrow \text{leg1} \mathbf{d}_{\text{B1}} \\ \mathbf{R}_{\text{leg1}} \rightarrow \text{beam} \mathbf{R}_{\text{Rm}} \rightarrow \text{leg1} \mathbf{!} \end{bmatrix} \right] = \mathbf{0} \quad (35) $$

$$ \sum_{j=1}^{5} S_{2j1}^{\perp} \left[ \mathbf{u}_2(q_2) - \begin{bmatrix} \mathbf{R}_{\text{leg2}} \rightarrow \text{beam} \mathbf{R}_{\text{Rm}} \rightarrow \text{leg2} \mathbf{d}_{\text{B2}} \\ \mathbf{R}_{\text{leg2}} \rightarrow \text{beam} \mathbf{R}_{\text{Rm}} \rightarrow \text{leg2} \mathbf{!} \end{bmatrix} \right] = \mathbf{0} \quad (36) $$

$$ \sum_{j=1}^{5} S_{3j1}^{\perp} \left[ \mathbf{u}_3(q_3) - \begin{bmatrix} \mathbf{R}_{\text{leg3}} \rightarrow \text{beam} \mathbf{R}_{\text{Rm}} \rightarrow \text{leg3} \mathbf{d}_{\text{B3}} \\ \mathbf{R}_{\text{leg3}} \rightarrow \text{beam} \mathbf{R}_{\text{Rm}} \rightarrow \text{leg3} \mathbf{!} \end{bmatrix} \right] = \mathbf{0} \quad (37) $$

$\mathbf{z}_1$ (resp. $\mathbf{z}_2$ and $\mathbf{z}_3$) is the unit vector along the direction of the ith prismatic joint, $i = 1, 2, 3$ as shown in Figure 8. The rotation matrix between the beam coordinate system and the leg coordinate system is expressed as:

$$ \mathbf{R}_{\text{beam}\rightarrow\text{leg}} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ -1 & 0 & 0 \end{bmatrix} \quad (38) $$

<figure><img src="image_12.png" /><figcaption>Figure 8. Leg and mobile platform coordinate systems and joint parametrization of Tripteor X7.</figcaption></figure>
The angles between the legs and the mobile platform are $\varphi_{24} = (\mathbf{z}_m, \mathbf{z}_1)$, $\varphi_{94} = (\mathbf{z}_m, \mathbf{z}_2)$ and $\varphi_{54} = (\mathbf{z}_m, \mathbf{z}_3)$ (Figure 8). These angles can be measured with a CAD model or computed with the geometric model of the Tripteor PKM described in [33]. With this notation, the rotation matrices between Legs 1, 2 and 3 and the mobile platform coordinate system $\mathcal{R}_m$ are:

$$ \mathbf{T}_{\text{leg1} \to \mathcal{R}_m} = \begin{bmatrix} \cos(\varphi_{24}) & 0 & \sin(\varphi_{24}) \\ 0 & 1 & 0 \\ -\sin(\varphi_{24}) & 0 & \cos(\varphi_{24}) \end{bmatrix} \qquad (39) $$

$$ \mathbf{T}_{\text{leg2} \to \mathcal{R}_m} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\varphi_{94}) & -\sin(\varphi_{94}) \\ 0 & \sin(\varphi_{94}) & \cos(\varphi_{94}) \end{bmatrix} \qquad (40) $$

$$ \mathbf{T}_{\text{leg3} \to \mathcal{R}_m} = \begin{bmatrix} \cos(\varphi_{54}) & 0 & \sin(\varphi_{54}) \\ 0 & 1 & 0 \\ -\sin(\varphi_{54}) & 0 & \cos(\varphi_{54}) \end{bmatrix} \qquad (41) $$

Consequently, in the case of Tripteor X7, by the application of a small-displacement screw relation at point $B_1$ and from Equations (35)–(37), nine constraint equations are obtained. Nine dependent parameters are chosen:

$$ \mathbf{u}_{\mathrm{D}}^{T} = (u_{2}(q_{2}), v_{2}(q_{2}), w_{2}(q_{2}), \theta_{\beta2}(q_{2}), u_{3}(q_{3}), v_{3}(q_{3}), w_{3}(q_{3}), \theta_{\alpha3}(q_{3}), \theta_{\gamma3}(q_{3})) \quad (42) $$

Thus, 16 independent parameters are considered:

$$ \mathbf{u}_{\mathrm{I}}^{T} = (\theta_{\beta1}(0), \theta_{\gamma1}(0), u_{1}(q_{1}), v_{1}(q_{1}), w_{1}(q_{1}), \theta_{\alpha1}(q_{1}), \theta_{\beta1}(q_{1}), \theta_{\gamma1}(q_{1}), \theta_{\alpha2}(0), \theta_{\beta2}(0), \theta_{\gamma2}(0), \theta_{\alpha2}(q_{2}), \theta_{\gamma2}(q_{2}), \theta_{\beta3}(0), \theta_{\gamma3}(0), \theta_{\beta3}(q_{3})) \quad (43) $$

The (9 × 9)-dimensional matrix $\mathbf{L}_D$ and the (9 × 16)-dimensional matrix $\mathbf{L}_I$ obtained from Equations (35)–(37), and the chosen vectors $\mathbf{u}_D$ and $\mathbf{u}_I$, are expressed in Equations (44) and (45).

$$ \mathbf{L}_{\mathrm{D}} = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & \sin(\varphi_{54}) & -\cos(\varphi_{54}) \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cos(\varphi_{54}) & \sin(\varphi_{54}) \\ 0 & 0 & 0 & -\sin(\varphi_{94}) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ -\sin(\varphi_{94}) & -\cos(\varphi_{94}) & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \cos(\varphi_{94}) & -\sin(\varphi_{94}) & 0 & -\cos(\varphi_{94}) & \mathbf{b}_{12} \cdot \mathbf{x}_m & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \sin(\varphi_{54}) & 0 & -\cos(\varphi_{54}) & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\cos(\varphi_{94}) & \mathbf{b}_{13} \cdot \mathbf{x}_m & \cos(\varphi_{54}) & 0 & \sin(\varphi_{54}) & 0 \end{bmatrix} \quad (44) $$

$$ \mathbf{L}_{\mathrm{I}} = \begin{bmatrix} 0 & 0 & 0 & -\sin(\varphi_{24}) & 0 & \cos(\varphi_{24}) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\cos(\varphi_{24}) & 0 & -\sin(\varphi_{24}) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\cos(\varphi_{24}) & 0 & -\sin(\varphi_{24}) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\cos(\varphi_{24}) & 0 & -\sin(\varphi_{24}) & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\sin(\varphi_{24}) & \cos(\varphi_{24}) & \cos(\varphi_{24}) & \mathbf{b}_{12} \cdot \mathbf{y}_m & 0 & \sin(\varphi_{24}) & \mathbf{b}_{12} \cdot \mathbf{y}_m & 0 & 0 \\ 0 & 0 & 1 & 0 & -\cos(\varphi_{24}) & \mathbf{b}_{12} \cdot \mathbf{x}_m & 0 & -\sin(\varphi_{24}) & \mathbf{b}_{12} \cdot \mathbf{x}_{\mathrm{MPS}} & 0 & 0 \\ 0 & 0 & -\cos(\varphi_{24}) & -\sin(\varphi_{24}) & -\sin(\varphi_{24}) & \mathbf{b}_{12} \cdot \mathbf{y}_m & 0 & \cos(\varphi_{24}) & \mathbf{b}_{12} \cdot \mathbf{y}_m & 0 & 0 & 0 \\ 0 & 0 & -\sin(\varphi_{24}) & \cos(\varphi_{24}) & 0 & 0 & 0 & 0 & 0 & -\sin(\varphi_{94}) & \mathbf{b}_{12} \cdot \mathbf{x}_m & 0 \\ 0 & 0 & 1 & 0 & -\cos(\varphi_{24}) & \mathbf{b}_{13} \cdot \mathbf{x}_m & 0 & -\sin(\varphi_{24}) & \mathbf{b}_{13} \cdot \mathbf{x}_m & 0 & 0 & 0 \\ 0 & 0 & -\cos(\varphi_{24}) & -\sin(\varphi_{24}) & 0 & 0 & 0 & 0 & 0 & -\sin(\varphi_{94}) & \mathbf{b}_{13} \cdot \mathbf{x}_m & 0 \end{bmatrix} \quad (45) $$

### 4.3. Cartesian Stiffness Matrix

As discussed in Section 2, Tripteor X7 is modeled as a structural assembly of three legs, which link the rigid fixed platform to the rigid mobile platform. The three legs are considered to be steel beams of rectangular cross-section of size consistent with Tripteor X7 legs (Figure 9).
<figure><img src="image_13.png" /><figcaption>**Figure 9.** CAD model of the Tripteor X7 PKM and cross-section of its first leg.</figcaption></figure>

To compute the Cartesian stiffness matrix of the robot at point $O_m$, the issue is to compute matrix $L_{XI}$ (Equation (29)). Finally, we can compute the stiffness map at a level $z = 1.18$ m for example (Figure 10).

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Principal stiffness coefficients of the Tripteor X7 along X, Y and Z axes in the XY-plane, $z = 1.18$ m.</figcaption></figure>

To plot Figure 11, the method was implemented in Matlab®. To emphasize the benefit of this proposed approach, the calculation was performed on a tablet PC with a I5-7300U processor. The time taken to compute this map is 128.72 s and is mainly due to the time needed to calculate the Inverse Kinematic Model (IKM) of Tripteor X7.

The main advantage of the method is the fast calculation of the natural frequencies and associated modes of over-constrained parallel robots.

### 4.4. Natural Frequencies

The proposed method is used for the calculation, with a low computational cost, of the first natural frequencies and associated natural modes of the Tripteor X7 PKM.

For this first computation, the leg lengths are $q_1 = 1.212$ m, $q_2 = 1.314$ m and $q_3 = 1.22$ m, and the mass of the mobile platform is considered to be 8 kg. The approach introduced in Section 3 is applied that gives the following estimations for the first natural frequencies in 0.4 s:
$$f_0 = \begin{pmatrix} 46.37 \text{ Hz} \\ 51.26 \text{ Hz} \\ 52.35 \text{ Hz} \\ 90.20 \text{ Hz} \\ 140.11 \text{ Hz} \\ 249.16 \text{ Hz} \\ 297.70 \text{ Hz} \\ 328.35 \text{ Hz} \\ 350.00 \text{ Hz} \\ 389.74 \text{ Hz} \end{pmatrix} \qquad (46)$$

The associated modes of the first three natural frequencies given in Equation (46) are schematized in Figure 11. The variation of first natural frequency value at a level z = 1.18 m is illustrated in Figure 12.

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Figure 11. Modal shapes associated with the first three natural frequencies of Triptor X7.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 12. Variation of the first natural frequency.</figcaption></figure>
## 4.5. Comparison with Experimental Results

In the literature, previous works studied the modal analysis of the Tripteor X7 PKM architecture using a shaker in three configurations [34] (Figure 13). The displacements of several points of the PKM architecture are measured with accelerometers located on the legs, and on the spindle. Nine DoF are considered for each leg, eight about *fi*-axis or *fl*-axis and one about *ff*-axis, and four DoF for the spindle (Figure 14). The obtained results show four natural frequencies close to 40, 50, 70 and 90 Hz.

<figure><img src="image_19.png" /><figcaption>Figure 13. First natural frequency measurement.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 14. Configuration of the Tripteor X7 for the experimental measurement.</figcaption></figure>

<table><thead><tr><th></th><th>Config. 1</th><th>Config. 2</th><th>Config. 3</th></tr></thead><tbody><tr><th>q1</th><td>1212 mm</td><td>1330 mm</td><td>921 mm</td></tr><tr><th>q2</th><td>1314 mm</td><td>1584 mm</td><td>1038 mm</td></tr><tr><th>q3</th><td>1220 mm</td><td>1528 mm</td><td>1167 mm</td></tr></tbody></table>

The same configurations 1, 2 and 3 were considered with the proposed beam approach, and the associated results are given in Table 1. Experimental and simulated results are close for the first mode. For the other modes, the simplification of the legs modeling and the difference of mass distribution between the model and the real machine tool can explain the gap between estimated and measured natural frequencies. Indeed, the modeling of the leg under the assumption that their cross-section area is constant affects the leg stiffness, inertia and the position of its gravity center. This impact is different regarding the leg lengths and so to the studied configuration. It should be noted that this method is dedicated to the embodiment design of parallel manipulators. Therefore, such an error is acceptable,
and the simplification of the leg shape allows the reduction of the calculation time of the elasto-dynamic model.

<table><caption>**Table 1.** Comparison of measured natural frequencies and calculated ones with the proposed method.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Config. 1</th><th colspan="3">Config. 2</th><th colspan="3">Config. 3</th></tr><tr><th>Beam Model</th><th>Measure</th><th>Error</th><th>Beam Model</th><th>Measure</th><th>Error</th><th>Beam Model</th><th>Measure</th><th>Error</th></tr></thead><tbody><tr><td>Mode 1</td><td>46.45 Hz</td><td>46.2 Hz</td><td>0.5%</td><td>34.84 Hz</td><td>39.6 Hz</td><td>12%</td><td>55.40 Hz</td><td>49.4 Hz</td><td>12.1%</td></tr><tr><td>Mode 2</td><td>50.98 Hz</td><td>55.3 Hz</td><td>7.8%</td><td>36.51 Hz</td><td>51.6 Hz</td><td>29.2%</td><td>66.66 Hz</td><td>50.2 Hz</td><td>32.8%</td></tr><tr><td>Mode 3</td><td>52.33 Hz</td><td>71.3 Hz</td><td>26.6%</td><td>44.46 Hz</td><td>58.4 Hz</td><td>23.9%</td><td>77.66 Hz</td><td>71.3 Hz</td><td>8.9%</td></tr></tbody></table>

## 5. Conclusions

This paper described a new methodology to compute the first natural frequencies of over-constrained PKMs based on the application of screw theory under the assumption of small displacements. Specifically designed to require few computing resources, the obtained model enables the estimation of stiffness maps, first natural frequencies and associated modes in a simple way without computing any Jacobian matrix. The proposed methodology was validated by comparing the obtained theoretical natural frequencies and associated modes of the Tripteor X7 to the measured experimental values. Errors between experimental and simulated results are less than 13% for the first mode and 33% for the second and third modes.

As a conclusion, the proposed methodology allows the mathematical expression of the simplified elasto-dynamic model of over-constrained parallel robots to reduce considerably the computation time of their natural frequencies. This low computation time allows the designer of over-constrained parallel robots to quickly estimate the natural frequencies of candidate robot architectures at the conceptual design stage and thus make the right choices of robot architecture with respect to required elasto-dynamic performance and a given task.

Later, variations in leg cross-sections will be implemented and joint stiffness will be added as an additional stiffness at beam extremities.

**Author Contributions:** Conceptualization, H.C., A.G., B.B. and S.C.; methodology, H.C. and S.C.; validation, H.C.; formal analysis, H.C., B.B. and S.C.; writing—original draft preparation, H.C.; writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This work was carried out within the Manufacturing 21 working group, which comprises 18 French research laboratories. The topics approached are: modeling of the manufacturing process, virtual machining, emergence of new manufacturing methods.

**Conflicts of Interest:** The authors declare no conflict of interest

## Abbreviations

The following abbreviations are used in this manuscript:

* $q_i$: ith joint variable

* $x$: Beam abscissa

* $\mathbf{u}(x)$: Displacement vector of the beam neutral axis

* $\mathbf{u}_{nod}$: Vector collecting the nodal values of the kinematic functions

* $u(x)$: Beam displacement along the neutral axis

* $v(x)$: Beam displacement along fi-axis

* $w(x)$: Beam displacement along fl-axis

* $\theta_{\alpha}(x)$: Beam section rotation about neutral axis

* $\theta_{\beta}(x)$: Beam section rotation about fi-axis

* $\theta_{\gamma}(x)$: Beam section rotation about fl-axis
<table><tr><td><b><i>f</i><sub>beam</sub>(<i>x</i>)</b></td><td>Beam force vector at a given abscissa</td></tr><tr><td><b><i>m</i><sub>beam</sub>(<i>x</i>)</b></td><td>Beam moment vector at a given abscissa</td></tr><tr><td><b><i>K</i><sub>e</sub></b></td><td>Beam stiffness matrix</td></tr><tr><td><b><i>M</i><sub>e</sub></b></td><td>Beam mass matrix</td></tr><tr><td><b><i>d</i><sub>B<sub>i</sub></sub></b></td><td>Vector of small displacements of point <i>B</i><sub>i</sub></td></tr><tr><td>!</td><td>Rotation vector of the mobile platform</td></tr><tr><td><b><i>u</i><sub>D</sub></b></td><td>Vector of dependent parameters</td></tr><tr><td><b><i>u</i><sub>I</sub></b></td><td>Vector of independent parameters</td></tr><tr><td><b><i>K̃</i></b></td><td>Global stiffness matrix</td></tr><tr><td><b><i>M̃</i></b></td><td>Global mass matrix</td></tr><tr><td><b><i>f</i><sub>0</sub></b></td><td>First natural frequency vector</td></tr></table>

## References

1. Patel, Y.; George, P. Parallel manipulators applications—A survey. *Mod. Mech. Eng.* **2012**, *2*, 57. [CrossRef]

2. Li, Y.; Ma, Y.; Liu, S.; Luo, Z.; Mei, J.; Huang, T.; Chetwynd, D. Integrated design of a 4-DOF high-speed pick-and-place parallel robot. *CIRP Ann.* **2014**, *63*, 185–188. [CrossRef]

3. Chanal, H.; Duc, E.; Ray, P. A study of the impact of machine tool structure on machining processes. *Int. J. Mach. Tools Manuf.* **2006**, *46*, 98–106. [CrossRef]

4. Terrier, M.; Dugas, A.; Hascoet, J.Y. Qualification of parallel kinematics machines in high-speed milling on free form surfaces. *Int. J. Mach. Tools Manuf.* **2004**, *44*, 865–877. [CrossRef]

5. Pritschow, G.; Eppler, C.; Garber, T. Influence of the dynamic stiffness on the accuracy of PKM. In Proceedings of the Chemnitz Parallel Kinematic Seminar, Chemnitz, Germany, 23–25 April 2002; pp. 313–333.

6. Pateloup, S.; Chanal, H.; Duc, E. Process parameter definition with respect to the behaviour of complex kinematic machine tools. *Int. J. Adv. Manuf. Technol.* **2013**, *69*, 1233–1248. [CrossRef]

7. Neumann, K. Exechon concept-parallel kinematic machines in research and practice. In Proceedings of the Chemnitz Parallel Kinematics Seminar, Chemnitz, Germany, 25–26 April 2006.

8. IFTOMM. IFTOMM Dictionary—3. Dynamics. Available online: http://www.iftomm-terminology.antonkb.nl/2057_1036/frames.html (accessed on 6 March 2022).

9. Liu, W.L.; Xu, Y.D.; Yao, J.T.; Zhao, Y.S. Methods for force analysis of overconstrained parallel mechanisms: A review. *Chin. J. Mech. Eng.* **2017**, *30*, 1460–1472. [CrossRef]

10. Mei, B.; Xie, F.; Liu, X.J.; Yang, C. Elasto-geometrical error modeling and compensation of a five-axis parallel machining robot. *Precis. Eng.* **2021**, *69*, 48–61. [CrossRef]

11. Benti, A.; Durieux, S.; Chanal, H. Mesure du Comportement Vibratoire D'une Machine-Outil à Structure Parallèle. 8ème Assises MUGV 2014. Available online: https://tel.archives-ouvertes.fr/tel-00451081/document (accessed on 9 March 2022).

12. Pierrot, F.; Pierrot, F. Modelling and design issues of a 3-axis parallel machine-tool. *Mech. Mach. Theory* **2002**, *11*, 1325–1345.

13. Yan, S.; Ong, S.; Nee, A. Stiffness analysis of parallelogram-type parallel manipulators using a strain energy method. *Robot. Comput.-Integr. Manuf.* **2016**, *37*, 13–22. [CrossRef]

14. Wang, Y.; Huang, T.; Zhao, X.; Mei, J.; Chetwynd, D.G.; Hu, S.J. Finite element analysis and comparison of two hybrid robots-the tricept and the trivariant. In Proceedings of the 2006 IEEE/RSJ International Conference on Intelligent Robots and Systems, Beijing, China, 9–15 October 2006; pp. 490–495.

15. Klimchik, A.; Pashkevich, A.; Chablat, D. CAD-based approach for identification of elasto-static parameters of robotic manipulators. *Finite Elem. Anal. Des.* **2013**, *75*, 19–30. [CrossRef]

16. Germain, C.; Briot, S.; Caro, S.; Wenger, P. Natural frequency computation of parallel robots. *J. Comput. Nonlinear Dyn.* **2015**, *10*, 021004. [CrossRef]

17. Zhang, J.; Zhao, Y.Q.; Jin, Y. Elastodynamic modeling and analysis for an Exechon parallel kinematic machine. *J. Manuf. Sci. Eng.* **2016**, *138*, 031011. [CrossRef]

18. Kermanian, A.; Kamali, A.; Taghvaeipour, A. Dynamic analysis of flexible parallel robots via enhanced co-rotational and rigid finite element formulations. *Mech. Mach. Theory* **2019**, *139*, 144–173. [CrossRef]

19. Klimchik, A.; Pashkevich, A.; Chablat, D. Fundamentals of manipulator stiffness modeling using matrix structural analysis. *Mech. Mach. Theory* **2019**, *133*, 365–394. [CrossRef]

20. Guo, J.; Zhao, Y.; Chen, B.; Zhang, G.; Xu, Y.; Yao, J. Force Analysis of the Overconstrained Mechanisms Based on Equivalent Stiffness Considering Limb Axial Deformation. *Res. Sq.* **2021**. Available online: https://www.researchsquare.com/article/rs-634709/v1 (accessed on 6 March 2022). [CrossRef]

21. Yang, C.; Li, Q.; Chen, Q.; Xu, L. Elastostatic stiffness modeling of overconstrained parallel manipulators. *Mech. Mach. Theory* **2018**, *122*, 58–74. [CrossRef]

22. Bi, Z. Kinetostatic modeling of Exechon parallel kinematic machine for stiffness analysis. *Int. J. Adv. Manuf. Technol.* **2014**, *71*, 325–335. [CrossRef]
23. Bi, Z.M.; Jin, Y. Kinematic modeling of Exechon parallel kinematic machine. *Robot. Comput.-Integr. Manuf.* **2011**, *27*, 186–193.
[CrossRef]

24. Amine, S.; Caro, S.; Wenger, P. Constraint and singularity analysis of the Exechon. *Appl. Mech. Mater.* **2012**, *162*, 141–150.
[CrossRef]

25. Selig, J.; Ding, X. A screw theory of timoshenko beams. *J. Appl. Mech.* **2009**, *76*, 031003. [CrossRef]

26. Bourdet, P.; Mathieu, L.; Lartigue, C.; Ballu, A. The concept of small displacement torsor in metrology. *Ser. Adv. Math. Appl. Sci.* **1996**, *40*, 110–122.

27. Aydin, Y.; Kucuk, S. Quaternion based inverse kinematics for industrial robot manipulators with Euler wrist. In Proceedings of the 2006 IEEE International Conference on Mechatronics, Budapest, Hungary, 3–5 July 2006; pp. 581–586.

28. Bonnemains, T.; Chanal, H.; Bouzgarrou, B.C.; Ray, P. Stiffness computation and identification of parallel kinematic machine tools. *J. Manuf. Sci. Eng.* **2009**, *131*, 041013. [CrossRef]

29. Klimchik, A.; Pashkevich, A.; Caro, S.; Chablat, D. Stiffness modeling of robotic-manipulators under auxiliary loadings. *International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Am. Soc. Mech. Eng.* **2012**, *45035*, 469–476.

30. Kucuk, S. Dexterous workspace optimization for a new hybrid parallel robot manipulator. *J. Mech. Robot.* **2018**, *10*, 064503.
[CrossRef]

31. Yoo, H.; Ryan, R.; Scott, R.A. Dynamic of flexible beams undergoing overall motions. *J. Sound Vib.* **1995**, *181*, 261–278. [CrossRef]

32. Kong, X.; Gosselin, C.M. *Type Synthesis of Parallel Mechanisms*; Springer: Berlin/Heidelberg, Germany, 2007; Volume 33.

33. Pateloup, S.; Chanal, H.; Duc, E. Geometric and kinematic modelling of a new Parallel Kinematic Machine Tool: The Tripteor X7 designed by PCI. *Adv. Mater. Res.* **2010**, *112*, 159–169. [CrossRef]

34. Bonnemains, T.; Chanal, H.; Bouzgarrou, B.C.; Ray, P. Dynamic analysis of the Tripteor X7: Model and experiments. *Robot. Comput.-Integr. Manuf.* **2010**, *29*, 180–191. [CrossRef]