R Jha

D Chablat
email: damien.chablat@irccyn.ec-nantes.fr

F Rouillier
email: fabrice.rouillier@inria.fr

G Moroz

AN ALGEBRAIC METHOD TO CHECK THE SINGULARITY-FREE PATHS FOR PARALLEL ROBOTS

Trajectory planning is a critical step while programming the parallel manipulators in a robotic cell. The main problem arises when there exists a singular configuration between the two poses of the end-effectors while discretizing the path with a classical approach. This paper presents an algebraic method to check the feasibility of any given trajectories in the workspace. The solutions of the polynomial equations associated with the trajectories are projected in the joint space using Gröbner based elimination methods and the remaining equations are expressed in a parametric form where the articular variables are functions of time t unlike any numerical or discretization method. These formal computations allow to write the Jacobian of the manipulator as a function of time and to check if its determinant can vanish between two poses. Another benefit of this approach is to use a largest workspace with a more complex shape than a cube, cylinder or sphere. For the Orthoglide, a three degrees of freedom parallel robot, three different trajectories are used to illustrate this method.

INTRODUCTION

One of the crucial steps in the trajectory planning is to check the singularity-free paths in the workspace for the parallel manipulators. It becomes a necessary protocol to validate the trajectory when the parallel robot is used in practical applications such as precise manufacturing operations. A trajectory verification problem is presented in [START_REF] Merlet | A generic trajectory verifier for the motion planning of parallel robots[END_REF] based on some validity criteria like whether the trajectory lies fully inside the workspace of the robot and is singularity-free. Singularity-free path planning for the Gough-Stewart platform with a given starting pose and a given ending pose has been addressed in [START_REF] Dasgupta | Singularity-free path planning for the Stewart platform manipulator[END_REF] using the clustering algorithm is presented in [START_REF] Dash | Singularityfree path planning of parallel manipulators using clustering algorithm and line geometry[END_REF]. An exact method and an approximate method are described in [START_REF] Bhattacharya | Comparison of an exact and an approximate method of singularity avoidance in platform type parallel manipulators[END_REF] to restructure a path close to the singularity locus into a path that avoids it while remaining close to the original path. Due to the geometry of the mechanism, the workspace may not cover fully the space of poses [START_REF] Dash | Singularityfree path planning of parallel manipulators using clustering algorithm and line geometry[END_REF], hence it is necessary to analyze the workspace of the manipulator.

The workspace of a parallel robot mainly depends on the actuated joint variables, the range of motion of the joints and the mechanical interferences between the bodies of the mechanism. There are different techniques based on geometric tools [START_REF] Gosselin | Determination of the workspace of 6-DOF parallel manipulators[END_REF][START_REF] Merlet | Geometrical determination of the workspace of a constrained parallel manipulator[END_REF], discretization [START_REF] Castelli | A fairly general algorithm to evaluate workspace characteristics of serial and parallel manipulators[END_REF][START_REF] Ilian | A new approach to orientation workspace analysis of 6-DOF parallel manipulators[END_REF][START_REF] Chablat | An Interval Analysis Based Study for the Design and the Comparison of 3-DOF Parallel Kinematic Machines[END_REF], and algebraic methods [START_REF] Zein | An exhaustive study of the workspace topologies of all 3R orthogonal manipulators with geometric simplifications[END_REF][START_REF] Ottaviano | Identification of the workspace boundary of a general 3-R manipulator[END_REF][START_REF] Chablat | Non-singular assembly mode changing trajectories in the workspace for the 3-RPS parallel robot[END_REF][START_REF] Chablat | Uniqueness domains and non singular assembly mode changing trajectories[END_REF] which are used to compute the workspace of a parallel robot. An algebraic method to solve the forward kinematics problem specifically applied to spatial parallel manipulators is described in [START_REF] Rolland | Certified solving of the forward kinematics problem with an exact algebraic method for the general parallel manipulator[END_REF]. The main advantage of the geometric approach is that, it establishes the nature of the boundary of the workspace [START_REF] Siciliano | Springer handbook of robotics[END_REF]. A procedure to automatically generate the kinematic model of parallel mechanisms which further used for singularity free path planning is reported in [START_REF] Lahouar | Singularity Free Path Planning for Parallel Robots[END_REF] .

An algorithm for computing singularity-free paths on 1 closed-chain manipulators is presented in [START_REF] Bohigas | Planning singularity-free paths on closed-chain manipulators[END_REF], also this method attempts to connect the given two non-singular configurations through a path that maintains a minimum clearance with respect to the singularity locus at all points. The main drawback of any numerical or discretization methods is that there might be a singular configuration between two poses of the end-effectors while discretizing the path. This paper illustrates a technique based on some algebraic methods to check the feasibility of any given trajectories in the workspace : it allows to write the Jacobian of the manipulator as a function of time and to check whether its determinant vanishes between two poses. Also, when the trajectory meets a singularity, its location can also be computed.

The outline of this paper is as follows. We first introduce the modelization and the basic algebraic tools for computing trajectories in the workspace for parallel robots. We then extended the process to check the singularity-free paths within the workspace and ensures the existence of a singular configurations between the two poses of the end-effector by analysing the Jacobian as a function of time. In later sections we give some detailed examples using three different trajectories for the Orthoglide, a three degrees of freedom parallel robot.

METHODOLOGY

To ensure the non existence of a singular configurations between two poses of the end-effector, it is necessary to express the Jacobian of the manipulator as a function of the time. The general procedure to check the feasibility of the trajectories in the workspace as well as a function to compute the projection of the trajectories in the joint space can be decomposed as follows :

• Defining the constraint equations, articular and pose variables associated with the parallel mechanism; • Computing the singularities and their projections in the workspace and in the joint space; • Computing the workspace and joint space boundaries; • Computing a parametric form of the trajectories in the Cartesian space; • Computing the projections of the trajectories in the joint space; • Computing the Jacobian of the manipulator as a function of the time.

By computing, we mean, by default, getting a full characterization as solutions of an exact system of algebraic equations. By abuse of language, computing a projection might thus mean computing the algebraic closure of the projection, for example by using classical elimination strategies based on Gröbner bases (which thus could introduce a subset of null measure made of spurious points).

Kinematics involves the study of the position, velocity, acceleration, and all higher order derivatives of the position variables (with respect to time or any other parameter/variables).

Hence, the study of the kinematics of manipulators refers to all the geometrical and time-based properties of the motion. For a translational manipulator, the set of (algebraic) relations that connects the input values ρ ρ ρ ρ ρ ρ to the output values X will be denoted by

F(ρ ρ ρ ρ ρ ρ, X) = 0 (1)
where ρ ρ ρ ρ ρ ρ is the set of all actuated joint variables and X is the set of all pose variables of the end-effector.

For the study of manipulators, two problems must be considered: the direct kinematic problem (DKP) and the inverse kinematic problem (IKP). Specifically, given a set of joint values, solving the DKP consists in computing the position and the orientation of the end-effector relatively to the base, or, equivalently to change the representation of the manipulator position from a joint space description into a Cartesian space description. Given the position and orientation of the end-effector of the manipulator, solving the IKP problem consists in computing all the possible sets of joint angles or parameters that could be used to attain this position and orientation, or, equivalently, mapping locations in three-dimensional Cartesian space to locations in the robot's internal joint space.

When considering an algebraic modelization, the direct kinematics has basically several solutions, which refers to the several poses of the end-effector for given values of the joint coordinates. It is therefore possible to assemble the manipulator in different ways, and these different configurations are known as the assembly modes of the manipulator [START_REF] Chablat | Working Modes and Aspects in Fully-Parallel Manipulator[END_REF]. Similarly, the multiple inverse kinematic solutions induce multiple postures for each leg of the manipulator and is termed as the working modes of the manipulator. An analytical relation can be given as ρ = γ(X) for IKP whereas X = β (ρ) for DKP. Differentiating Eq. ( 1) with respect to time leads to the velocity model:

A Ẋ + B ρ = 0 where A = ∂ F ∂ X , B = ∂ F ∂ ρ ρ ρ ρ ρ ρ (2)
The matrices A and B are respectively the direct-kinematics and the inverse-kinematics Jacobian matrices of the manipulator. These matrices are used for characterizing different kinds of singularities. The parallel singularities occur whenever det(A) = 0 and the serial singularities occur whenever det(B) = 0. The parallel singular configurations are located inside the workspace. They are particularly undesirable because the manipulator cannot resist to any forces and its control is lost.

Eliminating ρ in the system F(ρ, X) = 0, det(A) = 0 by means of a Gröbner basis computation for a suitable elimination ordering (see [START_REF] Cox | Using Algebraic Geometry[END_REF]) defines (the algebraic closure of) the projection ξ (X) of the parallel singularities in the workspace. In the same way, one can compute (the algebraic closure of) the projection ε(X) of the parallel singularities in the joint space. Both are then defined as the zero set of some system of algebraic equations and we assume that the considered robots are generic enough so that both are hypersurfaces.

The set of equations associated with the joint limits of the actuator χ(ρ) can also be projected, with the same elimination method, in the workspace µ(X) as in Eq. ( 3), where ρ min , ρ max are the minimum and maximum values of the articular variables: χ(ρ) and µ(X) are the crucial parameters in determining the number of assembly modes and working modes of the parallel manipulator.

χ(ρ) → µ(X) ∀ρ ∈ (ρ min , ρ max ] (3)
To use algebraic tools such as Gröbner bases, we must represent the trajectories in an algebraic form. A classical approach can be used to transform the trigonometric equations to algebraic ones as in [START_REF] Moroz | On the determination of cusp points of 3-RPR parallel manipulators[END_REF]. Equation ( 4) represents the trajectory in parametric form within the workspace as the function of time t.

X → φ (t) (4)
Ψ(X, ρ, t) in Eq. ( 5) is the system of equations which contains kinematic equations F(ρ, X) and the parametric equations of trajectory φ (t).

Ψ(X, ρ, t) = [φ (t) -X, F(ρ, X)] (5)
This system of equations is projected in the joint space as ϒ(ρ, t) as shown in Eq. ( 6) using Gröbner based elimination method.

Ψ(X, ρ, t) → ϒ(ρ, t) (6)
The parametric equations of the trajectory as a function of the time in the joint space can then be obtained by solving ρ ← ϒ(ρ, t).

