Julien Favier

Alfredo Pinelli

Ugo Piomelli

Control of the separated flow around an airfoil using a wavy leading edge inspired by humpback whale flippers

Keywords: Flow control, Biomimetics, Immersed Boundary, humpback whale flippers Contrôle d'écoulements, biomimétisme, Méthode des

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

Résumé

Contrôle du décollement autour d'un profil d'aile présentant un bord d'attaque ondulé inspiré des ailerons de la baleine à bosse. L'influence d'ondulations géométriques le long du bord d'attaque d'un profil d'aile est étudiée numériquement à faible nombre de Reynolds, dans une optique de contrôle passif du décollement et en se focalisant sur les mécanismes physiques mis en jeu. Inspiré des tubercules présents sur les ailerons des baleines à bosse, ce bord d'attaque ondulé est modelisé par une sinusoïde le long de l'envergure, dont la longueur d'onde et l'amplitude constituent les paramètres de contrôle. Une simulation numérique directe est effectuée sur un profil NACA0020 dans une configuration d'écoulement massivement décoll'e (α = 20 o ), avec et sans ondulation de bord d'attaque. Les frontières solides complexes engendrées par la variation des paramètres de l'ondulation géométrique sont traitées par la méthode des frontières immergées (IBM). Les noyaux des opérateurs d'interpolation/diffusion sont construits en utilisant la méthode RKPM [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finitevolume Navier-Stokes solvers[END_REF]. Une étude paramétrique permet d'extraire un jeu de paramètres d'ondulation qui amène à une modification significative de la topologie de l'écoulement décollé, qui se retrouve dominé par des tourbillons orientés vers l'aval et générés sur les cotés des protubérances de bord

Introduction

The humpback whale, Megaptera novaeangliae, is one of the most massive animals living in the oceans; at the same time, this whale is well known for astonishing maneuvering capabilities. One of the features of this mammal that affect its hydrodynamic performance is its very specific surface geometry, made of protuberances, especially on the leading edges of the pectoral flippers, called tubercles (Figure 1), which inspired the present work. Tubercles are the result of a long process of multiple-task natural optimization: they may have for instance thermo-regulation purposes in addition to their hydrodynamical role [START_REF] Stanway | Hydrodynamic effects of leading-edge tubercles on control surfaces and in flapping foil propulsion[END_REF]. Thus, the structural and shape parameters of the flipper may be suboptimal, as far as only hydrodynamics are concerned. Nonetheless, there is some interest in understanding the hydrodynamical features of tubercles, to adapt them to technological applications, within the renewed interest that biomimetics has recently received.

The potential benefits of tubercles on the aerodynamic performances of a wing airfoil were suggested by Bushnell & Moore [START_REF] Bushnell | Drag Reduction in Nature[END_REF] and later, from a more morphological point of view, by Fish & Battle [START_REF] Fish | Hydrodynamic design of the humpback whale flipper[END_REF], who showed that the geometrical properties of the tubercles can influence the aerodynamical performances of a wing. It was believed since the beginning that the tubercles would have an influence on stall, delaying the angle of attack at which lift decreases, and also reducing the form drag at high angles of attack.

Several studies were performed over the last two decades to assess the influence of the tubercles, both experimentally and numerically. The first numerical study was performed by Watts & Fish, using a numerical 3D panel method [START_REF] Watts | The influence of passive, leading edge tubercles on wing performance[END_REF]; they found a slight increase of lift coefficient (4.8%) and a 11% reduction of drag for a pre-stall angle of attack. The wind tunnel experiments of Miklosovic et al. [START_REF] Miklosovic | Fish Leading-edge tubercles delay stall on humpback whale (Megaptera novaeangliae) flippers[END_REF] exhibited drastic enhancements of lift for a post-stall angle of attack, also showing a delay of 40% of the stall angle. Both studies were performed at Reynolds numbers of the order of Re = 10 5 (lower than those observed in nature). More recent experiments confirmed these findings as far as the action of bumps on stall at similar Reynolds numbers is concerned [START_REF] Stanway | Hydrodynamic effects of leading-edge tubercles on control surfaces and in flapping foil propulsion[END_REF][START_REF] Johari | Effects of leading-edge protuberances on airfoil performance[END_REF]. In 2003, Patterson et al. [START_REF] Paterson | General-purpose parallel unsteady RANS CFD code for ship hydrodynamics[END_REF] solved the Unsteady Email address: Julien.favier@ciemat.es (Julien Favier). Reynolds-Averaged Navier-Stokes (URANS) equations, and showed that the tubercles can modify the shape on the recirculating zone and improve global performances. A simpler approach was explored by Van Nierop et al. [START_REF] Van Nierop | How bumps on whale flippers delay stall: An Aerodynamic Model[END_REF], who used the Kutta condition to estimate the circulation Γ in a potential flow model. Globally the numerical and experimental studies converge towards similar conclusions, but, as pointed out by in Ref. [START_REF] Van Nierop | How bumps on whale flippers delay stall: An Aerodynamic Model[END_REF], the underlying physical mechanisms and the associated fluid flow features still remain unclear, although there are speculations about the action of the bumps as vortex generators [START_REF] Fish | Passive and active flow Control by swimming fishes and mammals[END_REF], which create vortical structures similar from those emitted at the tip of a delta wing [START_REF] Stanway | Hydrodynamic effects of leading-edge tubercles on control surfaces and in flapping foil propulsion[END_REF]. The experimental study by Stanway [START_REF] Stanway | Hydrodynamic effects of leading-edge tubercles on control surfaces and in flapping foil propulsion[END_REF] also investigated, by means of PIV, the effect of the tubercles on the dynamic stall (as opposed to the static stall studied before) by flapping the wing. Different behaviours were observed in the static and dynamic cases, and the efficiency parameters for a fixed wing are not necessarily the same as those for a flapping one.

In this work we present a three-dimensional Direct Numerical Simulation of Navier-Stokes equations, focussing on the effect of the protuberances on the topology of the flow at low Reynolds number. Our goal is to obtain an insight of the underlying physical mechanism responsible of the hydrodynamic performance enhancements in the absence of transitional effects, or a broadband perturbation spectrum; the Reynolds number, based on the free stream velocity and the chord length, is, therefore, chosen to be Re = 800; this value is lower than that associated with transition to turbulence, but at this Reynolds number, the wake experiences three-dimensional effects issued from the first spanwise instabilities rising before transition to turbulence [START_REF] Bourguet | Capturing transition features around a wing by reduced-order modeling based on compressible Navier-Stokes equations[END_REF]. Although this value is much lower than the one observed in nature, it may allow us to obtain a first picture of the physical mechanism responsible of the performances enhancements. As the morphological parameters of the humpback whale flippers are similar to those of a wing section [START_REF] Fish | Hydrodynamic design of the humpback whale flipper[END_REF] with a thickness of approximately 20% of the chord length, we will perform our investigation on a static NACA0020 wing airfoil, whose chord can vary along the span to model the protuberances. The angle of attack is fixed to α = 20 o , a configuration corresponding to a massive boundary layer separation.

In the following, we will first describe briefly the numerical method that has been used. We will then present some results, including quantitative measures of the lift and drag modifications, as well as instantaneous flow visualizations. On the basis of the results obtained, we will propose a conceptual mechanism to explain the different flow topologies encountered in the fluid flow in the presence of different configurations of leading edge undulations. We will finally draw concluding remarks and recommendations for future works.

Problem formulation

Numerical Method

