Daniel Weisz-Patrault
email: weisz@lms.polytechnique.fr

Inverse Cauchy method with conformal mapping : application to latent flatness defect detection during rolling process

Keywords: Conformal Mapping, Inverse Problem, Residual stress, Cauchy problem

In this paper a semi-analytical inverse Cauchy problem is presented for a finite cylinder with a cylindrical hole under the surface called the detecting roll. Applications to latent flatness defect detection during rolling process are considered. A steel strip wraps at a known angle around the detecting roll. The presented inverse method aims at evaluating the residual stress profile of the strip. The inputs are the displacements and free surface traction at the hole surface and the outputs are the resultant forces per unit length in the contact between the strip and the roll, thus by equilibrium of the strip the residual stress profile is inferred. This paper is based on plane complex elasticity with conformal mapping techniques. Several 2D-computations are performed at several axial positions. This work is theoretical, the measurement device is not available and technical issues are not broached. Some numerical examples with synthetic data are provided in order to evaluate accuracy and noise sensitivity. Results encourage further investigations and technical developments.

Material coefficients for complex elasticity

Introduction

Inverse Cauchy problems refer to the class of inverse problems that take as inputs overdetermined conditions over a part of the boundary1 (typically displacement and surface traction for elastic problems) in order to infer the solution and in particular boundary conditions at the remaining part of the boundary. Mathematical studies (existence, uniqueness and stability of the solution) gathered for instance by [START_REF] Isakov | Inverse problems for partial differential equations[END_REF] give a framework for practical methods. For example, [START_REF] Delvare | An iterative boundary element method for cauchy inverse problems[END_REF] proposed an iterative boundary-element method for an inverse Cauchy problem related to the Laplace's equation. [START_REF] Marin | The method of fundamental solutions for the cauchy problem in two-dimensional linear elasticity[END_REF] proposed an inverse method based on fundamental solutions for the Cauchy problem in two dimensions related to isotropic elasticity, [START_REF] Andrieux | An energy error-based method for the resolution of the cauchy problem in 3d linear elasticity[END_REF] proposed an energy error-based method for 3D inverse Cauchy problems (both displacement and surface tractions are known for a part of the boundary and the remaining part should be determined), and [START_REF] Delvare | An iterative method for the cauchy problem in linear elasticity with fading regularization effect[END_REF] also proposed an iterative 3D inverse method for Cauchy problems. [START_REF] Liu | A modified collocation trefftz method for the inverse cauchy problem of laplace equation[END_REF] developed an inverse method for the Laplace equation using overdetermined data on a part of a circle and recovering data on the remaining part. Moreover, concept of topological sensitivity is used by [START_REF] Bellis | A fem-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF] for identification of cavities from overdetermined data. In addition, [START_REF] Liu | An analytical method for the inverse cauchy problem of laplace equation in a rectangular plate[END_REF] proposed an analytical method using Fredholm integral equation dedicated to the inverse Laplace equation in a plate.

