Emmanuel Abraham
email: emmanuel.abraham@u-bordeaux.fr

Harsono Cahyadi

Mathilde Brossard

Jérôme Degert

Eric Freysz

Takeshi Yasui

Development of a wavefront sensor for terahertz pulses

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Hartmann sensor was invented more than a century ago for optical beam analysis [START_REF] Hartmann | Bemerkungen über den Bau und die Justierung von Spektrographen[END_REF]. It makes it possible to locally measure the wavefront slopes of an incoming optical radiation by using a mask composed of a circular holes array placed just in front of a 2D sensor. It was improved in 1970 by Shack who replaced the mask by a lens array [START_REF] Shack | Production and use of a lenticular Hartmann screen[END_REF]. At present, in the visible spectral region, commercially available Shack-Hartmann sensors can perform wavefront characterization for adaptive optics and a wide variety of applications including astronomy, ophthalmology and microscopy. These sensors are widely used to measure optical aberrations such as astigmatism, coma, spherical aberration, etc [START_REF] Platt | History and principles of Shack-Hartmann wavefront sensing[END_REF]. In the terahertz (THz) spectral domain, it is still challenging to fully measure the spatial profile and wavefront of a THz beam due to the lack of effective THz cameras.

A few previous studies were oriented towards the measurement of THz wavefronts for both continuous wave and pulsed sources. In 2008, Bitzer et al. determined the beam profiles of THz pulses after passing through a hyperhemispherical silicon lens [START_REF] Bitzer | Beam-profiling and wavefront-sensing of THz pulses at the focus of a substrate-lens[END_REF]. They observed an asymmetric spatial-temporal field dynamic in the focus of the lens attributed to a distortion of the incident THz wavefront. However, they did not measure the THz wavefront but only the THz beam profile and their indirect and time-consuming method used a xy-scanning system in order to determine point-by-point the beam profile. In 2012, Cui et al. proposed the first THz wavefront measurement associated with a Hartmann sensor and indirect xy-scanning method to detect the radiation with a pyroelectric sensor [START_REF] Cui | Terahertz wavefronts measured using the Hartmann sensor principle[END_REF][START_REF] Cui | Beam and phase distributions of a terahertz quantum cascade wire laser[END_REF]. The same year, Richter and associates presented a similar arrangement with a 2D detection composed by a microbolometer camera [START_REF] Richter | Terahertz wavefront measurement with a Hartmann sensor[END_REF]. However, this method cannot be applied with short THz pulses produced by femtosecond laser pulses owing to the low sensitivity of the thermal incoherent detector. Another aperture scanning device was proposed in 2013 to characterize the THz beam propagation in a time-domain spectrometer [START_REF] Molloy | Characterization of Terahertz beam profile and propagation[END_REF]. Associated with a Hartmann mask, it has been used to produce 2D topographical image of the THz wavefront. From these few previous works, we can conclude that only two studies have been proposed for the wavefront reconstruction of THz pulses produced by femtosecond laser pulses. Both approaches were mainly limited by a 2D aperture scanning [START_REF] Bitzer | Beam-profiling and wavefront-sensing of THz pulses at the focus of a substrate-lens[END_REF][START_REF] Molloy | Characterization of Terahertz beam profile and propagation[END_REF].

For THz pulse detection, it is well-known that spatial-temporal profiles of THz transient fields can be directly determined through electro-optic (EO) sampling in a nonlinear crystal [START_REF] Wu | Free space electro-optic sampling of terahertz beams[END_REF]. By varying the time delay between the probe and THz pulses, EO sampling can record the spatial dependence of the full (amplitude and phase) THz electric field, E THz , in the crystal plane. Then, from the phase variation of E THz in this plane, it is numerically possible to determine the surface of equi-phase, which constitutes the definition of the wavefront. This approach is close to the principle of THz pulse time-domain holography and the concept of phase retrieval in THz imaging [START_REF] Petrov | Holography and phase retrieval in terahertz imaging[END_REF]. However, the method requires to propagate the THz waves, using the Rayleigh-Sommerfeld equation or the angular spectrum algorithm [START_REF] Goodman | Introduction to Fourier Optics[END_REF], in order to calculate the location of the THz wavefront at the crystal position.