In this work we use an Immersed Boundary Method (IBM) to enforce the presence and effects of the body in the fluid flow. With this approach, the grid is not required to conform to the body geometry, and the solid wall boundary conditions are satisfied on the body surface by using appropriate boundary forces [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finitevolume Navier-Stokes solvers[END_REF]. This allows to simulate a range of complex geometries on a cartesian grid, using fast and accurate solvers that conserve discretely the mechanical energy of the flow, and are, therefore, suitable to the simulation of transitional cases.

We use a three-dimensional incompressible Navier-Stokes solver to model the fluid flow. The equations are discretized on a staggered grid system [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. Both convective and diffusive fluxes are approximated by second-order centered differences. The fractional time-step method is used for the timeadvancement [START_REF] Chorin | Numerical solution of Navier-Stokes equations[END_REF][START_REF] Kim | Application of a fractional step method to incompressible Navier-Stokes equations[END_REF], in the form of a second-order semi-implicit pressure correction procedure [START_REF] Van Kan | A second-order accurate pressure-correction scheme for viscous incompressible flow[END_REF]. The alternating direction implicit method (ADI) is used for the temporal discretization of the diffusive terms, allowing to transform three-dimensional problem into three one-dimensional ones by an operator-splitting technique, while retaining the formal order of the scheme. The solution of the three-dimensional pressure Poisson equation, needed to carry out the correction step, is also reduced into a series of two-dimensional Helmoltz equations by performing a Fourier transform in the spanwise direction, z, and solving twodimensional problems for each wavenumber using a fast block-tridiagonal solver. The code parallelization relies upon the Message-Passing Interface (MPI) library and the domain-decomposition technique.

The numerical strategy used to impose the desired zero velocity boundary condition at the solid surface (which is a solid and rigid wing) is the following: we first obtain the predicted velocity, u * , explicitly, without the presence of the embedded boundary:

u * = u n -∆t N l (u n , u n-1 ) -Gφ n-1 + 1 Re L(u n ) , (1)
where u n is the divergence-free velocity field at time-step n, ∆t is the time step, N l is the discrete nonlinear operator, G and D are, respectively, the discrete gradient and divergence operators, L is the discrete Laplacian, φ is a projection variable (related to the pressure field). The operators include coefficients that are specific to the time scheme used in this study, a three-steps low-storage Runge Kutta.

Next, the presence of the solid geometry is imposed by using the IBM, via a process of interpolation and spreading [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF]: u * is interpolated onto the embedded geometry of the obstacle, Γ, which is discretized through a number of Lagrangian marker points with coordinates X k :

U * (X k , t n ) = I(u * ) (2)
At this stage, knowing the velocity U * (X k , t n ) at location of the Lagrangian markers, we are able to determine a distribution of singular forces that restore the desired zero velocity values U Γ (X k , t n ) on Γ:

F * (X k , t n ) = U Γ (X k , t n ) -U * (X k , t n ) ∆t . (3)
The singular surface force field given over Γ is then transformed by a spreading operator C into a volume force-field defined on the Cartesian mesh points x i,j,k surrounding Γ:

f * (x i,j,k , t n ) = C [F * (X k , t n )] . (4)
The predicted velocity is then re-calculated, using an implicit scheme for the viscous operator, adding the forces that accounts for the presence of the solid body:

u * -u n ∆t = -N l (u n , u n-1 ) -Gφ n-1 + 1 Re L(u * , u n ) + f * (5)
Finally, the algorithm completes the time step with the usual solution of the pressure Poisson equation and the consequent projection step:

Lφ = 1 ∆t Du * (6)
u n+1 = u * -∆tGφ n . (7)
The key elements of the present IBM are the transformations between the Eulerian mesh and the Lagrangian one, which are carried out through the interpolation and spreading operators, I and C. These two operators are built using a modified Reproducing-Kernel Particle Method (RKPM) [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finitevolume Navier-Stokes solvers[END_REF] that ensures a bijection between the interpolation and spreading operations, implying that the integral of the force is the same when computed in the Lagrangian or Eulerian frames. Important properties of the algorithm are the preservation of the global accuracy of the underlying differencing scheme, and the sharpness with which the interface is resolved. For further details the reader is referred to the work of Pinelli et al. [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finitevolume Navier-Stokes solvers[END_REF].

Numerical setup

The baseline mesh chosen for all the simulations is cartesian, uniform in the z-direction, and non-uniform in x (streamwise direction) and y (normal to the freestream direction). To ensure a fine clustering of the points near the airfoil, a stretching factor of 1.02 is introduced in both x and y-directions. The mesh is made of 6 × 10 6 points (512 × 386 × 32 in x, y and z directions, respectively). The computational box dimensions along the x, y and z coordinates are 38c × 17c × 2c (c being the average chord length), respectively. The chosen box size is large enough to avoid spurious effects due to the outer boundary conditions. In particular, at the inlet and outlet, the boundary conditions are set to uniform inflow and zero-velocity derivative. The top and bottom boundaries are planes of symmetry, and periodic conditions are used in the the spanwise direction. A spanwise cross section of the mesh in the region occupied by the airfoil is shown in Figure 2(a).

The baseline airfoil geometry has a thickness of 20 % of the chord, and follows the shape of a NACA0020 foil. To model the tubercles, a geometrical waviness on the leading edge is introduced, taking the form of a sinusoidal function along the span: f (t) = A * cos(2πz/λ * ), where A * = A/c is the non-dimensional amplitude, λ * = λ/c is the non-dimensional wavelength (Figure 2(b)). An optimization algorithm is used to distribute evenly the Lagrangian markers so that each Eulerian cell contains one marker (a condition required to obtain a reliable behaviour of the IBM [START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finitevolume Navier-Stokes solvers[END_REF]).

A validation study was first performed on a NACA0012 airfoil at a 20 o angle of attack to compare the present results with data obtained at a Reynolds number of Re = 800 by Bourguet et al. [START_REF] Bourguet | Capturing transition features around a wing by reduced-order modeling based on compressible Navier-Stokes equations[END_REF]. Typical values of mean lift and drag coefficients and Strouhal number (St = f c/U ∞ , where f is the frequency of vortex shedding) were used to assess the grid convergence of the simulation, and were found to be in good agreement with the mentioned study when using the mesh described above.

Results and discussion

To assess the hydrodynamic performance of the airfoil, the drag and lift coefficients defined by: were monitored. In [START_REF] Van Nierop | How bumps on whale flippers delay stall: An Aerodynamic Model[END_REF], ρ is the fluid density, U ∞ is the freestream velocity, s is the wing span, and F d and F l are the projections on the x-and y-axis of the aerodynamic force. In addition, three ratios based on these coefficients:

C d = F d (ρU 2 ∞ /2)s c ; C l = F l (ρU 2 ∞ /2)s c , (8)
R 1 = < C d > < C d,b > ; R 2 = < C l > < C l,b > ; R 3 = C ′ l C ′ l,b , (9)
are defined. In (9), < . > represents averaging over time (typically performed over 10 vortex shedding cycles), the subscript b refers to the baseline case (the one with a straight leading edge), and the lift fluctuations are defined by

C ′ l = max(C l -< C l >).
Improving the hydrodynamic performances of the airfoil translates into having

R 1 < 1, R 2 > 1 and R 3 < 1.
We performed a parametric study varying the amplitude A * from 0 (straight leading edge) to 0.1, and the wavelength λ * from 0.25 to 2.0. Note that A * = 0.1 corresponds to an amplitude of the bumps of 10% of the chord, and a wavelength of 2.0 corresponds to two periods of spanwise leading-edge waves per chord distance.

Contour plots of the ratios R 1 , R 2 and R 3 defined in ( 9) are shown in Figure 3 as functions of A * and λ * . Values of these ratios close to unity, which are observed for low amplitudes of the bumps and low wavelengths, indicate that there is no effect of the leading-edge waves on the flow in these geometries. The contours of R 1 (Figure 3(a)) show a maximum drag reduction at λ * ≃ 1 and A * ≃ 0.07, a case also associated to lift reduction (Figure 3(b)) and also to a drastic reduction of lift fluctuations (Figure 3(c)).

Next, we focus on the flow structure modifications that occur when fixing λ * = 1, a value common to all the minima of the ratios [START_REF] Fish | Passive and active flow Control by swimming fishes and mammals[END_REF], as shown in Figure 3, while varying the amplitudes A * . In Figure 4(a) we observe a decrease in drag (up to 35%), when increasing A * up to a minimum, which occurs for A * = 0.07. We also note, in Figure 4(b) the same decreasing trend of the mean lift coefficient. This may be due to the low value of Re for the present study, which is several orders of magnitude lower than the one considered in the experiments, where turbulent conditions are usually found. At λ * = 1 and A * = 0.07 the minima for the drag and lift coefficients are accompained by a complete cancellation of lift fluctuations, indicating the establishment of a quasi-steady wake and a strong weakening of the shedding. An almost-steady behaviour is found, in which the time fluctuations of the baseline Von Kàrmàn instability disappear. The spectrum of the lift coefficient, shown in Figure 5, indicates that as the undulation amplitude increases, the spectrum tends to flatten, and no frequency peak is recorded for the case corresponding to A * = 0.07, which confirms the trend mentioned above.

To gain more insight on the effects of these changes on the topology of the flow, Figure 6 shows isosurfaces of zero longitudinal velocity for the straight leading edge, compared with the case λ * = 1 and A * = 0.03. In these plots the zone inside the 3D envelope corresponds to negative longitudinal velocity, i.e., recirculating fluid. It appears that one of the effects of the bumps is to modify the topology of the separation zone along the spanwise direction. In Figure 6(b) this effect takes the form of two localized regions downstream of each bumps. Although the leading-edge bumps affect the shape of the recirculating zone, its global size is nearly the same for both cases.

In Figure 7 the recirculation zone for two other cases, including the one characterized by the minimum force coefficients (λ * = 1.0 and A * = 0.07) are shown. The major difference observed between the latter case and all the others is the fact that the boundary layer is partially attached at the spanwise location corresponding to the tip of the bumps (a behaviour that does not occur for other wavelengths). This partial reattachment was also observed by Weber et al. [START_REF] Weber | Computational evaluation of the performance of lifting surfaces with leading edge protuberances[END_REF]. Conversely, the global size of the recirculating zone is strongly increased, compared to the baseline case. In this particular regime, the wake is completely dominated by the presence of two strong streamwise rolls, as displayed in Figure 8; these rolls do not oscillate, as demonstrated by the corresponding spectrum of Figure 5. Also, as shown in Figure 4, these structures have an "anti-lift" action due to the recirculation zone created by the rolls, which pushes the fluid downwards on the suction side of the airfoil, thus counterbalancing the lift force. For other wavelengths these effects are still present but with a smaller impact. Figure 7(a) shows that for the same amplitude of the bumps but λ * = 0.5 the topology of the flow is quite different and the flow is fully separated along the span. This configuration also exhibits higher lift fluctuations, as shown in Figure 3(c).

For A * < 0.07, the streamwise rolls are still present but their intensity appears to be too weak to stabilize completely the separation region. Let us mention that in experimental studies at higher Reynolds numbers [START_REF] Johari | Effects of leading-edge protuberances on airfoil performance[END_REF][START_REF] Stanway | Hydrodynamic effects of leading-edge tubercles on control surfaces and in flapping foil propulsion[END_REF], PIV measurements show the trace of these streamwise rolls. Fish et al. [START_REF] Fish | The humpback whales flipper: Application of bio-inspired tubercle technology[END_REF] suggest a similar scenario claiming that the pairs of counter-rotating vortices generated on each side of the bump are generated through a mechanism similar to the one observed on a delta wing.

To shed more light on the flow dynamics, in Figure 8 we show instantaneous flow visualizations of two cases, one (A * = 0.07, λ * = 1) in which we observe a significant drag/lift reduction and another (A * = 0.03, λ * = 1) in which streamwise rollers are not strong enough to modify the wake. In both cases we show isosurfaces of the second invariant of the gradient of the velocity fluctuation Q ′ defined by:

Q ′ = - 1 2 ∂u ′ i ∂x j ∂u ′ j ∂x i , (10)
where u ′ i = u iu i z , and • z represents averaging in the spanwise direction. The second invariant of the velocity gradient tensor has been shown to be effective in visualizing the coherent structures in turbulence [START_REF] Hunt | Eddies, streams, and convergence zones in turbulent flows, studying turbulence using numerical simulation databases[END_REF][START_REF] Dubief | On coherent-vortex identification in turbulence[END_REF]. In the same figure, we also show the velocity fluctuation vectors in a horizontal plane near the leading edge, and the velocity fluctuations and vertical vorticity along a spanwise line in correspondance with the mean leading edge at an height of ≃ 0.15 c over the mean wing surface.

As expected, in both cases we observe that the velocity is lower near the crests and higher in the valleys. The fluctuating velocity forms an array of alternating jets pointing forwards and backwards. Correspondingly, regions of alternating shear (reflected by the behaviour of the normal vorticity fluctuactions) are formed, with maximum shear corresponding to inflection points in the velocity-fluctuation profile. Under those circumstances, it is expected that a Kevin-Helmholtz instability should produce rolls of vorticity pointing in the direction normal to the surface in the leading edge region. Those vortical tubes are visible in Figure 8, where they appear in form of vorticity tongues emanating from the near-wall region (in Figure 8(b) we have used arrows to indicate their positions). In particular, the vorticity appears to roll up into a vortex tube which is initially vertical, and then immediately tilted into the streamwise direction, as it is advected by the mean stream. This Kelvin-Helmholtz-driven scenario suggests a completely different mechanism for the generation of the streamwise rollers than the "delta wing" one proposed by other authors [START_REF] Stanway | Hydrodynamic effects of leading-edge tubercles on control surfaces and in flapping foil propulsion[END_REF][START_REF] Fish | The humpback whales flipper: Application of bio-inspired tubercle technology[END_REF] and is similar to conceptual models for the generation of quasi streamwise vortices from the destabilization of wall streaks in near wall turbulence [START_REF] Kawahara | Linear instability of a corrugated vortex sheet -a model for streak instability[END_REF].

Additionally, the action of the streamwise rolls generated explains the local reattachment of the boundary layer over the bumps seen on Figure 7(a), (i.e., in the case of A * = 0.07, λ * = 1). As indicated by the vectors of the fluctuating velocity field in Figure 8(b), the flow is accelerated downstream of the bumps, which is not the case in Figure 8(a). This is related to the presence of the counter-rotating streamwise vortices which tend to accelerate the fluid in between them by adding momentum, and thus confer to the fluid the extra amount of kinetic energy necessary to prevent local separation. This analysis, obtained in a massively separated flow configuration, is in good agreement with the one proposed recently in [START_REF] Fish | The humpback whales flipper: Application of bio-inspired tubercle technology[END_REF] to explain the delay of stall angle observed in the presence of leading edge bumps.

Conclusions

To assess the fundamental effects induced by the geometrical modulation of the leading edge of an infinite wing, a numerical investigation has been performed on a NACA0020 wing profile presenting a wavy leading edge, at low Reynolds number (Re = 800), and in a configuration of massively separated boundary layer (α = 20 o ). The waviness was introduced using a sinusoidal variation along the span of wavelength λ * and amplitude A * , and periodic boundary conditions along the same direction were used, for both fluid and wing geometry.

A parametric study on the effect of the wavelength and the amplitude of the leading-edge wave has been carried out. It has been found that a certain choice of the two parameters modify drastically the topology of the wake: λ ≃ 1 c, and an amplitude of the bump equals to 7 % of the chord. In this configuration, the shedding regime disappears, the flow is dominated by streamwise structures generated by the bumps and the boundary layer is partially attached to the wall in correspondance with the crests positions.

The influence of these streamwise structures on the boundary layer separation has been studied and a physical mechanism identified to explain the local reattachment observed downstream of the bumps. Moreover, a Kelvin-Helmholtz-like instability driven by the spanwise modulation of the streamwise velocity profile is proposed to be at the origin of the generation of the streamwise vortices which control the boundary layer separation.

In the future, higher Reynolds numbers, close to those observed in nature, and at least large enough to achieve fully turbulent regimes over the airfoil, need to be considered to explain the moderate increase in the lift coefficient observed in the experiments. Further investigations are also required to support the conceptual model proposed in this work to explain the generation of the streamwise rolls, and their subsequent effects on the boundary layer separation. Finally, it could be of interest to study the effects on heterogeneous bumps over the span, with different amplitudes along the span and decreasing wavelength near the tip of the wing, as illustrated on Figure 1.

Figure 1 .

 1 Figure 1. The Humpback whale jumping out of water. The pectoral flipper clearly exhibits tubercles on the leading edge.

Figure 2 .

 2 Figure 2. Model of the airfoil geometry. (a) Set of Lagrangian points plotted with the underlying Eulerian fluid mesh. (b) Undulated leading edge, showing the parameters of control λ and A. The span length is equal to 2c.

Figure 3 .

 3 Figure 3. Contours of the performance ratios varying the amplitude A * and the wavelength λ * . (a) R 1 ; (b) R 2 ; (c) R 3 .

Figure 4 .Figure 5 .

 45 Figure 4. (a) Drag coefficient and (b) lift coefficient for λ * = 1, varying the amplitude of the leading-edge bumps. Plain line: baseline case (A * = 0). Dashed line: A * = 0.03. Dotted line: A * = 0.05. Thick line: A * = 0.07.

Figure 6 .Figure 7 .

 67 Figure 6. Recirculating regions shown by isosurfaces of zero longitudinal velocity. (a) Baseline case; (b) A * = 0.03 and λ * = 1.

Figure 8 .

 8 Figure 8. Isosurfaces of Q ′ and velocity fluctuation vectors in the plane y = 0.05c (front view), and profiles of u ′ (solid line) and ωy (dashed line) along a line parallel to the leading edge. (a) A * = 0.07, λ * = 1; (b) A * = 0.03, λ * = 1.


Martin Weiler

Mitsuhiko Miyazaki

Shun-Ichi Géraldine Féraud

Claude Ishiuchi

Christophe Dedonder

Masaaki Jouvet

Fujii
email: mfujii@res.titech.ac.jp

Shun-Ichi Ishiuchi

Claude Dedonder

Christophe Jouvet
email: christophe.jouvet@univ-amu.fr

Masaaki Fujii

Unusual behavior in the first excited state lifetime of catechol

2

ABSTRACT: We are presenting vibrationally selective pump-probe measurements of the first electronic excited-state (ππ*) lifetime of jet-cooled neutral catechol (1,2-dihydroxybenzene). The lifetime of the 0-0 transition is very short (7 ps) as rationalized by the small ππ*/πσ* gap calculated. However the lifetimes implying higher out-of-plane vibrational levels are longer (~11 ps). This emphasizes the role of the out-of-plane vibration in the ππ*/πσ* coupling not only in its nature but also in the number of quanta.

TOC GRAPHICS

The 0-0 transition (S 1 <-S 0 ) of catechol has the shortest lifetime The excited-state dynamics of aromatic molecules containing an heteroatom (phenol being the prototype molecule) has received considerable attention since it was stated in 2002 that the excited-state dynamics is linked to an Hydrogen atom loss mediated by the tunneling of the H atom from the bound ππ* to the πσ* state dissociative along the OH coordinate. [START_REF] Sobolewski | Excited-state hydrogen detachment and hydrogen transfer driven by repulsive πσ* states: A new paradigm for nonradiative decay in aromatic biomolecules[END_REF] This statement, which was already based on clear experimental [START_REF] Grégoire | Has the Excited State Proton Transfer Ever Been Observed in Phenol-(NH 3 ) n Molecular Clusters?[END_REF][START_REF] Grégoire | Picosecond hydrogen transfer in the phenol-(NH 3 ) n=1-3 excited state[END_REF][START_REF] Pino | Intracluster hydrogen transfer followed by dissociation in the phenol-(NH 3 ) 3 excited state: PhOH(S 1 )-(NH 3 ) 3 →PhO • +(NH 4 )(NH 3 ) 2[END_REF] and theoretical work, [START_REF] Sobolewski | Photoinduced Electron and Proton Transfer in Phenol and Its Clusters with Water and Ammonia[END_REF] has been reinforced by many experimental evidences and calculations. [START_REF] Ashfold | πσ* excited states in molecular photochemistry[END_REF][START_REF] Dixon | Tunnelling under a conical intersection: Application to the product vibrational state distributions in the UV photodissociation of phenols[END_REF][START_REF] King | Vibrational energy redistribution in catechol during ultraviolet photolysis[END_REF][START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF][START_REF] Roberts | Direct Observation of Hydrogen Tunneling Dynamics in Photoexcited Phenol[END_REF][START_REF] Ishiuchi | Hydrogen transfer in photoexcited phenol/ammonia clusters by UV-IR-UV ion dip spectroscopy and ab initio molecular orbital calculations. I. Electronic transitions[END_REF][START_REF] Ishiuchi | Four-color hole burning spectra of phenol/ammonia 1:3 and 1:4 clusters[END_REF][START_REF] Vieuxmaire | Ab initio characterization of the conical intersections involved in the photochemistry of phenol[END_REF] In particular, a clear correlation between the ππ* and πσ* energy gap and the ππ* excited-state lifetime has been measured recently: [START_REF] Pino | Excited state hydrogen transfer dynamics in substituted phenols and their complexes with ammonia: ππ*-πσ* energy gap propensity and orthosubstitution effect[END_REF] the smaller the gap, the shorter the lifetime. It should be mentioned that in C s symmetry (to which ground state phenol and catechol belong) the ππ* (A') and πσ* (A") are not of the same symmetry, the former being A' and the latter A", so that they are not electronically coupled, the coupling occurring through out-of-plane vibrations. [START_REF] Dixon | Tunnelling under a conical intersection: Application to the product vibrational state distributions in the UV photodissociation of phenols[END_REF][START_REF] Vieuxmaire | Ab initio characterization of the conical intersections involved in the photochemistry of phenol[END_REF] The catechol (1,2-dihydroxybenzene) case is quite interesting since it contains two adjacent hydroxyl groups linked by a hydrogen bond. Some unexpected and complicated behavior has been recorded in the product energy distribution. [START_REF] King | Vibrational energy redistribution in catechol during ultraviolet photolysis[END_REF] In particular, the energy distribution in the product obtained in detecting the kinetic energy of the H atom is remarkably dependent on the initial S 1 out-of-plane vibration. [START_REF] Nix | High resolution photofragment translational spectroscopy studies of the near ultraviolet photolysis of phenol[END_REF] The excited-state vibrational spectroscopy of the molecule is quite complex and the interpretation has been subject of controversy. [START_REF] Bu Rgi | O-H torsional vibrations in the S 0 and S 1 states of catechol[END_REF][START_REF] Gerhards | Structure and vibrations of catechol and catechol H 2 O(D 2 O) in the S 0 and S 1 state[END_REF] The agreement is that the low vibrational modes observed are due to out-of-plane OH vibrations and that the excited state has lost the C s symmetry of the ground state.

A 12 ps excited-state lifetime has been measured with a vibrational excess energy in S 1 of ൎ1800 cm -1 . [START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF] This short lifetime as compared to that of phenol was assigned to the loss of planarity and to the small calculated ππ*/πσ* energy gap (0.29 eV). It was also measured that the S 1 lifetimes of resorcinol (1,3-dihydroxybenzene) and of hydroquinone (1,4-dihydroxybenzene) are in nanosecond range (> 1 ns and 0.43 ns respectively) [START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF] similarly to phenol. [START_REF] Pino | Excited state hydrogen transfer dynamics in substituted phenols and their complexes with ammonia: ππ*-πσ* energy gap propensity and orthosubstitution effect[END_REF] These isomers have a larger ππ*/πσ* energy gap (0.51 eV for resorcinol and 0.55 eV for hydroquinone) than catechol.

Very recently, ultrafast time-resolved velocity map ion imaging and time-resolved ionyield experiments allowed a very detailed analysis of the electronic state relaxation dynamics of photoexcited catechol. [START_REF] Chatterley | Manipulating dynamics with chemical structure: probing vibrationallyenhanced tunnelling in photoexcited catechol[END_REF] In this nice experiment, the excited S 1 lifetimes have been recorded on a very broad energy scale from S 1 (v = 0) to 5000 cm -1 excess energy. Two components have been observed in the decay: one in the 500-900 fs range and one in the 8-10 ps range. The first decay has been ascribed to the IVR in S 1 and the second one to the S 1 (ππ*)/S 2 (πσ*) H atom tunneling.

The observation of the fast component on the v = 0 level of S 1 is somewhat surprising but it is clear that owing to the laser bandwidth (around 500 cm -1 ), S 1 (v = 0) cannot be excited selectively.

The catechol is a chromophore involved in many neurotransmitters and the conformational landscapes of these catecholamines have been recently unraveled in gas phase (supersonic jet) through UV/UV and UV/IR hole burning spectroscopy. [START_REF] Ishiuchi | Conformational reduction of DOPA in the gas phase studied by laser desorption supersonic jet laser spectroscopy[END_REF][START_REF] Mitsuda | Evidence for Catechol Ring-Induced Conformational Restriction in Neurotransmitters[END_REF] It has been found that the number of isomers for neurotransmitters containing catechol is much less than for the other chromophores (phenol, resorcinol). [START_REF] Ishiuchi | Conformational reduction of DOPA in the gas phase studied by laser desorption supersonic jet laser spectroscopy[END_REF][START_REF] Mitsuda | Evidence for Catechol Ring-Induced Conformational Restriction in Neurotransmitters[END_REF] The possibility of a very short lifetime of some isomers preventing their detection through REMPI method has to be investigated and the measurements of the catechol lifetime is the first step.

We present picosecond pump-probe experiment performed on catechol to see if the number of vibrational quanta observed in the internal distribution of the H-loss product [START_REF] King | Vibrational energy redistribution in catechol during ultraviolet photolysis[END_REF] is also observed in the lifetime of the initially excited-state by measuring the excited lifetime of selectively excited vibrational levels.

The origin band at 35 646 cm -1 (280.54 nm) is in agreement with what was measured before (35 649 cm -1 ). [START_REF] Gerhards | Structure and vibrations of catechol and catechol H 2 O(D 2 O) in the S 0 and S 1 state[END_REF] A low frequency (113 cm -1 ) vibrational progression is observed and was previously characterized. [START_REF] Bu Rgi | O-H torsional vibrations in the S 0 and S 1 states of catechol[END_REF][START_REF] Gerhards | Structure and vibrations of catechol and catechol H 2 O(D 2 O) in the S 0 and S 1 state[END_REF] The lifetimes have been recorded for selected vibrational levels in this paper. The transient signals for the origin and v=1 of the torsional mode τ 1 8,17 at 113 cm -1 are presented in Figure 1 and the fitted lifetime values (from the convolution of a monoexponential decay with 3 ps Gaussian laser width) are gathered in Table 1. Except for the 0-0 transition, all the lifetimes are the same (11 ps) within the experimental errors. They are similar to the ones observed at high vibrational excess energy (ൎ1800 cm -1 ) [START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF] and by Chatterley et al. on a very broad energy scale. [START_REF] Chatterley | Manipulating dynamics with chemical structure: probing vibrationallyenhanced tunnelling in photoexcited catechol[END_REF] The S 1 (v = 0) transition has a significantly shorter lifetime (7 ps) than the excited levels which is not so common. No evidence of a shorter lifetime component could be found, but the laser temporal resolution (3 ps) may not allow such detection. Pump energy (cm -1 ) Pump energy above origin (cm -1 )

Assignment (from Gerhards et al. [START_REF] Gerhards | Structure and vibrations of catechol and catechol H 2 O(D 2 O) in the S 0 and S 1 state[END_REF] Table 1. Lifetimes of the origin and of different vibrational levels in S 1 of jet-cooled catechol.

In agreement with previous calculations, [START_REF] King | Vibrational energy redistribution in catechol during ultraviolet photolysis[END_REF][START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF][START_REF] Chatterley | Manipulating dynamics with chemical structure: probing vibrationallyenhanced tunnelling in photoexcited catechol[END_REF] the ππ*/πσ* energy gap at the ground state geometry is very small (0.07 eV) as compared to other phenolic systems. [START_REF] Pino | Excited state hydrogen transfer dynamics in substituted phenols and their complexes with ammonia: ππ*-πσ* energy gap propensity and orthosubstitution effect[END_REF] The calculations agree quite well with the experimental data (4.42 eV) for the adiabatic S 1 <-S 0 transition, calculated at 4.34 eV (including the zero-point energy correction). The optimized S 1 excited-state is not planar as already stated in previous papers, [START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF][START_REF] Bu Rgi | O-H torsional vibrations in the S 0 and S 1 states of catechol[END_REF][START_REF] Gerhards | Structure and vibrations of catechol and catechol H 2 O(D 2 O) in the S 0 and S 1 state[END_REF] and it is not purely of ππ* character but has a mixed ππ*-πσ* character. The out of plane deformation involves mainly the both OH (the oxygen are out of plane) but also the H nearby the OH. The carbon skeleton is also out of plane deformed.

(see supplementary material the geometries and the vibrations can be viualised using molden or molekel software). The deformation is quite important so that the Cs symmetry rules are not proper to describe the system and that σ and π nomenclature have to taken as reminiscent of the ground state symmetry but do not apply for the excited state. The same is true for the second excited state, which has a mixed πσ*-ππ* character.

The vibrational modes have been calculated, and the active vibrational modes cannot be assigned only to the OH out-of-plane bending coordinate since the carbon atoms of the aromatic part are strongly involved in the dynamics. The lowest active frequency mode can be seen in In the insert, the spectrum recorded with the picosecond laser shows that the spectral resolution is good enough to excite selectively v=0 and v=1 of the τ 1 mode (observed at 113 cm -1 ).

In the low energy part of the spectrum, the calculated lowest frequency τ 1 (137 cm -1 ) which is responsible for the progression is in reasonable agreement with the observed one (113 cm -1 ). The simulated Franck-Condon envelope is not badly reproduced. As usual the CC2 method tends to increase too much the out-of-plane deformation underestimating the Franck-Condon factor of the 0-0 transition. The barrier calculated for the "inversion coordinate" i.e., when going from the optimized non planar S 1 structure to the planar one, can be obtained by optimizing the ππ* state in C s symmetry. We found a barrier of 240 cm -1 which is probably overestimated since the calculated frequency is 15 % higher than the experimental one.

Due to the small ππ*/πσ* energy gap and to the absence of symmetry in the excited-state, one expects a very fast H loss through the ππ*/πσ* tunneling 1 i.e., a short lifetime as observed here and in other experiments. [START_REF] Livingstone | Time-resolved photoelectron imaging of excited state relaxation dynamics in phenol , catechol , resorcinol , and hydroquinone[END_REF][START_REF] Chatterley | Manipulating dynamics with chemical structure: probing vibrationallyenhanced tunnelling in photoexcited catechol[END_REF] The real unexpected result of the ps pump-probe experiment is coming from the fact that the 0-0 transition has a shorter lifetime than higher vibrational levels. This is not common and has not been observed in the phenol. We can tentatively explain this using a simplified scheme (Figure 3) similar to that of Gerhards et al. [START_REF] Gerhards | Structure and vibrations of catechol and catechol H 2 O(D 2 O) in the S 0 and S 1 state[END_REF] In this picture, the vibrational wavefunctions are presented in a potential along an out-ofplane torsion coordinate, assuming that the barrier is lying between the v=0 and v=1 excited vibronic levels, as calculated. The ππ*/πσ* coupling strongly depends on the out-of-plane deformation, the coupling being essentially zero at the planar geometry. [START_REF] Vieuxmaire | Ab initio characterization of the conical intersections involved in the photochemistry of phenol[END_REF] The v=0 level is below the barrier to planarity and thus the vibrational wavefunction is essentially 0 at the planar geometry. Because geometry at v=0 level is always bent, the coupling is quite large and the lifetime should be short. For other vibrational levels the probability of presence in the planar geometry is not zero and the average of the coupling along the vibration is smaller ("the system spends more time in the planar geometry") thus the lifetime should be longer. resonance and scanning the delay of the probe ps laser set at 310 nm, by steps of 0.08 ps.

Ab initio calculations have been performed with the TURBOMOLE program package, [START_REF] Ahlrichs | Electronic structure calculations on workstation computers: The program system turbomole[END_REF] making use of the resolution-of-the-identity (RI) approximation for the evaluation of the electronrepulsion integrals. [START_REF] Weigend | RI-MP2: optimized auxiliary basis sets and demonstration efficiency[END_REF] The equilibrium geometry of ground electronic state (S 0 ) has been determined at the RI-MP2 24 level. Excitation energies and equilibrium geometry of the lowest excited singlet state (S 1 ) have been determined at the RI-CC2 level. [START_REF] Christiansen | The 2nd-Order Approximate Coupled-Cluster Singles and Doubles Model CC2[END_REF][START_REF] Hattig | Geometry optimizations with the coupled-cluster model CC2 using the resolution identity approximation[END_REF] Calculations were performed with the correlation-consistent polarized valence triple-zeta aug-cc-pVTZ basis. [START_REF] Woon | Gaussian basis sets for use in correlated molecular calculations . III . The atoms aluminum through argon[END_REF] The vibrations in the ground and excited-states have been calculated and the vibronic spectra simulated using PGOPHER software [START_REF] Western | a Program for Simulating Rotational Structure[END_REF] for the Frank-Condon analysis.

  vibrationally selective pump probe, Hydrogen transfer, non adiabatic coupling, fast decay, aromatic molecules

Figure 1 .

 1 Figure 1.Pump-probe delay of catechol when exciting S 1 origin (higher panel, thin line) and the vibronic band 113 cm -1 above origin (lower panel, thin line). Their fit using a monoexponential decay convoluted by a Gaussian laser width of 3 ps is represented by thick lines.

Figure 2 .Figure 2 .

 22 Figure 2. (see supplementary material)

Figure 3 .

 3 Figure 3. artistic view of the wavefunction for the out-of-plane torsional coordinate. For v=0 the probability of presence is non-zero only for the bent geometry, which is not the case for v≥1, so that the ππ*/πσ* coupling is larger for v=0 than for v≥1, possibly explaining the shorter lifetime for the v=0 level adapted from reference 17.

  It would be interesting to see if theoretical calculations can reproduce this experimental result. The excited-state lifetimes of different vibronic levels in the first excited state of jetcooled catechol have been measured with ps pump-probe experiment. The lifetimes are very short (11 ps) and even shorter (7 ps) on the 0-0 transition. This can be understood in terms of variation of vibronic coupling between the bound ππ* and dissociative πσ* states due to the absence of C s symmetry in the excited state. The question of the evolution of this lifetime in substituted catechols and their conformers remains open.EXPERIMENTAL TECHNIQUES AND CALCULATIONSCold catechol is obtained in a molecular beam and investigated using nanosecond and picosecond laser systems, featuring pulse lengths of 5 ns and 3 ps, respectively. Details of the experimental setup, molecule generation, spectroscopic methods, and employed laser systems have been given elsewhere.[START_REF] Ishiuchi | IR signature of the photoionization-induced hydrophobic-->hydrophilic site switching in phenol-Ar n clusters[END_REF] Briefly, catechol in a supersonic expansion is seeded in 3 bar of Helium and ionized in the extraction region of a time-of-flight mass spectrometer via resonant two-color two-photon ionization (1 + 1' REMPI). The first UV photon excites the molecules into the S 1 electronic state, whereas a second UV photon is used for soft ionization with low excess energy. REMPI spectra are obtained by scanning ns UV laser through the vibronic resonances of the S 1 state. Lifetime measurement spectra are obtained by fixing the UV wavelength of the ps laser to a particular S 1

AUTHOR INFORMATION AUTHOR INFORMATION

Corresponding Authors

* Email: christophe.jouvet@univ-amu.fr * Email: mfujii@res.titech.ac.jp

Notes

The authors declare no competing financial interests.

ACKNOWLEDGMENT

This study was supported in part by a Grant-in-Aid for Scientific Research KAKENHI in the priority area (477) "Molecular Science for Supra Functional Systems" from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) Japan and the France-Japan Collaboration Program (SAKURA) France, and the Japan Society for Promotion of Science. MW thanks for Japan Society for Promotion of Science fellowship. We acknowledge the use of the computing facility cluster GMPCS of the LUMAT federation (FR LUMAT 2764)


Y F Zhao

N Perry

H Andriankaja

INFORMATION STRUCTURE FOR MANUFACTURING SUSTAINABILITY ASSESSMENT: STEP FOR LCA

Keywords: manufacturing process, information model, STEP

Environmental performances need product, processes and life cycle modeling and evaluations. There is direct need of assessment tools to monitor and estimate environmental impact generated by different types of manufacturing processes. Indeed, eco design and products optimizations can be done with the manufacturing process choices. But the manufacturing systems are very complex and they are driven by physicals laws that are heavy to model. This research proposes a manufacturing informatics framework for the assessment of manufacturing sustainability based on an EXPRESS information model developed to represent sustainability information. It is the first step of association of sustainability information with product design specification. In the next phase of research, investigation will be conducted to integrate sustainability information model and existing standardized product design model ISO 10303 AP 242.

INTRODUCTION

The manufacturing industry is often cited as the cause of many environmental and social problems, yet it is acknowledged as the main mechanism for change through economic growth (Baldwin 05). Industry is confronted with the challenge of designing sustainable products and manufacturing processes. Many sustainable development strategies have been proposed by government agencies and academia since the late 1990s. Research in sustainable manufacturing field can be categorized into the following four main themes [START_REF] Mayyas | Design for sustainability in automotive industry: A comprehensive review[END_REF][START_REF] Ashby | Life-cycle cost analysis: Aluminum versus steel in passenger cars[END_REF][START_REF] Barker | Organizing reuse managing the process of design for remanufacture (DFR)[END_REF][START_REF] Govetto | End-of-life vehicle recycling based on disassembly. Determining the Environmental Impact of Disposal[END_REF][START_REF] Perry | Improving design for recycling -Application to composites[END_REF][START_REF] Vijayaraghavan | Automated energy monitoring of machine tools[END_REF]): a) Life-Cycle Assessment (LCA), b) design-for-X principles and design for sustainability, c) end-of-life studies, and d) energy efficiency monitoring and studies. But it remains a lack of manufacturing system and sustainability information integration hampers the widespread adaptation of the best sustainability practices in the manufacturing industry.

INFORMATION MODELS OF PRODUCT DESIGN AND MANUFACTURING PROCESSES

In recent years, information technology has become increasingly important in the manufacturing enterprise. Effective information sharing and exchange among computer systems throughout a product's life cycle has been a critical issue (ITM 1995, Lee 1999). An information model is a representation of concepts, relationships, constrains, rules, and operations to specify data semantics for a chosen domain of discourse. Amongst these information models, the STEP and STEP-NC information models are the most advanced and standardized ones, developed by ISO committee to provide the basis for product design, machining standardization, and integration of part inspections with machining. The information is modeled in EXPRESS language. STEP Application Protocol (AP) 203 editions 1 & 2 (ISO10303 19941 & 2 (ISO10303 , 20071 & 2 (ISO10303 & 2009) ) (Configuration Controlled 3D Designs of Mechanical Parts and Assemblies) provides the data structures for the exchange of configuration-controlled 3D designs of mechanical parts and assemblies. AP 203 comes from ISO 10303 product data standard, and is not a data standard for configuration management of a product throughout its entire life cycle. The AP is centered on the design phase of mechanical parts and the high-level information entities.

PROPOSED MANUFACTURING INFORMATICS FRAMEWORK FOR THE INTEGRATION OF LCA AT DESIGN STAGE

In order to integrate LCA information into the product design stage, necessary information must be properly represented and associated to the product PLM information. In the proposed research, a case study was first conducted to examine what LCA information should be modeled. Composite parts give interesting examples on a simple piece such as a pedal crank developed with recycled carbon fibers for thermoset organic matrix composite as shown in Figure 3. During the design, product models have to integrate materials information such as matrix composition, reinforcement type (glass, carbon, aramid, and natural), their architecture (unidirectional, woven) and the structure composition (orientations of the different layers in the depth of the product). Other information like inserts or coating completes the bill of material. Currently in STEP AP 203, only very limited material information is modeled as shown in the following entities in table 1. To associate composite material information to a product design, the following entities were developed in this research. The new entities are written in bold-italic. Manufacturing data must also been included in order to ensure the product / material / process combined design in the case of composite parts. Promoting recycled carbon fibers eco-design [START_REF] Perry | Knowledge integration based on road mapping and conceptual framework approach to ease innovation management[END_REF] imposes to access to material data (from end of life scenario and properties) at the early stage of the design. It means developing specific data storage for material properties (depending to the recycling process) and properties models prediction (based on the different type of fibers) as illustrated on Figure 2. Figure3:Proposed information framework for integrating LCA in product design stage Currently in STEP AP 203, the limited material information, shown above in unbolt entities, is not associated to design geometry. The material information is defined only for annotation display. In order to fully integrate material information with design process, the association between product geometry and its material must be semantic. Today no integrated environment allows taking all these facets into consideration. And the next step is to compare, for an expected set of functions, the n-plet [START_REF] Mayyas | Design for sustainability in automotive industry: A comprehensive review[END_REF]) in order to optimize both technically, but also environmentally with specific focus such as material resources minimization (mix of less material and improved end of life possibilities), energy optimization and/or pollution reduction. These multi objective optimization needs to handle all these data, models, life stages, functional unit definition in a coherent environment and at the early stage of product development. STEP and STEP-NC is the most suitable information structure framework that needs to be enriched by the environmental aspects of the product/process data as shown in Figure 3. Eventually, with the proposed information framework, LCA related information is integrated to the design and process planning aspects of product development. The LCA knowledge and data accumulated throughout a product life cycle can be fed back to the product design stage to improve new product development and to reduce environmental impact.

FUTURE WORK AND CONCUSIONS

The research reported in this paper is the very beginning stage of developing an information framework to support the integration of LCA at the product design stage. A simple case study was conducted to examine what composite material information should be defined in association with product design information. The next steps of the proposed research are for the data structure level: i) identify the Product Life Phases needs in terms of environmental data versus the existing data into STEP, ii) identify a dynamic ontology of environmental data for product development, iii) structure the data and evaluation models to propose, usable simulation results for decision making, and iv) propose a STEP extension for End of Life support. For the data acquisition level: i) low level implementation (at the process level for the manufacturing phase) in order to enrich STEP with environmental concepts and start creating data base for process evaluation and optimization, ii) develop design approach and integrated tools specification based on STEP.