The workspace analysis allows the characterization of the workspace regions where the number of real solutions for the inverse kinematics is constant. A cylindrical algebraic decomposition (CAD) algorithm is used to compute the workspace of the robot in the projection space X with joint constraints χ(ρ) taken in account [START_REF] Chablat | Workspace and joint space analysis of the 3-RPS parallel robot[END_REF][START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Moroz | Sur la dcomposition relle et algbrique des systmes dpendant de paramtres[END_REF]. The three main steps involved in the analysis are:

• Computation of a subset of the joint space (resp. workspace)

where the number of solutions changes: the Discriminant Variety . • Description of the complementary of the discriminant variety in connected cells: the Generic Cylindrical Algebraic Decomposition (CAD). • Connection of the cells belonging to the same connected component in the counterpart of the discriminant variety: interval comparisons.

The joint space analysis predicts the feasible and nonfeasible combinations of the prismatic joint variables which are essential for the parallel robot control. The joint space analysis allows the characterization of the regions where the number of real solutions for the direct kinematic problem is constant. The joint space analysis is done using CAD, which gives the number of cells corresponding to the number of solutions in the joint space. χ(ρ) parameter significantly changes the number of assembly modes and working modes of the manipulator.

Substituting φ (t) in ξ (X) (in the polynomial defining ξ (X) ) we then get ξ (φ (t)), which defines the vanishing values of the Jacobian as a function of the time. This equation in cos(t), sin(t) can be turned into a zero-dimensional bivariate system and its solutions can thus be computed exactly, either by means of a rational parametrization or by means of isolating intervals with rational bounds that can be refined to any arbitrary precision 1 .

Whatever the chosen (exact) representation for a solution, it can easily be checked if it vanishes between the two poses of the end-effector. The solutions of ξ (φ (t))contain the singular points on the trajectory φ (t) and eventually few spurious points due to the projection of det(A) in Cartesian space. Spurious singular points can then be differentiated from real singular points by substituting φ (t) and the solutions of ϒ(ρ, t) in det(A).

METHOD VALIDATION

We propose to validate our approach (and related tools) by checking the feasibility of three different trajectories for the Othoglide. The first two trajectories are heart shaped planar parametric curves and the third trajectory is a parametric curve in three dimensions.

These trajectories are fully inside the workspace and are selected such that their projections in the joint space is defined by parametric equations containing trigonometric functions. In addition, we have shown that the singular and singular-free trajectories are well discriminated, as one of the trajectory cuts the singularity surface.

Manipulator Architecture and Kinematics

The manipulator under the study is an Orthoglide parallel robot with three degrees of freedom. The mechanism is driven by three actuated orthogonal prismatic joints ρ i , made of three parallelogram connected by revolute joints to the tool center point on one side and to the corresponding prismatic joint at another side. The assembly modes of these robots depends on the solutions of the DKP as shown in Fig. 1(a). The point P i represents the pose of corresponding robot. However more than one value of i for the point P i denote multiple solutions for the DKP. A i B i is equal to ρ i , where ρ i represents the prismatic joint variables whereas P represents the position vector of the tool center point.

The constraint equations for the Orthoglide with the actuated variables ρ = [ρ 1 , ρ 2 , ρ 3 ] and the pose variables X = [x, y, z] are: where l is the leg length of the manipulator and is equal to two for all the computations. The set of equations associated with the joint limits of the actuator χ(ρ), projected in the workspace µ(X), are given in Eq. ( 8). χ(ρ) plays an important role in determining the shape of the workspace and singularity surfaces. It also affects the number of solutions for the IKP i.e. the working modes associated with the manipulator. Figure 1(b) represents the workspace of the Orthoglide. A CAD algorithm is used to compute the workspace of the robot in the projection space (x, y, z), taking into account the joint constraints χ(ρ) . Without considering the joint limits, the Orthoglide admits two assembly modes and eight working modes [START_REF] Pashkevich | Kinematics and workspace analysis of a three-axis parallel manipulator the Orthoglide[END_REF].

F(ρ, X) : (x -ρ 1 ) 2 + y 2 + z 2 = l 2 x 2 + (y -ρ 2 ) 2 + z 2 = l 2 x 2 + y 2 + (z -ρ 3 ) 2 = l 2 (7)
χ(ρ) = [ρ 1 , ρ 2 , ρ 3 , -4 + ρ 1 , -4 + ρ 2 , -4 + ρ 3 ] χ(ρ) → µ(X) ρ i → x 2 + y 2 + z 2 -4 i = 1, 2, 3 -4 + ρ i → x 2 + y 2 + z 2 -8x + 12 -4 + ρ i → x 2 + y 2 + z 2 -8y + 12 -4 + ρ i → x 2 + y 2 + z 2 -8z + 12 (8)
In the Figure 2(b), the yellow region (resp. green region), corresponds to the region where the inverse kinematic model has two real solutions (resp. height).

The joint space analysis is done using CAD to find the regions where the direct kinematics problem admits a fixed number of real solutions. For example, in Fig. 1(d), the red cell corresponds to to the region where the DKP has two solutions.

There is a difference in the shape of the singularity surface, depending if we consider the joint constraints µ(X) or not (Fig. 2).

det(A) = -8ρ 1 ρ 2 ρ 3 + 8ρ 1 ρ 2 z + 8ρ 1 ρ 3 y + 8ρ 2 ρ 3 x det(A) → ε(ρ) det(A) → ξ (X) ε(ρ) = ρ 4 1 ρ 2 2 + ρ 4 1 ρ 2 3 + ρ 2 1 ρ 4 2 + 2ρ 2 1 ρ 2 2 ρ 2 3 + ρ 2 1 ρ 4 3 + ρ 4 2 ρ 2 3 + ρ 2 2 ρ 4 3 -16ρ 2 1 ρ 2 2 -16ρ 2 1 ρ 2 3 -16ρ 2 2 ρ 2 3 ( 9
)
In Eq. ( 9), det(A) is the parallel singularity of the Orthoglide and ξ (X) is the projection of det(A) in workspace. The mathematical expression for ξ (X) is not displayed in Eq. ( 9) due to the lack of space. Fig. 1(c) shows the projection ξ (X) which is plotted with µ(X) as one of the input parameter. The degree of this characteristic surface is 18 and it represents the singularities associated with the eight working modes.

Trajectory definition in the workspace

Trajectories 1 & 2 are heart shaped parametric curves which ( Fig. 3(a)). Eq. ( 10) and Eq. ( 11) are the mathematical definitions of Trajectory 1 and Trajectory 2, respectively. As these equations are trigonometric equations, it is necessary to represent them in an algebraic form. φ 1 (t) and φ 2 (t) are the parametric equations for Trajectory 1 and Trajectory 2, respectively, which are defined for t ∈ [-π, π]. From now, sin and cos are replaced by s and c respectively to reduce the size of the expressions. φ 1 (t) :

x = 8 7 s 3 (t)

y = 13 14 c(t) - 5 14 c(2t) - 1 10 c(3t) - 1 14 c(4t) z = 1 (10) φ 2 (t) : x = 4 5 s 3 (t) y = 13 20 c(t) - 1 4 c(2t) - 1 20 c(3t) - 1 20 c(4t) z = 1 ( 11
)
Trajectory 3 is a parametric curve in three dimensions. Fig. 3(b) represents the trajectory in the workspace of the manipulator. The parametric equations of the trajectory is defined by Eq. [START_REF] Zein | An exhaustive study of the workspace topologies of all 3R orthogonal manipulators with geometric simplifications[END_REF]. φ 3 (t) is the parametric equation for Trajectory 3 and is defined for t ∈ [0, 20]. In order to turn the system to an algebraic one, we add the following equations

φ 3 (t) : x = s(t) y = c(t) z = 1 20 t (12)
sin(t) = sin t cos(t) = cos t sin t 2 + cos t 2 = 1 ( 13
)
and we remark that this change of variables does not introduce any spurious solutions since it is bjective (t ∈ [-π, π]).

Projection in the joint space φ 1 (t), φ 2 (t) and φ 3 (t) are the trajectories which are defined in the workspace. To project these trajectories in the joint space, it is necessary to formulate a system of equations corresponding to each trajectory which also consists the kinematic equations of the manipulator. Ψ 1 , Ψ 2 and Ψ 3 are the corresponding systems of equations for the Trajectories 1, 2 and 3, respectively (see Eq. ( 14)).

Ψ 1 (X, ρ, t) = [φ 1 (t) -X, F(ρ, X)] Ψ 2 (X, ρ, t) = [φ 2 (t) -X, F(ρ, X)] Ψ 3 (X, ρ, t) = [φ 3 (t) -X, F(ρ, X)] ( 14
)
Each system of is projected in the joint space as ϒ 1 (ρ, t), ϒ 2 (ρ, t) and ϒ 3 (ρ, t) (see Eq. ( 15)). By solving ϒ 1 (ρ), ϒ 2 (ρ) and ϒ 3 (ρ) for ρ, we get the corresponding parametric equations for the trajectories in the joint space, as shown in Fig. 4 for the Trajectories 1&2 and in Fig. 5 for the Trajectory 3. Eqs. ( 16) -( 18) are used to plot all possible images of the trajectories in the joint space. All the computations are done without considering the joint limits. As the Orthoglide has eight working mode, there exists eight possible trajectories in the joint space ( Fig. 4 and Fig. 5 for Trajectories 1&2 and Trajectory 3, respectively).

Ψ 1 (X, ρ, t) → ϒ 1 (ρ, t) ∀t ∈ [-π, π] Ψ 2 (X, ρ, t) → ϒ 2 (ρ, t) ∀t ∈ [-π, π] Ψ 3 (X, ρ, t) → ϒ 3 (ρ, t) ∀t ∈ [0, 20] (15)
Solving ϒ 1 (ρ, t) and ϒ 2 (ρ, t) for ρ gives eight possible solutions ( Eq. ( 17) for Trajectory 1 and Eq. ( 18) for Trajectory 2), which inferred eight different possible images of Trajectory 1&2 in joint space. But from Fig. 4 one can see that only one trajectory lies inside the joint space boundary for the both trajectories in workspace.

ρ 1 = s(t) ± 1 20 400s 2 (t) -t 2 + 1200 ρ 2 = c(t) ± 1 20 400c 2 (t) -t 2 + 1200 ρ 3 = 1 20 t ± √ 3 ( 16
)
Figure 5 shows all the possible images of Trajectory 3 in the joint space. Solving ϒ 3 (ρ, t) for ρ gives eight possible solutions which (Eq. ( 16)). It can be seen in the Fig. 5 that only one trajectory (which is marked as 1), lies inside the joint space boundary.

Singularity analysis

To check if there exists any singular configuration between the two poses of the end-effector it is necessary to express the Jacobian of the manipulator as a function of the time or of some independent variable. The proposed algebraic method enables us to write the Jacobian of the manipulator as a function of the time and to check if its determinant vanishes between two poses.

ρ 1 = 8 7 s 3 (t)
By substituting the values of φ 1 (t), φ 2 (t) and φ 3 (t) from Eq. ( 10), Eq. ( 11) and Eq. ( 12) in µ(X), we will get µ(φ 1 (t)), µ(φ 2 (t)) and µ(φ 3 (t)) as the determinant of the Jacobian for Trajectory 1, Trajectory 2 and Trajectory 3 respectively (Eq. ( 19)). Due to the large expressions, the equations for µ 1 (t), µ 2 (t) and µ 3 (t) are not presented but their roots define the singular configurations : Figures 8,10 and 11 show the values of these functions when t varies.

µ 1 (t) = µ(φ 1 (t)) ∀t ∈ [-π, π] µ 2 (t) = µ(φ 2 (t)) ∀t ∈ [-π, π] µ 3 (t) = µ(φ 3 (t)) ∀t ∈ [0, 20] (19)
The number of solutions of Eq. ( 19) gives the total number of singular points on the corresponding trajectory. The solutions for µ 1 (t), µ 2 (t) and µ 3 (t) are shown in Eq. [START_REF] Moroz | On the determination of cusp points of 3-RPR parallel manipulators[END_REF]. From Eq. ( 20) we get that there exists four solutions for µ 1 (t) and zero solutions for µ 2 (t) and µ 3 (t), which confirms the presence of singular points on Trajectory 1 whereas Trajectory 2&3 are singularity-free trajectories. Note that numerical approximations are given for readability, but, in practice, isolating intervals with rational bounds with arbitrary width can be computed (see section METHODOLOGY or [25, chapter 8] for more details) so that the result is certified.

] ← µ 1 (t) = 0 ∀t ∈ [-π, π] t = [0.97, 1.51] ← det(A) = 0 ∀t ∈ [-π, π] t = [ ] ← µ 2 (t) = 0 ∀t ∈ [-π, π] t = [ ] ← µ 3 (t) = 0 ∀t ∈ [0, 20](20)
µ 1 (t) is the determinant of the Jacobian corresponding to the Trajectory 1. By solving a zero-dimensional system of equa- tions, it can be shown that Trajectory 1 cuts the singularity surface in four points s 1 , s 2 , s 3 and s 4 (see Fig. 8 and Eq. ( 20)). The image of these points in workspace is computed by substituting the values of Eq. ( 20) in Eq. [START_REF] Zein | An exhaustive study of the workspace topologies of all 3R orthogonal manipulators with geometric simplifications[END_REF]. Similarly, the image of these points in joint space is computed by substituting the values of Eq. (20) in Eq. [START_REF] Bohigas | Planning singularity-free paths on closed-chain manipulators[END_REF]. The obtained values are tabulated in Table 1. Figures 6 represents the Trajectory 1&2 along with singularity surface ξ (X), also, all the singular points s 1 , 2 , s 3 and s 4 are located on the singularity surface. s 1 and s 2 are the spurious singular points which were introduced by the projection of det(A) in Cartesian space. Substituting µ 1 (t) and φ 1 (t) in det(A) gives two solutions ( Eq. ( 20)). These two solutions s 3 and s 4 are the real singular points on Trajectory 1. The variation of det(A) along Trajectory 1 is shown in Fig. 9. The true singular postures s 3 and s 4 , belonging to the singularity surface ξ (X), which is associated with the working mode used by the Orthoglide prototype ( Fig. 7). µ 2 (t) and µ 3 (t) are the determinanst of the Jacobians corresponding to the Trajectories 2&3 respectively.