In this paper, derived from visible or infrared technology, we use a more standard method for THz wavefront measurement, which is based on the classical utilization of a Hartmann mask. This approach can appear less straightforward since it does not use the phase but the amplitude of E THz to measure the wavefront. However, it employs robust algorithms such as the modal reconstruction method to determine the final wavefront decomposed on the basis of the Zernike polynomials [START_REF] Southwell | Wave-front estimation from wave-front slope measurements[END_REF]. This calculation does not need any assumptions on the THz wave propagation. This constitutes the main advantage of the Hartmann method. The originality of our approach, compared to refs. [START_REF] Bitzer | Beam-profiling and wavefront-sensing of THz pulses at the focus of a substrate-lens[END_REF] and [START_REF] Molloy | Characterization of Terahertz beam profile and propagation[END_REF], concerns the use of 2D measurement system which does not need to scan any aperture through the THz beam cross section. For that purpose, we don't use a THz camera but an indirect detection of the THz beam thanks to a 2D EO imaging system [START_REF] Yasui | Real-time terahertz color scanner for moving objects[END_REF]. It means that the THz wavefront to analyze is transferred to a visible laser beam that can be more simply analyzed using a visible camera.

2

Combined with a Hartmann mask, this method has been applied to fully determine the frequency-resolved wavefront of THz pulses. Especially, the Zernike modal reconstruction method has been used to measure the wavefronts and optical aberrations of planar and convergent THz beams.

Experimental setup

The proposed method requires THz pulses with electric field strengths roughly larger than 1 kV/cm for effective 2D EO sampling. Therefore, it requires an amplified femtosecond laser source and applies well with THz generation by laser-induced air plasma, optical rectification in nonlinear crystals (ZnTe, LiNbO 3 , DAST) and large interdigitized photoconductive antennas. In our study, THz pulses were generated by optical rectification of amplified femtosecond laser pulses (800 nm, 1 mJ, 150 fs) in a 1 mm-thick ZnTe crystal with a clear aperture of 25 mm (Fig. 1(a)). The collimated THz beam is sent into the ZnTe crystal for collinear EO detection with a time-delayed laser probe pulse, reflected by the front face of the crystal. In the crystal, the spatial distribution of the broadband (0.1-2 THz) THz electric field modifies the intensity distribution of the probe beam that is finally captured by a 232x232 pixels CMOS camera after passing through an analyzer and a lens that images the crystal plane (magnification γ = -0.23). Similar optical arrangement was previously used for realtime 2D and 3D THz imaging [START_REF] Jewariya | Fast three-dimensional terahertz computed tomography using real-time line projection of intense terahertz pulse[END_REF]. For a given time delay between the THz and the probe pulses, the system is able to provide a 2D (x-horizontal and y-vertical directions) image at 800 nm corresponding to the distribution of the THz electric field at the crystal position. This image can be obtained with a 500 Hz acquisition rate thanks to the dynamic subtraction. However, signal-to-noise ratio higher than 100 usually requires an averaging time of 0.5 s. By changing the time delay between the THz and probe pulses, it is possible to record the temporal evolution of this THz electric field. Thanks to this indirect THz imaging, our system will be able to analyze the spatial distribution of the THz electric field at the crystal position.

The Hartmann mask consists of a 1 mm-thick metallic plate with 9 by 9 circular holes (1 mm diameter with 2 mm periodicity in the horizontal and vertical directions). The mask is placed 10 mm in front of the ZnTe crystal, as indicated by the dashed line in Fig. 1(a). This distance D corresponds to a good compromise between THz beam diffraction and the wavefront slope measurement, as already noticed in refs. [START_REF] Cui | Terahertz wavefronts measured using the Hartmann sensor principle[END_REF] and [START_REF] Richter | Terahertz wavefront measurement with a Hartmann sensor[END_REF]. To confirm this mask design, we calculated the intensity distribution of the THz electric field at the distance D = 10 mm behind the mask, assuming a Fresnel-Kirchhoff wave propagation of an incident monochromatic plane wave at 1 THz. Figure 1(b) shows the corresponding simulated THz intensity distribution, confirming that the mask design is appropriated since each hole produces a clearly projected imaging spot in the crystal plane.

Experimental results