Sustainable product design and manufacturing is the future of today's manufacturing industry. The lack of manufacturing system and sustainability information integration hampers the widespread adaptation of the best sustainability practices in the manufacturing industry. In order to provide product designer with comprehensive material and environmental related information at the design stage, a complete information structure must be developed to associate sustainability information with product design information and product manufacturing process. This research proposes to develop such an information framework based upon STEP and STEP-NC information models. An initial case study was conducted to identify composite material information that should be defined for product design.

Figure 1 :

 1 Figure 1: pedal crank made with recycled Carbon Fibers (developed at I2M) (Perry 2012)]

Figure2:

  Figure2: Recycling End of Life scenario and product/material properties evaluation tool.

ENTITY composite_matrix ABSTRACT SUPERTYTPE OF (ONEOF (mud, cement, polymers, metals, ceramics)); name

	(Anisotropic_material,	continuous_fibre_reinforced));
	Braided_assembly,	name:	STRING;
	Coating_layer,	direction:	fibre_directions;
	Isotropic_material,	END_ENTITY;
	Laminate_table,
	Part_view_definition,
	Substance_view_definition,
	Woven_assembly
	Composite_material);
	END_TYPE;
	ENTITY Material_identification;
	material_name : STRING;
	items : SET[1:?] OF material_item_select;
	END_ENTITY;	:	STRING;
		END_ENTITY;
	TYPE material_item_select = SELECT
	(Anisotropic_material,	ENTITY composite_resin
	Braided_assembly,	ABSTRACT SUPERTYPE OF (ONEOF
	Coating_layer,	(polyester_resin, vinylester_resin,
	Isotropic_material,	epoxy_resin, shape_memory_polymer_resin));
	Laminate_table,	name:	STRING;
	Part_view_definition,	material_property: STRING;
	Substance_view_definition,	END_ENTITY;
	Woven_assembly);
	END_TYPE;	ENTITY composite_reinforcement
		ABSTRACT SUPERTYOP OF (ONEOF (glass_fibre,
	ENTITY Composite_material_identification	carbon_fibres, aramid_fibres,
	SUBTYPE OF (Material_identification);	boron_fibres));
	DERIVE	name:	STRING;
	composite_material_name : STRING :=	architecture:	fibre_architecture;
		END_ENTITY;
	SELF\Material_identification.material_name;
	END_ENTITY;	ENTITY fibre_architecture;
		ABSTRACT SUPERTYPE OF (ONEOF
	TYPE material_item_select = SELECT	(short_fibre_reinforced,

ENTITY fibre_directions; SUPERTYPE OF (ONEOF (continuous_aligned, discontinuous_aligned, discontinuous_random, unidirectional, woven)); END_ENTITY;Table 1 :

 1 STEP AP203 entities enriched with composites material information


Delphine Mathias
email: delphine.mathias@gmail.com

Aaron Thode

Jan Straley

Russ Andrews

Olivier Le Bot

Cedric Gervaise

Jerome I Mars

Cédric Gervaise

Jérôme Mars

Range-depths tracking of multiple sperm whales over large

come

Introduction

In recent years, passive acoustic methods have become increasingly widespread for monitoring the general assessment of marine environments [START_REF] Dadouchi | Automated segmentation of linear time-frequency representation of marine mammal sounds[END_REF]- [START_REF] Anagnostou | Passive aquatic listener (PAL): An adoptive underwater acoustic recording system for the marine environment[END_REF] and expanding knowledge about marine mammals vocalization repertoire, distribution and habitat characterization. In the past decade, considerable efforts have been made for this purpose using a combination of ocean science, signal processing, statistics and computational (algorithms) science [START_REF] Andr | Listening to the Deep: Live monitoring of ocean noise and cetacean acoustic signals[END_REF].

This paper is concerned with the localization and tracking of sperm whales using a vertical array comprising only two hydrophones. Indeed, passive acoustic monitoring has become an important tool to study sperm whale behavior in the Gulf Of Alaska and their interaction with longlining fisheries [START_REF] Book | Detection, Classification and Localization of Marine Mammals using passive acoustics[END_REF]- [START_REF] Mathias | Acoustic and diving behavior of sperm whales (Physeter macrocephalus) during natural and depredation foraging in the Gulf of Alaska[END_REF]. Each click event generated by a sperm whale can arrive on a hydrophone via multiple ray paths. In this paper, the ray path that arrives first on a hydrophone will be called the primary path, and other ray path arrivals that arise from the same click event are called secondary paths, or multipath.

Most methods developed for localizing marine mammals use wide-baseline hydrophone arrays and the time-difference-of-arrival (TDOA) of a sound on pairs of hydrophones [START_REF] Mathias | Acoustic tracking of sperm whales in the Gulf of Alaska using a two-element vertical array and tags[END_REF]- [START_REF] Tiemann | Localization of marine mammals near Hawaii using an acoustic propagation model[END_REF]. Methods are often based on ray-trace acoustic propagation modeling and exploit multipath arrival information from recorded sperm whale clicks. The algorithm compares the arrival pattern from a sperm whale click to range and depth dependent modeled arrival patterns in order to estimate whale location. The technique can account for waveguide propagation physics like interaction with range-dependent bathymetry and ray refraction. Tieman et al. [START_REF] Glotin | Whale cocktail party: real-time multiple tracking and signal analyses[END_REF] managed to track a sperm whale in three dimensions using only one acoustic sensor and a model of the azimuthally dependent bathymetry.

When multiple whales are simultaneously clicking, the biggest challenge is to arrange clicks into separate click-trains corresponding to individual whales, and then classify clicks as primary paths and multipaths. In the past decade several authors proposed algorithms for separating multipaths from the primary click-train, either on single hydrophones [START_REF] Glotin | Whale cocktail party: real-time multiple tracking and signal analyses[END_REF] two-hydrophone arrays [START_REF] Tiemann | Three-dimensional localization of sperm whales using a single hydrophone[END_REF] or wide-baseline acoustic arrays [START_REF] Baggenstoss | Separation of sperm whale click-trains for multipath rejection[END_REF]. These algorithms exploit the slowly varying multipath structure of individual whales or the slowly varying features of clicks within a train (such as waveform, power). Recently, a few papers discussed how sparse coding can be an effective technique for solving the multiple-marine mammal tracking problem [START_REF] Bahl | Towards identification of sperm whales from their vocalizations[END_REF]- [START_REF] Doh | Dcomposition et analyse par codage parcimonieux des chants de ctacs[END_REF]. Sparse coding seems to be a promising alternative to usual time-frequency feature analysis.

To our knowledge the long-range tracking of multiple whales on a single deployment has not been performed yet. For many applications, the deployment of several hydrophones is impractical and too expensive. Here we discuss how a two-element vertical was used to track the range of multiple whales up to a 35 km range over a 3-day period and how the method could be automated using rhythmic properties of click-trains and sparse coding.

Semi-Automated tracking of multiple whales

A two-element vertical array deployed at the sound speed minimum was used to track sperm whales in the Gulf of Alaska between 15 and 17 August 2010. The vertical arrival angles and relative arrival times of multiple refracted and surface-reflected ray paths contain enough information for range-depth tracking without knowledge of the bottom bathymetry. A ray-tracing program was used to model the acoustic travel times from each candidate source location, using a measured sound speed profile. By comparing modeled and measured time lags and vertical angles, an ambiguity surface was created, displaying the best-fit whale position. A tagged sperm was tracked up to 35 km range under Beaufort 3 conditions, using satellite tag data to independently verify tracking estimates [START_REF] Mathias | Acoustic tracking of sperm whales in the Gulf of Alaska using a two-element vertical array and tags[END_REF]. The technique also permitted to measure the drift of multiple whales away from the vertical array. The method and results are described in detail in Mathias et al. [START_REF] Mathias | Acoustic tracking of sperm whales in the Gulf of Alaska using a two-element vertical array and tags[END_REF].

However we were not able to automate the tracking process in the presence of as much as six whales simultaneously vocalizing. Techniques described above to separate click-trains such as the crosscorrelation or a rhythmic analysis failed in our case, because of the high number of multipaths received at the hydrophones produced by whales at various ranges. 3 Towards better localization of sperm whales using rhythmic properties of click-trains and sparse coding

Relative time [s]

Localizing a sperm whale using a two-element vertical array requires measuring the relative arrival times of at least two ray paths (the primary path and a multipath) on both hydrophones [START_REF] Mathias | Acoustic tracking of sperm whales in the Gulf of Alaska using a two-element vertical array and tags[END_REF]. Therefore, we need to de-interleave click-trains and associate a primary path click-train with a multipath click-train. Two approaches seem promising for performing this task when three whales or more are vocalizing simultaneously. The first approach takes advantage of the two hydrophone array configuration and use the arrival angle information as an additional source of information for grouping sperm whale clicks into trains associated with a given whale and propagation path. Therefore on-going work includes combining the arrival angle information with an algorithm developed by Le Bot et al. [START_REF] Bot | Odontocete click train deinterleaving using a single hydrophone and rhythm analysis[END_REF] that uses the rhythmic properties of odontocete click-trains to separate interleaved click-trains. The algorithm only uses the time of arrival of each click and a complex-autocorrelation function to compute a histogram that exhibits peaks at inter-click intervals (ICI) corresponding to the interleaved click trains, while suppressing harmonics due to ICI multiples. This complex autocorrelation is calculated in a window sliding along the click train leading to a time-ICI representation, which is thresholded to detect the different interleaved click trains. This sequential search could use some complementary features such as the click arrival angle, level or its frequency content.

The second approach is based on sparse coding [START_REF] Glotin | Sparse coding for large scale bioacoustic similarity function improved by multiscale scattering[END_REF]- [START_REF] Paris | Physeter catodon localization by sparse coding[END_REF] and recent publications on the application of this technique on marine mammal sounds [START_REF] Bahl | Towards identification of sperm whales from their vocalizations[END_REF]- [START_REF] Doh | Dcomposition et analyse par codage parcimonieux des chants de ctacs[END_REF]. We propose that a sparse transform of the clicks in the time-frequency domain can help determine the stable components between multipaths belonging to an individual and a given propagation path. In order to reduce the signal dimension for more efficient computation, sets of Mel Frequency Cepstral Coefficients (MFCC) can be computed and a dictionary of features can be generated. Any click detected on the hydrophone can therefore be represented in this space of reduced dimension. The similary between each projected click can be computed using the cosine similarity measure for example.

FrequencyFigure 1 :

 1 Figure 1: Spectrogram of multipaths produced by two sperm whales on 15 August 2010 at 19:01:30 and recorded on vertical array top and bottom hydrophones

  [START_REF] Glotin | Sparse coding for large scale bioacoustic similarity function improved by multiscale scattering[END_REF] showed that this technique worked for tracking the sounds produced by the same minke whale during 30 minutes. It is also possible to work directly on the spectrogram to select areas of interest. A specific algorithm by Dadouchi et al.[START_REF] Bahl | Towards identification of sperm whales from their vocalizations[END_REF]-[START_REF] Doh | Dcomposition et analyse par codage parcimonieux des chants de ctacs[END_REF] has been developed to detect click and whistles. Based on a two-stage methodology, this algorithm estimates the instantaneous frequency law of non-linear frequency modulations under several constraints (high resolution estimation, ability to cope with multiple overlapping and/or close signals in the time-frequency plane). The first step of the methodology is applied on the square modulus of any linear time-frequency representation, and aims at detecting the time-frequency support of the signals of interest under probabilistic models. A Chisquared model is used to do the detection of time-frequency bins hosting signal, a Poisson model for the gathering of detected bins into regions of interest (RoIs). Once the RoIs are detected, a high resolution estimator using local polynomial frequency law estimation and phase continuity criteria is used to link local approximation to get a whole estimate of the instantaneous frequency law.

Acknowledgments

This work was supported by the North Pacific Research Board, the Alaska SeaLife Center, ONR, NOAA and ANR-12-ASTR-0021-03 "MER CALME"]

[Work supported by the North Pacific Research Board, the Alaska SeaLife Center, ONR, NOAA and ANR-12-ASTR-0021-03 "MER CALME"]


Adel Olabi
email: adel.olabi@ensam.eu

Mohamed Damak
email: m.damak@geomnia.eu

Richard Bearee

Olivier Gibaru

Stephane Leleu

Improving the Accuracy of Industrial Robots by offline Compensation of Joints Errors

The use of industrial robots in many fields of industry like prototyping, pre-machining and end milling is limited because of their poor accuracy. Robot joints are mainly responsible for this poor accuracy. The flexibility of robots joints and the kinematic errors in the transmission systems produce a significant error of position in the level of the end-effector. This paper presents these two types of joint errors. Identification methods are presented with experimental validation on a 6 axes industrial robot, STAUBLI RX 170 BH. An offline correction method used to improve the accuracy of this robot is validated experimentally.

I. INTRODUCTION

Industrial robots are usually used to realize industrial tasks like material handling, welding, cutting and spray painting. The mobility, flexibility and important work space of these robots allow using them in new fields of industry such as prototyping, cleaning and pre-machining of casts parts as well as end-machining of middle tolerance parts. These new applications require high level pose accuracy and to achieve a good path tracking. Unfortunately industrial robots are designed to have a good repeatability but not a good accuracy. Their repeatability ranges from 0.03 to 0.1mm for small and medium sized robots and can exceed 0.2 mm for big ones. Meanwhile the accuracy is often measured to be within several millimeters [START_REF] Damak | Vision robot based absolute accuracy measurement, calibration and uncertainty[END_REF]. This poor accuracy is caused by geometric factors, such as geometric parameters, joints offset errors and TCP definition, as well as by non-geometric factors such as compliancess, thermal effects gear, encoder resolution, gearboxes backlashes and kinematic errors.

Many fields of investigation are proposed to increase the accuracy of industrial robots like: robot calibration, process development and control system (see figure 1). Robot calibration improves the accuracy of positioning by reducing the deviation between the commanded pose and the real one.

The complete procedure of robot calibration basically consists of four stages: modeling, measurement, identification, and compensation [START_REF] Chen | Practical industrial robot zero offset calibration[END_REF]. In standard kinematic calibration, geometric errors are modeled and compensated; robot joints are assumed to be perfectly rigid [START_REF] Elatta | An overview of robot calibration[END_REF]. In non-kinematic calibration, robot joints compliances and other non geometric errors are concidered [START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF]. In [START_REF] Abele | Modeling and identification of an industrial robot for machining applications[END_REF], authors have worked on modeling the Cartesian compliance of an industrial robot according to its joints compliances in order to analyze the system's stiffness. Regarding the control field, a large number of works have been done on trajectory planning, feedback control, system compensation and feedforward control [6] [7]. The work in this paper is situated in the field of nonkinematic calibration, as shown in figure 1. We focus on two non-geometric errors: joints compliances and gearboxe kinematic errors. Different techniques could be used to compensate these errors. Typically, they can be classified in two categories:

• Model based method: in this method, each source of error must be identified and modeled by using a series of measurements of the real position of the end-effector. These models are integrated in a complex model of the robot to anticipate the positionning error in order to adjust the reference to correct this error [8] [9]. • Sensor based method: in this method the position error is measured by means of external position sensors. These sensors can be placed either in joint space or in Cartesian space. The measured error is injected then in the control loop to adjust the reference position in order to get the desired position of the end-effector [START_REF] Wang | Improving machining accuracy with robot deformation compensation[END_REF]. The sensor based method offers a good position accuracy. However its implementation on an existing robot is difficult. Many problems could be encountered like the positioning of sensors on the robot structure and their integration in the control loop besides the extra cost due to the integration of these sensors. The model based technique allows to realize correction of the robot position without any change in the robot structure or in the control loop. In this paper, the model based technique is used to compensate two joints major errors: compliance errors and gearboxes kinematic error. For each kind of error two phases are necessary: identification and modeling. The model of each error is integrated then in a trajectory planner, adapted for continous trajectories applications such as machining [11] [12], in order to realize an offline compensation of joints errors.

II. JOINTS COMPLAINCES

There are many sources of compliance in serial robots: base, joints and links, as well as an active compliance induced by the position control loop. Generally, the compliance of robot's links is negligible in comparison to joints compliance. Joints compliance consists of both axial and radial compliances. In this work, we take into account only the axial compliance, which is the dominating one. It is represented by an equivalent joint stiffness [START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]. Due to joints compliance, external forces applied on the end-effector induce elastic deformations in robot's structure. The figure 2 presents a a model of a robot joint with the positioning of the sensor before the kinematic chain. These elastic deformations are not corrected by the controller because of the positioning of sensors on the motor axis and not on the joint axis. Therefore, to improve the accuracy of robots these deformations must be compensated.

A. Identification of Joints stiffnesses

The identification methods can be classified in two categories: global and local.

1) Global method: This method is based on measurements realized in the end-effector level. An example of setting up of this method is shown in figure 3, where the external force F is measured by a force sensor. The end-effector position and the orientation deviation, due to F, is measured by an external 3D measurement system, like a laser tracker. Thus, the relationship between the force F and the displacement vector δX is given by:

F = K c δX = J(q) -T K J J(q) -1 δX, . (1)
where J(q) is the Jacobian matrix and

K J = diag [K j1 , K j2 , . . . , K j6 ] is the joint stiffness matrix.
A least square solution could be done to calculate the stiffness matrix K J .

The advantage of this method is the easiness of the installation. On the other side, it introduces many drawbacks:

• The results depend on the configuration used during the identification procedure, • The position measurement, realized in the end-effector level, integerats the deformations of the robot links, which are not perfectly stiff. (a) • The identification quality depends on the number of parameters to be identified at the same time.

F L Cable L q 1 Load(Kg) X Z W[N] (b) L Y F q 2 W[N]
To overcome these drawbacks we apply a local identification method, called axis by axis, where the stiffness of one joint is realized at each step.

2) local method: The only parameter to be identified in this approach is one joint stiffness. This stiffness is calculated by measuring the angular displacement resulted of an external torque applied on the joint. Our goal is to identify the stiffness of the six joints of our robot without disassembling it. This needs a special identification procedure to avoid the problem of stiffness coupling. This procedure changes according to the joint whose the stiffness is being identified. Where we chang the configuration of the robot for each joint stiffness identification, we measure the angular displacement by usinga Laser Tracker or an electronic level. We used a Laser Tracker for joints 1, 2 and an electronic level for joints 3,4,5 and 6.