This paper presents an inverse Cauchy method (related to Lamé's equations) that interprets displacement measurements on a free surface of a cylindrical hole drilled under the surface of a finite cylinder (defined as "detecting roll" in the following). This work is theoretical and no real measurements are provided since the measurement system is not available, thus technical issues are not broached. Synthetic data are extracted from Finite Element simulations instead, in order to test the accuracy of the method. The paper focuses on the inverse method itself that gives as outputs the resultant forces per unit length along the axial direction applied in the contact between the detecting roll and a steel strip that wraps at a known angle around it, as shown basically in Figure 1. Let Ω be the studied body (detecting roll). Ignoring flat edges of the latter, the boundary can be split into two parts ∂Ω = ∂Ω s ∪ ∂Ω i where ∂Ω s is the external surface and ∂Ω i the inner hole surface. The inverse problem solved in this paper is summarized as follows: finding displacement and stress fields (u, σ) such as:

on Ω,

           div σ = 0 σ = λtr(ϵ)1 + 2µϵ ϵ = 1 2 ( ∇u + ∇ T u ) on ∂Ω i , { u = u m σ.n = 0 (1)
where ϵ is the strain tensor, (λ, µ) the Lamé's coefficients, n the unit normal vector to the boundary and u m is the prescribed or measured displacements at the hole surface. Notations are listed in Table 1. As a matter of fact, the contact angle between the strip and the roll is quite small (around 15 • ) so the surface traction is rather sharp. Poor conditioning of this ill-posed problem with such a sharp loading makes impossible to evaluate accurately displacement and stress profiles at the surface of the roll, because the regularization level is very high (as detailed in Section 4). However, in this paper, the quantities of interest are not the complete displacement and stress fields (especially at the surface), but only the profile along the axial direction of the resultant forces per unit length of the contact pressure that takes place in the contact between the strip and the detecting roll. As detailed in following sections residual stress profile in the strip can be inferred from ( 2) and (3) on the basis of the presented inverse method. Short computation times are desired and a semi-analytical inverse method is proposed. In the field of semi-analytical inverse methods one can mention some works (elasticity and thermal problem) adapted for rolling process developed by [START_REF] Weisz-Patrault | A new sensor for the evaluation of contact stress by inverse analysis during steel strip rolling[END_REF], 2012a, 2013a) in 2D and Weisz-Patrault et al. (2013b, 2014b) in 3D. Moreover thermal experimental studies partially dedicated to the design of the measurement system have been proposed on this basis by Weisz-Patrault et al. (2012b); [START_REF] Legrand | Analysis of roll gap heat transfers in hot steel strip rolling through roll temperature sensors and heat transfer models[END_REF][START_REF] Legrand | Characterization of roll bite heat transfers in hot steel strip rolling and their influence on roll thermal fatigue degradation[END_REF]. However this paper deals with a more complicated shape (cylinder with a not centered hole). Nevertheless a semi-analytical solution is obtained by conformal mapping techniques. In the general field of inverse problems using conformal mapping techniques [START_REF] Kress | Inverse problems and conformal mapping[END_REF] proposed a unified review of some of his works related to the Laplace's equation where the geometry with an inclusion is unknown and determined by conformal mapping and Cauchy conditions (overdetermined conditions at the surface and simple condition in the inclusion). [START_REF] Hassen | Inverse impedance boundary problem via the conformal mapping method: the case of small impedances[END_REF] extended this work to the case of small impedances and [START_REF] Haddar | A conformal mapping method in inverse obstacle scattering[END_REF] extended recently the developed approach using conformal mapping techniques to inverse problems related to Helmholtz equation (inverse obstacle scattering for time-harmonic waves). This paper deals with Lamé's equation and the configuration is quite different, indeed the geometry is known and overdetermined conditions are prescribed at the inner hole surface, thus conformal mapping technique is used for simplifying the geometry and using some useful theorems about holomorphic functions as detailed in Section 3.

Detecing roll

It is demonstrated that several 2D inverse calculations using conformal mapping techniques applied at different axial positions are sufficient for the purpose of determining resultant forces per unit length distribution at the surface of the detecting roll. Thus, a Möbius transformation is used to map the real geometry E ζ into an annulus (transformed geometry E ξ ) as shown in Figure 6. Complex formalism proposed by [START_REF] Muskhelishvili | Some basic problems of the mathematical theory of elasticity[END_REF] is used and enables suitable introduction of conformal mapping. Thus an elegant semi-analytical solution to the presented inverse Cauchy problem is obtained. Weisz-Patrault et al. (2014a) published recently an extension in 3D of the classical elastic complex formalism using the 4D quaternionic algebra and the concept of monogenic functions (extension of holomorphy in the quaternionic algebra). Although Möbius transformations are still available in this extended formalism, the classical complex framework has been used for sake of simplicity and efficiency.

Applications

The inverse Cauchy problem introduced in Section 1 has a significant industrial application within the framework of rolling process. Indeed, flatness requirements for wider and thinner strips with harder grades tend to increase significantly in order to meet modern industrial needs. A general strategy of the steelmaking industry consists in selling less steel by reducing thickness and increasing material properties in order to remain competitive, for example for the automotive industry and packing. This becomes critical for rolling process of steel and aluminum. The difference between the incoming strip profile and the work roll deformed profile (due to bending and thermo-elastic deformations) results in heterogeneous plastic elongations in the strip that cause residual stresses. The complete mechanism is approached by numerical simulations that model the coupling between the thermo-elasto-visco-plastic behaviour of the strip and the thermo-elastic behavior of the roll and the post-bite buckling of the strip. For instance [START_REF] Abdelkhalek | Coupled approach for flatness prediction in cold rolling of thin strip[END_REF] proposed such a comprehensive model based on an older approach established by [START_REF] Hacquin | Modelisation thermo-mecanique tridimensionnelle du laminage: couplage bande-cylindres [3D thermomechanical modelling of rolling processes: coupling strip and rolls[END_REF].

The residual stress distribution on a strip section transversal to the rolling direction is called a flatness defect, which can be sufficiently compressive to cause local buckling. Manifest flatness defects consist in classical out of plane waves that can be observed in longitudinal, transverse or oblique directions considering which direction of the residual stress field is compressive. Thus, different types of flatness defects can be identified. However this paper focuses only on flatness defects where the residual stress profile consists in traction/compression according to the rolling direction e s . Longitudinal wrinkles (waves due to compressive residual stress according to the strip width direction e t ) are not broached in this contribution. Three major defects are taken into account: wavy edges or long edges (plastic elongations are more important at both edges), center buckle or long center (plastic elongations are more important at the center) and bad leveling (plastic elongations are more important at only one edge). Many optical methods quantify on-line manifest defects. [START_REF] Börchers | Topometric measurement of the flatness of rolled products -the system topplan[END_REF] and [START_REF] Paakkari | On-line measurement of large steel plates using Moiré topography[END_REF] used respectively fringe projection and moiré techniques. [START_REF] Garcia | Flatness measurement system based on a nonlinear optical triangulation technique[END_REF] and [START_REF] Molleda | Shape measurement of steel strips using a laser-based three-dimensional reconstruction technique[END_REF] used laser based techniques. [START_REF] Molleda | On-line flatness measurement in the steelmaking industry[END_REF] gives a significant review of flatness measurement systems and especially for optical methods.

However, during rolling process traction is applied to the strip. Thus, compressive parts of residual stresses can be compensated and defects remain not visible (no apparent buckling) on the production line. The flatness defect is called latent in this case, because the strip can therefore buckle after the traction is released (i.e. when the coil is unrolled). This paper is an attempt to detect in real time such latent flatness defects during the process. Optical methods are based on geometrical measurements and cannot measure directly residual stress profiles. Contact methods are suitable for such latent defects.

An intensive work has been done in this field and numerous patented industrial solutions are commonly used with a close-loop control that adapts the rolling parameters in real time when an estimation of flatness defects is measured. For all these solutions a detecting roll (also called tensiometer, shapemeter, stressometer etc...) is used to measure the residual stress profile when the strip wraps at a known angle around it, as shown in Figure 2. The initial principle proposed by [START_REF] Sivilotti | Strip flatness sensor[END_REF]; [START_REF] Miihlberg | Measurements up longitudinal stresses in moving metal bands and devices therefor[END_REF]; [START_REF] Berger | Apparatus for measuring stress distribution across the width of flexible strip[END_REF] was to divide the detecting roll into several slices or channels that measure independently the resultant force denoted by F j where j ∈ {1; ...; N c } is the channel index and N c is the number of channels. A transducer is fully embedded inside each slice. The delivered signal of each channel is proportional to the corresponding force, which is therefore directly inferred. Channel widths vary between 52 and 18 mm. Investigations and patented systems tend to limit the risk of scratching the strip such as the works of [START_REF] Kipping | Planarity measuring roller[END_REF]; [START_REF] Dahlberg | Seamless Stressometer Roll Eliminating the Strip marking and flatness measurement dilemma[END_REF]; [START_REF] Yang | Industrial shape detecting system of cold rolling strip[END_REF]. Basically these solutions consist in using a continuous detecting roll not divided into independent slices, several pressure sensors are embedded into the roll underneath the surface. Different kinds of transducers have been used: magnetostrictive sensor by [START_REF] Dahlberg | Seamless Stressometer Roll Eliminating the Strip marking and flatness measurement dilemma[END_REF]; [START_REF] Guodong | Shape meter of pressing magnetic type for cold rolling strip to applicating 4-roller reversible rolling mill[END_REF], piezzoelectric sensors by [START_REF] Berger | Deflection measuring roller[END_REF]; [START_REF] Bingqiang | A new type of contact shape meter and its industry application[END_REF] or Linear Variable Differential Transformer (LVDT) by [START_REF] Faure | Method of and a device for flatness detection[END_REF]. Signal interpretation being based on experimental calibration, sensors are sufficiently distant from one to another so that a punctual force at the surface of the detecting roll (with radial alignment with a sensor) does not activate the other sensors. A decorrelated signal is therefore obtained and the interpretation of experimental data is fairly easy, each sensor delivers a signal proportional to the resultant force considered on a fictive channel corresponding to the distance between two successive sensors. Thus the spatial resolution is also limited (around 20 mm). Since thin strips are considered, all authors assume for practical solution that shear stress in the strip is negligible compared with traction along the rolling direction so that one can consider strip channels independent and a trivial (string) equilibrium gives the resultant traction along the rolling direction T j (in the j-th strip channel) directly proportional to the measured resultant force F j as shown in Figure 3:

∀ j ∈ {1; ...; N c } , T j = F j 2 sin ( α 2 ) ( 2
)
where α is the wrapping angle known by geometrical considerations. Thus, there is equivalence between the resultant force profile at the surface of the detecting roll and the residual stress profile in the strip:

∀ j ∈ {1; ...; N c } , σ j ss = T j l j e (3)
where e is the strip thickness and l j the j-th strip channel width. Thus, contact methods that evaluate the resultant forces per unit length profile in the contact between the strip and the detecting roll enable the estimation of the residual stress profile in the strip according to the rolling direction e s .

F j α T j T j j-th strip channel

F j-1 T j-1 T j-1 F j+1 T j+1 T j+1
Figure 3: Simple equilibrium However, spatial resolution limitations lead to uncertainties near strip edges. Indeed, very local residual stresses are predicted by numerical simulations of the rolling process for instance in the model proposed by [START_REF] Abdelkhalek | Coupled approach for flatness prediction in cold rolling of thin strip[END_REF]) based on Hacquin (1996). Therefore there is a risk of undetected over-traction at strip edges leading to fracture or very local compression leading to buckling with short wave length (around 1 mm) as shown in principle in Figure 4. This paper aims at developing an alternative approach based on an inverse method. Current detecting rolls (tensiometers, stressometers etc...) measure directly at several axial positions the resultant forces (mean value over strip channel widths) by means of experimental calibration. The main idea of this paper is to develop an inverse method instead of direct calibration in order to overcome the spatial resolution limitation difficulty and to control the residual stress distribution at strip edges. To do so, the method interprets displacement measurements on the traction free surface of the axial hole made under the surface of the detecting roll, as shown in Figure 1. The idea is that a free hole surface emphasizes displacements and is accessible for measurements. Thus, displacements u m

x and u m y (m meaning measured), at the free hole surface are assumed to be accessible via a suitable measurement system. The technical solution is not developed here and the paper constitutes a motivation to adapt displacement measurement techniques, in order to measure displacements in a cylindrical hole.

As mentioned, this paper focuses on theoretical aspects only, technical issues are not broached and the measurement system is not available. Thus synthetic data extracted from Finite Element simulations are provided instead of real measurements. The aim is to encourage further technical developments of the measurement device. The choice of measuring displacements is motivated by the fact that very small devices exist, such as captive probes. Therefore the spatial resolution could be significantly improved. Moreover, capacitive probes can be used in harsh environment, which is suitable for real industrial conditions. Furthermore, vibrations that occur on the manufacturing line generate rigid body motions. Thus the measurement system should be fixed to the detecting roll so that independent vibrations are avoided and the same rigid body motions affect both the measurement system and the detecting roll in order to limit vibration effects.

The contact angle between the strip and the detecting roll can be controlled and kept constant with positioning rolls as presented in Figure 5. A backup roll also enables to avoid significant bending. The detecting roll has a rotation speed between 1 cycle per second and 10 cycles per second. In principle the inverse method could work at any time during the roll rotation, but in practice when the hole is not roughly aligned with the contact between the strip and the roll, the amplitude of displacements is very low and conditioning of the inverse Cauchy problem becomes critical. Displacement amplitude clearly reaches a maximum when the hole passes under the contact. Therefore this position is optimal for interpreting measurements with the inverse method. Thus, this inverse method is applied once for each cycle of the detecting roll. So the residual stress profile in the strip is not evaluated continuously but at several positions along the rolling direction e s according to the circumference of the detecting roll. Since flatness defects spread over a few meters along the rolling direction, a radius of 200 mm (circumference around 1.26 m) is a consistent design. Furthermore, since measurements are done only once during each cycle, an uneven rotation speed does not affect the inverse method (which would be the case if measurements were recorded all along the cycle).

Backup roll

Detecting roll

Positionning roll

Strip

Inverse method

The inverse method presented in this paper relies on two-dimensional elastic inverse calculations, applied several times at different axial positions. Each section is solved under plane strain assumption. This is justified by the fact that the detecting roll is much wider than the strip, thus even at strip edges the detecting roll is locally almost infinite, thus the plane strain assumption is well verified. The real geometry E ζ of each section is a disk of radius R s (s meaning surface) with a not centered hole of radius R i (i meaning internal). Let b denote the distance between the hole center and the roll center as shown in Figure 6 (left). For sake of simplicity, the inverse method is detailed without notation that specifies the axial position of the section, however the inverse method is applied several times at different axial positions and estimates the profile along the axial direction of resultant forces per unit length characterizing the contact between the strip and the roll or equivalently the profile of residual traction in the strip according to the rolling direction. Complex formalism due to [START_REF] Muskhelishvili | Some basic problems of the mathematical theory of elasticity[END_REF] is used, stress and displacement fields are related to two holomorphic potentials ϕ 0 and ψ 0

2 :            σ xx + σ yy = 2(ϕ ′ 0 (ζ) + ϕ ′ 0 (ζ)) σ yy -σ xx + 2iσ xy = 2(ζϕ ′′ 0 (ζ) + ψ ′ 0 (ζ)) 2µ(u x + iu y ) = κϕ 0 (ζ) -ζϕ ′ 0 (ζ) -ψ 0 (ζ) with κ =          3 -4ν plane strain 3 -ν 1 + ν plane stress (4)
Conditions of the inverse Cauchy problem are written at the hole surface ∂Ω i , which is thus described by a complex variable t i such as |t i -b| = R i . Let u m x (t i ) and u m y (t i ) denote the measured displacements around the hole surface according to e x and e y directions (m meaning measured). A complex quantity g 0 (t i ) = u m x (t i ) + iu m y (t i ) represents the displacement vector. Furthermore surface traction vanishes on ∂Ω i , therefore inverse Cauchy conditions are:

|t i -b| = R i , { ϕ 0 (t i ) + t i ϕ ′ 0 (t i ) + ψ 0 (t i ) = 0 κϕ 0 (t i ) -t i ϕ ′ 0 (t i ) -ψ 0 (t i ) = 2µg 0 (t i ) (5)
The general inverse Cauchy problem (1) reduces for each axial position (considered as a 2D problem) to ( 4) and ( 5). The first and second equations of ( 5) correspond respectively to free surface traction and measured displacements (inputs of the inverse method) at the hole surface. The difference between a classical well-posed direct problem and an inverse Cauchy problem is that for the direct problem displacement or surface traction (not both) are imposed on the whole boundary (roll surface and inner hole surface) and for the ill-posed inverse Cauchy problem both displacement and surface traction are imposed only at the inner hole surface. Thus, the inverse Cauchy problem reduces to seek two holomorphic potentials ϕ 0 and ψ 0 on E ζ that verify (5). Then displacement and stress fields are inferred from (4) especially at the surface of the roll, therefore the resultant forces per unit length is easily computed. The purpose of all following developments is to seek efficiently ϕ 0 and ψ 0 holomorphic on E ζ that verify (5). A well-known method useful for solving elastic problems with simple shapes is to expand both holomorphic potentials into power series. Indeed, all holomorphic functions admit an expansion into powers series on disks and annuli (with positive powers only on disks and with both positive and negative powers on annuli). A Möbius transformation ω such as ω -1 maps the real domain E ζ (disk with a not centered hole) into E ξ an annulus of unit external radius and internal radius r i < 1 (i meaning internal), is introduced:

r i E ζ ζ-plane E ξ ξ-plane ω θ i θ s τ i τ s b R s R i t i t s φ s φ i 1 e x e y
ω :          E ξ → E ζ ξ → ζ = ω(ξ) = R s ξ -a 1 -ξa ω -1 :          E ζ → E ξ ζ → ξ = ω -1 (ζ) = ζ + R s a R s + ζa (6)
2 Holomorphy of ϕ 0 and ψ 0 is guaranteed because there is no resultant force applied on the hole, but ϕ ′ 0 and ψ ′ 0 are always holomorphic as exposed by [START_REF] Lu | Complexe variable methods in plane elasticity[END_REF] Parameters a and r i are completely determined by the real geometry (i.e. R s , R i and b). Indeed it is well known that Möbius transformations map circles to circles, thus both conditions specifying that the circle of radius R s is transformed into the circle of radius 1 and the not centered circle of radius R i is transformed into the centered circle of radius r i completely determine r i and a.

Conformal mapping ω is holomorphic such as ω ′ (ξ) 0 (for all ξ ∈ E ξ ). Thus, ϕ and ψ defined by ( 7) are holomorphic on E ξ . Moreover the hole surface parameterized by t i such as |t i -b| = R i is transformed into the circle of radius r i parameterized by τ i such as |τ i | = r i . Thus, measured displacements g 0 (t i ) are transformed into g(τ i ) as defined by (8).

∀ξ ∈ E ξ , { ϕ(ξ) = ϕ 0 (ω(ξ)) ψ(ξ) = ψ 0 (ω(ξ)) ( 7
)
|τ i | = r i , g(τ i ) = g 0 (ω(τ i )) (8)
By introducing holomorphic functions ϕ and ψ into (5), inverse conditions expressed in the transformed geometry are obtained:

|τ i | = r i ,                ϕ(τ i ) + ω(τ i ) ω ′ (τ i ) ϕ ′ (τ i ) + ψ(τ i ) = 0 κϕ(τ i ) - ω(τ i ) ω ′ (τ i ) ϕ ′ (τ i ) -ψ(τ i ) = 2µg(τ i ) (9)
The inverse Cauchy problem is therefore reduced to seek two holomorphic potentials ϕ and ψ on the annulus E ξ such as ( 9) is verified. This simple transformed shape is suitable for finding a solution, because holomorphic functions can be expanded into a Laurent series, which is not the case in the disk with a not centered hole.

Hence, after basic manipulations on (9):

|τ i | = r i ,              ϕ(τ i ) = 2µ 1 + κ g(τ i ) ψ(τ i ) = -       ϕ(τ i ) + (1 -τ i a) 2 (r 2 i -τ i a) (τ i -r 2 i a)(1 -a 2 ) ϕ ′ (τ i )       (10)
The potential ϕ is expanded into a Laurent series (power series with both positive and negative powers), because it is holomorphic on the annulus E ξ :

ϕ(ξ) = ϕ 0 + +∞ ∑ k=1 ϕ k ξ k + +∞ ∑ k=1 ϕ * k ξ k (11)
An analogous form is sought for g(τ i ) so that coefficients ϕ k and ϕ * k can be identified by using (10). Measured displacements g(τ i ) can be expanded into a Fourier series. Let τ i = r i exp(iθ i ) where θ i denotes the angular position as shown in Figure 6. Measured displacements parameterized with θ i are abusively noted by g(θ i ), thus:

g(θ i ) = +∞ ∑ k=-∞ g k exp(ikθ i ) (12)
with the expression of the Fourier coefficients, which can be computed numerically with fast integration techniques using Fast Fourier Transform (fft) algorithms:

g k = 1 2π ∫ π -π g(θ i ) exp(-ikθ i )dθ i ( 13
)
Thus:

|τ i | = r i , g(τ i ) = g 0 + +∞ ∑ k=1 g k τ k i + +∞ ∑ k=1 g * k τ k i ( 14
)
Where:

               ∀k ∈ N g k = g k r k i ∀k ∈ N * g * k = g -k r k i ( 15
)
By re-writting the first equation of (10):

|τ i | = r i , ϕ 0 + +∞ ∑ k=1 ϕ k τ k i + ϕ * k τ k i = 2µ 1 + κ g 0 + 2µ 1 + κ +∞ ∑ k=1 g k τ k i + g * k τ k i (16)
Both equations are divided by 2iπ(τ i -ξ) and integrated over the circle of radius r i by considering ξ alternatively inside or outside the disk of radius r i :

1 2iπ

∫ |τ i |=r i τ k i τ i -ξ dτ i = { ξ k |ξ| < r i 0 |ξ| > r i and 1 2iπ ∫ |τ i |=r i 1 τ k i (τ i -ξ) dτ i =          0 |ξ| < r i - 1 ξ k |ξ| > r i (17)
Therefore ( 16) can be re-written as follows:

                   ϕ 0 + +∞ ∑ k=1 ϕ k ξ k = 2µ 1 + κ g 0 + 2µ 1 + κ +∞ ∑ k=1 g k ξ k |ξ| < r i +∞ ∑ k=1 ϕ * k ξ k = 2µ 1 + κ +∞ ∑ k=1 g * k ξ k |ξ| > r i ( 18
)
Hence by identification:

∀k ∈ N, ϕ k = 2µg k 1 + κ ∀k ∈ N * , ϕ * k = 2µg * k 1 + κ (19)
By re-writing the second equation of ( 10), it is sufficient to chose ψ(ξ) as the following holomorphic function on E ξ :

ψ(ξ) = -        ϕ 0 + +∞ ∑ k=1       ϕ k r 2k i ξ k + ϕ * k ξ k r 2k i       + (1 -ξa) 2 (r 2 i -ξa) (ξ -r 2 i a)(1 -a 2 ) +∞ ∑ k=1 [ kϕ k ξ k-1 -k ϕ * k ξ k+1 ]        (20)
A straightforward calculation gives:

ψ ′ (ξ) = +∞ ∑ k=1       kϕ k r 2k i ξ k+1 -kϕ * k ξ k-1 r 2k i       - (1 -ξa) 2 (r 2 i -ξa) (ξ -r 2 i a)(1 -a 2 ) +∞ ∑ k=1 [ k(k -1)ϕ k ξ k-2 + k(k + 1) ϕ * k ξ k+2 ] + (1 -aξ)(-2a 2 r 4 i -2a 2 ξ 2 + r 2 i (1 + a(ξ + a(-1 + 3aξ)))) (1 -a 2 )(ξ -ar 2 i ) 2 +∞ ∑ k=1 [ kϕ k ξ k-1 -k ϕ * k ξ k+1 ] (21)
Thus ϕ(ξ) and ψ(ξ) are known as a function of the Fourier coefficients of the displacements at the hole surface.

Stress and displacement fields in the whole detecting roll and especially at the surface are evaluated directly in the transformed geometry as follows:

                           σ xx + σ yy = 2       ϕ ′ (ξ) ω ′ (ξ) + ϕ ′ (ξ) ω ′ (ξ)       σ yy -σ xx + 2iσ xy = 2       ω(ξ) [ ω ′ (ξ) ] 2 ϕ ′′ (ξ) - ω(ξ)ω ′′ (ξ) [ ω ′ (ξ) ] 3 ϕ ′ (ξ) + ψ ′ (ξ) ω ′ (ξ)       2µ(u x + iu y ) = κϕ(ξ) - ω(ξ) ω ′ (ξ) ϕ ′ (ξ) -ψ(ξ) (22)
Only the resultant of surface traction is sought.

Points of the roll surface in the real geometry E ζ are denoted by t s such as |t s | = R s , hence t s = R s exp(iφ s ) where φ s is the angular position. Thus, the normal vector is n = cos(φ s )e x + sin(φ s )e y . Surface traction denoted by T = T x e x + T y e y is:

       T x = σ xx cos φ s + σ xy sin φ s = Re [ (σ xx -iσ xy ) exp(iφ s ) ] T y = σ xy cos φ s + σ yy sin φ s = Im [ (σ yy + iσ xy ) exp(iφ s ) ] (23)
In the transformed geometry, points of the roll surface are denoted by τ s such as |τ s | = 1, hence τ s = exp(iθ s ) where θ s is the angular position. Resultant forces are obtained by integrating surface traction. Since each section is self-equilibrated at the surface of the detecting roll (free surface traction being imposed at the hole surface), the integration over the whole circle [-π, π] gives zero. In order to obtain the resultant force corresponding to the contact between the strip and the roll, the integration is computed in the real geometry over [ -φ 0 , φ 0 ] where φ 0 is the half contact angle. With the parameters chosen in this paper and listed in Table 2 numerical values are φ 0 ≃ 6.78 • or 0.12 rad. This 13.56 • contact angle is ensured with the system presented in Figure 5. Since integration in the transformed geometry is wanted, one can introduce θ 0 corresponding to the half contact angle in the transformed geometry. Thus by using (6) and parameters listed in Table 2, θ 0 ≃ 1.92 rad or 110 • . Therefore, following expressions are obtained for the resultant forces:

F x = R s ∫ φ 0 -φ 0 T x dφ s F y = R s ∫ φ 0 -φ 0 T y dφ s (24)
Hence, by plugging ( 23) in ( 24):

F x = R s Re [∫ φ 0 -φ 0 (σ xx -iσ xy ) exp(iφ s )dφ s ] F y = R s Im [∫ φ 0 -φ 0 (σ yy + iσ xy ) exp(iφ s )dφ s ] (25)
Conformal mapping (6) gives:

exp(iφ s ) = exp(iθ s ) -a 1 -a exp(iθ s ) (26)
By differentiation of (26):

exp(iφ s )dφ s = (1 -a 2 ) exp(iθ s ) (1 -a exp(iθ s )) 2 dθ s (27)
Resultant forces are obtained as an integral over [-θ 0 , θ 0 ] in the transformed geometry:

                     F x = R s Re [∫ θ 0 -θ 0 (σ xx -iσ xy ) (1 -a 2 ) exp(iθ s ) (1 -a exp(iθ s )) 2 dθ s ] F y = R s Im [∫ θ 0 -θ 0 (σ yy + iσ xy ) (1 -a 2 ) exp(iθ s ) (1 -a exp(iθ s )) 2 dθ s ] (28)
The system ( 22) is evaluated at the roll surface (i.e. for ξ = τ s with |τ s | = 1):

|τ s | = 1,                  σ xx -iσ xy = ϕ ′ (τ s ) ω ′ (τ s ) + ϕ ′ (τ s ) ω ′ (τ s ) - ω(τ s ) [ω ′ (τ s )] 2 ϕ ′′ (τ s ) + ω(τ s )ω ′′ (τ s ) [ω ′ (τ s )] 3 ϕ ′ (τ s ) - ψ ′ (τ s ) ω ′ (τ s ) σ yy + iσ xy = ϕ ′ (τ s ) ω ′ (τ s ) + ϕ ′ (τ s ) ω ′ (τ s ) + ω(τ s ) [ω ′ (τ s )] 2 ϕ ′′ (τ s ) - ω(τ s )ω ′′ (τ s ) [ω ′ (τ s )] 3 ϕ ′ (τ s ) + ψ ′ (τ s ) ω ′ (τ s ) (29)
Finally by introducing

I x = Re [ R s (σ xx -iσ xy )(1 -a 2 )τ s /(1 -aτ s ) 2 ] and I y = Im [ R s (σ yy + iσ xy )(1 -a 2 )τ s /(1 -aτ s ) 2 ] : |τ s | = 1,                  I x = R s Re       (1 -a 2 )τ s (1 -τ s a) 2       ϕ ′ (τ s ) ω ′ (τ s ) + ϕ ′ (τ s ) ω ′ (τ s ) - ω(τ s ) [ω ′ (τ s )] 2 ϕ ′′ (τ s ) + ω(τ s )ω ′′ (τ s ) [ω ′ (τ s )] 3 ϕ ′ (τ s ) - ψ ′ (τ s ) ω ′ (τ s )             I y = R s Im       (1 -a 2 )τ s (1 -τ s a) 2       ϕ ′ (τ s ) ω ′ (τ s ) + ϕ ′ (τ s ) ω ′ (τ s ) + ω(τ s ) [ω ′ (τ s )] 2 ϕ ′′ (τ s ) - ω(τ s )ω ′′ (τ s ) [ω ′ (τ s )] 3 ϕ ′ (τ s ) + ψ ′ (τ s ) ω ′ (τ s )             (30)
By using ( 11) and ( 6) surface traction (30) becomes:

|τ s | = 1,                                                                                                I x = Re        +∞ ∑ k=1 ( kϕ k τ k-1 s -k ϕ * k τ k+1 s ) ( τ s + 2aτ s (1 -τ s a) 2 (1 -a 2 )(τ s -a) ) + ( kϕ k τ k+1 s -kϕ * k τ k-1 s ) τ s (τ s -a) 2 (1 -τ s a) 2 - ( k(k -1)ϕ k τ k-2 s + k(k + 1) ϕ * k τ k+2 s ) τ s (1 -τ s a) 3 (1 -a 2 )(τ s -a) -       kϕ k r 2k i τ k+1 s -kϕ * k τ k-1 s r 2k i       τ s + τ s (1 -τ s a) 2 (r 2 i -τ s a) (τ s -r 2 i a)(1 -a 2 ) ( k(k -1)ϕ k τ k-2 s + k(k + 1) ϕ * k τ k+2 s ) - τ s (1 -τ s a)(-2a 2 r 4 i -2a 2 τ 2 s + r 2 i (1 + a(τ s + a(-1 + 3aτ s )))) (1 -a 2 )(τ s -ar 2 i ) 2 ( kϕ k τ k-1 s -k ϕ * k τ k+1 s )      I y = Im        +∞ ∑ k=1 ( kϕ k τ k-1 s -k ϕ * k τ k+1 s ) ( τ s - 2aτ s (1 -τ s a) 2 (1 -a 2 )(τ s -a) ) + ( kϕ k τ k+1 s -kϕ * k τ k-1 s ) τ s (τ s -a) 2 (1 -τ s a) 2 + ( k(k -1)ϕ k τ k-2 s + k(k + 1) ϕ * k τ k+2 s ) τ s (1 -τ s a) 3 (1 -a 2 )(τ s -a) +       kϕ k r 2k i τ k+1 s -kϕ * k τ k-1 s r 2k i       τ s - τ s (1 -τ s a) 2 (r 2 i -τ s a) (τ s -r 2 i a)(1 -a 2 ) ( k(k -1)ϕ k τ k-2 s + k(k + 1) ϕ * k τ k+2 s ) + τ s (1 -τ s a)(-2a 2 r 4 i -2a 2 τ 2 s + r 2 i (1 + a(τ s + a(-1 + 3aτ s )))) (1 -a 2 )(τ s -ar 2 i ) 2 ( kϕ k τ k-1 s -k ϕ * k τ k+1 s )      (31)
Finally the resultant forces per unit length profile (along the axial direction z) is obtained by analytical integration:

F x = ∫ θ 0 -θ 0 I x (θ s )dθ s F y = ∫ θ 0 -θ 0 I y (θ s )dθ s (32)
4. Validation

Long edge defect

All inverse methods are ill-posed as defined among others by [START_REF] Isakov | Inverse problems for partial differential equations[END_REF]. Most of the time the stability condition is not fulfilled, that is to say that arbitrary small errors on the inputs of the inverse method lead to arbitrary large errors on the outputs. In Section 3, since analytical solution is developed, ill-posedness can be identified clearly. Indeed, 0 < r i < 1 so 1/r k i goes to infinity when k goes to infinity. Fourier coefficients (13) are computed from measured displacements with a given accuracy, so errors are amplified by 1/r k i in ( 15). An other source of instability is factors such as 1/(1a 2 ) in (31), because in practice |a| is very close to 1 as it can be seen in Table 2 and so 1/(1a 2 ) becomes very large. Thus, regularization that stabilizes the solution is needed. The presented inverse method uses expansions into series, therefore the regularization simply consists in truncating the sums. This paper aims at developing an inverse method for potential industrial purposes, and it has been discussed in Section 2 that only the resultant forces per unit length profile along the axial direction is needed and not the detailed pressure repartition in the contact. This weak information can be obtained with very few terms in the sum as shown in Figure 13 (convergence is reached with 3 terms). This is very convenient mainly because stability is guaranteed (keeping only very few terms in the sum is a very strong regularization) and the truncation number has not to be studied during the inversion procedure, no additional algorithm has to be developed. The criterion of truncation (3 terms) is the same for noisy inputs because in this case the regularization level is very high even without noise.

This section aims at evaluating numerically 3 the accuracy of the inverse Cauchy solution proposed in Section 3. Thus, the direct problem that provides the inputs for the inverse method is solved with FEM using Cast3m a free source software developed by CEA (2011). However, the purpose is not to build a comprehensive direct numerical simulation of the contact problem between the strip with a prescribed residual stress profile and the detecting roll. In this paper, the strip is not modeled in order to avoid long iterative non-linear contact problem. Prescribed resultant forces per unit length profiles (denoted by F p (z)) that represent the contact are imposed at the detecting roll surface so that linear FE simulations are obtained. Cubic elements are used, the mesh is presented in Figure 7 and model parameters are listed in Table 2. Boundary conditions for the direct FE simulation are: the centers of flat edges of the detecting roll are blocked and the prescribed contact pressure is applied at the roll surface with radial alignment with the hole center and the same pressure is applied at the diametrically opposite position in order to model the backup roll, as shown in Figure 5. This simple direct problem is sufficient to evaluate the quality of the presented inverse method. The influence of not modeling the strip is evaluated by considering two different contact pressures having the same resultant forces per unit length denoted by F p (z). It is shown that the method does not depend on the detailed contact repartition but only on the resultant forces per unit length F p (z). Since the inverse method enables to evaluate the resultant forces per unit length regardless of real local contact repartitions, the residual stress profile in the strip is available with (2) and (3) without modeling the contact and the strip. The relationship between the resultant forces per unit length and the residual stress profile in the strip is obtained by simple string equilibrium of the strip, which holds only if shear stress is negligible in the strip. This latter condition is satisfied because very thin strips are considered.

The two different contact pressures having the same resultant forces per unit length are denoted in the following by "quadratic conditions" and "constant conditions" because the local pressure repartitions are respectively quadratic or constant with respect to the circumferential direction as shown in Figure 8. Very similar results (evaluation of the resultant forces per unit length profile along the axial direction) are obtained for both conditions. Equations ( 2) and (3) give a relationship between the resultant forces per unit length and the residual stress profile in the strip. Thus, the prescribed resultant forces per unit length profile F p (z) (shown in Figure 9a) is set so that it corresponds to a 1 mm thick strip presenting a symmetric long edges defect as shown in Figure 9b. From the direct problem, synthetic data u m x and u m y are extracted at the hole surface and presented at several axial positions in Figure 10. In order to be consistent with a possible measurement system only a few points (around 10) are considered for the inputs for the inverse method as show in Figure 11. Then the inverse method presented in this paper is applied at each axial position of the mesh. Since symmetric pressures are applied with respect to the circumferential direction F y = 0, thus F x is compared with F p and the accuracy of the method is determined by the mean value of relative errors along the axial direction:

Err = 100 ∫ L -L [F(z) -F p (z)] 2 dz ∫ L -L [ F P (z) ] 2 dz (33)
Measurements are carried out practically with noise. In order to evaluate noise sensitivity, random numbers are added to the inputs. Capacitive probes have a very high resolution and can measure variations of displacements under 1 nm. Thus, three levels of noise are considered, with respective amplitude 0.01µm (which corresponds to 10 times the resolution of capacitive probes), 0.1 µm (100 times the resolution) and 0.3 µm (300 times the resolution). Examples of noisy inputs are given in Figure 12. Results without and with noise are presented in Figure 14 and quantified errors are listed in Table 3. Results are satisfying for the two first noise conditions, but the method fails for the third one (amplitude 0.3 µm). The measurement system should be designed such as the level of noise does not exceed 0.1 µm. The method converges with only 3 terms in the sums as demonstrated in Figure 13. This truncation number is the same for noisy inputs since it corresponds to a strong regularization.

Model parameters, listed in Table 2, have been chosen in order to be consistent with rolling process dimensions such as strip width, thickness, roll radius and length. The hole dimension has been set so that it seems technically possible to drill the hole under the roll surface. The inversion quality depends on the detecting roll design parameters. A convenient design is presented in this Section. Different designs are also studied in Section 5, where parameter b is varying from 148 mm to 138 mm. Sensitivity with respect to geometrical parameters are therefore highlighted.

0 -2 2 -3 -1 1 3 -2.5 -1.5 -0.5 0.5 1.5 2.5 0 -10 -8 -6 -4 -2 -11 -9 -7 -5 -3 -1 φ (rad) i z=-490 mm z=-455 mm z=-420 mm z=-385 mm u (μm) x m (a) Quadratic conditions u m x 0 -2 2 -3 -1 1 3 -2.5 -1.5 -0.5 0.5 1.5 2.5 0 -8 -6 -4 -2 -9 -7 -5 -3 -1 φ (rad)

Long center defect

A second numerical example is given where the prescribed resultant forces per unit length F p (z) correspond to a long center defect in the strip as presented in Figure 15. As mentioned the relationship between the resultant forces per unit length and the residual stress profile in the strip is given by ( 2) and (3). Both quadratic and constant local contact pressures are tested. Displacements extracted from the FE simulation at the hole surface are presented at several axial positions in Figure 16. Then, 10 points are taken into account in these displacements as presented in Figure 11. Measurement noise sensitivity is investigated alike in the previous section for long edges defects, thus three level of noise are considered (0.01 µm, 0.1 µm and 0.3 µm). Noisy inputs are similar to those presented in Figure 12. Finally, reconstructions of resultant forces per unit length along the axial direction are presented in Figure 17 for both quadratic and constant conditions, without and with noise on the inputs. Error estimates are given in Table 4. Similar results to the long edges defect conditions are obtained. Noise amplitudes under 0.1 µm do not compromise the method.

Parametric study

In this section a parametric study with respect to the position of the hole center b is presented in order to highlight geometrical constrains. The parameter b is varying from 148 mm (alike in Section 4) to 138 mm. Thus, several FE computations have been performed with geometrical and mechanical parameters listed in Table 5a and five different values of b listed in Table 5b. It is clear that the method quality is a compromise between the amplitude of displacements at the hole surface and the parameters r i and a as mentioned at the beginning of Section 4. When b decreases, r i and a also decrease. There is a regularization effect due to the fact that a becomes more and more different from 1, which reduces the instability due to factors 1/(1a 2 ). This latter regularization competes with the fact that the amplitude of displacements at the hole surface decreases which tends to degrade the integrals quality of (13) and increase errors amplitudes. Displacements at the hole surface are presented in Figure 18 at the axial position z = -385 mm for the five values of b. In the meanwhile r i decreases, thus errors are amplified by 1/r k i which becomes greater and greater when b decreases, however since the sums are truncated at only three terms this latter instability effect is not prevailing. Results for the same long edges defect as presented in Section 4 for quadratic conditions are presented in Figure 19 and errors are listed in Table 6, for both conditions without noise and with perturbed inputs with 0.1 µm amplitude. Accuracy decreases and noise sensitivity increases rapidly when b decreases. It is reasonable for this set of design parameters to fix b = 148 mm in order to ensure a good reconstruction quality.

Conclusion

A semi-analytical inverse method based on complex plane elasticity and conformal mapping techniques has been developed and tested by numerical examples. Application of this inverse method to latent flatness defect detection during rolling process has been considered (long edges and long center defects have been tested). Consistent quantity is evaluated (resultant forces per unit length profile along the axial position) in order to infer residual stress profile in steel strips. Good accuracy is observed and the added noise under 0.1 µm amplitude (which corresponds to 100 times the resolution of capacitive probes) does not compromise the method. This work encourages further technical developments for a feasible measurement system dedicated to the aimed application.

In this paper a Möbius transformation has been used, but more generally the same principles can be used for other geometries and other conformal mappings. The key point is the Riemann mapping theorem that ensures that one can transform any plane geometry simply or doubly-connected into a disk or an annulus respectively. Thus (9) arises for the inverse Cauchy problem where displacements and surface tractions are known on the internal or external surface. Therefore by adding both equations of (9) one holomorphic potential is directly identified, then the other one is obtained. Thus semi-analytical inverse solution can be established, alike the presented developments.

Acknowledgement

The author takes the opportunity to acknowledge Alain Ehrlacher (Laboratoire Navier), Nicolas Legrand (ArcelorMittal research and development) and Caio Pereira Santos (Ecole des Ponts ParisTech) for their ideas and fruitful discussions.

Figure 1 :

 1 Figure 1: Proposed detecting roll

  Figure 2: Detecting roll

Figure 4 :

 4 Figure 4: Spatial resolution problem

Figure 5 :

 5 Figure 5: Possible system

Figure 6 :

 6 Figure 6: Real and transformed geometries

Figure 7 :

 7 Figure 7: Mesh of the direct problem

Figure 9 :

 9 Figure 8: Prescribed contact pressure

Figure 10 :Figure 13 :

 1013 Figure 10: Displacements from FEM extracted at the hole surface for long edges defect

Figure 14 :

 14 Figure 14: Reconstruction of resultant forces per unit length for long edges defect

  long center defect in the strip

FigureFigure 16 :Figure 17 :

 1617 Figure 15: Prescribed long center defect

Figure 18 :

 18 Figure 18: Displacement at the hole surface for various b

Figure 19 :

 19 Figure 19: Reconstruction for long edges defect and different values of b

Table 1 :

 1 Nomenclature

		Strip
	e s	Rolling direction
	e t	Transverse direction
	e	Strip thickness
	N c	Number of strip channels
	l	Strip half-width
	l j	Width of channel j ∈ {1; ...; N c }
	T j	Tension along the rolling direction on the j-th channel
	σ j ss	Residual stress on the j-th channel
		Detecting roll
	x, y, z	Cartesian coordinates in the roll
	e x , e y	Directions in the plane of a roll section
	e z	Axial direction
	Ω	3D Elastic domain (detecting roll)
	∂Ω	Boundary ignoring flat edges
	∂Ω s	External curved surface
	∂Ω i	Inner hole surface
	σ	Stress tensor
	ϵ	Strain tensor
	u	Displacement vector
	u m x , u m y	Measured displacements
	F j	Resultant force on the j-th channel
	F x , F y Evaluated profile of resultant forces per unit length

Table 2

 2

		: Parameters
	R s	(mm)	200
	R i	(mm)	50
	b	(mm)	148
	L	(mm)	700
	a	(-)	-0.920357
	r i	(-)	0.783857
	E	(MPa)	210 000
	ν	(-)	0.3
	α = 2φ 0	( • )	13.56
	l	(mm)	500
	e	(mm)	1

Table 3 :

 3 Error estimate for long edges defect

	Noise amplitude Quadratic conditions Constant conditions
	Without noise	4.12 %	4.5%
	0.01 µm	4.3 %	4.6 %
	0.1 µm	5.7 %	8.35 %
	0.3 µm	29.7 %	34.6 %

Table 4 :

 4 Error estimate for long center defect

	Noise amplitude Quadratic conditions Constant conditions
	Without noise	0.5 %	0.5%
	0.01 µm	0.52 %	0.52 %
	0.1 µm	1.1 %	0.9 %
	0.3 µm	6.7 %	7.15 %

Table 5 :

 5 Parametric study

	(a) Fixed parameters		(b) Studied parameters
	R s R i L E ν α = 2φ 0 l	(mm) (mm) (mm) (MPa) 210 000 200 50 700 (-) 0.3 ( • ) 13.56 (mm) 500	b (mm) 148 147 146 145 138	a -0.920357 0.783857 r i -0.902601 0.742748 -0.887552 0.710052 -0.874199 0.682676 -0.820466 0.587304
		Table 6: Error estimate for parametric study
		b (mm) Without noise Noise amplitude 0.1 µm
		148	4.12 %	5.7 %
		147	4.2 %	8.6 %
		146	4.6 %	9.2 %
		145	5.4 %	14.5 %
		138	9.9 %	53 %

Preprint submitted to International Journal of Solids and StructuresDecember 8, 2014

If the problem is an Ordinary Differential Equation or a Partial Differential Equation of order k, derivatives of order 0, 1, • • • , k -1 are involved in the Cauchy conditions

No experiment has been done, the measurement system being not developed at the current state. Thus, measured quantities are replaced by synthetic data obtained by numerical simulation.


Leonard Pashkevich

Fr

Anatol Fatiha

Rançois Rédi

A Klimch

A Pashkev

S Garni

S Caro

F Leona

G Abb

J Qin

Lcfc Descrip

Modèle Simp Modèle

Modèle Modèle

G A Corousso

In r Corous Nom BOUJDAINE

ent de simula rocédé d'usin rocédé FSW

INTRODUCTION

L'objectif de cette tâche est de détailler les modèles utilisés en simulation et les résultats obtenus en simulation aussi bien pour le procédé d'usinage que pour le procédé FSW. De cette manière, on peut se rendre compte du comportement du robot durant le procédé et ainsi prévoir les meilleurs réglages mais également de limiter les coûts des essais qui seront réalisés dans la tâche 5.

Ce livrable contient donc dans une première partie les détails sur la modélisation des flexibilités des robots essentiellement par le biais d'un modèle de flexibilités localisées. Les flexibilités sont exprimées aussi bien dans l'espace cartésien des robots que dans l'espace articulaire et en tenant compte des éventuels couplages.

Une seconde partie traite du modèle dynamique utilisé dans le simulateur et de l'environnement de simulation. Un travail non négligeable a été de modéliser également le contrôleur de robot de l'entreprise Kuka.

Les procédés d'usinage et de soudage FSW sont eux-mêmes modélisés par des modèles simples mais traduisant bien la réalité du comportement. Ces modèles sont issus des livrables 2.1 à 2.3 et du livrable 3.4.

Ce rapport se termine par des résultats obtenus en simulation en usinage d'une part avec un robot Kuka KR 270 à l'Iut de Nantes et d'autre part en soudage FSW avec un robot Kuka KR 500 MT2 à l'institut de Soudure de Goin. The dynamic behaviour of the robot under the loading F caused by technological process can be described as

C C C + + = M δt C δt K δt F && & (1)
where

C
M is 6 6 × mass matrix that represents the global behaviour of the robot in terms of natural frequencies, C C is 6 6 × damping matrix, C K is 6 6 × Cartesian stiffness matrix of the robot under the external loading F , , δt δt & and δt && are dynamic displacement, velocity and acceleration of the tool end-point in a current moment respectively [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF].

In general, the cutting force F c has a nonlinear nature and depends on many factors such as cutting conditions, properties of workpiece material and tool cutting part, etc. [START_REF] Ritou | A new versatile in-process monitoring system for milling[END_REF]. But, for given tool/workpiece combination, the force F c could be approximated as a function of an uncut chip thickness h, which represents the desired thickness to cut at each instant of machining.

Hence, to reduce errors caused by cutting forces in the robotic-based machining it is required to obtain an accurate elasto-static model of robot and elasto-dynamic model of machining process. These problems are addressed in the following sections taking into account some particularities of the considered application (robotic-based milling).

Modèle avec raideurs localisées

Modèle de flexibilités dans l'espace cartésien

Elasto-static model. Elasto-static model of a serial robot is usually defined by its Cartesian stiffness matrix, which should be computed in the neighborhood of loaded configuration. Let us propose numerical technique for computing static equilibrium configuration for a general type of serial manipulator. Such manipulator may be approximated as a set of rigid links and virtual joints, which take into account elasto-static properties (Figure 1). Since the link weight of serial robots is not negligible, it is reasonable to decompose it into two parts (based on the link mass centre) and apply them to the both ends of the link. All this loadings will be aggregated

... n ⎡ ⎤ = ⎣ ⎦ J J J , T T T 1 ... n ⎡ ⎤ = ⎣ ⎦ G G G
. Finally, taking into account the virtual spring

reaction θ θ = ⋅ τ K θ , where (
)

1 n θ θ θ
,..., diag = K K K , the desired static equilibrium equations can be presented as

(G)T (F)T θ θ θ ⋅ + ⋅ = ⋅ J G J F K θ (6)
To obtain a relation between the external loading F and internal coordinates of the kinematic chain θ corresponding to the static equilibrium, equations (6) should be solved either for different given values of F or for different given values of t . Let us solve the static equilibrium equations with respect to the manipulator configuration θ and the external loading F for given end-effector position ( ) = t g θ and the function of auxiliary-loadings ( )

G θ ( ) ( ) (G)T (F)T θ θ θ ; ; ⋅ = + = = K θ J G J F t g θ G G θ (7)
where the unknown variables are ( )

, θ F . θ θ - - - - + + + + = ⋅ ⋅ - + - = ⋅ + ⋅ F J K J t g θ J θ J K J G θ K J G J F (10)
The proposed algorithm allows us to compute the static equilibrium configuration for the serial robot under external loadings applied to any point of the manipulator and the loading from the technological process.

Stiffness matrix.

In order to obtain the Cartesian stiffness matrix, let us linearize the force-deflection relation in the neighborhood of the equilibrium. Following this approach, two equilibriums that correspond to the manipulator state variables ( , , ) F θ t and ( δ , δ , δ ) + + + F F θ θ t t should be considered simultaneously. Here, notations δF , δt define small increments of the external loading and relevant displacement of the end-point. Finally, the static equilibrium equations may be written as

δ δ δ δ δ δ δ + = + ⋅ + = + ⋅ + + + ⋅ + t t g θ θ K θ θ J J G G J J F F (12)
where θ , , , , t F G K θ are assumed to be known.

After linearization of the function ( )

g θ in the neighbourhood of the loaded equilibrium, the system (11), ( 12) is reduced to equations

(F) θ (G) (G) (F) (F) θ θ θ θ θ δ δ δ δ δ δ δ = ⋅ = + + + t J θ K θ J G J G J F J F (13)
which defines the desired linear relations between δt and δF . In this system, small variations of Jacobians may be expressed via the second order derivatives

(F) (F) θ θ θ δ δ = ⋅ J H θ , (G) (G) θ θ θ δ δ = ⋅ J H θ ,
Corousso (G) 2 θθ 1 2 T j j j n = = ∂ ∂ ∑ H g G θ , (F θθ 2 ) 2 T = ∂ ∂ H g F θ .
Also, the auxiliary loading G may be computed via the first order derivatives as

δ δ = ∂ ∂ ⋅ G G θ θ
Further, let us introduce additional notation

(F) (G) (G)T θθ θθ θθ θ = + + ⋅∂ ∂ H H H J G θ
, which allows us to present system (13) in the form

(F) θ (F)T θ θ θ θ δ δ δ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ = ⋅ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ -+ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 0 J t F 0 θ J K H (14)
So, the desired Cartesian stiffness matrices C K can be computed as ( )

1 (F) 1 (F)T C θ θ θ θ θ ( ) - - = - K J K H J (15)
Below, this expression will be used for computing of the elasto-static deflections of the robotic manipulator.

Mass matrix. To evaluate dynamic behaviour of the robot under the loading, in addition to the Cartesian stiffness matrix C K it is required to define the mass matrix C M . Comprehensive analysis and definition of this matrix have been proposed in [START_REF] Briot | Reduced elastodynamic modelling of parallel robots for the computation of their natural frequencies[END_REF]. Here, let us summarise the main results that will be used further in the error compensation technique. M is defined in the joint coordinates it can be recomputed into the Cartesian coordinates associated with the tool end-point using the Jacobian matrix θ J (which depend on the robot configuration q and computed with respect to virtual joint coordinates θ ) using following expression

θ C θ θ T = M J M J (17)
Thus, using expressions (16) and (17) it is possible to compute the mass matrix La première étude des problèmes liés à la présence de transmissions flexibles dans les robots manipulateurs remonte à (Sweet and Good, 1985), (Good et al., 1985), avec les premiers résultats expérimentaux sur le bras GE P-50. Des considérations mécaniques impliquées dans la conception des bras robotiques et dans l'évaluation de leurs éléments peuvent être trouvées dans (Rivin, 1988). Une analyse bibliographique sur la dynamique des manipulateurs flexibles peut être trouvée dans (Dwivedy and Eberhard, 2006). Du point de vue de la modélisation des manipulateurs flexibles, on peut supposer que la flexibilité est localisée au niveau des articulations et/ou distribuée le long des corps. En effet, ces deux types de flexibilité peuvent être présents en même temps. Les flexibilités aux articulations sont souvent présentées sur des robots industriels lorsque certains éléments de transmission de mouvement/réduction sont utilisés : les courroies, les longues tiges, les câbles, les réducteurs harmonic drive ou les trains d'engrenages épicycloïdaux. Ces composants permettent de placer les actionneurs au plus près de la base du robot, d'améliorer l'efficacité dynamique, ou de garantir des taux de réduction élevés. Par contre, sous des contraintes de forces importantes, ces composants deviennent intrinsèquement flexibles (Siciliano and Khatib, 2008).

Le modèle à articulations flexibles suppose que les élasticités sont concentrées dans la chaine de transmission mécanique entre les moteurs et les corps rigides du robot (Makarov, 2013). Les robots étudiés ici sont des robots séries à chaines ouvertes. Ces robots peuvent être considérés avec des articulations flexibles et une chaine cinématique ouverte à n + 1 corps rigides. La base et les n corps sont reliés par n articulations et actionnés par n moteurs électriques. De plus, toutes les articulations sont considérées comme flexibles (Siciliano and Khatib, 2008). De Luca a présenté en détail la modélisation de robots avec des éléments flexibles dans (De Luca and Book, 2008). Une étude sur les sources de flexibilités des robots séries a été faite par Dumas (Dumas, 2011). Les résultats montrent que les flexibilités localisées aux articulations sont les sources principales des flexibilités des robots séries. Quand il y a des réducteurs, on peut placer la déformation avant ou après le réducteur. L'élasticité au niveau de la transmission peut être modélisée par un ressort avec une raideur de torsion sur les axes de rotation et avec une raideur de traction/compression pour les axes en translation.

Les flexibilités localisées aux articulations peuvent s'expliquer aisément si on suppose qu'elles sont dues à la torsion de l'arbre de transmission des réducteurs. Si on adopte alors le modèle de la Figure 1, on obtient l'équation suivante :

) ( ) ( q K q N K l v - = - = Γ θ θ (18)

Coro

Avec Γ l'inverse de vecteur des Le mod Où :

• J(q) (6 • T e = [ dans • Γ m = • Γ = [Γ • q = [ articu • θ = [θ des m • D(q) (6 et bo • I a = d • I = [I 1 usso Γ le vecteu e la matrice s angles artic dèle dynamiq 6×6) est la ma [Fx Fy Fz Cx le repère d'o [Γ m1 Γ m2 Γ m3 Γ 1 Γ 2 Γ 3 Γ 4 Γ 5 [q 1 q 2 q 3 q 4 ulaires, θ 1 θ 2 θ 3 θ 4 θ 5 moteurs, 6×6) est la ma ornée, diag([I a1 I a2 I 1 I 2 I 3 I 4 I 5 I 6 ]
que complet q D ) ( trice jacobie x Cy Cz] T est outil (R t ) , 3 Γ m4 Γ m5 Γ m6 Γ 6 ] T est le v q 5 q 6 ] T , q & θ 6 ] T , θ & et θ & atrice d'inert I a3 I a4 I a5 I a6 ]
] T est la matr Pro L ésultats en s n robot avec es à la sortie ts de réduct Flexibilité loc • Fm et Fs : les matrices de frottement de Coulomb aux niveaux des moteurs des articulations en Nm,

t du robot ut q H q ( - Γ = & & θ & & m a I Γ = ) ( B q F fr & = ) ( B q F fm & =
fr q F q q ( ) , - & T v F N - Γ - (sig F q B s s & + (si F q B m m & + I k t m = Γ ère outil (R t )
T T q J q ) ( ) - & ) (θ & fm F )) (q gn & )) (q ign &
• N la matrice (6x6) des rapports de réduction et N v son inverse.

Tous les paramètres de ce modèle ont été identifiés (voir livrable 4.2, projet ANR Corousso) pour les robots Kuka KR270-2F et KR500-2MT.

Les valeurs des flexibilités (en Nm/rd) et des rapports de réduction sont données dans les tableaux cidessous.

K ij K 1 K 2 K 3 K 4 K 5 K 6 KR270-2F
3.8 ×10 6 6.6 ×10 6 3.9 ×10 6 5.6 ×10 5 6.6 ×10 5 4.7 ×10 Ici, est le vecteur des angles après réduction et avant la flexibilité des réducteurs. de dimension 6 6 est la matrice des rapports de réduction. Elle est non diagonale à cause du couplage entre axe 4, 5 et 6.

Environnement de simulation

Le simulateur a été développé dans un environnement de travail Matlab/Simulink©. Il est composé d'un fichier Simulink décrivant l'ensemble des équations de comportement du robot, du procédé et de la commande (y compris les capteurs et les boucles de retour en effort). La Figure 2

Corousso

Les détails des explications du simulateur sont donnés dans la thèse Qin, 2013. Un résumé succinct extrait de ce mémoire est donné ci-dessous. Les blocs en bleu représentent la partie de la commande du robot, en rouge la régulation de force, en gris le générateur de trajectoire, en vert le système simulant le comportement dynamique des moteurs, et en orange la partie correspondant à la dynamique des axes du robot, ensuite en magenta la partie liée aux flexibilités, et en violet le modèle du procédé. Dans ce simulateur, les modèles du robot et du procédé peuvent être remplacés par ceux d'un autre robot ou procédé.

Les fonctions des blocs dans le simulateur sont :

• Générateur de trajectoire : transfert les trajectoires désirées ( )

d d d d a v p T =
(positions, vitesses et accélérations désirées, respectivement) du repère cartésien aux articulaires ( )

d d d q q q & & & ;
• P contrôleur : régulateur proportionnel d'asservissement de la position ;

• PI contrôleur : régulateur proportionnel-intégral d'asservissement de la vitesse ;

• DNL : Dynamique Non-Linéaire ;

• Modèle Dynamique du Moteur : modèle dynamique des moteurs ;

• Réducteurs : avec les rapports de transmission et intégral ;

• Flexibilité : modèle de la flexibilité aux articulations ;

• Modèle Dynamique du Robot : modèle dynamique du robot ;

• Force Procédé : modèle d'effort du procédé ;

• Régulateur de force : régulateur proportionnel-intégral des efforts.

Modèle du procédé d'usinage

Let us obtain the model of the cutting force which depends on the relative position of the tool with respect to the workpiece at each instant of machining. As follows from previous works (Brissaud, 2008), for the known chip thickness h, the cutting force F c can be expresses as:

( ) ( ) 2 0 , 0 1 s s c s p h h r h h F h k a h h h + = ≥ + (24)
where p a is a depth of cut,

0 1 r k k ∞ =
< depends on the parameters k ∞ , k 0 that define the so called stiffness of the cutting process for large and small chip thickness h respectively (Figure 3) and h s is a specific chip thickness, which depends on the current state of the tool cutting edge. The parameters k 0 , h s , r are evaluated experimentally for a given combination of tool/working material. To take into account the possible loss of contact between the tool and the workpiece, expression (24) should be supplement by the case of 0 h < as

( ) 0, if 0 c F h h = < ( 25
)
For the multi-edge tool the machining surface is formed by means of several edges simultaneously. The number of working edges varies during machining and depends on the width of cut. For this reason, the total force F c of such interaction is a superposition of forces F c,i generated by each tool edge i, which are currently in the contact with the workpiece. Besides, the contact force F c,i can be decomposed by its radial F r,i and tangential F t,i components (Figure 4). In accordance with Merchant's model [START_REF] Merchant | Mechanics of metal cutting process. I-Orthogonal cutting and type 2 chip[END_REF], the tcomponent of cutting force F t,i can be computed with the equation ( 24). The r-component F r,i is related with F t,i by following expression (Laporte, 2009)

, , r i r t i F k F = ( 26
)
where the ratio factor k r depends on the given tool/workpiece characteristics.

It should be mentioned that in robotic machining it is more suitable to operate with forces expressed in the robot tool frame {x,y,z}. Then, the corresponding components F x , F y (Figure 4) of the cutting force F c can be expressed as follows

, , 1 1 , , 1 1 cos sin sin cos z z z z n n x r i i t i i i i n n y r i i t i i i i F F F F F F ϕ ϕ ϕ ϕ = = = = = - + = + ∑ ∑ ∑ ∑ (27)
where n z is the number of currently working cutting edges, φ i is the angular position of the i-th cutting edge (the cutting force in z direction F z is negligible here). So, the vector of external loading of the robot due to the machining process can be composed in the frame {x,y,z} using the defined components F x , F y as F=[F x ,F y ,0,0,0,0] T .

Projet It should be stressed that the cutting force components F r,i , F t,i mentioned in equation ( 24),(26) are computed for the given chip thickness h i , which should be also evaluated. Let us define model for h i using mechanical approach. Then the chip thickness h i removed by i-th tooth depends on the angular position φ i of this tooth and it can be evaluated using to the geometrical distance between the position of the given tooth i and the current machining profile (Figure 4). It should be mentioned, that the main issue here is to follow the current relative position between the i-th tooth and the working material or to define whether the i-th tooth is involved in cutting for given instant of process. Because of the robot dynamic behaviour and the regenerative mechanism of surface formation [START_REF] Tlusty | Basic non-linearity in machining chatter[END_REF] this problem cannot be solved directly using kinematic relations. In this case it is reasonable to introduce a special rectangular grid, which decomposes the workpiece area into segments and allows tracking the tool/workpiece interaction and the formation of the machining profile (Figure 5). , N w is the number of nodes) of the grid can be marked as "1" or "0": "1" corresponds to nodes situated in the workpiece area with material (rose nodes in Figure 6), "0" corresponds to nodes situated in workpiece area that was cut away (white nodes in Figure 6).

In order to define the number of currently cut nodes by the i-th tooth, the previous instant of machining process should be considered. Let us define A i as an amount of working material that is currently cut away by the i-th tooth (Figure 6). So, if node j marked as "1" is located inside the marked sector (green nodes in Figure 6), it changes to "0" and A i is increasing by x y s s Δ Δ . Analysing all potential nodes and computing A i , the chip thickness h i , removed at given instant of the process by the i-th tooth, can be estimated by ,

i i i h A R α = Δ 1, z i N =
. The angle Δφ i determines the current angular position of the i-th tooth regarding to its position at the instant τ-Δτ and referred to the position of TCP at τ-Δτ.

Figure 6 : Evaluating the tool/workpiece intersection A i and computing the corresponding chip thickness h i

Described mechanism of chip formation and the machining force model (24) allow computing the dynamic behaviour of the robotic machining process where models of robot inertia and stiffness are discussed in the section 3 of the paper. The detailed algorithm that is used in numerical analysis is presented in 0, where the analysis of the robot dynamics is performed in the tool frame with respect to the dynamic displacement of the tool δt dyn fixed on the robot end-effector around its position on the trajectory.

Modèle du procédé FSW

Le principe du FSW consiste à apporter une énergie d'origine mécanique par l'action d'un outil à l'interface entre les pièces à souder. La liaison se crée de proche en proche à l'état solide. Le procédé FSW a été mis en oeuvre, en premier lieu, dans les industries du transport, qui ont recours à l'aluminium pour alléger les structures mécaniques. En simulation, on utilise des modèles présentés dans le rapport Livrable 2.3 : Modèle exploitable pour la définition de la commande du robot. Les forces du soudage dépendent des paramètres physiques : matériaux des pièces à souder, épaisseur, forme et taille de l'outil, etc. Elles dépendent également de paramètres dynamiques : position et angle relatifs entre l'outil et la pièce, la vitesse avance, la vitesse de la rotation, etc.

Les paramètres opératoires avec lesquels on contrôle la phase de soudage sont :

• La vitesse d'avance de l'outil ou vitesse d'avance .

• L'effort axial appliqué sur l'outil suivant son axe de rotation.

• La fréquence de rotation .

• L'angle d'inclinaison de l'outil dénommé angle de tilt. Les paramètres opératoires lors de la phase de soudage changent en fonction de l'épaisseur à souder, le matériau et la géométrie de l'outil. L'ensemble conditionne l'apport d'énergie, le flux de matière, la formation du cordon, les propriétés mécaniques de l'assemblage et les efforts générés [START_REF] Mishra | Friction stir welding and processing[END_REF]].

Modèle :

Selon les conditions dans notre étude, les forces extérieures pendant le procédé FSW peuvent être modélisées dans le simulateur comme suit :

F A . V . . F F A . V . N

Corousso

RÉSULTATS DE SIMULATION EN USINAGE

Exemples avec un modèle d'usinage « Grid based Model »

CASE A: RIGID TOOL FIXATION. The objective of this case study is to understand the mechanics of the tool/workpiece interaction, while any dynamic aspect related to the robot compliance is excluded from the analysis. In that case the applied feed rate and spindle rotational speed totally determine the position {x i , y i },

1, z i N =
of all tool cutting edges referred to robotic base frame at each instant of machining τ

( )( ) sin , cos 2 1, 1 , i f i i i z i z x v R y R N i i N τ φ φ φ τ π = + = = Ω + - = (28)
The orientation α i of the i-th tooth velocity v i can be defined by the feed rate v f and the spindle rotational speed Ω at each instant as ( ) ( )

atan2 i f x y v +v v α Ω Ω = with sin , cos x i y i v R v R ϕ ϕ Ω Ω = Ω =Ω . The angle α i
provides computing the chip thickness h i as the displacement of TCP corresponding to the one tooth period 2 z N π Ω and referenced to the i-th tooth, when it is situated inside the working material. If the given tooth is located outside the working material, the corresponding chip thickness h i is equal to zero. The following expressions allow evaluating h i for all possible positions of the i-th tooth on its path while machining 0, sin , 2 , 1, 2 sin , 2

i z i i i i f z f z i i f z x R h x R x v N i N v N x v N α π π α π ⎧ < ⎪ = ≤ < Ω = ⎨ ⎪ Ω ≥ Ω ⎩ (29)
The advantage of the presented algorithm of computing the chip thickness is that different phases of tool/workpiece interaction illustrated in Figure 9 can be identified.

It should be mentioned that the phase of tool approaching to the workpiece corresponding to the zero machining force is not considered here. For the remaining phases a detailed analysis is presented below:

The phase of tool engaging into the working material (phase I) corresponds to the variable contact area between the tool and the workpiece. The TCP during this phase is located always outside the workpiece. In fact, the phase I can be divided into two sub phases: Phase Ia: In the beginning of milling operation a small area of workpiece is affected by the machining process. This fact and presence of two types of motion (feed, spindle rotation) form the case, when only one tooth can participate in cutting at the same time (Figure 9-a). Such behaviour produces intermittent machining forces F x and F y with the frequency 2 z N π Ω Hz. The sub phase 1 is very limited in time and its duration depends on the feed rate, the spindle rotational speed and the number of teeth N z . For example, if v f =4m/min, Ω=10 4 rpm, N z =4 the duration of phase Ia is only 0.04sec.

Phase Ib: It corresponds to a case, when several teeth can participate in cutting at the same time, but the TCP does not reach the workpiece border (Figure 9-b). As a result an oscillatory periodic behaviour in machining forces F x and F y is observed. But, because of different number of currently working teeth, the force patterns are not homogenous.

The phase of machining with fully engaged tool (Phase II) starts when TCP reaches the workpiece border (Figure 9-c). In that case always the same number of teeth (n c =2 for the tool with N z =4) is working at every instant of cutting process. It produces harmonic periodic machining forces

Corousso

So, even with the analysis of tool motion during machining without any dynamic aspect, the oscillatory behaviour of machining forces can be detected for whole process (Figure 10). The high frequency of such oscillation (for example, Ω=10 4 rpm, N z =4 give frequency of 667Hz) does not affect the motion of robot but it can be crucial for robot control system and should be considered in design of robotic machining process.

CASE B: TOOL FIXATION WITH COMPLIANCE IN X DIRECTION.

In contrast to Case A, the dynamic aspect of tool motion associated with the robot compliance is considered here. Thus, at each instant of machining, the position x TCP is defined as a superposition of tool displacement x f due to feed and a dynamic displacement δx due to compliance of the fixation:

TCP f x x x δ = + . The first component f f x v τ = is known at each instant of
process while the second one depends on the current position of the tool regarding to the machining profile. In this case, the dynamic displacement δx can be obtained by reducing the equation (1) to a one-dimensional problem and by introducing the damping related to the machining process and robot control algorithms

x x x M x C x K x F δ δ δ + + = && & ( 30
)
where M is the equivalent mass of the tool fixation, K x and C x are its stiffness and damping respectively. The damping 2

x x C K M ζ =
is related to the damping factor ζ, which can be estimated experimentally. The machining force F x acting in x direction, depends on the chip thickness h which in its turn depends on the relative position between teeth of the tool and the workpiece (i.e. δx) at each instant of machining.

Next, the position of the i-th tooth in the robot base frame {x,y} can be easily determined:

sin , cos , 1, z i T C P i i i x x R y R i N ϕ ϕ = + = = .
Comparative analysis of this position with respect to the current machining profile defines the chip thickness h i removing by i-th tooth. But, the main issue here is to define whether i-th tooth participates in cutting for given instant of process. For this reason, it is proposed to create a mesh on the workpiece, where each node j ( , N w is the number of nodes) can be filled with "1" or "0": "1" corresponds to nodes situated in the workpiece area with material, "0" corresponds to nodes situated in workpiece area that was cut away.

In order to define the number of currently cut nodes by the i-th tooth, the previous instant of machining process should be considered. Let us define A i as an amount of working material that is currently cut away by the i-th tooth (Figure 11). So, if node j filled with "1" is located inside the sector, it changes to "0" and A i is increasing by x y s s Δ Δ (Δs x , Δs y are node steps in x and y directions respectively). Analysing all potential nodes and computing A i , the chip thickness h i , removed at given instant of process by the i-th tooth, can be estimated by ,

i i i h A R α = Δ 1, z i N =
. The angle Δα i determines the current angular position of the i-th tooth regarding to its position at the instant τ-Δτ (Δτ is the time step) and referred to the position of TCP at τ-Δτ.

Figure 11 : Evaluating the tool/workpiece intersection A i and computing the corresponding chip thickness h i

Here, in contrast to the previous case, the dynamic aspect of the tool motion allows:

Estimating of deviation in tool motion from the desired one because of the robot compliance in the feed direction. It should be mentioned that this deviation affects the Cartesian stiffness of robot but does not influence the machining profile quality. For example, following parameters M=100 kg, K x = 5 3 10 × N/m, ζ=0.05, N z =4, Ω=10 4 rpm, v f =4 m/min provide deviation in the feed direction of 0.13mm (Figure 12). But it is not essential for this application.

Detecting of vibratory behaviour of tool motion while it is engaging into the workpiece. In some cases the low frequency of such motion can excite robot natural frequency, destabilize machining operation and may even damage the tool or/and workpiece. For example, the milling process with following parameters M=100 kg, K x = 5 3 10 × N/m, ζ=0.05, N z =4, Ω=10 4 rpm, v f =4 m/min generates vibration of f 1 =8.8Hz from the beginning of cutting process in addition to the frequency 667Hz related to the spindle rotation. More details on influence of the tool fixation parameters (M, K x ) on the dynamics of its motion during machining are presented in

Corousso

Suitable robot configurations to perform given machining operation could be defined. But, it should be mentioned that this one dimensional equivalent model presented here allows analysis of machining process dynamics in the feed direction only. In order to evaluate behaviour of the tool motion more closely to the real machining operation, this model should be extended. ,

M r,xx =M

f f x f f y x v y v τ τ = =
. The dynamic displacements δx, δy could be obtained from the equation (1) which, in this case, is reduced to

x y F x x x F y y y δ δ δ δ δ δ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ + + = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ r r Μ C Κ && & && & (31)
where M r (2×2) is the equivalent mass matrix of the tool fixation, the matrices K r (2×2) and C (2×2), where

Corousso

In order to illustrate the advantages of this two dimensional model and its ability to detect some phenomena (that are not visible in Cases A, B) the robotic milling process is simulated for KUKA KR 270 robot with the following parameters 5 5 100 0 3 10 0 550 0 , kg, , , , 0 100 0 550 0 3 10 sec

N k g m ⎡ ⎤ × ⎡ ⎤ ⎡ ⎤ = = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ × ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ r r M K C .
Here, the equivalent mass matrix M r is computed in accordance with the method presented in the Section 2 of this paper. The stiffness K r is the structural stiffness of the robot, referred to its end-effector. It should be noted that in practice the non-diagonal elements in these matrices are non-zero. But, with this simplified case it is possible to identify qualitatively the dynamic nature of the tool behaviour in the direction (y) orthogonal to the feed. Simulation results corresponding to this case study are presented in Fig. 9.

It should be also stressed that the tool dynamic behaviour in the feed direction (x) is similar to the results, which were obtained in the Case B. The displacement in y-direction has an essential dynamic component during the phase of tool engagement (Phase I) into the workpiece and becomes constant while machining with fully engaged tool. The corresponding frequency f 1y =7.3Hz is comparable with the frequency f 1x =8.2Hz of the tool dynamic displacement in feed direction. More details on the tool motion during milling process regarding to the parameters of the tool fixation are presented in Table 2 (it covers the range of values for M r,xx , M r,yy , K r,xx , K r,yy computed for different configurations of the robot KUKA KR 270 using the methodology presented in Section 2)

It should be noted that changing the robot configuration affects the tool dynamics in y direction which is crucial regarding to the quality of the final product. Hence, considering the tool displacement orthogonal to the feed direction is essential and the Case C gives more realistic results comparing Cases A and B. As it is shown in 0, the deviation (0.31 -0.56mm) in machining profile from the desired one has a vibratory behaviour during the phase of tool engagement into the workpiece. Thus it cannot be suppressed by straightforward compensation methods and other compensation techniques should be proposed.

Exemple avec un modèle d'usinage à effort spécifique de coupe

Le modèle d'effort de coupe choisi est adapté à une opération de détourage. L'outil est une fraise cylindrique revêtue de grains abrasifs. L'outil est solidaire de la partie terminale de l'outil et son axe de rotation est perpendiculaire au plan du travail. On suppose que la trajectoire du bout d'outil reste dans le plan , . Les forces sont exercées au point par la pièce sur l'outil qui est fixé au 6 ème corps du robot. Les forces peuvent être calculées par une modélisation des efforts de coupe. En première approximation, on peut considérer que la force de coupe est proportionnelle à la section de matière enlevée et à l'effort spécifique de coupe . Lorsque le robot est soumis à un effort, les raideurs en torsion non infinies entraînent un déplacement de celui-ci et donc une modification de la section de matière enlevée comme le montre la Figure 14. On obtient l'effort de coupe par : avec la profondeur de coupe et représentant la largeur de matière enlevée. L'effort spécifique de coupe dépend de la matière, de l'outil et de la vitesse d'avance qui est imposée par le robot. Dans le cas du détourage considéré sur la Figure 14, la vitesse d'avance est définie par qui dépend à l'instant t des déformations du robot dans la direction . La modélisation de la coupe des matériaux composites n'a pas fait l'objet de travaux aussi approfondis que la coupe des métaux. Par analogie avec les modèles de coupe considérant un cisaillement primaire, l'effort spécifique peut être exprimé par :

  ............... .................. .................... .................... .................. .................. .................... .................... .................... .................. .................... .................... .................. .................... .G, S.C, F A, J.Q .................. .................. .................... .................... .................. .................. .................... .................... .................... .................. .................... .................... .................. .................... .................. .................. en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4, S.G, S.C, F.L, G.A., J.Q. Corousso Projet COROUSSO Livrable n°4.3 Résultats en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4

  Projet COROUSSO Livrable n°4.3 Résultats en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4Modèle simplifié du comportement flexible du robot

  , where i G is the loading applied to the i-th node-point. Besides, it is assumed that the external loading F (caused by the interaction of the tool and the workpiece) is applied to the robot endeffector.Projet COROUSSO Livrable n°4.3 Résultats en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4

  where Projet COROUSSO Livrable n°4.3 Résultats en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4, S.G, S.C, F.L, G.A, J.Q

  en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4, S.G, S.C, F.L, G.A., J.Q.

  de flexibilités dans l'espace articulaire La flexibilité des robots est souvent bien analysée sur les robots légers, par exemple, le Kuka-DLR LWR, le Barrett WAM, le Bras ASSIST, etc. Les robots manipulateurs lourds comme les robots utilisés dans cette étude sont souvent considérés comme rigides. Mais pour une opération qui demande des efforts très importants, la déformation à cause de la flexibilité ne peut pas être négligée, si l'on désire obtenir de bonnes performances de fabrication. De nombreuses recherches ont été faites sur la modélisation des robots. Les modèles les plus souvent utilisés pour les robots manipulateurs sont décrits dans les ouvrages (Craig et al., 1986), (Sciavicco and Siciliano, 2000) et (Spong et al., 2005). Les principes fondamentaux et les méthodes qui sont utilisées pour modéliser, concevoir et piloter un système robotique ainsi que la planification et la commande d'un robot réel sont discutés dans (Siciliano and Khatib, 2008).

  Figure 1 :

Γ

  en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4.3 15/09/2014 indice A Page 12/34 Auteurs : A.K, A.P, S.G, S.C, F.L, G.A., J.Q.

  des effets centrifuges, de Coriolis et de gravité et le système de compensation de gravité. Dans notre modèle, on tient donc également compte de l'effet du compensateur de gravité,• Bm et Bs : les matrices des coefficients de frottement visqueux au niveau des arbres moteurs et aux articulations en Nm/[rad/s],

  Figur

Figure 3 :

 3 Figure 3 : Fractional cutting force model F c (h)

Figure 4 :Figure 5 :

 45 Figure 4 : Forces of tool/workpiece interaction

  Projet COROUSSO Livrable n°4.3 Résultats en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4.3 15/09/2014 indice A Page 17/34 Corousso Here, Steps Δs x , Δs y between grid nodes are constant and depend on the tool geometry, cutting condition and time discretization Δτ. Each node j ( 1, w j N =

Figure 7 :

 7 Figure 7 : Algorithm for numerical simulation of robotic machining process dynamics

Figure 8 :

 8 Figure 8 : Algorithm for numerical simulation of robotic machining process dynamics

Figure 9 :

 9 Figure 9 : Different phases of tool/workpiece interaction (case N z =4) and corresponding machining forces; (a) -tool entry into the workpiece, only one tooth can be in contact with workpiece at the same time, (b) -increasing tool engagement into the workpiece, several teeth could be in contact with workpiece at the same time, (c) -slot machining with fully engaged tool.

  F x and F y with the frequency 2 , S.G, S.C, F.L, G.A., J.Q.

Figure 10 :

 10 Figure 10 : Machining force patterns with average forces referenced to the frame {x, y}.

  , S.G, S.C, F.L, G.A., J.Q.

Figure 12 :

 12 Figure 12 : Machining force patterns and TCP dynamic displacement in case of 1DOF model (M=100 kg, K x =

Figure 13 :

 13 Figure 13 : Machining force patterns and TCP dynamic displacements in case of 2DOF model (M r,xx =M r,yy =100 kg, K r,xx =K r,yy =

  Figure14

  Figure 16 : Co

Corousso

Valeurs des raideurs articulaires des robots utilisés.

										5
	KR500 MT-2	6.61 ×10 6	7.16 ×10 6	3.08 ×10 6	-	-		-
	Tableau 1 : N ij N 11 N 22	N 33	N 44	N 55	N 66	N 54	N 64	N 65
	KR270-2F	-219.00 -286.00		-	-	-	157.552	1.002	2.020	2.520
			270.177	222.158	221.765
	KR500 MT-	469.375 469.375		-	-	-	164.570 -1.0964 -1.5836	1.5311
	2		504.770	260.619	255.977

Tableau 2 : Valeurs des coefficients de réduction des transmissions mécaniques des robots utilisés.

  . F

			Projet COROUSSO		ANR-10-SEGI-003-LI 4.3
			Livrable n°4.3
			Résultats en simulations de l'utilisation
		d'un robot avec flexibilités pour l'usinage et	15/09/2014
				le soudage		indice A
						Page 20/34
	Par identification, on peut calculer les valeurs des coefficients du modèle avec la méthode des moindres
	carrés. On choisit les valeurs suivantes pour les simulations :	0,54; α	0,938 ; β	0,647; μ
	0,021 et	13 . 10 ; α	0,596 ; β	0,37; μ	0,02
	Évolution des variables principales en FSW
	• Identification des phénomènes limites aux frontières du domaine de fonctionnement
	• Évolution des composantes du torseur des interactions mécaniques outil/matière dans le domaine de
	fonctionnement
	• Modélisation statistique des composantes du torseur des efforts en fonction des principaux
	paramètres de conduite du procédé

Auteurs : A.K, A.P, S.G, S.C, F.L, G.A., J.Q.

Table 1 ,

 1 which covers the range of values for M, K x computed for different configurations of the robot KUKA KR 270. M, kg K x , N/m τ s , sec |x s |, mm PO, % f 1 , Hz

	100	0.05 10 ×	6	1.2	0.80	52	3.5
	100	0.30 10 ×		6	0.6	0.13	30	8.8
	100	0.60 10 ×		6	0.5	0.07	23 12.4
	100	1.00 10 ×	6	0.4	0.04	17 15.6
	100	2.00 10 ×	6	0.4	0.02	11 22.6
	150	2.00 10 ×	6	0.5	0.02	14 18.3
	200	2.00 10 ×	6	0.5	0.02	18 15.6

Table 1 : Influence of tool fixation parameters on the tool motion; τ s is the settling time, PO is the overshoot, f 1 is the first frequency of the tool dynamic displacement in feed direction ζ=0.05, N z =4.10 4 rpm, v f =4 m/min

 1 As it can be observed from the table, changing the fixation parameters (i.e. the robot configuration) influences low frequencies (about 10 -20 Hz) of the tool motion.

	Projet COROUSSO	ANR-10-SEGI-003-LI 4.3
	Livrable n°4.3
	Résultats en simulations de l'utilisation
	d'un robot avec flexibilités pour l'usinage et	15/09/2014
	le soudage	indice A
		Page 25/34

Table 2 : Influence of the tool fixation parameters on tool dynamic behavior; f 1x , f 1y are first frequencies of the tool dynamic displacement in x and y directions respectively; N z =4. 10 4 rpm, v f =4 m/min CASE C: TOOL FIXATION WITH COMPLIANCE IN X AND Y DIRECTIONS.

 2 In this case, a dynamic aspect of tool motion in feed direction (x) and orthogonal to it (y) is considered. Then, similarly to the Case B, at each instant of machining process, the position of TCP is determined by ,

			K r,xx =K r,yy ,	|x s |,	|y s |,	f 1x ,	f 1y ,
		r,yy , kg	N/m			mm	mm	Hz	Hz
		100	0.05 10 ×	6	0.80 2.61	3.2	2.2
		100	0.30 10 ×		6	0.13 0.43	8.2	7.3
		100	0.60 10 ×		6	0.07 0.22 12.4	11.4
		100	1.00 10 ×	6	0.04 0.13 15.6	14.6
		100	2.00 10 ×	6	0.02 0.06 22.4	21.5
		150	2.00 10 ×	6	0.02 0.06 18.3	17.4
		200	2.00 10 ×	6	0.02 0.06 15.6	15.6
								TCP x	f = + x	δ	TCP x y	f = + y	δ	y	, where
	,	,

  F x and F y are the machining forces in x and y directions.

	, i j C experimentally), In contrast to the previous case, the position of the i-th tooth at each process instant t includes dynamic , , , 2 i j i j i j ζ = K M characterize the fixation stiffness and damping respectively (which can be estimated components in both directions: sin , i T C P i x x R ϕ = + cos , i T C P i y y R ϕ = + 1, z i N = . The algorithm of computing Livrable n°4.3 Résultats en simulations de l'utilisation d'un robot avec flexibilités pour l'usinage et le soudage ANR-10-SEGI-003-LI 4.3 15/09/2014 indice A the chip thickness h Projet COROUSSO Page 26/34

i for given position of tooth {x i , y i } is similar to Case B. Auteurs : A.K, A.P, S.G, S.C, F.L, G.A., J.Q.

Corousso

CONCLUSION

Dans cette tâche, des modèles simplifiés et avec flexibilités localisées sont proposés pour traduire au mieux les déformations des robots manipulateurs utilisés dans les procédés de production. Les modèles de déformations peuvent être utilisés dans l'espace articulaire ou dans l'espace cartésien. Ces modèles sont ensuite intégrés dans une simulation qui comprend également la modélisation dynamique des robots, des actionneurs et des transmissions mécaniques et de la commande.