Before inserting the Hartmann mask, the THz pulse is directly sent onto the ZnTe crystal for EO sampling with the femtosecond laser probe pulse. The THz beam incident on the crystal is characterized by an expected planar wavefront. Then, we inserted a plan-convex THz lens (f' = 100 mm focal length) at e = 65 mm before the crystal (not represented in Fig. 1). This second configuration provides a converging spherical THz beam at the crystal position with a theoretical ' 35 f e -= mm radius of curvature.

Figure 2(a) (resp. 2(c)) shows the transversal xy-image of the collimated planar (resp. converging spherical) THz beam for a 0 ps time delay between the THz and probe pulses (corresponding to the main negative peak of the THz signal in the center of the image). Videos are available as supplementary materials to see the complete evolution of the transversal THz beam profile for different time delays. For both configurations, Fig. 2(b) and 2(d) show the temporal evolution of the THz pulse in the vertical y-direction. These data make it possible to qualitatively assume that the THz wavefront is planar (resp. spherical and convergent), even if no quantitative optical aberrations can be measured at this stage. However, the radius of curvature of the spherical wavefront can be estimated as follows.

From the temporal waveform of the THz electric field, as presented in Fig. 3, we can measure a time delay Δt = 0.51 ps between the main negative (blue color in Fig. 2(c)) and positive peaks (red color in Fig. 2(c)) of the THz pulse. Then, we measure the distance r = 3.13 mm between the center of the main blue circle and the radius of the circular red ring in Fig. 2(c), which corresponds to the intersection of the crystal plane with two spherical wavefronts, one for the negative peak and one for the positive peak (dashed line in Fig. 2(d)). Assuming a perfect spherical wavefront, the radius of curvature is given by:

2 r R c t = Δ ( 1
)
where c is the speed of light. Our measurement gives R = 32 mm which is close to the 35 mm theoretical value. The difference can be explained by optical aberrations that have not been taken into account in this first simple measurement. We then inserted the Hartmann mask at the distance D = 10 mm before the crystal in order to perform a complete quantitative evaluation of both the planar and converging wavefronts. First, one has to record a reference image indicating the locations of the mask holes in the plane of the ZnTe crystal, in the absence of the THz beam. This reference image is presented in Fig. 4(a) and can be used for any further THz wavefront measurement. The strong direct reflection from the ZnTe crystal has been removed by slightly tilting this crystal without any visible deviation of the spot positions in the image. Then, Fig. 4(b) (resp. 4(c)) shows the locations of the mask holes in the plane of the crystal in the presence of the planar (resp. converging) THz beam. Both images have been calculated at 1.2 THz after FFT of the temporal data. However, the inset of Fig. 3, which represents the THz spectrum transmitted by the ZnTe crystal and Hartmann mask, indicates that our analysis can be extended from 0 to 1.5 THz owing to the large spectral bandwidth of the THz source. Since no significant difference has been observed by changing the frequency of the wavefront analysis from 0 to 1.5 THz, we limited our demonstration to 1.2 THz in this paper.

Wavefront reconstruction and discussion

The calculation of the wavefront was performed using Igor Pro and ImageJ softwares. For the three images presented in Fig. 4, taking into account the 7 by 7 central spots, the program determines the coordinates of the centroids of each spot which are directly related to the local slopes of the THz wavefront. It means that the wavefront is reconstructed for a beam diameter of 12 mm. The final calculation of the wavefront is performed by using a least-squares method in which the wavefront W(x,y) is represented by a linear combination of the Zernike polynomials [START_REF] Southwell | Wave-front estimation from wave-front slope measurements[END_REF]:

( , ) ( , ) i i i W x y a Z x y =  (2)
where a i represents the amplitudes of the polynomial Z i , each polynomial corresponding to one specific aberration. The calculated wavefront of the planar (resp. converging) beam is presented in Fig. 5(a) (resp. 5(b)). The colored scale indicates the local wavefront deformation in µm (accuracy and dynamic range of the measurement will be discussed later). The planar wavefront deformation extends from -50 to + 50 µm (-λ/5 to + λ/5 at 1.2 THz) with a visible tilt along the horizontal x-axis. The converging wavefront exhibits a clear spherical shape extending from -500 to + 450 µm (-2λ to + 1.8λ at 1.2 THz), with slightly visible additional aberrations. A more precise representation of the aberrations can be observed by representing the amplitudes of the 14 first Zernike polynomials (Fig. 6(a) and (b)). The x-tilt amplitude of the planar wavefront is a 2 = 18 µm (λ/14 at 1.2 THz). Then, the predominant aberration is the astigmatism along the x-and y-axis (a 5 = -12 µm or -λ/21 at 1.2 THz). For the converging spherical wavefront, it is noteworthy that the same x-tilt is present, which means that the lens is properly aligned since it does not provide any additional wavefront tilt. A large defocus of a 4 = 272 µm (1.1λ at 1.2 THz) is also measured together with the presence of spherical aberration (a 12 = -25 µm or -λ/10 at 1.2 THz). Assuming a perfect spherical wavefront (i.e. a i = 0 except for i = 4), the wavefront W(r) writes:

2 2 4 ( ) ( ) 2 3 2 W r a R αρ ρ = = (3)
where 2 2 1 / 2 ( ) r x y = + and α = r/ρ is a normalization factor since the Zernike polynomials are defined for a maximum circular aperture with ρ = 1. In our measurement, with α = 8.48 corresponding to the analysis of the 7 by 7 central spots in Fig. 4, one get R = 38 mm, which is close to the theoretical value f'-e + D = 45 mm. The 15% deviation can be attributed to the imperfect spherical wavefront, i.e. the previously noticed presence of x-tilt and spherical aberration.

The first order x-tilt and y-tilt terms are not actually true optical aberrations, as they do not represent a model curvature in the wavefront. Similarly, the second order defocus is directly induced by the presence of the converging lens and is generally not considered for the final analysis of the optical aberrations. We thus plotted again the reconstructed THz wavefronts of the planar and converging spherical THz beam after fixing to zero the x-tilt, y-tilt and defocus amplitudes (Fig. 5(c) and(d)). At present, the planar wavefront is mainly distorted by the presence of 0° and 90° astigmatism (-12 µm). The converging spherical wavefront is now principally distorted by the primary spherical aberration (-25 µm), even if vertical secondary astigmatism (9 µm) and vertical quadrafoil (-9 µm) are present. Table 1 summarizes the main characteristics of the THz wavefront sensor. Let us first recall that it is appropriate to analyze THz pulses with relatively high peak power since it requires a 2D EO sampling detection with a large aperture ZnTe crystal. The working spectral bandwidth of the sensor extends from 0.1 to 4 THz, mainly limited by the THz absorption of the ZnTe crystal. In this paper, we demonstrated the wavefront analysis for a single frequency (1.2 THz). However, the analysis can be extended for any other frequency as long as the EO signal can be detected. It is important to note that if the THz emission spans over a larger spectral range (e.g. THz pulses generated from two color laser-induced air plasma), a 0.1-8 THz sensitivity could be obtained for instance with a GaP crystal, with the limitation of the reduced efficiency of the EO effect.

To evaluate the performance of our Hartmann wavefront sensor, we have to consider the characteristics of the mask, the distance D = 10 mm between the mask and the ZnTe crystal and the magnification γ = -0.23 of the imaging lens between the crystal and the camera (Fig. 1(b)). We first calculate the absolute accuracy of the wavefront sensor. It corresponds to the minimum detectable wavefront slope that can be measured by the system. It depends on D, γ, the pixel size of the camera, the precision of the centroid determination and the signal-tonoise ratio of the sensor. We consider that the minimum slope corresponds to a spot deviation of 1 pixel (20 µm) on the camera sensor. As a consequence, the minimum spot deviation on the ZnTe crystal is (λ/35 at 1 THz). This parameter is comparable to that of a commercial Hartmann sensor operating in the infrared range. It could be reduced by increasing the magnification factor γ or the distance D. This measurement also indicates that a Zernike coefficient with an amplitude smaller than 8.7 µm is not significant since it is smaller than the wavefront measurement accuracy Δz min . This limit is indicated in Fig. 6 by the vertical dashed lines. It shows that, for the planar wavefront we analyzed in this paper, only a 2 (x-tilt) and a 5 (0/90° astigmatism) coefficients are really measured, other Zernike coefficients being too small to be detected by our system. For the spherical wavefront, only a 2 (x-tilt), a 4 (defocus), a 12 (spherical aberration), a 13 (vertical secondary astigmatism) and a 14 (vertical quadrafoil) are significant. We can also evaluate the tilt dynamic range of the sensor which represents the maximum wavefront slope that can be measured. It is generally defined with a spot displacement Δx max = d/2 where d = 2 mm is the distance between two adjacent holes on the Hartmann mask. The tilt dynamic range is therefore max max 5.7

x D θ = Δ = °. Considering the 20 mm clear aperture of the crystal, this will provide a maximum wavefront tilt of 1 mm (3.3λ at 1 THz). This tilt is similar to that obtained with a commercial system operating in the visible even if in this case it corresponds to a dynamic range of nearly 1000λ. The dynamic range of our sensor can be increased by using a hole mask with either a larger spot separation or a smaller product D γ .

However, increasing the dynamic range will also decrease the number of Zernike coefficients available to represent the wavefront. Conversely, increasing the dynamic range by shortening D γ will also decrease the absolute efficiency (larger minimum wavefront tilt).

Conclusion

We demonstrated the combination of a Hartmann mask, a ZnTe EO crystal and a CMOS camera for the quantitative wavefront determination of THz pulses. The proposed method is achromatic and can be applied for all frequencies carried by the THz pulses, limited by the THz absorption of the detection crystal. Further optimization of the technique includes the use of a Shack-Hartmann mask constituted by a THz matrix of lenses in order to improve the accuracy of the centroid spot determination. Future applications of the sensor will be the precise THz wavefront analysis of intense THz sources. It is well known that THz beams generated from two color laser-induced air plasma or optical rectification of femtosecond laser pulses in LiNbO 3 crystal can suffer from severe optical aberrations due to the complex nonlinear optical processes involved into the THz generation [START_REF] Klarskov | Experimental three-dimensional beam profiling and modeling of a terahertz generated from a two-color air plasma[END_REF]. We believe that our THz wavefront sensor could provide an effective tool to control the THz emission from these sources. Another interesting application concerns the association of the Hartmann sensor with deformable mirrors in order to open the route to THz adaptive optics. This will provide a real advance for time-domain (imaging) spectrometers which require a perfect focalization of the THz beam or any other THz devices sensitive to wavefront distortions.

Figure 1 (

 1 c) illustrates the deviations of the hole positions in the plane of the ZnTe crystal, induced by a THz wavefront distortion. It also shows the strategic role of the crystal as an intermediate plane between the mask and the CMOS camera. The displacement Δx of one hole centroid in the plane of the crystal is clearly converted into a displacement x x δ γ = Δ at the position of the CMOS sensor. Similar deviations Δy and δy can be measured along the perpendicular y-direction.

Fig. 1 .

 1 Fig. 1. (a) Experimental setup of the THz wavefront sensor. (b) Monochromatic THz plane wave (1 THz) passing through the Hartmann mask (1 mm diameter with 2 mm periodicity). Simulation of the intensity distribution at the crystal position (D = 10 mm behind the mask). (c) Schematic representation of THz wavefront distortion and spot displacement. ZnTe crystal plane is imaged onto the CMOS sensor with magnification factor γ = -0.23. D = 10 mm, d = 2 mm.

2

 2

Fig. 2 .

 2 Fig. 2. (a) xy-image of a collimated planar THz beam for a 0 ps time delay. (b) Corresponding yt-image. (c) xy-image of a converging spherical THz beam for a 0 ps time delay. (d) Corresponding yt-image. Red color is for a positive electric field and blue color is for a negative one. Associated supplementary materials: Visualization 1 corresponding to the temporal evolution of (a) and Visualization 2 to the temporal evolution of (c).

Fig. 3 .

 3 Fig. 3. Temporal waveform of the THz electric field. Inset: corresponding THz spectrum.

Fig. 4 .

 4 Fig. 4. (a) Reference locations of the mask holes in the plane of the ZnTe crystal. (b) Locations of the projected light spots in the plane of the crystal in the presence of the planar collimated THz beam. (c) Locations of the projected light spots in the plane of the crystal in the presence of the converging spherical THz beam.

Fig. 5 .

 5 Fig. 5. Reconstructed THz wavefront. (a) Planar beam. (b) Converging spherical beam. (c) Planar beam without tilt and defocus. (d) Converging spherical beam without tilt and defocus.

Fig. 6 .

 6 Fig. 6. Amplitudes a i of the Zernike polynomials. (a) Planar THz wavefront. (b) THz converging spherical THz wavefront. The vertical dashed lines show the wavefront measurement accuracy, as calculated in the paper and indicated in Table 1. The amplitude of the defocus term a 4 = 272 µm is out of range of the scale in Fig. 6(b) for a better visualization of the other Zernike terms.

  Considering the d/2 = 1 mm diameter of the mask holes, this will provide a local wavefront measurement accuracy

Table 1 . Main characteristics of the THz waveform sensor

 1

	Aperture dimension	25x25 mm 2
	Hartmann mask	9x9 holes, 1 mm diameter,
		2 mm periodicity
	Working frequency range with ZnTe	0.1 -4 THz
	detection crystal
	Wavefront measurement accuracy	0.5 (λ/25 @ 1 THz)
	Tilt dynamic range	5.7 (3.3 λ @ 1 THz)

Acknowledgments

The authors thank Junji Tamatani for the machining of the metallic Hartmann mask. They also thank the financial support from CNRS (PICS 6507).


Julien Prades

Franck Jourdan

Olivier Company

Sébastien Krut

François Pierrot

Dynamics effects on natural frequencies in modal analysis of PKMs

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Parallel robots are said to be faster and stiffer than their serial counterparts. They benefit from continuous technical improvments and their mechanical performances are more impressive, as we can see in [START_REF] Corbel | ŞTowards 100G with PKM. Is actuation redundancy a good solution for pick-and-place?[END_REF] (robot LIRMM 100 g). Inertia effects due to speed and acceleration are more significant on the mecanism. Duringa path tracking sequency, positions, speeds and accelerations are calculated. Actuated torques can be deduced in respect with dynamical model. Then the command set up all the previous datas and the closed loop of control behave as "spring" in regards to actuations. Errors on following, speed and acceleration of the Mobile Platform will induce error on inertia forces and on the command. This paper introduce a methode to allow a quatification of the vibrations and will bring to light their effect on natural frequencies of the robots. Inertia forces of the robot are linked with geometry weight and acceleration. During apick and place strategy, the goal is to move from a point A to a point B in the quickest way possible. This induces high inertia forces and leads the system to vibrate. In order to calculate the vibration frequencies, the mass and stiffness matrices must been known for the robot industriel. These matrices depend on the position. Andreas Müller [START_REF] Müller | ŞStiffness control of redundantly actuated parallel manipulators[END_REF] shown that the stiffness matrix can be modified by external forces, command forces and speed. In this paper, the influence of acceleration on stiffness of the robots and therefore on their vibration frequencies will be explain. This paper is organized as follows. In the next section, we describe the family of redundant parallel manipulators we are focusing on. Then, a simplified dynamic model based on jacobian matrices is derived. In section IV, the proposed analysis methodology is presented, based on the linearization of the dynamic model. The following section gives preliminary simulation results for two manipulators and provides some comments.

II. PKMS CONCERNED WITH THE STUDY

This study is focused on a specific type of robots, light ones. Moreover we are considering fast ones with parallel kinematic structure that are redundantly in actuation. As an example, we can consider the Par2 in [START_REF] Baradat | Par2: a Spatial Mechanism for Fast Planar, 2-dof, Pick-and-Place Applications[END_REF]. The robot has, u actuated joints and its Mobile Platform (MP) has v degrees of freedom . The structure has also u kinematic chains between its base frame and the MP with u ≥ v due tue the actuation redundancy. The objective of this paper is to determinate the impact of the dynamics on robot eigen frequencies for such parallel robots.

In this paper we are considering robots with the following architecture:

• Kinematic chains are composed by two solids and three joints • Actuators can be linear motors or rotative motors, they are modelized by the L 1i joints which can be revolute or primatic • L 2i et L 3i joints can be revolute in the case of a planar mechanism or spherical (or universal) for spatial mechanisms The robot generic joint graph is shown in figure 7. According to this description, we can mention that this family of robots includes R4 Par2 and Dual-V see [START_REF] Corbel | ŞTowards 100G with PKM. Is actuation redundancy a good solution for pick-and-place?[END_REF] and [START_REF] Baradat | Par2: a Spatial Mechanism for Fast Planar, 2-dof, Pick-and-Place Applications[END_REF] and [START_REF] Van Der Wijk | Design and experimental evaluation of a dynamically balanced redundant planar 4-RRR parallel manipulator[END_REF] be done separately as we will see in section III of this paper -The robot dynamics equation can be written using only the Jacobian matrix of the robot (the same as the one used in kinematic study) This is the key to have simple equations.

III. PRELEMINARY SIMPLIFIED DYNAMIC MODEL

Hypothesis 1 allows to write the Dynamic Model related to Jacobian Matrices. Adding the following hypothesis : frictions is neglected, solids keep their shapes and the links are perfect. Preliminary calculus will start with the writing of Inverse Geometric Model(IGM) before establishing robot Dynamic equations. Robot kinematic study gives IGM which links operational coordinates vector X of MP to actuated joint coordinates vector q.

J q q = J x Ẋ

(1)

As J q is a full rank square matrix (away from singular configurations), for our family of robots we can express:

q = J m Ẋ (2)
with:

J m = J q -1 J x (3)
Actuator dynamic equation is obtained by isolating the actuator moving parts. It is written depending on actuators inertia and bars inertia contribution m, forces applied by bars F b and actuators torque τ.

τ + J T q F b = m q (4)
MP dynamic equation is achieved by isolating it. It is written depending on platform inertia and bars inertia contribution M, forces applied by bars F b and external forces F ext :

F ext -J T x F b = M Ẍ (5)
By derivating the IGM with respect to time (eq. 2), connection between the two dynamic equations 4 and 5 is obtained :

q = J m Ẍ + J m Ẋ (6)
Which give us the complete robot dynamic model.

F ext = (M + J T m m J m ) Ẍ + J T m m J m Ẋ -J T m τ (7)
where:

• J q , J x and J m matrices are written based on the following variables X and q • J m matrix is written based on the following variables X, Ẋ and q, q For a sake of clarity, dependence of associated variables won't be specified, except in case where numericals values will have to be calculated. This point will be important in section IV where matrices numerical evaluation and derivation are necessary. Then we will have to keep in mind that q = f (X) and q = f (X, Ẋ).

IV. PROPOSED ANALYSIS METHODOLOGY

Robot dynamic equation is non linear. To extract informations of phenomenons which interest us, meaning actuating and moving effects over natural frequencies, we apply first order Taylor formula near to functional point :

P f = (X 0 , Ẋ0 , Ẍ0 , τ 0 , F 0 ) (8)
It should be noted that in this study we take into account torques in the functional point, as cable robot studies. To simplify the writing, we use function h as :

h = (M + J T m m J m ) Ẍ + J T m m J m Ẋ -J T m τ -F ext (9)
Taylor equation of h gives :

∂ h ∂ Ẍ P f ∆ Ẍ + ∂ h ∂ Ẋ P f ∆ Ẋ + ∂ h ∂ X P f ∆X + ∂ h ∂ τ P f ∆τ + ∂ h ∂ F ext P f ∆F ext = 0 (10)
To go any further, we need informations about τ, which we can get thanks to the control law. In general, the simplest axis control is a proportional derivative controler (PDC). The PDC will be set in actuated joint coordinates to avoid singularity problems which can occur. PDC will be integrated in the equation. PDC law depends on joint position error ε, joint speed error ε and torques that don't cause movements of the MP. These torques allow to solve many problems like mechanical clearance and other issues [START_REF] Müller | Internal Preload Control of redundantly actuated Parallel Manipulators Ű Backlash avoiding Control[END_REF]. Control law equation is expressed:

τ(ε, ε) = τ 0 -K D ε -K P ε (11)
By studying the variations and supposing that control law stay constant near to the functional point P f , we have :

∆τ = -K D ∆ q -K P ∆q (12)
Using kinematic equations, we can write equation 13 depending on X :

∆τ = -K D J m ∆ Ẋ -K P J m ∆X (13)
Looking at the equation terms, we notice that we get a second order differential equations system in ∆X with constant coefficients depending on P f . To simplify the writing, we will use matrices M, C, K depending on inertia, Jacobian matrix and control law.

Then, we get the following equation:

∆F ext = M ∆ Ẍ + C ∆ Ẋ + K ∆X (14)
Matrices M, C and K are now expressed:

M(X 0 ) = M + J T m (X 0 ) m J m (X 0 ) (15)
Matrices C and K which are too complex are decomposed as follows:

-C v , kinematic damping matrix depending on speed and position.

C v ( Ẋ0 , X 0 ) = ∂ J T m m J m Ẋ ∂ Ẋ Ẋ0 ,X 0 (16)
-C d , command law damping matrix depending on derivating command law and position.

C d (X 0 ) = J T m (X 0 ) K d J m (X 0 ) (17)
-K i , dynamic stiffness matrix depending on acceleration and position.

K i ( Ẍ0 , X 0 ) = ∂ J T m m J m Ẍ ∂ X Ẍ0 ,X 0 (18)
-K v , kinematic stiffness matrix depending on speed and position.

K v ( Ẋ0 , X 0 ) = ∂ J T m m J m Ẋ ∂ X Ẍ0 ,X 0 (19)
-K g , geometrical stiffness matrix depending on torques and position.

K g (X 0 , τ 0 ) = - ∂ J T m τ ∂ X X 0 ,τ 0 (20)
-K a , proportional command stiffness matrix depending on proportional command law and position.

K a (X 0 ) = J T m (X 0 ) K q J m (X 0 ) (21)
K g and K a are the same matrices that in paper [START_REF] Kock | A parallel x-y manipulator with actuation redundancy for high-speed and active-stiffness applications[END_REF].

The different partial derivatives have to be done, taking into account that q = f (X). We have therefore, for each partial derivatives of w :

∂ w(X, q(X)) ∂ X P f = ∂ w(X, q 0 ) ∂ X P f + ∂ w(X 0 , q) ∂ q P f J m (22)
The equation 14 becomes :

∆F ext = M ∆ Ẍ + C v + C d ) ∆ Ẋ + K i + K v + K g + K a ∆X (23
) These mathematical developments show the impact of the functional point P f over the different exposed matrices. These equations are able to describe the instant vibration behaviour of robot when it follows a trajectory, typically when there is a shock.