a) Joints 1,2: As shown in figure 4, the same configuration was selected for the stiffness identification of joints 1 and 2. An external force F is applied on the robot link to produce a torque on the joint axis. An angular displacement δq, is resulted of this torque; it is measured by the Laser Tracker. This operation is done for different external forces, then the stiffness is calculated by applying a linear interpolation model as follows: F = Kδq. The same strategy, force/displacement, is applied for joints 3,4,5 and 6. However, in the cas of joint 3 identification, the application of an external force on the robot causes an angular displacement in joints 2 and 3 at the same time, as shown in figure 5. For this reason two electronic levels are necessary: one is used as a reference to measure the angular displacement of axis 2 according to the vertical direction and another one to measure the combined angular displacements of joints 2 and 3. The positions of joints 2 and 3 are measured before applying the external torque, the relative angular position between these two joints before deformation is given by δq = q 3q 2 . After applying the force F, the new values of the levels are q2 = q 2 + δq 2 and q3 = q 3 + δq 2 + δq 3 , where δq 2 and δq 3 are the angular displacement of the joints. The difference after deformation is δq = q3 -q2 . Thereby, the angular displacement of joint 3 is given by: δqδq = (q 3 -q2 ) -(q 3q 2 ) , = (q 3 + δq 2 + δq 3 + δq 2 ) -(q 2 + δq 2 ) , = δq 3 .

