Yann Marco-Luc

Chevalier-Gilles Regnier

Arnaud Poitou

Induced Crystallization and Orientation of Poly( ethylene terephthalate) during Uniaxial and Biaxial Elongation

Stretching PET at a high strain rate above the glass transition temperature has a positive effect on the strength of the material. In a recent paper[!], we presented the influence of stretch and blow molding parameters on the properties of the final product, especially on the crystallinity induced by stretching. In this paper, we focus on the effects of loading, temperature, elongation and strain rate on macromolecular orientation and crystallization kinetics. We present experimental results from uniaxial and biaxial elongation tests carried out on injected PET specimens. To minimize the effect of quiescent crystallization, specimens are quickly heated with infrared lamps before the test and temperature is regulated during the test. Both uniaxial and biaxial tests are analyzed using a cross correlation technique[ 2 J that compares a picture used as reference and the picture of the deformed specimen. This technique allows us to determine all strain components at each point of the specimen, even when the strain field is not homogeneous. In a second part, we present measurements of macromolecular orientation and crystallinity ratio performed after each test. The infrared dichroism technique is used to determine the orientation of the microscopic morphology of PET before and after the testing. DSC measurements and density measurements are carried out to calculate the crystallinity ratio. Influences of strain rate, temperature and strain path sequence are evaluated in order to build a database for recent models of induced crystallization[ 3 J,[ 4 J,[SJ.

Introduction

What we present here is part of a global study of the injection blow-molding process of PET (PolyEthylene Terephthalate) bottles. PET is widely implanted in the plastic bottles market: the consumption of PET bottles increased from 2,8 million tons in 1995 to 5,3 million tons in the year 2000. The optimization of blowing process needs a good knowledge of the different stages of the process and in particular, the associated morphology evolution during stretching and blowing.

Blow molding machines make PET bottles from injected preforms. These preforms are heated up and blown up through bi-orientational -longitudinal then radial-stretching. Like most polymers PET has a low heating conductivity. Heating techniques using convection or conduction not only require a long heating time, but also lead to microstructure heterogeneity between the skin and the core of material. An alternative is radiation heating with infrared waves: solution used industrially. It leaves one face of the specimen available for image acquisition with a CCD camera.

Adjustable parameters

Figure 1: Infra red heating apparatus The heating apparatus is composed of three lamps. Each one is fed independently and the distances between them are adjusted in such a way as to heat our sample uniformly (Fig. 1 ). Since we want to have a deep penetration and a maximum absorption whatever the heating power required, we allow an adjustment range of between 65 and 100 % of the maximum voltage. The maximum variation during the test is of ±1 o compared to the imposed temperature. The deformation of the specimen during the test is captured by a CCD camera. To determine the displacement field of one image with respect to a reference image, we consider the image of a square region which from now on will be called the 'zone of interest' (ZOI). The aim of the correlation method is to match the zone of interest in both images. The displacement of one ZOI with respect to the other is reflected by the two-dimensional shift of the intensity signal which has been digitized by the CCD camera. To estimate a shift between two signals, a standard approach is to utilize a correlation function. The theoretical aspects of correlation are developed in a previous paper from Chevalier et al. [ 21 . The procedure is implemented in Matlab ™. The typical results obtained after an uniaxial test is presented in Fig. 2. The precision of the method is at least on the order of 2/1 00 pixel and the minimum detectable displacement is also in the order of 2/100 pixel.

Zone of interest

Pattern size

Patterns centres +

Uniaxial and biaxial tension tests on PET samples

In this section we present the tension tests carried out to analyze the effects of temperature, elongation speed and loading type (uniaxial, equibiaxial or sequenced biaxial elongations) on mechanical and microstructural characteristics. Uniaxial tension tests have been carried out at different temperature (from 80 to l10°C) for a wide range of tension speeds (from 0.01 to 1 s-1 ).

Each specimen is stretched from 50 up to 150 mm (A= 3, where A is the elongation) and then strain is maintained while stress relaxes. Typical results are presented in this section.

Uniaxial tension test on PET at 1 oooc

Uniaxial tension tests on PET at 25 mm/s strain rate more or less constant. The second stage is the relaxation step in which the strain remains constant. The strain field is continuously acquired by the CCD camera. We present typical tension and relaxation results and particularly the strain hardening index (SHI) and relaxation characteristics (9, crr) with respect to the strain rate and the temperature imposed during the tension test. In Fig. 3 Cauchy stress is plotted versus elongation A (i.e. A=1 +En), which shows the typical behavior where strain hardening occurs. For all tests performed above the glass transition temperature, tensile load decreases while the temperature increases. It has been admitted that strain hardening effect was due to induced crystallization but some authors[ 6 J,( 7 l point out in recent papers that macromolecular orientation should explain this strain hardening effect and that crystallization only appears during the second step of the test. As the temperature of the test increases, this strain hardening effect decreases. For example, it does not occur during the test at 1 05°C. To quantify this effect we define SHI index:

10 ~ .------------------------------, !---• 3mm/s ~ i--• 8mmfs ;/ !--- / / ' / !---• 25mrnls ,~/~"'" / !---'• 100mm/s ~:~ ~,
SHI = cr(E = 2.5) cr(E = 0.25) (1)
In Fig. 4 this index is plotted versus the strain rate for various test temperature. It can be observed that an important jump appears at the glass transition temperature. Indexes for 80°C and 85°C are higher then all others. Above 90°C the increase in temperature is associated with a decreasing index value. Moreover the strain rate has a small influence on index value. crr by minimizing the distance between experimental curves and the analytical expression given by the following simple viscoelastic model: The calculation was done for each couple of speed and temperature. Figure 5 shows the relaxation time values versus the strain rate at various temperatures. It clearly shows that near the glass transition threshold, relaxation times are low and the influence of the strain rate is minimal. At higher temperatures, the relaxation time increases and the influence of speed is greater. For a given temperature, the relaxation time decreases with the strain rate.

a(t) = CTo-(ao-err { 1-exp( -~)) (2)
The biaxial tests presented here have been carried out on a triaxial testing machine Astree (Fig. 6). This electro-hydraulic testing machine has six servohydraulic actuators. This machine has been developed by LMT-Cachan and Schenck AG, Darmstadt, Germany. A hydraulic power station generates a maximum flow of 330 limn. Closed loop control for each actuator is provided by a digital controller. The controller monitors and provides signal conditioning for each load actuator position channel. Each axis (X1, X2 and X3) has its own dedicated strain channel for signal conditioning and control. Strain input signals can be generated through a variety of strain measuring devices such as strain gages, extensometers, etc. Computer test control and data acquisition are performed by an object-oriented programming software, Lab VIEW®.

Grip

Infra red heating lamps As in the uniaxial tension test, both tests show a strain hardening effect and the combined influence of temperature and strain rate. Figure 9 presents the results of the equi-biaxial test. One can observe that the influence of speed and temperature is the same as in the uniaxial tension tests. In superposition the first step of sequential biaxial tension is plotted. It is worth noting that the elongation in the central region is lower than the final elongation in the same direction during biaxial testing. This is obviously due to poor boundary conditions in the fixed grips. We can also observe that for the same elongation, the stress is lower. This is a classical result for plane strain in regard of equibiaxial strain even if some shearing appears from boundary conditions. BiaxiaiTension Tests

~"----'"'---------=:----1• Equibiaxial 1oo•c 20 mm/s 1,2 ~ Plane strain 1 oo•c 40 mm/s 0 Equibiaxial 1oo•c 40 mm/s 1,4 1,6 1,8 2 2,2 2,4 2,6 2,8
Elongation in central region of the specimen L IS the diagonal length and e the thickness of the specimen. With the incompressibility hypothesis, we can determine e from elongation ' A. On strain-stress curves (Fig. 11) we can compare biaxial and uniaxial behaviors. This confirms that for identical conditions (same elongation, speed and temperature), uniaxial response is almost twice as weak as the biaxial one . were tested. The curves in Fig. 12 show the crystallinity ratio evolution versus time. Induction times, which one classically finds for the curves of thermal crystallization, are very small here.

This is certainly due to the low but nonzero initial crystallinity ratio in our injected specimen. The identification of k and n parameters in Avrami's model shows that with the usual value of n=3 model does not exactly fit on the whole curve (Fig. 12). In particular, we can see that experimental data present a lack of crystallization kinetics when the ultimate crystallization is nearly obtained. This effect is traditionally explained by secondary interlamelar crystallization In addition Table 1 shows the important influence of temperature on parameter k: the higher the heating temperature, the faster crystallization spreads. Hieber As a partial conclusion, this study of thermal crystallization of PET allows us to validate the experimental protocol we used by comparing our values with those in the bibliography. It clearly appears that thermal regulation using convection heating is too long to assure that no quiescent crystallization superpose flow induced crystallization during elongation tests. Our choice, infrared heating, appears to be quick enough to solve this problem. In addition, the comparison of the experimental curves with simple models of the Avrami type or Nakamura, confirms the failure of these models to take the secondary crystallization into account. Lastly, the study of the ultimate crystallinity makes it possible to prolong Verhoyen study for temperature near the glass transition.

(see for example [ 9 J,[ 1 0J,[' 1 l).
Infrared spectroscopy is based on a selective absorption of monochromatic radiation by a portion of molecule. Jasse and Koenig[ 14 l indicate that 1340 cm-1 corresponds to 'trans' conformation of CH 2 group, which can easily crystallize. 875 cm-1 and 1020 cm-1 peaks are associated with benzene ring as well as 1420 cm-1 . For that last peak, Cole et al[ 1 SJ precise that its intensity is neither sensible to crystallization, nor orientation. Infrared dichroism measurements are carried out to evaluate specimen's orientation given via the harmonic functions P2:

3 < cos 2 8 > -1 2 f = = .?z ( < cos 8 >) 2 (5)
The chosen thickness for our PET samples is 1 Of..Lm. Results have been obtained for an almost amorphous specimen, a slightly crystallized specimen (1 00°C for 90 min), a strongly crystallized specimen (130°C for 3 hours). If we focus on the dependence of intensity of 875, 975, 1020 and 1340 cm-1 peaks with respect to crystallinity (Fig. 14), it appears that 1340 peak is very sensitive to thermal crystallization. One can suppose the others will represent amorphous phase or average value between crystallized phase and amorphous phase orientation. 'Z 0,50 -......r-1020 .

.;

• 0,25 ----1 340 .!

.. .. 0,00

1 0 - .! 5 1 0 2 0 2 0 3 0 3 0 ' Q •0,25 -0,50 -------•-•----•-. . -•-.. -•--"-•-•--•--- -0,50

Deoth luml

Depth from specimen akin ( pm)

Figure 15: Macromolecular orientation in thermally crystallized PET specimen (left) and non stretched amorphous specimen (right)

Density and infrared dichroism measurements have been done on the uniaxialy and biaxialy stretched specimens in order to quantify the effect of tension conditions on induced crystallization and morphology changes. The typical result obtained from uniaxial tests (Fig. 16) is that strain rate as an increasing effect on crystallization kinetic as temperature as the opposite effect from Tg to higher values. These results are to be compared with the one related to relaxation time: as tension test temperature growths, relaxation time decreases which is representative of a macromolecular mobility. This mobility is reduced if crystallization occurs during elongation and changes the material morphology. For high temperature (20°C above Tg), relaxation is easier and elongation of the specimen is not enough to change the material morphology. Very small induced crystallinity is observed and typical relaxation time is equal to amorphous phase relaxation time. (left) Influence of thermal conditions on induced crystallization for a given tension test elongation speed (8 mm/s) (right) Influence of strain rate conditions on induced crystallization for a given test temperature (90°C)

Figure 17 shows crystallinity measurements done on samples cut off the central zone of biaxial elongation tests. Since thickness is different from uniaxial tension test specimen, we also made density measurements on non deformed specimens to confirm the low initial level of crystallinity before elongation. The first chart (top left) shows the increasing effect of strain during an equibiaxial test on induced crystallinity for identical values of strain rate, temperature. The elongation ratio corresponding to high strain is about 3 since the one corresponding to low strain is about 2. As in uniaxial elongation, if we admit that crystallinity ratio is correlated with strain hardening, it appears that during the first stage of elongation no consequent morphology change can be detected. Strain hardening and crystallinity increasing appears during the second stage of elongation: typically above "A = 2. For high strain, top-right chart in Fig. 17 shows the influence of strain rate. As speed increases, final crystallinity is higher. Of course this correlation saturates for very high tension speed since ultimate crystallinity ratio is an upper bound. In low-left chart, the effect of temperature is illustrated on two sequential biaxial tests. This result has also been observed on equibiaxial tests: as temperature increases, induced crystallinity decreases. The last (low-right) chart gives a comparison between equibiaxial and sequential biaxial tests. This final observation has to be confirmed, it is in fact very difficult to manage the same final strain field for both kinds of test. It seems that initial orientation during the first stage of the sequential biaxial test helps induced crystallization during the second stage. The sequence leads to an higher final crystallinity ratio. This confirms a posteriori the temperature homogeneity obtained by the infrared heating apparatus. Orientation can thus be represented by only one variable f for the all specimen description, since then one can study the influence of strain rate and temperature on this variable. Let us notice that orientation values obtained with 875 and 1020 peaks are quite similar, that confirms those two peaks represent the same chain element vibration. For specimens stretched at 1 05°C quasi no orientation is observed but those stretched at 90°C have an obvious orientation in the elongation direction for the same strain rate. The influence of temperature on orientation appears to follow the same rules than for crystallization ratio: the effect is more important for temperature near Tg and quickly decreases when temperature rises up to 100° or more. For the same temperature (90°C) the influence of strain rate is also similar to crystallinity ratio evolution: as elongation speed increases, as orientation is important. --.

I Tension conditions

;

./! 0,00 -- = .

. Figure 18: Influence of uniaxial testing conditions on macromolecular orientation. On top plots one can see that from 90° to 105° the final orientation decreases nearly to 0. On left plots one can see a similar influence due to decreasing of tension speed.

Once again, the conjugate effect of temperature and speed is highlighted: orientation appears for high strain rate and for low temperature (till it remains over Tg). Orientation values obtained with 1340 and 975 peaks are higher than the one calculated with 1020 or 875 cm-1 peaks. This can be explained with the partial induced crystallization of the specimen since 1340 cm-1 peak is very sensible to crystallinity ratio. For a given crystallization ratio Xc, the average orientation is obtained by relation:

From Eq. 6 it is easy to extract the orientation /c of the crystalline phase. Let's examine the 90oc specimen stretched at 25 rnrnls: the density measured showed an induced crystallization ratio up to 35%. According to Cole and al., if we make the assumption that 1020 and 875 cm• 1 peaks represent the amorphous phase orientation Cfa = 0.4) and that 1340 em•' is representative of the average orientation (f = 0.6), it appear that the crystalline phase is nearly completely oriented in the stretching direction Uc = 0,97). This result is confirmed by the orientation value given by the 975 peak which is very high also and in the same range than fc• The natural conclusion is that this peak should be representative of the crystalline phase.

Central Zone 875 970 1020 1340 " 0,7 + ------1

• @ 0,5 +-------L--'--~---j .2 0,3 +-------------1 . § 0,1 ~~5::Jiiii.i~====~~ ! § -0,1 + ij; •0,3 -1 ------------j ; § • 0,5 -'-----------~ Arm Zone (Direction 1) 875 970 1020 Cl ~ 0,9 +===========! ~ 0 , 7 -1 ---- .g 0,5 -1 ---- u ~ 0 , 3 -1 ---- g 0,1 p llil:::;= ~ -0,1 +-------1 ~ -0,3 +-------1
•.:: shows that the crystalline phase is oriented (f= 0.8) in the elongation direction 1 in the arm zone of the lower left chart. As well the arm zone of the upper right chart shows that crystalline phase is oriented in the elongation direction 2 (f= -0.3). Since the biaxial test does not lead to an homogeneous strain field, the orientation is very similar to the one measured after uniaxial tests.

0 -0,5 .L.....----~--------__,_,.t-----1 100°C
Both arms are uniaxialy stretched. One can see that the 970 cm-1 peak is the most sensitive to the small difference of nominal strain in 1 and 2 directions. This comparison remains qualitative because strain rate and temperature conditions are not identical from a test to another.

Nevertheless, when strain is higher in direction 1 it appears that the orientation function is positive. Orientation values become negative when strain in direction 2 becomes higher than strain in direction 1. For identical strain rate, testing temperature also influences the effect of a difference between direction 1 and 2 elongations on orientation function. For temperature near Tg (90°C) a small elongation difference leads to an important value of orientation function. The same elongation difference has a smaller effect when testing temperature is higher (1 00°C for example). This is one more strain rate-temperature opposite influence.

Conclusions

Blow molding process generates major mechanical properties modifications linked to thermo mechanical history of the material illustrated by induced crystallization and macromolecular orientation. A complete range of uniaxial elongation tests has been carried out and correlation between strain induced crystallization and orientation with the mechanical behavior of PET during the test has been presented.

The effects of tension speed and temperature during the test on final morphology have been highlighted. Biaxial elongation tests have been achieved using infrared heating apparatus and image correlation technique for strain field calculation during the test. Different complex path of biaxial can be managed: equibiaxial and sequential biaxial elongation tests have been presented.

Morphology measurements on biaxialy stretched samples have been presented in terms of crystallinity measured by density and orientation function measured by infrared dichroism.

First results have been discussed in regard with mechanical properties. Further work in this field will be done to improve strain induced crystallinity modeling.

Figure 2 :

 2 Figure 2: Digital image intercorrelation technique used to measure displacement field (typical uniaxial result)

Figure 3 :

 3 Figure 3: Typical strain-stress curves for PET behavior during tension tests at constant elongation speed

Figure 4 :

 4 Figure 4: Strain hardening index for given tension test temperature and speed conditions

Figure 5 :

 5 Figure 5: Asymptotic stress crr and relaxation time 8 for different tension test temperatures at a given speed

Figure 6 :Figure 6 5 A. = 3 Figure 7 :

 66537 Figure 6: Biaxial testing aparatus, triaxial testing machine Astree, biaxial specimens, IR heating apparatus

Figure 8 :

 8 Figure 8: Sequential biaxial loading, Grip displacement u, and Uz, Load F1 and F2, Elongation in the central region of specimen J...1 and J...2

Figure 9 :

 9 Figure 9: Strain-stress curves obtained from equibiaxial and plane strain (i.e., first step of biaxial sequential test) tension tests

Figure 10 :Figure 11 :

 1011 Figure 10: Typical inter-correlation result for the equibiaxial test

Figure 12 :

 12 Figure 12: Experimental and analytical crystallization kinetics for some temperatures

4 )Figure 13 :

 413 Figure 13: Ultimate crystallinity ratio dependence on temperature

Figure 14 :

 14 Figure 14: Influence of thermal crystallization: influence on peak absorbance

~

  Figure 16:(left) Influence of thermal conditions on induced crystallization for a given tension test elongation speed (8 mm/s) (right) Influence of strain rate conditions on induced crystallization for a given test temperature (90°C)

10 Figure 17 :

 1017 Figure 17: Influence of elongation, strain-rate, temperature and strain sequence on induced crystallization.

J

  : so•c -25 mm/s I l Tension conditions : 105• • 25 mm/s

3 ~ -0, 5 Figure 19 :

 3519 Figure 19: Orientation function for different positions in the non homogeneous strain field

Table 1 :

 1 k values for our PET transformed in Ak and compared with Hieber's values[tl]

	Temperature (°C)	k	A.k (min) A (min)
	100°C	0,0000006 119	125
	105°C	0,00001	46	40
	l10°C	0,00007	24	25
	l15°C	0,0002	17	10

Acknowledgement

To Eastmann who provides PET for injection.

To Pierre and Laura for English corrections.


Lennart Rubbert
email: lennart.rubbert@insa-strasbourg.fr

Isabelle Charpentier

Simon Henein

Pierre Renaud

Higher-order continuation method for the rigid-body kinematic design of compliant mechanisms

Keywords: compliant mechanism, mechanical design, continuation method, bifurcation diagram, 3-US mechanism, singularity, workspace analysis

à 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.

Introduction

The advantages of compliant mechanisms for precision applications are well known [START_REF] Howell | Handbook of Compliant Mechanisms[END_REF]. Their design can be efficiently performed using the so-called rigid body replacement (RBR) method [START_REF] Howell | Compliant Mechanisms[END_REF]. The synthesis is then achieved by relying on a rigid-body mechanism, taking into account the properties of compliant joints by spring models. Given their stiffness properties, parallel mechanisms are in this context of particular interest for the development of compliant structures. Several studies are focused on the selection of appropriate parallel architectures for the design of compliant mechanisms [START_REF] Awtar | Constraint-Based Design of Parallel Kinematic XY Flexure Mechanisms[END_REF][START_REF] Tian | Design and dynamics of a 3-DOF flexure-based parallel mechanism for micro/nano manipulation[END_REF][START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF], which is a difficult design issue for two reasons. Firstly, a great number of parallel mechanisms exist, including solutions with less than six degrees-of-freedom (DOF) that have the advantage of simple mechanical design [START_REF] Fang | Structure synthesis of a class of 4-DOF and 5-DOF parallel manipulators with identical limb structures[END_REF] and large workspace [START_REF] Merlet | Parallel robots[END_REF]. During the synthesis of a compliant mechanism using the RBR method, one difficulty is to select an adequate parallel mechanism, which requires comparing numerous designs. Secondly, compliant mechanisms usually have a limited workspace, and their synthesis therefore consists in determining the most adequate configuration in the rigid-body mechanism workspace. The determination of a relevant mechanism configuration is therefore also an issue, since parallel mechanisms can have a complex workspace, and kinematic properties can largely vary within that workspace with presence of kinematic singularities. Efficient kinematic analysis with determination of the presence of singularities is needed.

Several studies have been focused on the determination of mechanism workspace. Monte Carlo method [START_REF] Rasetgar | Generation of manipulator workspace boundary geometry using the monte carlo method and interactive computer graphics[END_REF], branch-and-prune algorithm [START_REF] Bohigas | A complete method for workspace boundary determination on general structure manipulators[END_REF] as well as geometrical approaches [START_REF] Ahmad | Analysis of kinematic singularities for robot manipulators in Cartesian coordinate parameters[END_REF][START_REF] Lai | A new method for the singularity analysis of simple six-link manipulators[END_REF] have been proposed. These methods exploit a priori knowledge of the mechanism properties to reduce computation time, which constitutes a serious drawback during a mechanism synthesis when several designs have to be compared. It has been demonstrated that the workspace boundaries corresponds to singularity locii [START_REF] Litvin | Application of theorem of implicit function system existence for analysis and synthesis of linkages[END_REF]. Haug et al. [START_REF] Haug | Numerical algorithms for mapping boundaries of manipulator workspaces[END_REF] and recently Hentz et al. [START_REF] Hentz | Higher-order continuation for the determination of robot workspace boundaries[END_REF][START_REF] Boehler | Definition and Computation of Tensegrity Mechanism Workspace[END_REF] have consequently developed continuation-based methods to determine the workspace outer envelope by following onedimension paths on it. Any mechanism can be analyzed after simply deriving the geometrical closed-loop equations, and with a limited computation time. This constitutes very interesting features for a systematic evaluation of mechanism, but no information is obtained on accessible configurations within the workspace. In this paper, we introduce an evaluation method based on the same theoretical background but with a complementary approach giving information on the boundary and the inner properties of the workspace.

Kinematic analysis to study the geometrical model and properties of robots have been developed using different methods including analytical [START_REF] Khalil | General solution for the dynamic modeling of parallel robots[END_REF][START_REF] Briot | Self-Motions of General 3-RPR Planar Parallel Robots[END_REF] and geometrical approaches [START_REF] Yu | Geometric approach to the accuracy analysis of a class of 3-DOF planar parallel robots[END_REF][START_REF] Merlet | Parallel robots[END_REF]. Both have their limitations since obtaining an analytical model is not always possible and the geometrical approach efficiency strongly depends on the mechanism and the user expertise. Moreover, considering the analysis of robot singular configurations, analytical approaches generally fail. For these particular situations, Grassman-Cayley algebra and screw theory can be used to prove the existence of singularity and to identify the nature of singular movements [START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF]. However, these methods lack generality and require advanced knowledge on the mechanism as well as designer expertise to determine the mechanism kinematic properties. This is an issue, since singularities are either explicitly avoided or used for the synthesis of compliant structures as detailed in [START_REF] Rubbert | Using singularities of parallel manipulators to enhance the rigid-body replacement design method of compliant mechanisms[END_REF].

A first main contribution of the present paper is a method identifying the presence of singularities for a mechanism along any given set of end-effector configurations, at no extra computation time and with no loss of generality. A second main contribution is the ability that offers the method to determine workspace and kinematic properties of mechanisms that can be employed using only basic knowledge of a mechanism. It builds on the higher-order continuation method [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée: l'approche Diamant[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] for the continuous computation of joint angular positions along any parameterized end-effector path. The results are presented in the form of so-called bifurcation diagrams displaying joint positions as functions of the end-effector position. A third main contribution is that precise information on the mechanism kinematic properties can be extracted, helping the designer in choosing a mechanism, its configuration, the actuated joints or the joints to be suppressed in the particular situations of singularities.

The method is exemplified by a study of the kinematic behavior of the 3-US mechanism. Composed of three identical legs connected to the end-effector (Fig. 1), this mechanism is of simple architecture and can be of great interest for compliant mechanisms. Each leg is linked to the base with a universal joint (U) and to the end-effector with a spherical joint (S). The mechanism is also interesting because of the motions it can exhibit when the endeffector lies in the plane of the base (Fig. 1, left). As outlined in [START_REF] Rubbert | Using singularities of parallel manipulators to enhance the rigid-body replacement design method of compliant mechanisms[END_REF][START_REF] Qin | Forward displacement analysis of two foldable 3-US parallel mechanisms[END_REF][START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF], the end-effector then exhibits two rotations around axes in the base plane and a translation in the direction perpendicular to the plane. Such motions have been successfully used in the surgical context for a compliant device [START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF]. Beyond this use, very little information exists on the properties of the mechanism. It is therefore a case study that is suited for the introduction of the proposed approach.

The paper is structured as follows: the 3-US kinematic modeling is introduced in Section 2 to outline the available information on the mechanism. This sets the general formulation used in the method. The higher-order continuation method is introduced in Section 3: it is the core of the general analysis. We propose to use the Diamanlab software [START_REF] Charpentier | Diamanlab -an interactive taylor-based continuation tool in matlab[END_REF] and describe its main underlying principles. The analysis of the mechanism properties, including the workspace boundaries and interesting features such as screw motion and self-motion are discussed in Section 4. A summary and an outlook are proposed in Section 5.

Case study description: the 3-US mechanism

Introduction and parameterization

We consider the parameterization of the 3-US mechanism as given in Fig. 2. It is designed with three identical legs equally distributed along a circle of radius R and center O. Each leg i, i ∈ {1, 2, 3}, is composed of a universal joint rigidly connected to a spherical joint. The universal joint center A i is located at the intersection of the axes of its two revolute joints. The axes of the first revolute joint of the universal joints intersect at O. The base frame is chosen with O as center, vector x aligned with OA 1 and vector y in the base plane. For symmetry reasons, the points A i are positioned with angles ∠A 1 OA 2 = ∠A 2 OA 3 = ∠A 3 OA 1 = 2π/3. This joint arrangement is of particular interest as it allows for the deployment of the structure [START_REF] Qin | Forward displacement analysis of two foldable 3-US parallel mechanisms[END_REF]. The single plane containing O and the points A i is designated as the base plane in the following.

The center of the spherical joint of leg i is denoted by B i . The centroid of the triangular platform B 1 B 2 B 3 is denoted by G. The configuration of the spherical joint centers is such that GB i = OA i , i ∈ {1, 2, 3}, so that the mechanism architecture is defined only by two geometrical parameters: the radius R = OA i = GB i and the angle α = ∠A i OB i , i ∈ {1, 2, 3}.

The angular positions of the first and second revolute joints of the universal joints are respectively defined by θ 1i and θ 2i , i ∈ {1, 2, 3}. In the planar configuration depicted in Fig. 2, the angular positions are defined such that θ 1i = 0 and θ 2i = 0, i ∈ {1, 2, 3}. From a technological point of view, actuating the spherical joints is difficult. As the 3-US mechanism has three DOF, we consider that either the first or the second revolute joints are actuated. Both situations and their consequences are numerically studied in section 4. The input and intermediate coordinates are consolidated in the state vector designated by u = (θ 11 , θ 21 , θ 12 , θ 22 , θ 13 , θ 23 ).

The platform motions are rotational and translational in most of the cases [START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF]. Since the mechanism has three DOF, one single set of three variables among the input and output coordinates is sufficient to constrain the mechanism. The coordinates of the platform centroid G is chosen as output coordinate set in order to define the mechanism's configuration and is denoted by X.

Kinematic model

By definition, the output coordinate set X satisfies:

X = OG = 1 3 3 i=1 OB i , (1)
and can be related to the input and intermediate coordinates by expressing the vectors OB i , i ∈ {1, 2, 3}, as:

OB i =             (R -L 2 cos (θ 2i )) cos 2π 3 (i -1) -(L 1 cos (θ 1i ) -L 2 sin (θ 2i ) sin (θ 1i )) sin 2π 3 (i -1) (L 1 cos (θ 1i ) -L 2 sin (θ 2i ) sin (θ 1i )) cos 2π 3 (i -1) + (R -L 2 cos (θ 2i )) sin 2π 3 (i -1) L 1 sin (θ 1i ) + L 2 sin (θ 2i ) cos (θ 1i )             , (2)
where

L 1 = R sin (α), L 2 = 2R sin (α/2) 2 .
Combining equations ( 1) and ( 2) yields three constraint equations linking the output coordinates X to the state vector u.

Besides, the mobile platform rigidly connects the points B i , that is:

B 1 B 2 = B 2 B 3 = B 1 B 3 = √ 3R. (3)
Expanding equation (3) using the vectors expressions from equations (2) gives three additional constraint equations. The resulting system of six constraint equations for six unknowns is written explicitly in the code presented in Appendix A.

The Inverse Kinematic Model (IKM) and the Direct Kinematic Model (DKM) are generally derived from constraint equations to determine the input-output relation of a mechanism, i.e. the geometrical relation between the end-effector and the actuator positions. Specifically the IKM allows computing the actuated joint coordinates for given output coordinates X. Reciprocally, the DKM allows computing X from the knowledge of the actuated joint coordinates. Even if the mathematical formulation of the constraint equations seems affordable, to the best of our knowledge, the derivation of the IKM or DKM for this mechanism is not achievable in an analytical manner in the general case. Their numerical computation could be considered with for instance a classical Newton-Raphson method, but such use remains difficult in the vicinity of singular configurations. This constitutes a strong limitation since the vicinity of singular configurations is particularly interesting for the design of compliant mechanisms: indeed it corresponds to particular kinematic behavior which can be used notably for motion amplification purpose for example [START_REF] Bachta | The Cardiolock Project: Design of an Active Cardiac Stabilizer for Cardiac Surgery[END_REF][START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF]. As outlined in the following, the proposed higher-order continuation method is particularly relevant and accurate in such situations.

Mechanism analysis using higher-order continuation method

The workspace exploration of a mechanism can be performed along parameterized paths. For instance, one of the three output coordinates related to the end-effector centroid may be considered as a varying parameter, hereafter denoted λ, while the other two are kept constant.

One deduces from the constraint equations ( 2) and ( 3) that the mechanism analysis can be formulated as the general nonlinear parametric problem linking the input and intermediate coordinates u ∈ R n (n = 6 here)

R(u, λ) = 0, (4)
where the parameter λ ∈ I R is to be varied along the path and R is a the vector-valued function. The solutions for the under-determined system (4) define one-dimensional continua, hereafter called branches of solutions. Continuation methods and related software, see [START_REF] Govaerts | Interactive continuation tools[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] enable to compute solution branches. Most of them rely on first order predictor-corrector principles described in [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF][START_REF] Allgower | Numerical Continuation Methods: an Introduction[END_REF][START_REF] Doedel | Numerical analysis and control of bifurcation problems (i) bifurcation in finite dimensions[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] and approximate the solution branches as a collection of points satisfying (4). Under analyticity assumptions, continuous solution branches may be actually computed [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF].

Higher-order continuation method

Under analyticity assumptions, the solution for nonlinear problems may indeed be carried out through a Taylorbased numerical method and higher-order Automatic Differentiation (AD) [START_REF] Chang | Solving ordinary differential equations using Taylor series[END_REF][START_REF] Guillaume | Computation of high order derivatives in optimal shape design[END_REF][START_REF] Berz | COSY INFINITY and its applications in nonlinear dynamics[END_REF][START_REF] Nedialkov | Solving differential-algebraic equations by Taylor series (I): Computing Taylor coefficients[END_REF]. In a nutshell, AD is a set of techniques that enable to augment a given computer code with derivative computations, with an obvious advantage in terms of ease of implementation and generality since derivative computations, Jacobians and higher-order terms are performed in an automatic manner [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF][START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF][START_REF] Charpentier | Fast higher-order derivative tensors with Rapsodia[END_REF].

The Diamant approach [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée: l'approche Diamant[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] is a higher-order continuation method that automates the Asymptotic Numerical Method discussed in [START_REF] Thompson | The nonlinear perturbation analysis of discrete structural systems[END_REF][START_REF] Cochelin | Méthode asymptotique numérique[END_REF]. Diamant is of particular interest for mechanism analysis since it may be used from the constraint equations only. The related Matlab software is named Diamanlab [START_REF] Charpentier | Diamanlab -an interactive taylor-based continuation tool in matlab[END_REF][START_REF] Charpentier | Diamanlab v1.0: An interactive Taylor-based continuation tool in MATLAB, Software protected by Agence de Protection des Programmes[END_REF].

Following [START_REF] Cochelin | Méthode asymptotique numérique[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF], the solution for ( 4) is performed considering the nonlinear problem:

R(u(a), λ(a)) = 0, (5)
a = u (a) -u(0), ∂u ∂a (0) + (λ(a) -λ(0)) ∂λ ∂a (0) , ( 6
)
where a is a measure of the pseudo arc-length (equation ( 6)) [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF] along the branch of solutions (u(a), λ(a)) issued from the point (u(0), λ(0)). This branch is locally approximated by means of two Taylor series:

(u(a), λ(a)) = K k=0 a k k! ∂ k u ∂a k (0), K k=0 a k k! ∂ k λ ∂a k (0) = K k=0 a k u k , K k=0 a k λ k , (7)
where u k and λ k are the Taylor coefficients of u(0) and λ(0) at order k, and K is the truncation order. Thanks to analyticity, the generic problem [START_REF] Rubbert | Design and Optimization for a Cardiac Stabilizer Based on Planar Parallel Compliant Mechanisms[END_REF] written in terms of the Taylor coefficients

R k of R R(u(a), λ(a)) = K k=0 a k R k = 0, (8)
allows to deduce that:

R k = 0, for k = 1 to K. ( 9
)
The Taylor coefficient R k (k = 1, .., K) may be split into [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II: Homotopy and complex nonlinear eigenvalue method[END_REF]

R k = {R 1 }(u k , λ k ) T + {R k|u k =0,λ k =0 } = 0, ( 10
)
where the Jacobian {R 1 } of R is the same over the orders. The higher-order term {R k|u k =0,λ k =0 } contains the nonlinear contributions that are independent from the unknowns u k and λ k . The higher-order terms and the solutions of linear system [START_REF] Yu | Geometric approach to the accuracy analysis of a class of 3-DOF planar parallel robots[END_REF] are computed each in turn following:

{R 1|u 1 =Id,λ 1 =1 }(u k , λ k ) T = -{R k|u k =0,λ k =0 }, for k = 1 to K, a = (u k -u 0 )u 1 + (λ k -λ 0 )λ 1 . (11)
Once the series are computed, the a posteriori error estimate [START_REF] Rubbert | Using singularities of parallel manipulators to enhance the rigid-body replacement design method of compliant mechanisms[END_REF] provides an adaptive step size control based on a numerical radius of convergence a * of the series [START_REF] Cochelin | Méthode asymptotique numérique[END_REF]

a * = ε ||U 1 || ||U K || 1 K-1 , ( 12
)
where ε is a small parameter usually chosen equal to 10 -6 . It also enables to evaluate the interest for a correction. Truncated Taylor series allow for the construction of accurate continuous solution branches and contain information for the detection of bifurcation and instability points [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF]. Bifurcation diagrams plot continuous solution branches comprising possible bifurcation points. These bifurcation points correspond to the intersection of branches of solutions. They may be detected either by monitoring a change of sign for the Jacobian determinant along the solution branch [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF] or by implementing higher-order bifurcation indicators [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF][START_REF] Guevel | Automatic detection and branch switching methods for steady bifurcation in fluid mechanics[END_REF] from Taylor series information. For the sake of generality, see subsection 3.2, Diamanlab (Matlab implementation of Diamant [START_REF] Charpentier | Diamanlab -an interactive taylor-based continuation tool in matlab[END_REF][START_REF] Charpentier | Diamanlab v1.0: An interactive Taylor-based continuation tool in MATLAB, Software protected by Agence de Protection des Programmes[END_REF]) implements the detection and branch switching processes proposed in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF].

Key points for the automation of the solution (5) are the generality of both the high-level nonlinear solver [START_REF] Yu | Geometric approach to the accuracy analysis of a class of 3-DOF planar parallel robots[END_REF] and AD for Taylor coefficient computations. In other words, only the constraint equations are necessary for the mechanism analysis.

Implementation with Diamanlab

Diamanlab [START_REF] Charpentier | Diamanlab v1.0: An interactive Taylor-based continuation tool in MATLAB, Software protected by Agence de Protection des Programmes[END_REF] builds on the general high-level nonlinear solver [START_REF] Yu | Geometric approach to the accuracy analysis of a class of 3-DOF planar parallel robots[END_REF] for the continuous solution of the general nonlinear problem [START_REF] Rubbert | Design and Optimization for a Cardiac Stabilizer Based on Planar Parallel Compliant Mechanisms[END_REF]. It is therefore of particular interest for mechanism analysis near to singular configurations.

From a user point of view, a nonlinear problem satisfying the analyticity assumption and written as in ( 5) may be solved by deriving a user class from the general system class of Diamanlab and by implementing the mechanism equations ( 1)-( 2) and (3) as a user method R. Additional user-defined graphics may be coded within the user class. In the following, each exploration is performed with a truncation order of K = 20, an threshold of ε = 1.10 -6 in ( 12) and a Newton-Raphson threshold of 1.10 -6 . Diamanlab moreover provides the detection of bifurcations [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF], exploited here to detect the mechanism singularities. A graphical user interface allows for an interactive management of the continuation.

Diamanlab is available for academic research. The software includes a benchmark [START_REF] Charpentier | Diamanlab -an interactive taylor-based continuation tool in matlab[END_REF] of classical nonlinear problems as a basis for the implementation of user-defined problems as well as the 3-US mechanism problem. The interested readers may contact the second author for a copy of the Diamanlab software including the 3-US example, see Appendix A.

Analysis of the 3-US mechanism properties for compliant mechanism design

The benefits for the design of compliant mechanisms coming from the use of bifurcation diagrams are illustrated considering two meaningful singular situations. The screw motion of the platform obtained when its centroid is displaced along the axis (O, z) is first presented in subsection 4.2 to emphasize the ability of the proposed method to detect parallel and serial singularities. The kinematic behavior is similar to the one exhibited by the Wren mechanism [START_REF] Vertechy | Static and stiffness analyses of a class of over-constrained parallel manipulators with legs of type us and ups[END_REF][START_REF] Andersen | Principles and preliminary concepts for compliant mechanically reactive armor[END_REF]. Second a self-motion due to architectural singularity, with a circular trajectory which has not been presented in the literature before, is described in subsection 4.3. Finally, in subsection 4.4 the workspace boundary is obtained from the bifurcation plots local and global extremum and displayed.

Preliminary observations

The radius R and the arc length α of the legs are the two geometrical parameters that define the mechanism behavior. The radius R only affects linearly the size of the workspace. It is therefore arbitrarily set to R = 10. The angle α influences the volume and the shape of the workspace as well as the kinematic behavior. The angle α is chosen equal to π 2 in order to avoid leg overlapping, which would strongly reduce the practical interest of the mechanism. Cylindrical coordinates are adopted with O as origin and (O, z) as axis in accordance with the architecture of the mechanism. The position of the centroid G is denoted by (ρ G , φ G , z G ) with ρ G the radius, φ G the angular position and z G the height. Taking into account the symmetry of the mechanism with respect to the (O, x, y) plane, the bifurcation diagrams are plotted for z G ≥ 0 only.

Screw motion, bifurcation and kinematic singularities

The proposed analysis is based on the information that can be extracted from bifurcation diagrams. These diagrams are representations of the variable set u(λ) that plot the continuous evolution of u along the variation of λ. The initial point (u(0), λ(0)) is computed from an initial guess using a Newton-Raphson algorithm. Then the continuation takes place and the solution branches of the nonlinear problem under study are computed as sets of Taylor-based solution branches stored in checkpoints [START_REF] Griewank | Algorithm 799: Revolve: An implementation of checkpoint for the reverse or adjoint mode of computational differentiation[END_REF][START_REF] Charpentier | Checkpointing schemes of adjoint codes: Application to the meteorological model Meso-NH[END_REF]. These checkpoints contain information about the initial point (u(0), λ(0)), the series and its range of validity, and one (or several) oriented tangent vector that indicates the continuation direction. A branch of solutions is constituted by the evolution of the septuplet (u, λ) = (θ 11 , θ 21 , θ 12 , θ 22 , θ 13 , θ 23 , λ). A projected bifurcation diagram is a 2D representation of one or several branches of solutions.

The bifurcation diagram presented in Fig. 3a plots the branches of solutions when the position of G varies along the axis (O, z). The end-effector motion path is set by (ρ G , φ G , z G ) = (0, 0, λ). For the sake of readability, curves for (θ 12 , θ 22 ) and respectively (θ 13 , θ 23 ), are shifted by 2π and 4π. From Fig. 3a, it is possible to compute the configuration of all the mechanism components. The representations of the mechanism for situations A to E (Fig. 3a) are given in Fig. 3c. Square symbols in Fig. 3a indicate values of θ 1i and θ 2i , i ∈ {1, 2, 3} that belong to the same septuplet. Arrows indicate how the joint values vary with λ for this situation.

For this situation, analytical relationships between the end-effector position and the joint values can be established after a specific analytical derivation, as detailed in Appendix B. The corresponding representation of the bifurcation diagram is presented in Fig. 3b for sake of comparison. The accuracy of the continuation-based analysis is confirmed: Fig. 3b corresponds to Fig. 3a after inclusion of the symmetry with respect to z G = 0.

A bifurcation point is identified when the end-effector is in the base plane, which happens for the two situations (θ 1i , θ 2i ) = (0, 0) and (θ 1i , θ 2i ) = (π, 0), i ∈ {1, 2, 3}. They are indicated by red circles in Fig. 3a. When leaving the configuration A or E (Fig. 3c) to reach non-planar configurations, it is possible to get two distinct pairs of values of θ 1i and θ 2i , i ∈ {1, 2, 3} for one pose of the end-effector. This indicates the presence of two configurations for the legs for a single end-effector pose. This can be observed in situations B, C and D in Fig. 3c where the legs 1 and 2 are in one configuration and the leg 3 is in the other configuration. The branches describing the two different leg configurations do not exhibit same shape, nor the same slope at the origin z G = 0. In other words, the V-profile of the curves representing the evolution of θ 2i , i ∈ {1, 2, 3}, differs from the variation of θ 1i , i ∈ {1, 2, 3}. One can see here that the proposed evaluation method allows the designer to easily apprehend such differences, without any linearization of the mechanism kinematics for local analysis as usually performed. It can be for instance beneficial to choose as active joints the ones with lowest slopes at the origin since the actuator ranges of motion will be lowered for a given end-effector displacement.

When the joints θ 1i , i ∈ {1, 2, 3} are actuated, we observe at points A and E that the curves θ 1i (z G ) present a vertical tangent: the platform can move locally without any variation of θ 1i , i ∈ {1, 2, 3}. This signifies that the mechanism undergoes a parallel singularity, which is in accordance with the results obtained in a previous work [START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF]. This means that we can visualize the presence of a parallel singularity from the bifurcation diagram. The detection is possible if the path used to build the diagram corresponds to the instantaneous motion gained by the mechanism end-effector, which is here the case with an out-of-plane translation [START_REF] Rubbert | A planar RRP compliant mechanism based on the singularity analysis of a 3-US parallel mechanism[END_REF]. Singular configurations are interesting as they can be candidate configurations for singularity-based compliant mechanism design as the one developed in [START_REF] Rubbert | Using singularities of parallel manipulators to enhance the rigid-body replacement design method of compliant mechanisms[END_REF]. In addition, in the vicinity of such a singular configuration, a high amplification ratio can be obtained for θ 1i when the branches of solutions are those presenting vertical tangents. This is of particular interest when piezoelectric actuators requiring mechanical amplification are needed, as it is the case in many compliant mechanical amplifiers.

When the end-effector is at the maximum height (configuration C in Fig. 3a and Fig. 3b), a so-called turning point is reached, i.e. an extreme point. The absence of infinitesimal displacement δz G about that configuration implies that the end-effector does not exhibit any translational motion. No end-effector displacement is obtained while the variables describing the leg joints vary. Within the same exploration we here identify a serial singularity corresponding to a workspace limit. Such situation can also be employed for a compliant mechanism as in its vicinity it can be used as a way to increase the mechanism resolution.

Multiple bifurcations, architectural singularity and self-motion

The bifurcation diagram presented in Fig 4a plots the branches of solutions when the centroid G is driven along a circle of radius R at height z G = R, that is by considering the path (ρ G , φ G , z G ) = (10, λ, 10). The diagram is composed of two curves for each of the six components (θ 11 , θ 21 , θ 12 , θ 22 , θ 13 , θ 23 ) that are part of the state vector u. In Fig. 4a, the plots of (θ 12 , θ 22 ) and (θ 13 , θ 23 ) are again shifted for the sake of clarity. For each joint variable, it is possible to switch between two paths at bifurcation points represented by red circles. The bifurcations do not occur for a single value of the angle φ G but with 2π/3 shifts: six bifurcations are detected along the end-effector path.

!" !"#$ % !"#$ !" $!" &!" '!" $ % $ & ' ( )% )$ )& )' ✁ z= Bifurcation diagram at point (✂,✄)=(0 ,0mm), Rayon=100mm, =90 * ☎ )) ☎ $) ☎ )$ +$✆ ☎ $$ +$✆ ☎ ), +&✆ ☎ $, +&✆ B D C B A E A E D C -! -!/2 0 !/2 ! 2! 4! 6! =" B D B D A E E A A C C A E E B B D C C D z G = λ [rad] (a) Bifurcation diagram. -π/2 0 π/2 π -16 -14 -12 -10 -8 -6 -4 -
One solution branch corresponds to θ 1i = π/2, i ∈ {1, 2, 3} (vertical curves on Fig. 4a) with simultaneous linear variation of the angles θ 2i , i ∈ {1, 2, 3}. As outlined by Fig. 4b, the end-effector exhibits a circular movement such that (O, z) remains normal to the end-effector plane. Since the first joint angles are constant, it means the end-effector can move in this situation even if the first joints on the base are maintained in their initial position. We can here observe the existence of an architectural singularity [START_REF] Ma | Architecture singularities of platform manipulators[END_REF], also called a self-motion. The continuation method allows us to observe this particular motion using only the closed-loop equations. In addition, we can see that it is also possible for the end-effector to describe a motion of circular shape without being permanently in singularity (see configurations E-F on Fig. 4b).

The three different joint evolutions observed during the platform motion have basically constant, linear or sinusoidal shapes. Once again, the bifurcation diagram gives the designer immediate and precise information on the joint angle evolutions so he can pick-up the most suitable configurations of the mechanism for the design of a compliant mechanism. Here, it is possible to obtain end-effector displacement with no velocity of θ 1i . This means the joints θ 1i could be removed and thus the compliant mechanism architecture could be greatly simplified without modification of the desired end-effector motion. At the same time, linear or sinusoidal evolutions of the other joint angles can be exploited. The designer can make his decision depending on the most suitable input-output relationship for the application. In particular, there is a factor two in the total range of motion on the joints (see Fig. 4a) between the two situations for one gyration of the platform, which may be a determinant indicator.

Finally, one can observe that the bifurcation points allow the legs to switch from one profile to the other, for example from a sinusoidal profile to a constant profile when considering the evolution of θ 1i . If the joints θ 1i are actuated and the mechanism passes through a bifurcation point, there is a chance to lose control of the platform. In this situation, the bifurcation points are probably configurations to be avoided by the designer.

Workspace estimation

The reachable workspace is the set of all positions that can be reached by the point G for at least one orientation of the mobile platform. The shape of the workspace can help the designer to choose a design configuration regarding the integration constraints. It is here obtained by extracting the local and global boundaries for each bifurcation diagram obtained for a discrete map of positions (φ G , ρ G ) with a continuous exploration of z G . We therefore obtain a semi-discrete volume approximating the maximal workspace.

Workspace estimation can be simplified by taking into account mechanism symmetries. The mechanism is invariant with respect to a rotation of 2π 3 radians around the axis (O, z). The workspace of each leg end is also symmetric with respect to the plane (O, OA i , z) (Fig. 1). As said before, the (O, x, y) plane is a symmetry plane for the mechanism. The workspace can thus be computed from one twelfth defined by 0 ≤ φ G ≤ π 3 and z G ≥ 0 and completed by circular pattern repetition and symmetry.

In order to efficiently build the semi-discrete volume, an adaptive technique inspired from [START_REF] Charpentier | Exact derivatives of the pair-correlation function of simple liquids using the linear tangent method: application to self-consistent integral equations[END_REF] is adopted. It is based on a Newton-Raphson (NR) method to find an initial approximation at a given (φ G , ρ G ). Then we associate Diamanlab for the computation of branches of solutions along the path λ = z G to an iterative management of the 2D grid construction. This construction encloses an angular do loop in the range [0, π 3 ] with a step of π/90 radians (2 • ), and an adaptive loop performed on the radii ρ. This runs from ρ = 1 to ρ = L with a step of δρ = 1. The checkpoints computed for the current (φ G , ρ G ) are used as initial guesses for the NR method (until the NR convergence). If the NR method fails, the step δρ is divided by two and the NR method is run again. The loop on ρ stops when δρ < 0.01. This allows the building of the upper part of Fig. 5. The lower part is built using the same process for φ G ∈ [0, π 3 ] and ρ G running from 10 down to 0. The lower checkpoints of the solution branches computed for ρ G = 10 are used as initial approximations. The surfaces that delimit the reachable workspace2 are represented in Fig. 6. The reachable workspace is mainly composed of two parts, which can be seen on Fig. 5a as upper and lower volumes. The upper volume can be described as the volume existing between a portion of the sphere of radius L (superimposed on Fig. 6) for z G ∈ [L 1 ; L] = [10; 10 √ 2] and a surface that is going through the center O. The funnel shape along the vertical z-axis indicate restricted displacement along x-and y-axis which can help the designer to make a decision of considering or not this workspace area. The lower volume is a set of end-effector positions that can be reached with flipped configurations, a priori complex to integrate and probably of lesser interest for the designer. Further analysis shows that this lower volume exists only for α ∈ [π/3, π/3 + π]. If α is not in this interval, only the upper volume remains.

0 pi/2 pi 2pi 3pi 7pi/2 4pi 5pi 6pi 13pi/2 7pi 8pi 9pi 0 pi/3 2pi/3 pi 4pi/3 5pi/3 2pi ! (rad) "/2 " 2" 4" 13"/2 # G =$ 3" 7"/2 5" 6" 7" 8" 9" "/3 2"/3 " 4"/3 5"/3 2" 0 0 A A A A A A B B B B B B C C C C C C E E E E E E D D D D D D F F F F F F ! ! 11 ! 21 ! 12 +2& ! 22 +2& ! 13 +4& ! 23 +4& 3
The shape of the workspace gives the designer the available volume and in addition it indicates the existence of workspace areas that are close but non-connected through a direct end-effector path. Once again, such information on the mechanism has been obtained using only basic knowledge on the mechanism, namely the closed-loop equations.

Conclusions and outlook

This paper presents the use of higher-order continuation method for the preliminary design of compliant mechanisms using the rigid-body replacement method. The method is general. It can be applied to any mechanism since only the geometrical closed-loop equations are required. It results in continuous bifurcation diagrams revealing the kinematic input-output relationships as well as the presence of parallel singularities, serial singularities, architectural singularities and workspace boundaries. The advantage of the method is fourfold:

1. The bifurcation diagrams give the input-output coordinate relationships for a given end-effector path. This is helpful in selecting the most appropriate mechanism working configuration. 2. The different possible joint velocity profiles for a given end-effector path are made visible. This is helpful in selecting the joints to be actuated. Those with non-linear behavior can be interesting for particular actuators such as piezoelectric actuators requiring mechanical amplification. 3. Singular configurations are identified which may be considered as candidates for singularity-based compliant mechanism design [START_REF] Rubbert | Using singularities of parallel manipulators to enhance the rigid-body replacement design method of compliant mechanisms[END_REF]. Some of the ones identified for the 3-US mechanism have never been described in the literature, to our best knowledge. 4. The workspace boundary defines the exact volume within which the end-effector can move. The possible presence of disconnected portions of the workspace which, despite mutual proximity, are not be connected by admissible trajectories is detected.

A 3-US mechanism has been analyzed in order to illustrate the relevance of the method for the design of compliant mechanisms. The proposed analysis clearly shows the mechanism's properties and interests: a large reachable workspace and particular kinematic behavior as screw motion and platform gyration. These latter are respectively illustrated with their bifurcation diagrams and commented.

Future works are to consider sensitivity analysis, already provided by Diamant [START_REF] Charpentier | Sensitivity of solutions computed through the asymptotic numerical method[END_REF][START_REF] Charpentier | Sensitivity computations in higher order continuation methods[END_REF], to give the designer additional information on, for instance, the sensitivity of kinematic properties to manufacturing errors. This point is currently under investigation. The derivative of z G with respect to θ 2 is not defined for θ 2 = 0, but a power series in θ 2 = 0 indicates that: Four values of θ 1 can be obtained as functions of z G :

             lim θ 2 →0 -
θ 1 ∈           

Figure 1 :

 1 Figure 1: In-plane (left) and out-of-plane (right) configurations of the 3-US mechanism.

Figure 2 :

 2 Figure 2: Parameterization of the 3-US mechanism.

  Poses at labelled points A to E. The planar configuration with null joint angles is plotted in green.

Figure 3 :

 3 Figure 3: Numerical and analytical bifurcation diagrams, with associated mechanism configurations, for (φ G , ρ G ) = (0, 0).

4 6 7 (

 7 a) Bifurcation diagram with labelled points. Poses at labelled points. The planar configuration with null joint angles is plotted in green.

Figure 4 :

 4 Figure 4: Bifurcation diagram and associated mechanism configurations for (ρ G , z G ) = (10, 10).

Figure 5 :

 5 Figure 5: 3D views of 1/12th of the reachable workspace (α = π/2; R = 10; φ G ∈ [π/3; 2π/3]).

Figure 6 :

 6 Figure 6: 3D view of the reachable workspace (α = π/2; R = 10).

  Maximal workspace of a leg for R = 10 and α = π/2. (b) Cylindrical coordinates of B.

Figure B. 7 :

 7 Figure B.7: Parameterization of the leg end B and its workspace.

θ 1 and θ 2 (R

 2 Fig. B.7b) are equal to: -L 2 cos (θ 2 ) L 1 cos (θ 1 ) -L 2 sin (θ 2 ) sin (θ 1 ) L 1 sin (θ 1 ) + L 2 sin (θ 2 ) cos (θ 1 )

2R sin α 2 .θ 2 = arccos 1 - z 2 G L 2 , 2 GL

 22222 This confirms the tangent discontinuity observed on the analytical bifurcation diagram of Fig. 3b. Analytical relationship between the altitude z G and the angle θ 1 The angle θ 2 can be expressed as a function of z G from equation (B.5): (B.7) with L = 2R sin α 2 . Using (B.7) in (B.1), one deduces: z G = L 1 sin (θ 1 ) + L 2 sin arccos 1 -z 2 cos (θ 1 ) = L 1 sin (θ 1 ) + L 2 1 -1 -

3D figures are provided as supplementary material.

Appendix A. Diamanlab implementation for the 3-US residual equations function Residual = R ( obj , U ) global Ck % Input geometrical data R =10; % mechanism radius alpha =90* pi /180; % angle AiOBi expressed in radius xw =0; yw =0; % fixed abcissa and ordinate for the centroid X of the platform % The Diamanlab unknown U =( X , zw ) is split as follows zw = U (7); % Height of the centroid X of the platform , continuation parameter X (1:6)= U (1:6); % output coordinates X if isa (U , ' Taylor ') Residual = Taylor ( get (U , ' order ') , zeros ( get ( obj , ' neq ' ) ,1)); else Residual = zeros ( get ( obj , ' neq ' ) ,1); end % Local geometrical parameters L =2* R * sin ( alpha /2); % leg chord alpha1 = alpha ; alpha2 = alpha ; alpha3 = alpha ; % angles AiOBi ( all equal , here ) % projection of the ith cord on the x -and y_axis , respectively L1x = R * sin ( alpha1 ); L1y =2* R *( sin ( alpha1 /2))^2; % independent from zw L2x = R * sin ( alpha2 ); L2y =2* R *( sin ( alpha2 /2))^2; L3x = R * sin ( alpha3 ); L3y =2* R *( sin ( alpha3 /2))^2; % Equation ( 2

Parameterization

The three legs behave identically during the mechanism screw motion. Only the first leg is therefore considered and we will not designate in the following the leg number by a subscript index. The cartesian coordinates of point B as a function of the angles 

, L 1 = R sin (α), L 2 = 2R sin (α/2) 2 . This allows us to build the second part of the bifurcation diagram in Fig. 3b