V. VALIDATION CASE STUDY

A. Type Analysis

In this section, only numerical simulations will be presented. Two cases of robots natural frequencies are shown: -1 st case : for a given trajectory, during each interval, position and torques static functional point lowest natural frequency is calculated in static case. -2 nd case : for a given trajectory, during each interval, position, speed, acceleration and torques dynamic functional point lowest natural frequency is calculated in dynamic case. We chose to work with redundant PKM easy to handle, but with Jacobian matrix which depends on the position. All studied robots structures can be modelized in a plane: PRR-4 and Dual-V manipulators, are both three degree of freedom robots (v = 3). MP motions are the two translations of the ( x, y) plane plus the rotation about z. Both manipulators have four actuation degrees (u = 4):

• four translations for PRR-4

• four rotations for Dual-V PRR-4 has been chosen and modelized to observe the impact of prestress and movement around one parallel singularity. Dual-V has been chosen because future experiments will be possible on an existing demonstrator. More details about these robots will be given in their respective subsection. To perform the calculations, the following trajectory has been chosen: a translation along the axis (O, y) beginning from y d = 0.1m to y f = -0.1m by blocking x translation and z rotation.

M ∆ Ẍ + K i + K v + K g + K a ∆X = 0 (24)
Natural frequencies are obtained by finding eigenvalues of the following matrix.