(

) 2
The strategy is used then to measure the stiffness of joints 4,5 and 6. Figure 6 shows the positioning of the electronic levels during the identification procedure. The table I presents the values of the joints stiffnesses of our robot, founded by following the previous procedure.

III. JOINTS KINEMATIC ERROR

The movement transmission in a robot joints is not ideal. High ratio gearbox induces a kinematic error, noted q err , which is given by the difference between the joint anglular position q out and the motor one q in , scaled by the ideal gear ratio [START_REF] Tuttle | Understanding and modeling the behavior of a harmonic drive gear transmission[END_REF] as follows:

q err = q in gear ratio -q out . ( 3
)
This kinematic error could limit the use of industrial robots for high accuracy demanding applications. This error can reach The three first joints of our robot are equipped with two stages of reduction system as shown in figure 7.

• The first level consists of two helical gears.

• The second level is a special reducer made by Staubli, similar to harmonic drives.

Thus, the kinematic error of each joint is composed of the defaults of the components of these two stages. Again, the objective is to measure the kinematic errors without disassembling the robot. Thereby, the identification is done by using measurements of the positions of robot links.

c) Identification protocol : To measure the kinematic error of a robot joint, we perform an angular displacement of this joint. Then the real position is measured by an external measurement system and compared to the rotary encoder commanded value. To avoid any dynamic effects, this operation is done by taking measurement in static positions. Thus, the identification of the kinematic error is done as follows:

• Measuring of the joint angular positions by an external system, for example a Laser Tracker, in a specified interval. • Comparing these measures q r with the joint reference q theo to calculate the position error Δq = q rq theo . This error cumulates the stiffness error q stif f and the kinematic error q gear . • Calculating the stiffness error by using robot parameters and then extract the kinematic error. The Laser Tracker target is fixed on the robot link. This link is actuated with a specified angular step. The link motion is stopped after each step in order to measure automatically the position after a short stabilization time. It has to be noted that the laser tracker measures cartesian positions, so these positions must be converted to angular positions. d) Joint 2: In this work the joint 2 of the robot is considered to validate this identification method. The selected configuration is shown in figure 8. The robot joint is rotated from q 2 = -50 • to +50 • (where q 2 =0 is the vertical position) with a step of 0.05 • . The waiting time after each step is one second before measuring the real position.

These measures are then converted to angular positions and compared to joint references for each step. The calculated position error is shown in figure 9. As mentioned before, this error cumulates:

• The stiffness error dq stif f , about many 1/100 of degree, is resulting from the joint deformation due to the residual torque applied on this joint. This residual torque is resulting from the difference between gravity and the compensation torques. In serial-type manipulators gravity effects are eliminated by systems of compensations [START_REF] Wongratanaphisan | Gravity compensation of spatial two-dof serial manipulators[END_REF]. For our robot, the gravity compensation is done in the level of joint 2 by a technique of spring suspension. . The residual torque is calculated by using our robot data in the test interval and it is shown in figue 10. This stiffness error is calculated then after having estimated the joint stiffness. • The kinematic error dq gear , about many 1/1000 of degree. The kinematic error δq gear for joint 2 is calculated by subtracting the position error δq from the stiffness error δq stif f . This kinematic error is show in figure 11.