The variation of µ 2 (t) along Trajectory 2 for t ∈ [-π, π] is shown in Fig. 10. Figure 11 represents the variation of µ 3 (t) for t ∈ [0, 20] along Trajectory 3. All the values given in Table 1 are associated with the trajectories which is marked as 1 in Fig. 4. These values are obtained by solving a zero-dimensional system of polynomial equations using RootFinding:-Isolate so that the related values for sin t and cos t are certified.

In Table 1, s 1 and s 2 are the spurious singular points and s 3 and s 4 are the real singular points on the Trajectory 1. Trajectory 2&3 are the singularity-free trajectories.

Conclusions

This paper reports the use of algebraic methods to check the feasibility of given trajectories in the workspace. The general method uses the projection of the polynomial equations associated with the trajectories in the joint space using a Gröbner based elimination strategy. There is a significant change in the shape of workspace and the singularity surfaces due to the joint constraints. The proposed method enables us to project the joint constraints in the workspace of the manipulator, which further helps to analyse the projection of the singularities in the workspace with the joint constraints. Such computations ensure the existence of a singular configuration between two poses of the endeffector unlike other numerical or discretization methods. This paper highlights that the singularity analysis should be done in the cross product of the workspace and joint space for parallel robots with several assembly and working modes. The single analysis of the singularities in a projection space can introduce spurious singularities. In fact, the algebraic tools does not allow to distinguish the parallel singularities according to the working mode.

FIGURE 1 .

 1 FIGURE 1. Architecture of an Orthoglide including two assembly modes (a) Workspace plot with χ(ρ) joint constraints (b) and Parallel singularity det(A) projected in the workspace as ξ (X) (c) and joint space plot with χ(ρ) joint limits (d)

FIGURE 2 .

 2 FIGURE 2. A comparison between the parallel singularity surface for an Orthoglide computed with the joint limits µ(X) (a) without µ(X) (b)

FIGURE 3 .

 3 FIGURE 3. Position of Trajectories 1&2 (a) Trajectory 3 (b) in the workspace of an Orthoglide

FIGURE 4 .

 4 FIGURE 4. A pictorial representation of the mapping of trajectories from workspace to joint space. Eight different pairs of trajectories, ϒ 1 (ρ, t) & ϒ 2 (ρ, t) in joint space are the image of corresponding φ 1 (t) & φ 2 (t). These eight different trajectories are associated with the eight working modes of the Orthoglide. Only one trajectory lies inside the joint space boundary due to to the joint constraints.

FIGURE 5 .

 5 FIGURE 5. Mapping of the Trajectory 3 from workspace to joint space. Eight different possible solutions of ϒ 3 (ρ, t) are marked in joint space which are the image of φ 3 (t). There is only one feasible trajectory (marked as 1) lies inside the joint space boundary.

FIGURE 6 .

 6 FIGURE 6. Representation of the trajectories with singularity surface. Trajectory 1 cuts the parallel singularity surface ξ (X) in four points s 1 , s 2 , s 3 and s 4 in the workspace of the Orthoglide. Also, it can be seen that Trajectory 2 lies inside ξ (X) as µ 2 (t) = 0 ∀t ∈ [-π, π].

FIGURE 7 .

 7 FIGURE 7. Representation of the real singular points s 3 and s 4 along with the singularity surface ξ (X), which is associated with one working mode.

FIGURE 8 .

 8 FIGURE 8. Variation of µ 1 (t) along Trajectory 1. s 1 , s 2 , s 3 and s 4 represents the four solutions of µ 1 (t) = 0 ∀t ∈ [-π, π].

FIGURE 9 .

 9 FIGURE 9. Variation of det(A) along the Trajectory 1. s 3 and s 4 represents the two solutions of det(A) = 0 ∀t ∈ [-π, π].

FIGURE 10 .

 10 FIGURE 10. Variation of µ 2 (t) along the Trajectory 2. There are no solutions for µ 2 (t) = 0 ∀t ∈ [-π, π].

FIGURE 11 .

 11 FIGURE 11. Variation of µ 3 (t) along the Trajectory 3. There are no solutions for µ 3 (t) = 0 ∀t ∈ [-π, π].

  ± 1 35 1996s 6 (t) -400s 8 (t) + 560s 6 (t)c(t) -100s 4 (t)c(t) -2009s 4 (t) + 2190c 3 (t) -1379c 2 (t) -1320c(t) + 3988 (t) -16s 8 (t) + 16s 6 (t)c(t) + 12s 4 (t)c(t) -85s 4 (t) + 96c 3 (t) -66c 2 (t) -60c(t) + 321 (t) -16c 7 (t) 7 + 52c 6 (t) + 60c 5 (t) -145c 4 (t) + 24c 3 (t) + 132c 2 (t) -32c(t) (18)

	ρ 2 = -	4 7	c 4 (t) -	2 5	c 3 (t) -	1 7	c 2 (t) +	43 35	c(t) +	2 7	±	1 7	64c 6 (t) -192c 4 (t) + 192c 2 (t) + 83
	ρ 3 = 1 ±	1 35	3200 -400c 8 (t) -560c 7 (t) + 1204c 6 (t) + 1580c 5 (t) -3221c 4 (t) + 710c 3 (t) + 3051c 2 (t) -860c(t)	(17)
		ρ 1 = 76s 6 ρ 2 = -4 5 s 3 (t) ± 1 10 2 5 c 4 (t) 1 5 c 3 (t) -	1 10	2 (t) +	4 5	c(t) +	1 5	±	1 5	16c 6 (t) -48c 4 (t) + 48c 2 (t) + 59
		ρ 3 = 1 ±	1 10	332 -16c 8

TABLE 1 .

 1 Singular postures on Trajectory 1

		Workspace			Joint space
	S	x	y	z	ρ 1	ρ 2	ρ 3
	S 1 -1.13 0.35 1.00 0.55 1.66 2.60
	S 2 -0.65 0.80 1.00 0.88 2.40 2.71
	S 3	0.65	0.80 1.00 2.18 2.40 2.71
	S 4	1.13	0.35 1.00 2.83 1.66 2.60

for example the functions Groebner[RationalUnivariateRepresentation] and RootFinding[Isolate] in Maple software

ACKNOWLEDGMENT

The work presented in this paper was partially funded by the Erasmus Mundus project "India4EU II".


Jafar Al-Muhammad
email: jafar.al-muhammad@irstea.fr

Séverine Tomas

Nassim Ait-Mouheb

Muriel Amielh

Fabien Anselmet

µ-PIV characterization of the flow in a milli-labyrinth-channel used in micro-irrigation