-1 st case : static case

M -1 K g + K a ( 25
)
K g is calculated with static torques only with redundant PKM if not K g = 0 -2 nd case : dynamic case.

M -1 K i + K v + K g + K a ( 26
)
K g is calculated with dynamic torques (redundant PKM could have J T m null space torques)

The only thing still needed is to formulate each robot matrices.

B. Dual-V Case

In this section the Dual-V case will be studied and robot design is shown in figure V-B. Actuation degrees are motors rotation at A i and operational coordinates vector is movements and rotation at P 5 . For more details about Dual-V see reference [START_REF] Van Der Wijk | Design and experimental evaluation of a dynamically balanced redundant planar 4-RRR parallel manipulator[END_REF]. We use the same Jacobian matrices. Geometry, inertia and control parameters are shown in table II.

The actuated coordinates vector is :

q =     θ 11 θ 21 θ 31 θ 41     (29)
The operational coordinates vector is :

X =   x 5 y 5 θ 5   (30)
In Figure 6, first natural mode frequencies of Dual-V will be calculated every millisecond. The blue curve shows static case frequencies and the green curve shows dynamic case frequencies. When comparing both curves, we can see that acceleration decreases the frequency and that speed inscreases the frequency. Informations given by this model could help to control robots because it allows to evaluate moving robot natural frequency.