The basic component of kinematic error for harmonic drive is due to a misalignment during assembly operation. The frequency of this component could be calculated according to drive input rotation. A typical signature of a drive error consists mainly of a fundamental harmonic occurring at a frequency of two cycles by input rotation [START_REF] Ghorbel | On the kinematic error in harmonic drive gears[END_REF]. The fundamental harmonic of the measured kinematic error according to the input rotation of the second level joint 2 transmission system confirms this explanation as shown in figure 11.

IV. OFFLINE TRAJECTORY CORRECTION

In this work we apply an offline strategy to eliminate the joints errors presented previously. This strategy is based on anticipating joints errors and correcting them by adjusting the position references of each joints as follows:

Q ref = Q theo -(ΔQ stif f + ΔQ gear ), (4)
where ΔQ stif f is a vector of angular displacements due to joints deformations. This vector is calculated by using an elasto-static model based on geometric parameters, joints stiffness, components weights and gravity centre of these components. ΔQ gear is a vector of joints kinematic errors. This correction strategy is applied in the trajectory planning stage, developed in previous works [START_REF] Olabi | Enhanced trajectory planning for machining with industrial six-axis robots[END_REF], where the joints references are modified according to the calculated errors.

A. Experimental tests

For the sake of simplicity without losing generality we only apply the offline correction strategy on joint 2 . The configuration selected for these tests is presented in figure 12.

The tests are realized without any charge on the end-effector. In these conditions joint 2 errors are dominating.

The desired path is represented by a logarithmic spiral, given by R = a.e bθ where a = 24 (mm), b = 0.15 (mm) and θ = 4 * π (rad). The spiral is placed in XY plan with Z robot =200 mm. A Laser Tracker is used to measure the trajectory of the end-effector. To highlight the contribution of the correction of each error, we compensated firstly, the stiffness error then the kinematic error compensation is added.

e) Correction tests for quasi-static trajectories: To avoid any dynamic effects due to the inertia of the robot components, a low speed was selected, S=0.03 ms -1 and an acceleration of γ=0.1 ms -2 . The same trajectory was measured three times: without taking into account joint 2 errors then, with correction of joint stiffness error and finally by considering both, stiffness and kinematic errors.

The figure 13 shows the three trajectories measured during these tests. To quantify the correction of position of the endeffector in Z direction we use, by analogy with the roughness parameters in machining, the following parameters:

• The average value R a = 1 n n k=1 |z k | • The maximum height of the profile R z , given by the difference between the maximum peak height z max and the maximum valley depth -z min . The table II shows the values of these parameters for each test. The figure 13 and the table highlight the effectiveness of the correction realized in these tests. The correction of the stiffness error reduces Ra up to 77% and Rz up to 27%. The correction of the kinematic error reduces Ra up to 32% and

V. CONCLUSION

In this paper, two types of joints errors in industrial robots, joints compliances and kinematic errors, were presented. A local identification method, called axis by axis, is used in this work to identify the joints stiffness of a six axes industrial robot. An identification protocol to measure the kinematic error of joint 2 is presented in this work as well. Then an offline correction method is developed to improve the accuracy of industrial robots. This method is validated experimentally by measuring the trajectory of the end effector with and without correction to highlight its effectiveness.

Fig. 1 .

 1 Fig. 1. functional chain of industrial robots

Fig. 2 .

 2 Fig. 2. model of robot joint wiht compliance

Fig. 3 .

 3 Fig.3. Set up for the golabal identification method[START_REF] Alici | Enhanced stiffness modeling, identification and characterization for robot manipulators[END_REF]

Fig. 4 .

 4 Fig. 4. Stiffness identification of joints 1 and 2

Fig. 5 .

 5 Fig. 5. Stiffness identification of joints 3

5 Fig. 6 .Fig. 7 .

 567 Fig. 6. Levels positionning for the identification of joints 3,4 and 5

Fig. 8 .

 8 Fig. 8. Robot configuration for the measurement of joint 2 error

Fig. 9 .

 9 Fig. 9. Position error of joint 2

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Residual torque applied on joint 2

Fig. 12 .

 12 Fig. 12. Robot configuration during correction tests

Fig. 13 .

 13 Fig. 13. measured trajectories without and with correction

TABLE I

 I

		STIFFNESS OF ROBOT JOINTS
	Joint N •	1	2	3	4	5	6
	K[N.m/rad].10 -6	0.204	0.85 0.57	0.49 0.12 0.005

b)

Joints 3,[START_REF] Gong | Nongeometric error identification and compensation for robotic system by inverse calibration[END_REF]5 and 6:

ACKNOWLEDGMENT

The authors would like to thank STAUBLI Robotics for providing the necessary information and Mr X.Helle for his help with the experimental validations.