The baffle-fitted labyrinth-channel is largely u sed in m icro-irrigation system s. The existing baffles, w hich p lay an im p ortant role for generating p ressu re losses and ensu re the flow regu lation on the irrigation netw ork, p rod u ce vorticities w here the velocity is low or zero. These vorticities favor the d ep osition of p articles or other biochem ical d evelop m ent cau sing em itter clogging. Flow top ology characterization in this labyrinth-channel m u st be d escribed to analyze em itter clogging sensibility w hich d rastically red u ces its p erform ance. This characterization can be p erform ed exp erim entally u sing the m icro-p article im age velocim etry (µ-PIV) techniqu e. In this stu d y, µ-PIV exp erim ents are p erform ed to visu alize the flow in tenp attern rep eating baffles u sed in m icro irrigation em itter. Cross section is 1m m 2 . Reynold s nu m ber varies from 400 to 800. This investigation allow s analyzing the flow regim e. It help s to fix the baffle nu m ber from w here the flow is d evelop ed . Finally, the µ-PIV and m od eling resu lts com p arisons are p resented . The global objective of this stu d y is to id entify and analyze the sensitive areas in ord er to red u ce them thanks to geom etry op tim ization. ρ p Particle d ensity (kg/ m 3 ) μ Dynam ic viscosity (Pa.s)

List of symbols

Introduction

Micro-irrigation is a technique characterized by low -velocity w ater flow s. The w ater d rops near the plants through em itters. This type of irrigation im proves efficiency by red ucing evaporation, d rift, runoff and d eep percolation losses w hen com pared w ith the other techniques such as sprinkler irrigation. In this technique, em itters are the m ost im portant and critical com ponents [START_REF] Karm Eli | Classification and flow regim e analysis of d rippers[END_REF].

They generate low flow rate betw een 0.5 and 8 l/ h for pressure betw een 0.5 and 4 bars. Em itter flow rate increases w ith static pressure in a lateral pipe accord ing to a pow er-law (Karm eli 1977) [START_REF] Karm Eli | Classification and flow regim e analysis of d rippers[END_REF]:

(

w here q is the emitter flow rate , K the proportionality coefficient that characterizes each em itter, P the pressure head and x the em itter d ischarge exponent. The value of this exponent d epend s on em itter conception. Ind eed , m anufacturers try to d esign em itters w hose flow rate is not d irectly d epend ent on the pressure head (x < 0.5). To reach this goal, they introd uce som e elem ents in the em itter such as a labyrinth -channel w hich generates local pressure head losses. The m ain d raw back of this technique is clogging. There are a few stud ies w here labyrinth-channel w as investigated experim entally. Som e researchers em ploy num erical m od eling to stud y the flow in the labyrinth-channel. The num erical stud ies emphasize that em itter perform ances are related to hyd raulic conditions; [START_REF] Berkow | H yd raulic perform ance of subsurface w astew ater d rip system s. In on-site w astew ater treatm ent. Proceed ings of the N inth N ational Sym posium on Individ ual and Sm all Com m unity Sew age System s[END_REF], [START_REF] Li | CFD and d igital particle tracking to assess flow characteristics in the labyrinth flow path of a d rip irrigation em itter[END_REF], [START_REF] Paulau | H yd raulic flow behavior through an in-line em itter labyrinth using CFD techniques[END_REF]. Clogging occurs in low -velocity regions and vortex as it is experim entally found [START_REF] Jun | N um erical investigation of the clogging m echanism in labyrinth channel of the em itter[END_REF], [START_REF] Zhengying | Flow behaviour analysis and experim ental investigation for emitter m icro-channels[END_REF]. When d esigning em itters, clogging can be prevented or at least significantly red uced by d ecreasing, as m uch as possible, vortices size. Thus, anti-clogging properties and consequently m icro-irrigation efficiency are strongly related to hyd raulic perform ance. Previous num erical stud ies [START_REF] Jun | N um erical investigation of the clogging m echanism in labyrinth channel of the em itter[END_REF] and [START_REF] Li | CFD and d igital particle tracking to assess flow characteristics in the labyrinth flow path of a d rip irrigation em itter[END_REF] analyze flow in several geom etries in ord er to choose the best d esign or to optim ize an existing geom etry [START_REF] Yan | N um erical stud y on flow property in d entate path of d rip emitters[END_REF] and [START_REF] Wei | Stud y on hyd raulic perform ance of d rip em itters by com putational fluid d ynam ics[END_REF]. Depend ing on the m od el used [START_REF] Al-Muham M Ad | Mod eling a w eak turbulent flow in a narrow and w avy channel: case of m icro irrigation[END_REF], flow is strongly d ifferent. It is thus necessary to valid ate these results. In this w ork, we propose to overcom e w ith this challenge in perform ing m icro particle im age velocimetry (µ-PIV) experim ent to analyze the flow w ithin such a labyrinth-channel. This analysis is focused on the vortices. The velocity field s and profiles are presented along the labyrinth-channel in the directions x, z. The results obtained by µ-PIV are com pared to those calculated , num erically, by [START_REF] Al-Muham M Ad | Mod eling a w eak turbulent flow in a narrow and w avy channel: case of m icro irrigation[END_REF] in ord er to valid ate the num erical hypothesis and choose the best turbulence m od el ad apted to the labyrinth-channel flow .

Materials and method

Hydraulic pilot

The labyrinth-channel

Ten-pattern labyrinth-channel has been fabricated by m achining and a second plan plate is used to seal the prototype. The m aterial is transparent PMMA. Sealing is assured in pressing the tw o PMMA plates betw een tw o 5m m thick steel plates. This geom etry is the m ost w id espread in m icro-irrigation. That is w hy this first stud y has been perform ed on it. The labyrinth-channel geom etry is show n in Fig. 7. The labyrinth-channel unit length, maxim um w idth and d epth are 3 mm, 2.7 mm and 1m m respectively. The cross-section is equal to 1m m 2 .

Optic components

Optical schem e or µ-PIV equipm ent is com posed by a laser source w hich is a Litron Nd-YAG Laser of 135 MJ, d oubled in frequency (532 nm). A HiSence 4M cam era w ith the resolution of 2048 × 2048 pixels equipped w ith a Canon MP-E 65mm f/ 2.8 lens w as m ounted inside of the labyrinthchannel.

Experimental setting

µ-PIV a. Time betw een pulses and frequencies

There are several parameters to adjust before starting the µ-PIV experim ents. Firstly, tim e betw een tw o laser pulses. This tim e is d efined by flow velocity. The flow in labyrinth-channel is com plex. Accord ing to [START_REF] Al-Muham M Ad | Mod eling a w eak turbulent flow in a narrow and w avy channel: case of m icro irrigation[END_REF], it consists in tw o d ifferent regions: one is the m ain flow and the other is the vortices zones. The m odeling is carried out by k-ε m od els and RSM. The flow is supposed to be turbulent. The velocity profile is plotted in Fig. 3. This profile corresponds to the m od eling results by the standard k-ε m odel for a flow rate of 48m l/ m in. We chose tw o velocities corresponding to the peaks for estim ating . is calculated in ord er that a particle m oves a quarter of interrogation w ind ow , w hich is equal to 64 pixels (0.93 m m ). This step is repeated to calculate the tim e betw een tw o laser pulses for the three flow rates im posed d uring the experim ent (Tab 1). Tests are carried out to d eterm ine w hich tim e w ill be better for µ-PIV experim ents. The tests began from the largest tim e of each flow rate, Tab.1 (w hich correspond s to vortex velocity) to the sm allest tim e (w hich correspond s to the m ain flow velocity). It has been observed that the velocity profiles m easured do not evolve anym ore w hen the time is low er than that of the m ain flow (for exam ple 17µs for flow rate of 48m l/ m in). Therefore, the sm allest tim e has been taken. Secondly, frequency betw een tw o im age couples is calculated by taking into account that images have to be non-correlated , i.e. the particle w hich has a low velocity leaves the im age. We thus choose to fix . Third ly, im age num ber is chosen to ensure convergence criteria tow ard m ean velocity and tw o fluctuating velocity second -ord er mom ents com ponents. In these experim ents, we analyze 250 image series. The resid ue values obtained (Fig. 4) are 10 -2 and 10 -5 , for respectively m ean and fluctuating velocities.

b. Seeding particles

During the experim ent, the particles concentration is ad ded to the w ater until there are four particles by 64 × 64 pixels interrogation w ind ow . This concentration w as 10 ml solution charged w ith fluorescent particles ad d ed to 150 ml pure w ater. The particles ad d ed are polystyrene covered w ith red fluorescence. Particle d ensity is nearly equal to the w ater d ensity (1.05 kg/ m 3 ). Typical d ye w ave lengths used in red fluorescent particles are λ abs =530 nm and λ em =607 nm. These particles have a d iam eter of 5.2 µm . The choice of particle d iam eter is m ad e to ensure that it is inferior to the Kolm ogorov length scale. The Kolm ogorov length scale is calculated (Eq. 2), accord ing to the results of [START_REF] Al-Muham M Ad | Mod eling a w eak turbulent flow in a narrow and w avy channel: case of m icro irrigation[END_REF], w here the d issipation rates ε for the flow rates 24 m l/ m in (Re=400) and 48 m l/ min (Re=800) w ere 8 and 100 m 2 /s 3 respectively.

(

) 2
The Kolm ogorov length scale varies from 9 µm to 13 µm for these flow rates based on the num erical results. Stokes num ber (Eq. 3) is also a param eter d eterm ining the particles choice. This num ber is d efined as the ratio of the characteristic tim e of a particle (Eq. 4) to a characteristic tim e of the flow or the Kolm ogorov tim e scale (Eq. 5).

(3) Stokes num ber for the sam e flow rates tend s tow ard zero. This m eans that the particle behaves like a tracer. The particle seed ing image is show n in Fig. 5.

Experiment post-treatment

The µ-PIV system analysis softw are used is comm ercial softw are provided by Dantec (Dynamics Stud io). During the treatm ent, a m ask is used to elim inate the area outsid e of the labyrinth-channel flow . Then, the m ean im age is calculated for 250 im ages. An arithm etic im age is then d ed uced . Finally, an adaptive correlation is applied to the arithm etic im ages w ith final interrogation area size of 64 × 64 pixels and overlap of 50%. As an option, high accuracy subpixel refinem ent is chosen. The adaptive correlation d ata are exported to Matlab in ord er to continue and illustrate the treatm ent results (velocity field s and profiles, fluctuating velocity second -ord er m om ents, turbulence intensity).

Results and discussion

The µ-PIV has been perform ed on the ten patterns labyrinth-channel of em itter. Three flow rates have been im posed , namely, 48 m l/ m in, 36 m l/ in and 24 m l/ m in. Based on the m ean velocity and inlet cross-section, the Reynold s num bers are 400, 600 and 800 respectively.

Inlet flow

The labyrinth-channel prototype is d esigned taking into account that the flow has to be established before the entrance in the baffle curve. Velocity profiles, before the baffle curve, are plotted to d eterm ine w hether the flow w as d eveloped or not. These profiles are show n Fig. 8 for Re=800 and Re=400. The velocity profiles do not evolve. Therefore, the flow is consid ered to be established. The parabolic shape of velocity profile for Re= 400 ensures this fact. Then, d im ensionless velocities in line a (Fig. 7) for the tw o Reynold s num bers are plotted in Fig. 9. The d im ensionless velocities have the sam e tendency near the w all. The d ifference is im portant in the flow core. At Re=400, the parabolic behavior confirm s that the flow is laminar w hereas at Re=800, the curve is m ore hom ogenized w hich can be attributed to a transitional regim e or an influence of the d ow nstream turn. The Reynold s stresses are calculated. The fluctuating velocity second -ord er m om ents show that the turbulence is w eak; about 10 -4 for and . 10 -2 for (for Re=800 and Re=400) w hich is im portant com paring w ith the other term s ( and ) Fig. 10 and Fig. 11 The turbulence intensity, at the inlet, is calculated from the follow ing formula:

w here:

As only tw o com ponents of the fluctuating velocity are m easured by 2C µ-PIV, the third com ponent has to be estim ated from the first tw o. Com paring the tw o components m easured by µ-PIV, the fluctuating velocity is greater than . The ratios are 20 (30) for Re=800 (Re=400). Therefore, the fluctuating velocity perpend icular to the m ain flow axis is 20-30 tim es low er than in the m ain flow d irection. Therefore, it w as d ecid ed to calculate the turbulent intensity by assuming axisym m etric flow [START_REF] Richard | PIV stud y of flow through porous structure using refractive ind ex m atching[END_REF], i.e. the tw o com ponents perpend icular to the m ain flow direction are consid ered equal. For Re=800 (Re=400), the turbulence intensity is large near the w all about 15-35 % (20-60%). N evertheless, it is very sm all, 3-3.5% in the channel m id dle (1.5-2.5%).

Flow characterizations in the labyrinth-channel

The flow stream lines and velocity fields obtained for the fifth baffle are show n in Fig. 12 and Fig. 13 for respectively Re=800 and 400. The flow in the labyrinth-channel is com posed of a main flow w here the velocity is high (reaching 2.8 and 1.2 m / s) and vortex zones at very low velocity. The velocity in the vortex zone varies from 0.3 to 0.6 m / s for Re=800 (respectively 0.1 to 0.3 m / s for Re=400). The vortex form is not changed w hen increasing the flow rate. The vortex center m oves a little tow ard s the w all for Re=400. Fluctuating velocity second -ord er m om ent fields , and are also calculated and plotted for the fifth baffle for Re=800 and Re=400 (Fig. 14 and Fig. 15). It can be observed that the Reynold s stresses are three tim es higher for Re=800 in com parison w ith Re=400. For both, and are m axim um in the m ain flow near the separation zone. Reynolds stresses is close to zero at the vortex zones, thereby facilitating the particle deposition and em itter clogging.

Axial development of the flow

The vortex zone, as m entioned above, has the sam e form for all baffles, for both cases Re=800 and Re=400. H ow ever, the question is how the flow d evelops along the ten baffles. In ord er to answ er to this question, the velocity profiles are plotted at the beginning (third baffle), the m id d le (fifth baffle) and the end of labyrinth-channel (ninth baffle) (Fig 16). At Re=800, in the m ain flow , velocity m agnitud e is the same for all positions but in the sw irl zone, som e d ifferences are observed (34% in the vortex zone). For Re=400, a greater d ispersion is observed (40% in the vortex zone). One of the m easurem ent errors is related to the flow m eter. In fact, tests are perform ed to calculate this error. The volum e of w ater is collected for a period of tim e im posed , and then the flow rate is calculated. Com paring w ith the flow m eter values, the errors w ere 4% (8%) for Re=800 (respectively Re=400). This partially explains the d ispersion. No conclusion about establishm ent can be thus em phasized . Therefore, we focus on the velocity evolution on the first three baffles plotting the velocity com ponents profiles v x , v y and v m od ule on the first three baffles in ord er to und erstand how the flow changes from one baffle to another (Fig. 17 and Fig. 18). v x is negative in the vortex zone and positive in the m ain flow w hen v y is positive in the vortex zone and negative in the m ain flow . The flow velocity reaches the m axim um value at y=1.95m m from the low w all. In the first baffle, the m ain flow and vortex zone are not w ell clear. It is d ifficult to d istinguish them . Flow structure is m ore or less established from the third baffle but sm all evolution is observed until the tenth baffle. Regard ing dim ensionless velocity for Re=800 and Re=600 (Fig. 19), it appears that the flow is d eveloped in the third baffle. velocity second -ord er mom ents profiles confirm that third baffle is the baffle w here the flow is d eveloped . It can be conclud ed that the flow properties do not change after the third baffle, for the d ifferent flow rates.

Vertical evolution of the flow

As m entioned above, the labyrinth-channel is fabricated by m achining and thus this surface is not sm ooth. The surface roughness has not yet been characterized . The velocity profiles are plotted in Fig. 22, for the third baffle at z=[0.25; 0.5; 0.75] mm above the rough w all.. Velocity profiles for 0.25m m and 0.75m m are superim posed ; the flow is sym m etric about the z axis for this range. At this distance from the w all, the roughness has no effect on velocity.

Comparison of µ-PIV results w ith CFD modeling

The flow m od eling has been perform ed by com m ercial softw are (AN SYS/ Fluent) [START_REF] Yan | N um erical stud y on flow property in d entate path of d rip emitters[END_REF], several m od els w ere tested such as k-ε and RSM m od els. It m ust be noted that k-ε model assumes that the flow is isotropic.

is calculated within the third baffle to know w hether the flow is isotropic or not; show n in Fig. 23. In ord er to better distinguish contour values, the high scale is fixed to be 2.5. H ow ever, the m axim um value of reaches 1000 near the w all. It is d ue to the low value of . The red zones highlight the zone w here is greater than 2.5. It confirms that the flow is not isotropic. Therefore, the PIV results w ill be com pared w ith those of RSM m od el only.

Fig. 23

ratio in the third baffle, Re=800. Stream lines obtained w ith µ-PIV are com pared w ith the m od eling by RSM m od el in Fig. 24 for Re=800. This m od el d oes not pred ict the vortex zone form as µ-PIV: vortex center is not at the sam e position and velocity m agnitud e is greater in the µ-PIV results (1.6 tim es greater) Fig. 25. Regard ing the Reynold s stresses (Fig. 26) is also und er-pred icted in the m ain flow by RSM m od el. N evertheless, is the sam e ord er and is negative in the m ain flow for both cases (µ-PIV experim ents and RSM m od el results). The RSM m od el d oes not pred ict the flow in the labyrinth-channel. Therefore, the choice of another m od el or an option mod ification of this m od el is ad vised for the follow ing m od eling.

Conclusion

The m ain objective of this w ork w as to better und erstand and characterize the flow inside labyrinth. This und erstand ing helps to d esign the labyrinth-channel in ord er to avoid the em itter clogging. The flow w ithin the labyrinth-channel experim entally is investigated by µ-PIV on tenpattern repeating of em itter. The velocity field s and profiles show that the flow , w hich is com posed of a m ain flow and vortex zones, is d eveloped from the third baffle and that the flow is sym m etric about the z axis. The flow is not isotropic. Experim ental and num erical results pred icted m ore or less the sam e structure but intensities are greater from µ-PIV d ata. Further investigations are needed in ord er to predict this flow structure. Experim ents using sm aller particles in ord er to d escribe m ore precisely the flow are in construction. Afterw ard , particle tracking velocim etry (PTV) m ethod w ill be d eveloped to analyze the particle path and agglom eration processes, in interaction w ith such a com plex flow . Such a stud y is im portant in ord er to optim ize em itter conception to limit particle d eposit insid e sw irling zones.

  u le (m / s) u m Mean flow velocity at the inlet (m / s) x Exp onent u ' The root-m ean-squ are of the tu rbu lent velocity flu ctu ations (m / s) , Reynold s stresses v x Velocity com p onent in the d irection x (m / s) v y Velocity com p onent in the d irection y (m / s) x, y, z Cartesian coord inates ∆t Tim e betw een p u lses (s) λ abs , λ em Absorp tion and em ission length (m ) η Kolm ogorov length scale (m ) ν Kinem atic viscosity (m 2 / s) ε Dissip ation rate of tu rbu lent kinetic energy ( m 2 /s 3 ) τ p Characteristic tim e of a p article(s) τ k Kolm ogorov tim e scale (s)

Fig. 1

 1 Fig. 1 Prototyp e su p p ort

Fig. 2

 2 Fig.2General schem e of exp erim ents and exp erim ental set-up.

Fig. 3 Tab. 1 Fig. 4

 314 Fig. 3 Velocity p rofile for the flow m od eled in 3D by the RSM m od el, flow rate=48m l/ m in

Fig. 5

 5 Fig. 5 Particles seed ing

Fig. 6

 6 Fig. 6 velocity p rofile in the m id d le of line a in the d irection z 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4

Fig. 7

 7 Fig. 7 Labyrinth-channel

  images record ed are treated by Dantec softw are (Dynam ics Stud io). The m ean and fluctuating velocities are calculated and traced by Matlab softw are.

Fig. 8

 8 Velocity p rofiles for variou s y p ositions at the inlet.

Fig. 9

 9 Fig. 9 Dim ensionless velocity p rofiles in line a.

  .

Fig. 10

 10 Fig. 10 Flu ctu ating velocity second -ord er field s in line a, Re=800.

Fig. 11

 11 Fig. 11 Flu ctu ating velocity second -ord er field s in line a, Re=400.

Fig. 12

 12 Fig. 12 Stream lines and velocity field s in the fifth baffle, Re=800.

Fig. 13

 13 Fig. 13 Stream lines and velocity field s in the fifth baffle, Re=400.

Fig. 14

 14 Fig. 14 Flu ctu ating velocity second -ord er m om ents field s in the fifth baffle, Re=800.

Fig. 15

 15 Fig. 15 Flu ctu ating velocity second -ord er m om ents field s in the fifth baffle, Re=400.

Fig. 16

 16 Fig. 16 Velocity p rofiles in 3rd , 5 th and 9 th baffles for Re=800 and Re=400.

Fig. 17

 17 Fig. 17Velocity p rofiles for the first three baffles Re=800, at the lines 1,2 and 3 in Fig.7

Fig. 18

 18 Fig. 18Velocity p rofiles for the first three baffles Re=600, at the lines 1,2 and 3 in Fig.7

Fig. 19

 19 Fig. 19 Dim ensionless velocity p rofiles Re=600, Re=800.

Fig. 20

 20 Fig. 20Flu ctu ating velocity second -ord er m om ents p rofiles, Re=800.

Fig. 21

 21 Fig. 21 Flu ctu ating velocity second -ord er m om ents p rofiles, Re=400.

Fig. 22

 22 Fig. 22 Velocity p rofiles, Re=800.

Fig. 24

 24 Fig. 24 Velocity stream lines in the second baffle, Re=800.

Fig. 25

 25 Fig. 25 µ-PIV velocity p rofiles com p aring to RSM m od el at the line 3, Re=800.

Fig. 26

 26 Fig. 26 µ-PIV Reynold s stresses p rofiles com p aring to RSM m od el at the line 3, Re=800.

  Experimental and numerical velocity, 48 ml/min z=0.5 mm, PIV z=0.5 mm, RSM model z=0.75 mm, PIV z=0.75 mm, RSM model

Acknow ledgement:

This w ork has been funded by the Seventh Fram ew ork Program m e FP7-KBBE-2012-6: Water4Crops.


S Caro

D Chablat
email: damien.chablat@irccyn.ec-nantes.fr

P Lemoine

P Wenger

KINEMATIC ANALYSIS AND TRAJECTORY PLANNING OF THE ORTHOGLIDE 5-AXIS

The subject of this paper is about the kinematic analysis and the trajectory planning of the Orthoglide 5-axis. The Orthoglide 5-axis a five degrees of freedom parallel kinematic machine developed at IRCCyN and is made up of a hybrid architecture, namely, a three degrees of freedom translational parallel manipulator mounted in series with a two degrees of freedom parallel spherical wrist. The simpler the kinematic modeling of the Orthoglide 5-axis, the higher the maximum frequency of its control loop. Indeed, the control loop of a parallel kinematic machine should be computed with a high frequency, i.e., higher than 1.5 MHz, in order the manipulator to be able to reach high speed motions with a good accuracy. Accordingly, the direct and inverse kinematic models of the Orthoglide 5-axis, its inverse kinematic Jacobian matrix and the first derivative of the latter with respect to time are expressed in this paper. It appears that the kinematic model of the manipulator under study can be written in a quadratic form due to the hybrid architecture of the Orthoglide 5-axis. As illustrative examples, the profiles of the actuated joint angles (lengths), velocities and accelerations that are used in the control loop of the robot are traced for two test trajectories.

INTRODUCTION

Parallel kinematics machines become more and more popular in industrial applications.This growing attention is inspired by their essential advantages over serial manipulators that have already reached the dynamic performance limits. In contrast, parallel manipulators are claimed to offer better accuracy, lower mass/inertia properties, and higher structural stiffness (i.e. stiffness-to-mass ratio) [START_REF] Merlet | Parallel robots[END_REF]. These features are induced by their specific kinematic architecture, which resists to the error accumulation in kinematic chains and allows convenient actuators location close the manipulator base. Besides, the links work in parallel against the external force/torque, eliminating the cantilevertype loading and increasing the manipulator stiffness.

Unlike the Variax proposed by Gidding & Lewis in Chicago in 1994, the delta robot invented by Clavel [START_REF] Clavel | DELTA, a Fast Robot with Parallel Geometry[END_REF] has known a great success for pick and place applications. One reason for this success is the simplicity of the kinematic and dynamic models compared to the models of the Gough-Stewart platform [START_REF] Merlet | Parallel robots[END_REF]. Indeed, the performance of parallel robots may vary considerably within their workspace, which is often small compared to the volume occupied by the machine. It is noteworthy that the inverse kinematics of a parallel manipulator is usually easy to calculate when the actuated joints are prismatic joints as the corresponding equations to be solved are quadratic. However, the inverse Jacobian matrix of such manipulators is more difficult to express and its computing time higher. Several five degrees of freedom (dof) parallel manipulators have been synthesized in the literature the last few decades. However, their complexity make them difficult to build and use in general. Moreover, the use of fully parallel manipulators leads to robots with five limbs whose mutual collisions or geometric constraints reduce the workspace size. The Tripteron is one of the simplest translational parallel robot with three degrees of freedom that can be found in the literature [START_REF] Gosselin | Cartesian parallel manipulators[END_REF]. However, this architecture is not suitable for machining operations since its legs are subjected to buckling.

As a consequence, a five dof hybrid machine, named Orthoglide 5-axis, has been developed at IRCCyN. This machine is composed of three dof translational parallel manipulator, named Orthoglide 3-axis, mounted in series with a two dof spherical parallel manipulator, named Agile Eye 2-axis. The Orthoglide 3-axis has the advantages of both serial and parallel kinematic architectures such as regular workspace, homogeneous performances, good dynamic performances and stiffness. The interesting features of the Orthoglide 3-axis are large regular dextrous workspace, uniform kinetostatic performances, good compactness [START_REF] Pashkevich | Design Strategies for the Geometric Synthesis of Orthoglide-type Mechanisms[END_REF] and high stiffness [START_REF] Pashkevich | Stiffness Analysis of 3-d.o.f. Overconstrained Translational Parallel Manipulators[END_REF]. Besides, the translational and rotational motions of the end-effector (tool) are partially decoupled with the hybrid architecture of the Orthoglide 5-axis.

This paper is organized as follows. The next section deals with the kinematic modeling of the Orthoglide 5-axis. Then, some trajectories are generated using a simplified computed torque control loop and tested experimentally. Finally, some conclusions and future work are presented.

THE ORTHOGLIDE 5-AXIS A hybrid architecture

Figure 1 depicts a CAD modeling of the Orthoglide 5-axis and figure 2 shows a semi industrial prototype of the Orthoglide 5-axis located at IRCCyN. The Orthoglide 5-axis is a hybrid parallel kinematics machine (PKM) composed of a 3-dof translational parallel manipulator, the Orthoglide 3-axis, mounted in series with two dof parallel spherical manipulator, the Agile Eye 2axis. The Agile Eye 2-axis is spherical wrist developed at Laval University [START_REF] Gosselin | The agile eye: a high performance three-degree-of-freedom camera-orienting device[END_REF]. The architecture of the Orthoglide 5-axis was presented in [START_REF] Ur-Rehman | Kinematic and Dynamic Analyses of the Orthoglide 5-axis[END_REF] as well as in [START_REF] Chablat | Device for the movement and orientation of an object in space and use thereof in rapid machining[END_REF].

Orthoglide 3-Axis

The Orthoglide 3-axis is composed of three identical legs. Each leg is made up of a prismatic joint, a revolute joint, a parallelogram joint and another revolute joint. The first joint, i.e., the prismatic joint of each leg, is actuated and the end-effector is attached to the other end of each leg. Hence, the Orthoglide 3-axis is a PKM with movable foot points and constant chain lengths.

The kinematics of the Orthoglide 3-axis was defined in [START_REF] Wenger | Kinematic Analysis of a New Parallel Machine Tool: the Orthoglide[END_REF]. Let ρ ρ ρ ρ ρ ρ = [ρ 1 ρ 2 ρ 3 ] t denote the vector of linear joint variables and p = [x y z] t denote the Cartesian coordinate vector of the position of the end-effector. The loop closure of the Orthoglide 3-axis leads to the following three constraint equations:

(x -ρ 1 ) 2 + y 2 + z 2 = l 2 1
(1)

x 2 + (y -ρ 2 ) 2 + z 2 = l 2 2
(2)

x 2 + y 2 + (z -ρ 3 ) 2 = l 2 3 (3)
Therefore, the inverse Jacobian matrix J -1 O of the Orthoglide 3axis takes the form:

J -1 O =       1 - y ρ 1 -x - z ρ 1 -x - x ρ 2 -y 1 - z ρ 2 -y - x ρ 3 -z - y ρ 3 -z 1       (4)
It appears that the Orthoglide 3-axis can have up to two assembly modes, i.e., two solutions to its direct geometric model, and up to eight working modes, i.e., eight solutions to its inverse geometric model. Moreover, those solutions can be easily obtained by solving simple quadratic equations [START_REF] Pashkevich | Kinematics and Workspace Analysis of a Three-Axis Parallel Manipulator: the Orthoglide[END_REF]. The lengths of the parallelogram joints and the joint limits were obtained by using the method presented in [START_REF] Chablat | Architecture Optimization, of a 3-DOF Parallel Mechanism for Machining Applications, The Orthoglide[END_REF] in order the manipulator to get a cubeshaped translational workspace of 500 mm edge.

Agile Eye 2-axis

Figure 3 illustrates a spindle mounted on the two dof spherical parallel manipulator, the latter being mounted in series on the Orthoglide 3-axis in order to get the Orthoglide 5-axis. The Agile Eye 2-axis, namely the two dof spherical parallel manipulator, consists of a closed kinematic chain composed of five components: the proximal 1 link, the proximal 2 link, the distal link, the terminal link and the base. These five links are connected by means of revolute joints, the two revolute joints connected to the base being actuated. Let us notice that all revolute joints axes intersect. The shape of the proximal and distal links were determined in [START_REF] Chablat | A six degree of freedom Haptic device based on the Orthoglide and a hybrid Agile Eye[END_REF] in order the wrist to have a high stiffness.

The kinematics of the Agile eye was defined in [START_REF] Gosselin | The agile eye: a high performance three-degree-of-freedom camera-orienting device[END_REF]. Note that the orientation workspace of the Orthoglide 5-axis is limited to ±45 degrees due to mechanical constraints. For machining operations, the orientation of the tool is defined by two rotation angles α and β rotating respectively around the x and y axes. A IBAG spindle HFK 95.1 is mounted on the spherical wrist. Its power is equal to 1 KW and its maximum speed is equal to 42000 rpm. The distance l between the tool tip and the geometric and rotation center of the wrist should be higher than 72 mm (l > 72mm).

The geometric, kinematic and dynamic models of the Agile Eye 2-axis are given in [START_REF] Caron | Analyse et conception d'un manipulateur parallèle sphérique deux degrès de liberté pour l'orientation d'une camŕa[END_REF]. Here, its inverse Jacobian matrix is recalled in order to compute the velocity of the geometric center of the wrist with respect to the tool tip and the corresponding actuated prismatic joint rates of the Orthoglide 3-axis.

The Agile Eye admits two assembly modes that are easily obtained and discriminated by solving a quadratic equation. The Cartesian coordinate vector p tt of the tool tip is expressed as a function of angles α and β as follows:

p tt =   -sin (β ) sin (α) cos (β ) -cos (α) cos (β )   (5)
The actuated joint angles θ 1 and θ 2 are expressed as a function of angles α and β as follows:

θ 1 = -arctan -sin(α) cos(β ) cos(α) cos(β ) (6)
θ 2 = arctan sin(β ) cos(α) cos(β ) (7)
The inverse kinematic Jacobian matrix, J -1 A , of the parallel spherical manipulator is expressed as:

J -1 A =   1 0 Sθ 2 Sθ 1 Cβ Sβ Sθ 2 +Cθ 1 CβCθ 2 SβCθ 1 Sθ 2 +Cθ 2 Cβ Sβ Sθ 2 +Cθ 1 CβCθ 2   (8)
where

Cθ i = cos(θ i ), Sθ i = sin(θ i ), Cα = cos(α), Sα = sin(α), Cβ = cos(β ), Sβ = sin(β ) and θ1 θ2 = J -1 A α β (9)
The relation between the Cartesian velocities of the geometric center of the spherical wrist and the output angle rates is given as a function of the tool length l by:

  ẋ ẏ ż   = J -1 C α β (10)
with

J -1 C =        ylCαCβ -x + ρ 1 + zlSαCβ -x + ρ 1 lCβ - ylSαSβ -x + ρ 1 + zlCαSβ -x + ρ 1 -lCαCβ + zlSαCβ -y + ρ 2 - xlCβ -y + ρ 2 + lSαSβ + zlCαSβ -y + ρ 2 ylCαCβ -z + ρ 3 -lSαCβ - xlCβ -z + ρ 3 - ylSαSβ -z + ρ 3 -lCαSβ        (11)

Coupling two parallel robots

The assembly of two parallel robots in series leads to a simple kinematic model of the Orthoglide 5-axis. Indeed, the direct kinematics of the latter is obtained by solving first the direct kinematics of the Orthoglide 3-axis in order to get the Cartesian coordinates of the geometric center of the spherical wrist and then by calculating the direct kinematics of the spherical wrist in order to express the pose of the tool. About the inverse kinematics of the Orthoglide 5-axis, the inverse of the spherical wrist is first solved, then the inverse kinematics of the Orthoglide 3-axis is solved. It is noteworthy that a translational motion of the geometric center of the wrist does not change the orientation of the tool. However, a rotation of the tool about its end leads to translational displacement of the geometric center of the spherical wrist and as a result to actuated prismatic joint displacements. The inverse kinematic modeling of the Orthoglide 5-axis is expressed as:

q = J -1 t (12)
where q = θ1 θ2 ρ1 ρ1 ρ1 T and t = α β ẋ ẏ ż T . We can describe the inverse Jacobian matrix as three sub-matrices coming from the both parallel robots.

J -1 = J -1 A 0 2×3 J -1 O J -1 C J -1 O (13)
as is written in 14. Note that the first derivative of J -1 with respect to time can be easily obtained by using a symbolic computation software. However, its expression is too lengthy to be displayed in the paper.

TRAJECTORY GENERATION

The trajectory planning approach described in [START_REF] Khalil | Modeling, Identification and Control of Robots[END_REF] has been used for the trajectory generation of the Orthoglide 5-axis. Note that the trajectories are defined in the task space, i.e., the robot workspace, whereas the control loop requires data in the joint space. According the design rules presented in [START_REF] Chablat | Architecture Optimization, of a 3-DOF Parallel Mechanism for Machining Applications, The Orthoglide[END_REF], the boundaries of the Orthoglide 3-axis workspace are defined in the Cartesian space. Therefore, the direct kinematics is computed at the frequency of the control loop in order to check whether the moving-platform of the robot is within its Cartesian workspace at anytime or not.

Here, the system dynamics is decoupled in order to consider five separate actuators moving an equivalent mass and analyze the stability of the system. In the robot control scheme, the dynamic effects have been taken into account in order to improve the classical PID control scheme. Therefore, the torque Γ R associated to the rotational motions and the torque Γ T associated to the translational motions are expressed as follows:

Γ R = J θi + K P R (θ d i -θ i ) + K D R ( θ d i -θi ) + K I R t t0 (θ d i -θ i ) Γ T = M ρ j + K P T (ρ d j -ρ j ) + K D T ( ρ p j -ρ j ) + K I T t t0 (ρ d j -ρ j ) for i = 1..2, j = 1..3, J = 0.2772 Kg.m 2 , M = 91.6278 Kg, θ d i
and θ d i denote the prescribed actuated joint angles and rates. θ i and θi denote the actual joint angles and rates. The values of coefficients J and M mainly depend on the motor inertia, the gears and the axis of the ball screws. Note that the effect of the inertia and mass of the elements in motion is somehow the corresponding inertia and mass divided by the square of the gearhead ratio. As a result, for the rotational and translational parts, we obtain:

ω = 49 rad/s K P R = K P T = 3ω 2 = 19200 K V R = K V T = 3ω = 240 K I R = K I T = ω 3 = 512000 ( 15
)
where ω is a function of the torque constant, the dielectric constant, the motor efficiency and its inertia. Figure 4 shows the

K P  Robot K I K D q q + + + q d q d q d e e 

Control hardware and control loop

The control hardware of the Orthoglide 5-axis is a 1103 DSPACE card [START_REF] Dspace | dSPACE -DS1103 PPC Controller Board[END_REF] with a 933 MHz PowerPC. The actuator positions are acquired with a frequency equal to 9 kHz and a 200 Hz

J -1 =              1 0 0 0 0 Sθ 2 SαCβ Sβ Sθ 2 +CαCβCθ 2 CαSβ Sθ 2 +Cθ 2 Cβ Sβ Sθ 2 +CαCβCθ 2 0 0 0 ylCαCβ -x + ρ 1 + zlSαCβ -x + ρ 1 lCβ - ylSαSβ -x + ρ 1 + zlCαSβ -x + ρ 1 1 - y -x + ρ 1 - z -x + ρ 1 -lCαCβ + zlSαCβ -y + ρ 2 - xlCβ -y + ρ 2 + lSαSβ + zlCαSβ -y + ρ 2 - x -y + ρ 2 1 - z -y + ρ 2 ylCαCβ -z + ρ 3 -lSαCβ - xlCβ -z + ρ 3 - ylSαSβ -z + ρ 3 -lCαSβ - x -z + ρ 3 - y -z + ρ 3 1              (14)
low pass filter is used to compute the actuator velocities. The robot motions are controlled thanks to a sub program working at 1.5 KHz. This program manages the control loop, the input/output relation and the security system. For instance, if the positional error of the end-effector is larger than 3 mm, the machine will shutdown. Another level of security is on the actuators with electrical sensors stroke end. If enabled, neither the PC nor the Dspace card operates, it amounts to the cable security.

The robot trajectories are calculated using Matlab and sent to the Dspace card via an optical fiber. In the Dspace card, the inverse Jacobian matrix is used when the robot is close to its workspace boundaries in order to know whether the Cartesian displacements of the end-effector lead to exceeding joint limits or not. A pneumatic compensator is mounted along the vertical axis in order to reduce the gravity effect. This compensator continuously generates an effort equivalent to the weight of the mobile platform, but in the opposite direction. The pressure in the cylinders is controlled by a purely mechanical system.

Two types of motor are used for the actuated joints.Two Harmonic Drive FFA-20-80-H-E2048 are used to actuate the two-dof parallel spherical wrist and three Parvex Brushless NX430EAF coupled with Kinetic TDU 200 ball screws are used to actuate the three-dof translational parallel manipulator. As a consequence, the Orthoglide 5-axis can reach linear velocities up to 1 m/s and a linear accelerations up to 1 G.

Note that the maximum angular velocity of the wrist motors is equal to 3.27 rad/s and their maximum angular accelerations is equal to 270 rad/s 2 . The maximum linear velocity of the prismatic joints is equal to 1.2 m/s due to mechanical constraints and their maximum acceleration is equal to 13 m/s 2 due to the value of the maximum continuous torques.

Trajectories

Three types of trajectories can be realized by the Orthoglide 5-axis:

1. The first trajectory is a linear interpolation between two postures defined by two angles α and β and by three Cartesian coordinates x, y and z for which we can specify a velocity smaller than the maximum velocity, namely, a ratio of the maximum velocity. 2. The second trajectory is a circle defined by a center point, a radius and a percentage of the maximum speed. The orientation of the spindle is changing linearly around the circle. 3. The last trajectory is described by a G-code file. So far, only G00 and G01 functions have been implemented. The trajectory planner aims to keep the velocity between two postures constant. Some corners were added in order to smooth the path and avoid robot breaks. The radii of corner curvature are defined in order to use acceleration capacity of the machine optimally and to minimize robot speed reduction.

For all trajectories, a linear trajectory is automatically added between the current pose of the robot and its the starting point of the trajectory. Another linear trajectory between the final point of the trajectory and the current pose in order to close the loop. For safety reason, the actual speed is limited to 10% of the maximal speed during all trajectories.

For the first two types of trajectory, a fifth degree polynomial equation is used to define the position r, the velocity ṙ and the acceleration r along a linear interpolation:

r(t) = 10 t t f 3 -15 t t f 4 + 6 t t f 5 (16) ṙ(t) = 30 t 2 t 3 f -60 t 3 t 4 f + 30 t 4 t 5 f (17) r(t) = 60 t t f 3 -180 t 2 t f 4 + 120 t 3 t f 5 (18)
A computed torque control in the joint space is used. Therefore, the trajectory to be followed is projected into the robot joint space by using the inverse geometric model of the Orthoglide 5axis, the inverse kinematic Jacobian matrix J -1 and its first time derivative J-1 , in order to know whether the maximum motor velocities and accelerations are reached or not.

Linear trajectory

The computing traveling time depends on the total distance between the two robot poses. The minimum traveling time t f is defined as:

t f = max |D T | k V T , |D R | k V R (19)
where |D T | is the norm of the Cartesian motion, |D R | is the norm of the angular motion, k V T and k V R are the maximum linear and angular velocities. We can notice that these parameters are defined in the Cartesian workspace and will be larger if the maximum joint velocity and acceleration are exceeded. Let P 1 and P 2 be the vectors defining the first pose and second pose of the robot end-effector, respectively. Any intermediate pose P is defined as:

P(t) = P 1 + (P 2 -P 1 )r(t) (20)
The velocity vector V(t) is expressed as:

V(t) = ṙ(t)(P 2 -P 1 ) (21)
The acceleration vector A(t) takes the form:

A(t) = r(t)(P 2 -P 1 ) (22)
The corresponding joint coordinate vector is obtained by:

q(t) = f(P(t)) (23)
with f denotes the inverse geometric model of the Orthglide 5axis. The joint velocities are computed with the inverse kinematic Jacobian matrix J -1 as:

q(t) = J -1 V(t) (24)
The joint acceleration by using the inverse Jacobian matrix and its first derivation J-1

q(t) = J -1 A(t) + J-1 V (t) (25)
Table 1 defines a trajectory where the maximum velocity and acceleration are requested between pose P 2 and pose P 3 (1.2 m/s and 1.2 G).

Figures 5, 6 and7 show the position, velocity and acceleration profiles between P 1 to P 2 and P 2 to P 3 in the Cartesian space, respectively. Note that the trajectory between poses P 3 and P 4 is not depicted as it amounts to the one between poses P 1 and P 2 .

We can observe that the maximum velocity is reached between poses P 3 and P 4 . It means that the projection in the joint space of the maximum speed is not greater than the maximum velocity or maximum acceleration able to be produced by the ball-screws.

Figures 8 and9 show the joint positions and joint velocities

Circular trajectory

Table 2 describes the path associated with a circular trajectory, the geometric center of the circle being point P C and its radius is equal to R. P 2 and P 3 are added before and after the circular trajectory, respectively. These poses are defined by the radius of the circle, the two angles η min and η max define both the ends of the arc of circle. Thus, from pose P 1 to pose P 2 and from pose P 3 to pose P 4 , a linear trajectory is generated. The orientation of the circular path is characterized by two angles α 1 and β 1 about the z and y axes, respectively. A 2 and B 2 , (A 3 and B 3 , resp.) define the orientation of the tool tip for poses P 2 and P 3 . For the illustrative example, R = 30mm and A 2 = 20, B 2 = 0, A 3 = 0 and The local coordinate vector X l of the end-effector is defined as:

X l =   x l y l z l   =   R cos(η min + ∆ηr(t)) R sin(η min + ∆ηr(t)) 0   (26)
Its global coordinate vector X g is defined as:

X g = P C + R(α 1 )R(β 1 )X l ( 27
)
The pose of the tool is then defined as:

P(t) =       A 2 + (A 3 -A 2 )r(t) B 2 + (B 3 -B 2 )r(t) x C + x g y C + y g z C + z g       (28)
The local velocity is:

Ẋl =   ẋl ẏl żl   =   -R sin(η min + ∆ηr(t))ṙ(t)∆η R cos(η min + ∆ηr(t))ṙ(t)∆η 0   (29)
And its velocity Ẋg in the global reference frame is

Ẋg = R(α 1 )R(β 1 ) Vl ( 30
)
The velocity is then defined as

V(t) =       (A 3 -A 2 )ṙ(t) (B 3 -B 2 )ṙ(t) ẋg ẏg żg       (31)
Finally, the local acceleration is

Ẍl =   ẍl ÿl zl   =   -R((ṙ(t)∆η) 2 cos(γ) + r(t)∆η sin(γ) R(-(ṙ(t)∆η) 2 sin(γ) + r(t)∆η cos(γ) 0   (32)
with γ = η min + ∆ηr(t). And the acceleration Ẍg in the global reference frame is

Ẍg = R(α 1 )R(β 1 ) Ȧl ( 33
)
The acceleration is then defined as

A(t) =       (A 3 -A 2 )r(t) (B 3 -B 2 )r(t) ẍg ÿg zg       (34)
Figures 13 and14 show the joint positions and joint velocities for the complete circular trajectory. During the circular trajectory, the maximum joint acceleration is reached and led to a maximum speed reduction. In that case, the execution time of the complete trajectory is increased. Then, the coefficient: γ V which is the ratio of the maximum joint velocity to the speed of the associated motor and γ A which is the ratio of the maximum joint acceleration to the acceleration of the associated motor are computed. As a result, the traveling time becomes larger and the maximum Cartesian velocity decreases up to 0.6 m/s. Figure 17 gives the execution time for both trajectories. We can notice that the time is stable and close to 10 µs. This time includes the computation of the velocity, the direct kinematics and othe control loop. As we have a real-time operating system, we have also data to send to the computer in which we have the Control Desk software able to make data acquisitions. Figure 18 shows the control signal sent to the drive. The variation is between -1 and 1, which amount to the maximum continuous torques that each motor may produce. For an acceleration equal to 1.3 G, is 80% of the maximum motor torque for the 3 rd and 4 th motors. As a consequence, we expect to increase the efficiency of the control loop later on by using a part of peak torques

t f = t f max(γ V , √ γ A ) (35
t[s] 0 1 2 3 4 5 6 7 8
Δθ1, Δθ2

[rad]

×10 -3 -4 -2 0 2 4 Δθ1 Δθ2 t[s] 0 1 4 5 6 7
Δρ2, Δρ3

CONCLUSIONS

This paper dealt with the kinematic modeling of the Orthoglide 5-axis and its use for a computed torque control of the semi-industrial prototype. From the semi-industrial prototype developed at IRCCyN, we studied the use of this model for generating movements. We have shown that unlike most of the parallel kinematics machines, the kinematics of the Orthoglide 5-axis The assembly of two parallel robots leads to an inverse kinematic matrix composed of three blocks, the inverse kinematic Jacobian matrix of each robot and block corresponding to the coupling between the two parallel robots. For performing computed torque control, we have generated a set combining positions, velocities and accelerations of the actuators that are sent to the control loop. Two types of trajectories were studied and validated experimentally. We have shown that for linear motions of up to 1 m / s, the tracking error in the actuator position is lower than one millimeter. Experiments are ongoing to refine the estimate of the masses and inertia in motion and to identify the friction parameters.

FIGURE 1 .FIGURE 2 .

 12 FIGURE 1. Digital mock-up of the Orthoglide 5-Axis

FIGURE 3 .

 3 FIGURE 3. Spindle mounted on the two dof spherical wrist

FIGURE 4 .

 4 FIGURE 4. Control scheme of the Orthoglide 5-Axis

FIGURE 5 .P 3 FIGURE 6 .

 536 FIGURE 5. Wrist angle and position of the end-effector along the linear trajectory

P 2 to P 3 FIGURE 7 .FIGURE

 37 FIGURE 7. Wrist angle acceleration and acceleration of the endeffector along the linear trajectory

FIGURE 9 . 2 .P 1 P C P 4 αB 3 =

 92143 FIGURE 9. Joint velocities along the complete linear trajectory

P 2 to P 3 FIGURE 12 .

 312 FIGURE 12. Wrist angle acceleration and acceleration of the endeffector along the circular trajectory

FIGURE 13 .FIGURE 14 .

 1314 FIGURE 13. Joint values along the complete circular trajectory

FIGURE 15 .FIGURE 16 .

 1516 FIGURE 15. Joint angles and position errors along the complete linear trajectory

FIGURE 17 .FIGURE 18 .

 1718 FIGURE 17. Execution time for (a) the linear trajectory and (b) the circular trajectory

TABLE 1 .

 1 Set of Orthoglide 5-axis end-effector poses used to define three linear trajectories

									P 1	P 2		P 3		P 4
						α[degres]		0	20		0		0
						β [degres]		0	0			20		0
							x[mm]		0	140 -240		0
							y[mm]		0	130 -230		0
							z[mm]		-72 60 -180 -72
			0.4							0.4
	α, β[rad]	0.2 0.3					α β	α, β[rad]	0.2 0.3							α β
			0.1							0.1
			0	0	0.5	1	t [s] 1.5	2	2.5	3	0	0	0.1	0.2	0.3	t [s]	0.4	0.5	0.6	0.7
		0.15							0.2
	x, y, z[m]	0 0.05 0.1					x y z	x, y, z[m]	-0.1 0.1 0							x y z
		-0.05							-0.2
			0 -0.1	0.5	1	t [s] 1.5	2	2.5	3	0 -0.3	0.1	0.2	0.3	t [s]	0.4	0.5	0.6	0.7
	(a) P 1 to P 2					(b) P 2 to P 3

  )

				0.4								0.4
	α, β[rad]	0.2 0.3					α β	α, β[rad]	0.2 0.3	α β
				0.1								0.1
					0	0	0.2	0.4	t [s] 0.6	0.8	1	1.2			0	0	0.1	0.2	t [s] 0.3	0.4	0.5	0.6
			0.05							0.04
	x, y, z[m]	-0.05 0					x y z	x, y, z[m]	0 0.02	x y z
				0 -0.1	0.2	0.4	t [s] 0.6	0.8	1	1.2	0 -0.02	0.1	0.2	t [s] 0.3	0.4	0.5	0.6
	(a) P 1 to P 2				(b) P 2 to P 3
	FIGURE 10. Wrist angle and position of the end-effector along the
	circular trajectory
				0.8									2
	α, β[rad/s]		0 0.2 0.4 0.6					α β	α, β[rad/s]	-1 0 1	α β
				0 -0.2	0.2	0.4	t [s] 0.6	0.8	1	1.2			0 -2	0.1	0.2	t [s] 0.3	0.4	0.5	0.6
				0.2									1
	ẋ, ẏ, ż[m/s]	0.05 0.1 0.15					ẋ ẏ ż	ẋ, ẏ, ż[m/s]	-0.5 0 0.5	ẋ ẏ ż
					0	0	0.2	0.4	t [s] 0.6		1	1.2			0 -1	0.1	0.2	t [s] 0.3	0.4	0.5	0.6
	(a) P 1 to P 2				(b) P 2 to P 3
	FIGURE 11. Wrist angle velocity and velocity of the end-effector
	along the circular trajectory
					4
											α
		[rad/s 2 ]	0 2						β
		α, β	-2							α, β
					-4	0	0.2	0.4	t [s] 0.6	0.8	1	1.2				0	0.1	0.2	t [s] 0.3	0.4	0.5	0.6
				0.5
											ẍ
	ẍ, ÿ, z ] [m/s 2			0						ÿ z
				-0.5
						0	0.2	0.4	t [s] 0.6	0.8	1	1.2
	(a) P 1 to P 2

ACKNOWLEDGMENT

The work presented in this paper was partially funded by the Region "Pays de la Loire", the European projects NEXT ("Next Generation of Productions System") and by the FEDER Robotic Project of IRCCyN.