Fig. 1 .

 1 Fig. 1. Robots type generic joint graph

Figures 2 ,

 2 3 and 4 show position, speed and acceleration with respect to time according to

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Position over time

Fig. 5 .

 5 Fig. 5. Dual-V PKM

Fig. 6 .

 6 Fig. 6. First mode natural frequency of Dual-V

Table I .

 I

	parametres	Values
	Jerk max	20000 m/s 3
	Acceleration max	200 m/s 2
	Vitesse max	10 m/s
	Distance max	0.2 m
	TABLE I
	TRAJECTORY PARAMETERS

C. PRR-4 Case

Robot dimensions have been chosen in order to have a parallel singularity in the center of the workspace to see the wanted effect whatever working conditions are. Actuated joints are the prismatic ones with the x direction at point A i . This actuation is done by linear motors (see Figure ). The robot kinematic scheme is given on Figure 7. Geometric, inertia and command parameters are shown in Table III The actuated coordinates vector is :

The operational coordinates vector is :

Robot Jacobian matrix will be found, based on bars speed (1 to 4) by using the following property for a rigid body (points A and B are supposed to belong to this body):

Here, we chose u i for bar i director vector from A i to B i .

In Figure 8, first natural mode frequencies of PRR-4 will be calculated every millisecond. The blue curve shows static case frequencies and the green curve shows dynamic case frequencies. When comparing both curves, we see the same kind of frequency behaviour. We also note that in static case, first mode frequency become null when y = 0 (parallel singularity) but it does not nullify in the dynamic case. This method shows how we could go through singularities of this kind. In this paper, we only have mathematical and numerical developments which show the impact of movements on robots natural frequency. Acceleration is harmful for first mode frenquency and speed has positive effect. The experimental confirmation of this model is in process with the Dual-V robot and another paper will show results of experiment. The next stage would be the improvement of the robots command by using informations given by this model.