Tiantian Xu
email: tiantian.xu@isir.upmc.fr

Gilgueng Hwang
email: gilgueng.hwang@lpn.cnrs.fr

Nicolas Andreff
email: nicolas.andreff@femto-st.fr

Stéphane Régnier
email: stephane.regnier@upmc.fr

Planar Path Following of 3D Steering Scaled-up Helical Microswimmers

Keywords: helical swimmer, magnetic actuation, 3D steering, path following, visual servoing, 3D pose estimation

Helical microswimmers capable of propulsion at low Reynolds numbers have great potential for numerous applications. Several kinds of artificial magnetic actuated helical microswimmers have been designed by researchers. However, they are primarily open-loop controlled. This paper aims to investigate methods of closed-loop control of a magnetic actuated helical swimmer at low Reynolds number by using visual feedback. For many in-vitro applications, helical swimmers should pass through a defined path, for example along channels, without any prerequisite on the velocity profile along the path. Therefore, the main objective of this paper is to achieve a velocity-independent planar path following task. Since the planar path following is based on 3D steering control of the helical swimmer, a 3D pose estimation of a helical swimmer is introduced based on the realtime visual tracking with a stereo vision system. The contribution of this paper is in two parts: the 3D steering of a helical swimmer is demonstrated by visual servo control; and the path following of a straight line with visual servo control is achieved, then compared to open-loop control. We further expect that with this visual servo control method, the helical swimmers will be able to follow reference paths at the microscale.

I. INTRODUCTION Helical microswimmers capable of propulsion at low

Reynolds numbers have been proposed for numerous applications, ranging from in vitro tasks on labs-on-chips (e.g. transporting and sorting micro objects; mechanical components micro assembly [START_REF] Chaillet | Microrobotics for Micromanipulation[END_REF]...) to in vivo applications for minimally invasive medicine (e.g. targeted drug delivery; brachytherapy; hyperthermia [START_REF] Nelson | Microrobots for minimally invasive medicine[END_REF]...), due to their micro sizes and accessibility to tiny and clustered environments. However, at the microscale, the fluid becomes extremely viscous and the Reynolds number of the microswimmers decreases dramatically. The corkscrew type rotating propulsion used by some natural micro-organism, for example E. coli bacteria, is one of the swimming technique suitable for swimming at low Reynolds numbers [START_REF] Purcell | Life at low reynolds number[END_REF], [START_REF] Purcell | The efficiency of propulsion by a rotating flagellum[END_REF]. Artificial helical microswimmers, driven by a rotating magnetic field, are inspired by E. coli bacteria [START_REF] Berg | Bacteria swim by rotating their flagellar filaments[END_REF]. They are able to convert rotary motion to a linear motion.

For more than ten years, researchers have developed several kinds of magnetic actuated helical swimmers or crawlers from several millimeters [START_REF] Honda | Micro swimming mechanisms propelled by external magnetic fields[END_REF], [START_REF] Ishiyama | Spiral-type micro-machine for medical applications[END_REF], [START_REF] Xu | Modeling and swimming property characterizations of scaled-up helical microswimmers[END_REF] down to a few micrometers [START_REF] Bell | Fabrication and characterization of three-dimensional ingaas/gaas nanosprings[END_REF], [START_REF] Zhang | Artificial bacterial flagella: Fabrication and magnetic control[END_REF], [START_REF] Hwang | Electro-osmotic propulsion of helical nanobelt swimmers[END_REF], [START_REF] Pooyath | Conformal cytocompatible ferrite coatings facilitate the realization of a nano-voyager in human blood[END_REF]. Ghosh et al. proved that helical swimmers with an overall length of 2 µm could follow a curved trajectory, for example "R@H", under the action of a preprogrammed magnetic field [START_REF] Ghosh | Controlled propulsion of artificial magnetic nanostructured propellers[END_REF]. Mahoney et al. showed an open-loop controlled helical swimmer with an overall length of 5 mm, which enabled a U-turn trajectory in the vertical plane [START_REF] Mahoney | Velocity control with gravity compensation for magnetic helical microswimmers[END_REF]. Tottori et al. fabricated helical swimmers ranging from 4 µm to 64.5 µm, which enabled 3D motion and could carry a sphere as well as raise it up [START_REF] Tottori | Magnetic helical micromachines: Fabrication, controlled swimming, and cargo transport[END_REF]. As far as we know, the existing helical swimmers are primarily open-loop controlled, mainly due to the challenges of visual feedback and tracking through micro and nanoscopic imaging.

A closed-loop control of helical swimmers is necessary for precise motion in presence of perturbations generated by thermal noise [START_REF] Ghosh | Velocity fluctuations in helical propulsion: How small can a propeller be[END_REF] or drifting due to boundary effects [START_REF] Peyer | Non-ideal swimming of artificial bacterial flagella near a surface[END_REF]. The aim of this paper is to control a magnetic actuated helical swimmer using visual feedback. Visual servo controls show good results with industrial robots [START_REF] Gangloff | High-speed visual servoing of a 6-dof manipulator using multivariable predictive control[END_REF], [START_REF] Dahmouche | Dynamic visual servoing from sequential regions of interest acquisi-tion[END_REF] or unmanned vehicles [START_REF] Courbon | Vision-based navigation of unmanned aerial vehicles[END_REF], but have not been applied to helical swimmers. A Scaled-up Helical Microswimmer (SHM) at the millimeter scale will be used, in order to eliminate the optical tracking problems. The fluid motion around the helical swimmer is only characterized by the Reynolds number. Therefore, the swimming properties of helical swimmers should be the same if they swim at the same Reynolds numbers, which makes a dimensionless study possible [START_REF] Xu | Modeling and swimming property characterizations of scaled-up helical microswimmers[END_REF]. At a larger scale, to emulate the environment of microscale swimmers swimming at low Reynolds numbers, one can use a more viscous liquid, such as glycerol. The SHMs facilitate fabrication and can be observed with regular optical devices. They are therefore good candidates for preliminary visual servoing experiments. Once the visual servo control is validated on the SHMs, it can be applied to helical microswimmers by using suitable optical devices with a higher resolution in future works. The SHMs can be scaled down to the microscale with the recently developed microswimmer micro-fabrication processes, such as the self-scrolling of a bilayer [START_REF] Bell | Fabrication and characterization of three-dimensional ingaas/gaas nanosprings[END_REF], the glancing angle deposition method [START_REF] Ghosh | Controlled propulsion of artificial magnetic nanostructured propellers[END_REF], and microscale 3D lithography [START_REF] Tottori | Magnetic helical micromachines: Fabrication, controlled swimming, and cargo transport[END_REF].

For most in-vitro applications of helical swimmers, such as labs-on-chips, they should pass through a defined path, for example pass through micro channels, with a constant altitude without any prerequisite on the velocity profile along the path. Therefore, the main objective of this paper is to establish a velocity-independent control law of the helical system to follow a planar reference path based on 3D steering. The steering is in the 3D space, because the helical swimmer has an inclination angle in order to decrease the friction with the substrate. In this paper, Section II depicts the visual servo control architecture of the path following based on its 3D steering. Section III presents the modelling of the magnetic actuated helical swimmer, including the linear propulsion and the 3D steering. Section IV details the kinematic equations of the planar path following task of a helical swimmer, and the control law design. Section V presents the visual servo control for 3D steering based on the orientation error of the helical swimmer. Section VI presents the overview of the experimental system, including the magnetic actuation system, stereo vision system, and the SHM at low Reynolds number. Section VII introduces a 3D pose estimation including the position and the orientation of the SHM based on its real-time tracked barycenter and axis. The results of the 3D steering and the planar path following are given in Section VIII. Finally, Section IX concludes the paper.

II. CONTROL ARCHITECTURE

A planar reference path, as shown in Fig. 1, is defined to be followed by a helical swimmer. A helical swimmer can be considered as a unicycle-like nonholonomic system, because one helical swimmer cannot generate a lateral motion voluntarily. As explained in [START_REF] Bloch | Nonholonomic dynamics[END_REF], "Nonholonomic systems are, roughly speaking, mechanical systems with constraints on their velocity that are not derivable from position constraints". All the lateral motions applied on a helical swimmer generated by thermal noise [START_REF] Ghosh | Velocity fluctuations in helical propulsion: How small can a propeller be[END_REF], or drifting due to boundary effect [START_REF] Peyer | Non-ideal swimming of artificial bacterial flagella near a surface[END_REF] can be considered as a perturbation, and should be corrected. Compensation of thermal noises is both a motivation of the closedloop control and a major technical challenge. Concerning the former, it is indeed very delicate to follow a time-dependent trajectory (e.g. to follow a virtual robot) in presence of thermal noise because thermal noise will yield time delays with the virtual robot on the trajectory and thus, the geometrical path will not be achieved. As for the latter, we have not investigated yet where the frontier lies between efficient propulsion at low Reynolds number and predominance of thermal noise.

Moreover, the velocity amplitude of a helical swimmer can not be reliably achieved in many cases. The first reason for that is the well-known cut-off frequency of the velocity amplitude [START_REF] Honda | Micro swimming mechanisms propelled by external magnetic fields[END_REF], which imposes a strong saturation of the available amplitude. As we have shown in [START_REF] Xu | Modeling and swimming property characterizations of scaled-up helical microswimmers[END_REF], this cut-off frequency depends on many uncontrolled factors (microfabrication quality and environmental viscosity, for the main ones). The second reason is related to the potential stiction of the helix onto the substrate, which creates a distortion in velocity amplitude. The third reason is related to the viscosity at low scales (through the Reynolds number). Indeed, at low scales, we are not sure that the fluid is really homogeneous (inherently or due to the presence of micro or nano-dirt in a non totally clean room). This is especially true when dealing with biological fluids (blood, cell culture) where the medium is definitely heterogeneous or, at least, where the cells represent mobile obstacles that need to be pushed aside to allow for the helical swimmer to follow the expected trajectory.

Disturbances have several origins (thermal noise, friction, nonhomogeneous fluid, cut-off frequency, whatever). These disturbances split in two effects. The first one is a lateral deviation with respect to the desired path and the second one is a delay with respect to the time-varying reference on the trajectory. Namely, the helical swimmer can effectively be on the desired trajectory but it is late. In such a case, trajectory tracking algorithms tend to resume the delay by "cutting the curves" (Fig. 1), whereas path planning algorithms consider that there is no error since the path is time-independent. Note that both kinds of algorithms tend to reject the lateral deviation. Our proposition is thus to merge non-holonomous path following algorithms with the specificity of the helical swimmer actuation. To do so, we chose to decompose the problem in two tasks: the first one is to keep a constant altitude and the other is to follow a path in the horizontal plane. The planar path following of a helical swimmer is based on its steering similar to a unicycle, as presented by a blockdiagram in Fig. 2. The notation will be detailed in Section IV and VII. The steering is in the 3D space, because the helical swimmer has an inclination angle to maintain its altitude as well as to decrease friction with the substrate. The realtime position and orientation of the helical swimmer can be obtained by a 3D pose estimation method based on two images taken by a stereo vision system, which will be introduced in Section VII. The distance and orientation error between the helical swimmer and the reference path can be then obtained and used to deduce the target orientation of the helical swimmer, which will be detailed in Section IV. The control of the magnetic field for the 3D steering based on the orientation error will be detailed in Section V.

III. MODELLING OF MAGNETIC ACTUATED HELICAL

SWIMMERS

A. Magnetic force and torque

All magnetized objects are subjected to both force and torque within an externally imposed magnetic field. The magnetic force and torque induced on a magnetized object are expressed as follows [START_REF] Jiles | Introduction to magnetism and magnetic materials[END_REF]:

f m = V m (M • )B dV m (1)
τ m = V m M × B dV m (2)
where V m is the volume of the magnetized object, B is the flux density of the applied field (T ), M is the magnetization of the object (A/m ).

In this paper, the helical swimmer is actuated by a rotating magnetic field generated by a 3D Helmholtz coils system. The magnetic field can be considered as uniform. Therefore, the helical swimmer is driven by magnetic torque only.

B. Inclination angle and direction angle

In order to describe more clearly the movement of the helical swimmer, we define herein two angles: the inclination angle θ i and the direction angle θ d as shown in Fig. 3 . The inclination angle θ i is defined as the angle between the axis of the helical swimmer and the horizontal plane yOz. This inclination angle is very important to control the height of the helical swimmer. The propulsive force projected on the upward direction can be increased by increasing the inclination angle or rotation frequency. If the propulsive force created by the helical swimmer projected on the upward direction is higher than its own gravity, the helical swimmer will go upwards. Otherwise, the helical swimmer will go downwards. Helical swimmers are kept inclined in order to decrease the surface friction as well.

The direction angle θ d is defined as the angle between the z axis and the projection of the axis of the helical swimmer on the horizontal plane. By controlling the direction angle of the helical swimmer, its forward direction projected on the horizontal plane can be changed.

The normal vector of the helical swimmer axis is noted as n, which is defined as the orientation of the SHM, and can be expressed as:

n =   n x n y n z   =   sin θ i cos θ i sin θ d cos θ i cos θ d   (3)

C. Modelling of linear advance

Purcell proved in his paper that the nonfluidic applied force ( f a ) and nonfluidic applied torque (τ a ) on a helical swimmer along its axis are linearly related to its propulsion velocity (v a ) and angular speed (ω a ) along its axis with four principal quantities described by a symmetric propulsion matrix [START_REF] Purcell | Life at low reynolds number[END_REF]:

f a τ a = a b b c v a ω a (4)
The parameters a, b, and c encapsulate the geometric and environmental properties of the helical swimmer. The nonfluidic applied forces are gravity, buoyancy, friction and magnetic forces, while the only nonfluidic torque is the magnetic torque.

The advance velocity v h is defined as the linear velocity projected on the horizontal plane, which can be expressed as:

v h = v a cos θ i (5)
where v a is the propulsion velocity along its axis, θ i is the inclination angle.

IV. PLANAR PATH FOLLOWING

A. Notations

The notations of the planar path following of the helical swimmer are depicted in Fig. 4: Fig. 4: Helical swimmer modelling of the planar path following.

• C (s) is the predefined reference path with respect to the curvilinear coordinate s of the curve. The curve should be twice differentiable in R 2 [START_REF] Cherubini | A position-based visual servoing scheme for following paths with nonholonomic mobile robots[END_REF], [START_REF] Cherubini | An image-based visual servoing scheme for following paths with nonholonomic mobile robots[END_REF]. • θ d is the real-time direction angle of the SHM.

• v h is the advance velocity of the SHM on the horizontal plane.

• G and G are the positions of the SHM's barycentre at time t and t + dt. • s is the curvilinear coordinate of C (s). S and S are the nearest points on the path to the SHM at time t and t +dt. • d = SG is the distance between the SHM and the path, which denotes the lateral error.

• - → i s is the tangent of path C at S. - → j s is perpendicular to - → i s .
• c denotes the curvature of the path C , and R denotes the radius of curvature at S, with R = 1/c. • θ c is the angle between the current tangent of the path and the z axis, which can be considered as the direction angle of the path.

• θ e = θ dθ c defines the error of the direction angle of the SHM to the path. • θ * d is the calculated target direction angle of the SHM for the control to achieve the path.

B. Kinematic equations

The aim of the path following control is to minimize the lateral and angular errors. Following [START_REF] Laumond | La robotique mobile. Traité IC[END_REF], [START_REF] Thuilot | Automatic guidance of a farm tractor relying on a single cp-dgps[END_REF], we propose here to express the state coordinates with respect to the path C . More precisely, the position G and the direction angle θ d of the SHM are described by the curvilinear coordinate s, the distance error d, and the direction angle error θ e . The state vector is then written as:

q =   s d θ e   (6)
The control variables of the system are the linear advance velocity and the angular velocity of the steering projected on the horizontal plane:

u = v h Ω x (7)
According to a simple geometry relationship, the state space model of the SHM moving on the horizontal plane is expressed as [START_REF] Laumond | La robotique mobile. Traité IC[END_REF], [START_REF] Thuilot | Automatic guidance of a farm tractor relying on a single cp-dgps[END_REF]:

         ṡ = v h cos θ e 1 + d/R ḋ = v h sin θ e θe = Ω x - v h cos θ e R + d (8)

C. Control law designed with chained form

The control objective is to ensure the convergence of the SHM to the reference path C . Therefore, the state variables d and θ e are expected to be brought and kept at zero. d = 0 means that the barycenter of the SHM belongs to path C . θ e = 0 ensures that the SHM will follow the path C .

The state space model of the SHM ( 8) is clearly nonlinear. Samson has established in control theory that models of mobile robots can be converted into almost linear models, namely chained forms, in an exact way [START_REF] Samson | Control of chained systems application to path following and time-varying point-stabilization of mobile robots[END_REF]. Such an approach is attractive and used by lots of researchers for car-like nonholonomic mobile robots [START_REF] Laumond | La robotique mobile. Traité IC[END_REF], [START_REF] Thuilot | Automatic guidance of a farm tractor relying on a single cp-dgps[END_REF], since it allows us to use linear systems theory while still relying upon the actual nonlinear model. This approach is followed hereafter.

The model of the SHM (8) can be converted into chained form dedicated to systems with two inputs and limited to dimension three. The general chained system of dimension three is written as:

   ẋ1 = m 1 ẋ2 = x 3 m 1 ẋ3 = m 2 (9)
where X = (x 1 , x 2 , x 3 ) T and M = (m 1 , m 2 ) T are respectively the state and control vectors. Following [START_REF] Thuilot | Automatic guidance of a farm tractor relying on a single cp-dgps[END_REF], the model of the SHM is converted into chained form by using the change of the coordinates:

x 1 = s ( 10
)
m 1 = ẋ1 = v h cos θ e 1 + dc (11)
x 2 = d (12) ẋ2 = v h sin θ e = x 3 m 1 ( 13
)
Therefore, from ( 11) and ( 13), the last state variable x 3 is given by:

x 3 = (1 + dc) tan θ e (14)
Finally, considering ( 9) and ( 8), the last control variable m 2 is given by:

m 2 = ẋ3 = cv h sin θ e tan θ e + dc ds v h cos θ e 1 + dc d tan θ e + 1 + dc cos 2 θ e Ω x - v h cos θ e c 1 + dc (15)
These transformations are invertible as long as d = -

1 c , v h = 0 and tan θ e = π 2 [π].
From a practical point of view, the SHM respects these conditions once properly initialized. Now that the model of the SHM is converted to a linear model, the control techniques for linear systems can be applied. As shown in [START_REF] Thuilot | Automatic guidance of a farm tractor relying on a single cp-dgps[END_REF], the simple expression for the virtual control law is:

m 2 m 1 = -k t x 3 -k d x 2 (k t , k d ) ∈ R +2 (16)
where k t and k d are the control parameters and are strictly positive, which implies that both x 2 and x 3 converge to 0. In view of ( 12) and ( 14), d and θ e converge to 0 as well. The path following is therefore achieved. Moreover, since the evolution of the error is driven by x 1 = s, the gains (k t , k d ) impose a settling distance instead of a settling time. Consequently, the helical swimmer trajectory will be identical if the initial error is the same. The guidance performances will be velocity independent, which are suitable for low Reynolds numbers environments.

The control law on the angular velocity of the steering of the SHM is expressed as follows, by reporting ( 11) and ( 15) into [START_REF] Ghosh | Velocity fluctuations in helical propulsion: How small can a propeller be[END_REF]:

Ω x = - cos 3 θ e 1 + dc k d d + k t (1 + dc) tan θ e + dc ds d tan θ e +c tan 2 θ e (1 + dc) + v h cos θ e R + d (17)
where v h is the advance velocity projected on the horizontal plane, which is usually constant, and depends on the frequency of the rotating field and the inclination angle. In the case where the reference path is a straight line, the control of the angular speed of the steering projected on the horizontal plane can be simplified as:

Ω x = -cos 3 θ e (k d d + k t tan θ e ) (18)
The target direction angle can be calculated by the real-time direction angle and the angular velocity of the steering of the SHM, which can be expressed as:

θ * d (t) = θ d (t) + Ω x dt (19)
The target orientation of the SHM n * can be calculated with the target direction angle θ * d and target inclination angle θ * i with (3). The target inclination angle is fixed for many tasks.

The path following problem is now converted into a 3D steering problem of the SHM from the real-time orientation n to a target orientation n * . The control of the magnetic field which enables the steering of the SHM from n to n * will be studied in the next section.

V. 3D STEERING CONTROL

A. Modelling of steering

The magnetization of the helical swimmer is noted as M, as shown in Fig. 3. The magnetic torque τ exerted on the helical swimmer can be expressed as :

τ = M × B ( 20
)
The magnetic field B can be decomposed in the directions perpendicular to its axis B ⊥n , and colinear to its axis B n . As the magnetization of the helical swimmer is always perpendicular to its own axis n, the magnetic torque can also be decomposed in the two directions:

τ = M × B ⊥n τ n + M × B n τ ⊥n (21)
The torque τ n colinear to n makes the helical swimmer rotate around its current axis. The torque τ ⊥n perpendicular to n makes the helical swimmer rotate around M × n. It means that the axis of the helical swimmer changes the direction, which enables the 3D steering of the helical swimmer.

B. Self-rotation control

The magnetic field B ⊥n rotating around the axis n can be expressed as:

B ⊥n = B 0 cos(2π f t) ũ + B 0 cos(2π f t)ṽ ( 22
)
where B 0 is the magnetic flux density at the center of the Helmholtz coils, f is the rotation frequency, and ( ũ, ṽ) are the basis vectors of the plane orthogonal to the axis n:

ũ = u u = 1 u   n y -n x 0   (23) ṽ = v v = 1 v   n x n z n y n z -n 2 x -n 2 y   (24)
with n x = 0 or n y = 0. Remark that if n x = 0 and n y = 0, we can chose ũ = (1, 0, 0) T and ṽ = (0, 1, 0) T .

C. Steering control

The magnetic field B n , which contributes to the steering, depends on the orientation error of the SHM. The orientation error is defined as the sine of the angle between the realtime orientation and the target orientation: sin(n, n * ), which is equivalent to n × n * .

Let us define a control proportional to this error. As shown in Fig. 5, the direction of B n depends on the projection of B ⊥n on the target orientation n * . If this projection is positive B ⊥n • n * > 0, the magnetic field B n should be in the direction of -n. If this projection is negative, B n should be in the same direction as n. Therefore, the control law of B n can be expressed as:

B n = -sign(B ⊥n • n * ) • λ n × n * n (25)
with λ the control parameter.

In conclusion, the control law of the magnetic field B which enables the 3D steering of helical swimmers can be expressed as:

B = B ⊥n + B n (26)
by replacing B ⊥n by ( 22) and B n by [START_REF] Laumond | La robotique mobile. Traité IC[END_REF]. The input of the system is the target orientation n * . The real-time orientation of the SHM n can be constructed from the real-time tracking on the images, which will be introduced in Section VII.

VI. SYSTEM OVERVIEW

A. Magnetic actuation and stereo vision system

The SHM in this paper is actuated by a 3D Helmholtz coils system. Three Helmholtz coil pairs are orthogonally arranged to generate a uniform rotating field in the 3D space. The 3D world coordinate system R w = (O, x, y, z) is defined in Fig. 6a: x, y, z are respectively the axis of the small, medium, and big coils. The maximum magnetic flux density is designed as 10 mT. Each pair is driven by a Maxon ADS 50/5 4-Q-DC servoamplifier, capable of 5 A continuous current. Analog communication between the PC and the amplifiers is accomplished with a Sensoray 626 Analog and Digital I/O card. The magnetic fields generated by the three coil pairs were measured and calibrated by a gaussmeter Hirst GM08. The stereo vision system consists of two firewire cameras, as shown in Fig. 6b: one is on the top (Pike F032B) of, the other is beside the Helmholtz coils system (Guppy Pro F032). Due to the limited side opening, an endoscope (Bipol) with diameter of 2.7 mm is mounted on the camera instead of a simple lens.

B. Helical swimmer at low Reynolds number

We designed a Scaled-up Helical Swimmer (SHM) with a magnetic head, shown in Fig. 7a, and as presented in [START_REF] Xu | Characterization of three-dimensional steering for helical swimmers[END_REF]. The SHM is painted in black in order to have a good image contrast for tracking. It has an overall length of 17.8 mm with 3.5 turns, and a diameter of 1 mm. The width of the helical tail is 1.4 mm, and the thickness is 0.3 mm. A permanent magnet is planted in the head. The magnet is magnetized in the direction of the diameter, so that the SHM will rotate around its own axis in a rotating magnetic field.

Helical swimmers with a magnetic head show a cut-off frequency [START_REF] Xu | Modeling and swimming property characterizations of scaled-up helical microswimmers[END_REF]: below the cut-off frequency, they rotate in sync with the external rotating magnetic field; beyond this cut-off frequency, their rotation frequency decreases. The cut-off frequency depends on the Reynolds number. In the following experiments, the SHM is put in pure glycerol liquid. , and a human spermatozoon moves with Re ≈ 0.01 [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF]. The cut-off frequency of the used SHM in pure glycerol is measured as 5 Hz. The rotation frequency of the magnetic field is always below the cut-off frequency in the following experiments, so that the SHM is supposed to always rotate synchronously with the magnetic field.

C. Real-time tracking on the image

A robust real-time tracking of the SHM is crucial for visual servoing tasks, which is realized by ViSP [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]. The image of the SHM is considered as a blob (vpDot class in ViSP) defined by a set of connected pixels with the same gray level. The image of the SHM is shown in Fig. 7b, and tracked blob is shown in Fig. 7c and presented in green.

The barycentre and the axis of a helical swimmer can be calculated from image moments [START_REF] Horn | Robot Vision[END_REF]. The barycenter of the SHM is presented in Fig. 7c by a red cross. The axis of a SHM in the image can be defined by the direction of the axis of least inertia. The axis of the SHM is presented in Fig. 7c by a red line. The 3D pose of the SHM can be estimated based on the barycenter and the axis tracked from the images, which will be presented in the next section.

VII. REAL-TIME 3D POSE ESTIMATION

The 3D pose of the SHM is described by its orientation and its position, which are respectively defined as the normal vector of the real-time axis and barycenter of the SHM in the 3D space. The methods to estimate the orientation and the position of the SHM, using the axes and the barycenters tracked from the two images taken with the top and side cameras, will be introduced hereafter. The calibration of the cameras are detailed in [START_REF] Xu | Propulsion Characteristics and Visual Servo Control of Scaled-up Helical Microswimmers[END_REF].

A. Real-time orientation estimation

The axis of the SHM in the 3D space is the intersection of the two planes P t and P s , which are defined by the projection centers of the cameras and the axis of the SHM n, as shown in Fig. 8. It can then be calculated by the vector product of the two normal vectors h t and h s of the two planes [START_REF] Andreff | Visual servoing from lines[END_REF]. The point P is a point on the axis of the SHM, with ct P = ( ct X, ct Y, ct Z) T in the top camera coordinate system. As P is in the plane P t , we can write:

h tx ct X + h ty ct Y + h tz ct Z = 0 (27)
with ct h t = (h tx , h ty , h tz ) T the normal vector h t in the top camera coordinate system. Then,

h tx ct X ct Z + h ty ct Y ct Z + h tz = 0 (28)
The point p t is the image of the point P,

with im p t = ct X ct Z , ct Y ct Z
T in the metric image plane of the top camera by definition, and im its homogeneous coordinates. The relationship [START_REF] Xu | Characterization of three-dimensional steering for helical swimmers[END_REF] defining the equation of the line in the image supporting the helical swimmer can be written under a matrix form:

ct h T t imt pt = 0 ( 29
)
with imt pt obtained by measuring the point p t in the pixel image plane of the top camera. Then, with h t = 1, we can obtain ct h t . The unit normal vector h s of the plane P s in the side camera coordinate system cs h s can be obtained by the same method. Therefore, the axis of the SHM in the top camera coordinate system can be expressed

ct n = ct h t × ct h s = ct h t × ( ct R cs cs h s ) (30)
where ct R cs is the rotation matrix from the side camera coordinate system to the top camera coordinate system. Then, the axis of the SHM in 3D world coordinate system can be written as:

n = w R ct ct n (31)
where w R ct is the rotation matrix from the top camera coordinate system to the world coordinate system. The orientation of the SHM is then obtained.

B. Real-time position estimation

The real-time position of the SHM is estimated by triangulation [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Based on the projective model of the camera, the images of the barycenter of the SHM on the images planes are linearly dependent with it. Therefore, their cross products are zero, with the transformation between the two camera coordinate systems, we can write:

imt gt × ( ct R cs cs G) + imt gt × ct t cs = 0 ims gs × cs G = 0 (32)
where ct R cs and ct t cs are respectively the rotation and translation matrix from the side camera to the top camera coordinate system. The system can be solved with the Singular Value Decomposition (SVD) method to get the point G in the side camera coordinate system cs G. Then, the barycenter G of the SHM in the 3D world coordinate system can be expressed as:

w G = w R cs cs G + w t cs (33)
where w R cs and w t cs are respectively the rotation and translation matrix from the side camera to the world coordinate system. The position of the barycenter of the SHM is then obtained.

Once the real-time orientation and position of the SHM are obtained, the control loop is closed, as shown in Fig. 2. According to our best knowledge, the control law of the magnetic field actuating a helical swimmer to follow a reference path with visual feedback is proposed for the first time.

VIII. RESULTS

A. Control parameter of the 3D steering

First of all, the control parameter of the 3D steering should be determined. The magnetic field is controlled by the output voltage of the PCI card U. As U is proportional to B, U can be written in the same form as [START_REF] Laumond | La robotique mobile. Traité IC[END_REF] with a control parameter of λ U . According to the maximum capacity of the Helmholtz setup, the control parameter λ U is limited to 3. The performances of the 3D steering of the SHM with different control parameters are shown in Fig. 9, with an initial position of 0/40 (direction angle / inclination angle) to a target position of 60/40. The rotation frequency is 2 Hz. Note that an oscillation of the calculated real-time angle of the SHM is observed, which is due to the asymmetry of the helical geometry. The control parameter is chosen as λ U = 2, because of the response time and reasonable overshooting, as shown in Fig. 9a. An impulse of the inclination angle is observed in Fig. 9b, although the target inclination angle is the same as the initial inclination angle. The inclination angle increases more with a higher control parameter, because the magnetic torque exerted on the SHM is higher.

The trajectory of the SHM for the 3D steering from 0/40 to 60/40 is traced on a unit sphere, as shown in Fig. 10. The movements are described in the reference frame relative to the center of the head of SHM. Static errors between the target orientation and the final orientation can be observed on the figure. As these errors are the same for different parameters, which are about -2 • (3%), we presume that these errors are due to the calibration error of the system. Once the control parameter is determined, visual servo control for 3D steering of the SHM will be tested with more command angles.

B. Visual servo control for 3D steering

The 3D steering with zero angle firstly validated. There is no initial error, the SHM stabilizes at the initial orientation. The control of the 3D steering is sensitve to small angle control, which is validated with an initial error of 5 • . The steering of direction angles are then tested every 60 • in a turn, with an inclination angle of 40 • . The errors in direction angle during the steering are shown in Fig. 11. The convergence of the errors states that the steering is achieved. The trajectories are traced on the unit sphere, as shown in Fig. 12. Note that when steering at negative direction angles, such as -60 • , the inclination decreases and then increases during the control, which is contrary to steering at a positive direction angle. This is due to the direction of the self-rotation of the SHM. Hereafter, the 3D steering control is achieved with the same target direction angle and different target inclination angles. The trajectories are traced in Fig. 13. The parameter k t is for the direction angle error, and the parameter k d is for the distance error. For the same k t , the travelled path of the SHM is the shortest with the least k d ; and the SHM reaches the reference path the fastest with the highest k d . Therefore, the choice of the control parameters depends on the requirement of the system: reaching the reference quickly or travelling less.

An image sequence of the SHM during the visual servo control for the path following of a straight line is shown in Fig. 15. The reference path is presented by a blue line. The SHM reached the reference path, then went along this reference path. In fact, the SHM presents drifting in the y axis due to the boundary effect [START_REF] Peyer | Non-ideal swimming of artificial bacterial flagella near a surface[END_REF]. As the control scheme takes into account both the orientation error and the position error, the closed-loop control corrects this distance error with a negative target direction angle. Therefore, a propulsive force is generated in the -y axis to counterbalance this drifting effect. An open-loop control for advancing along the line of y = 0 axis will be compared in the next section.

D. Open-loop vs closed-loop

The trajectory of the SHM for an open-loop control of following a straight line of y = 0 is traced on Fig. 16, as well as a trajectory for a visual servo control of following the same reference path without initial distance error. A perturbation of the force in the direction of y on the SHM is observed during the open-loop control, due to the friction with the substrate. In case of the closed-loop control by visual servoing, the SHM returned onto the reference line despite the perturbation. We define that the SHM is on the path while the distance between the SHM and the path is < 0.5 mm. Then, for the visual servo path following, the SHM is always on the path. As for the open-loop control, the SHM is off the path from about 15 s and never goes back onto the path. Therefore, visual servo control of SHM shows more precise motion than a simple open-loop control. According to our best knowledge, the closed-loop control for the path following of helical swimmers is realized for the first time.

IX. CONCLUSION AND PERSPECTIVES

In this paper, we proposed a control law of the rotating magnetic field for planar path following of helical swimmers. The helical swimmer is millimeter-scaled in a viscous liquid to simulate the low Reynolds number environment. A 3D pose estimation including the position and orientation of helical swimmers is based on the real-time tracking on the images taken by a stereo vision system. The results of the experiments show that the helical swimmer enables closedloop 3D steering by visual servoing for the first time according to our knowledge. The planar path following of a straight line was also achieved. The control parameters define the trajectories of the helical swimmer: reaching the reference path quickly or traveling less. In future works, we will study the 3D path following of helical swimmers. Then, we will scale the helical swimmers down to the microscale, and apply the same control method with suitable optical devices. We expect that helical swimmers will enable closed-loop 3D motion in the future.

Fig. 1 :

 1 Fig. 1: Comparison between path following and trajectory tracking algorithms in presence of strong disturbances on the helical swimmer velocity (e.g. friction on substrate, nonhomogeneous fluid).

Fig. 2 :

 2 Fig. 2: Block-diagram of the visual servo control for planar path following based on 3D steering of helical swimmers.

Fig. 3 :

 3 Fig. 3: Magnetic helical swimmer steering model. The inclination angle and direction angle of steering are respectively θ i and θ d .

Fig. 5 :

 5 Fig. 5: The direction of the magnetic field B n depends on the projection of B ⊥n on the target orientation of the SHM n * : (a) along n * , and (b) along -n * .

Fig. 6 :

 6 Fig. 6: (a) Three orthogonal Helmholtz coil pairs CAD design. (b) Experimental setup: Helmholtz coils, top camera, side camera with endoscope.

  The glycerol density measured by a portable density meter (DMA 35) is 1.26 × 10 3 kg/m 3 . The viscosity of glycerol measured by a falling ball viscometer (Brookfield KF10) is

Fig. 7 :

 7 Fig. 7: (a) CAD design of the SHM with a magnetic head. (b) The fabricated SHM. (c) The SHM is real-time tracked by ViSP [29]. The red cross is the barycentre and the red line is the axis.

Fig. 8 :

 8 Fig. 8: Determination of the axis of the SHM in the world coordinate system.

Fig. 9 :

 9 Fig. 9: Closed-loop responses of the real-time (a) direction angle (b) inclination angle of the SHM to the 3D steering with different control parameters.

Fig. 10 :

 10 Fig. 10: Trajectories of the SHM for 3D steering control with different control parameters.

Fig. 11 :

 11 Fig. 11: Error of the orientation during the 3D steering control.

Fig. 12 :

 12 Fig. 12: Trajectories of the SHM for the 3D steering control with different direction command angles.

Fig. 13 :

 13 Fig. 13: Trajectories of the SHM for the 3D steering control with different inclination command angles.The visual servo control for the 3D steering of the SHM in the 3D space is validated. According to our best knowledge, the closed-loop control of the 3D steering of helical swimmers is realized for the first time. The closed-loop 3D steering control of helical swimmers is crucial to realize a controlled 3D motion.

Fig. 16 :

 16 Fig. 16: Comparison of the trajectories of the open-loop control and visual servo control for following a straight line of y = 0.

Fig. 15 :

 15 Fig. 15: Image sequence of the SHM during the visual servo control for path following of a straight line with the control parameters of k t = 0.1 and k d = 0.8 taken by the top camera, and the side camera (insert at the top-left corner). The image sequence demonstrates that thanks to the control scheme, the barycenter of the SHM can reach the reference path and then follow it.

ACKNOWLEDGMENT

We acknowledge funding from Émergence-UPMC-2012 research program, the Franche-Comté region, and ACTION, the French ANR Labex no. ANR-11-LABX-01-01. The authors are also grateful to Sinan Haliyo, Sylvain Pledel, Christophe Grand, Thomas Seon, Tianyi Li and Jean-ochin Abrahamians for providing technical supports on the expermental setup. Gilgueng Hwang obtained both M.Sc. and Ph.D degrees of Electrical Engineering from University of Tokyo in 2005 and 2008 respectively. Part of this Ph.D work was done by the collaboration with Prof. Bradley J. Nelson's multiscale robotics group in the Institute of Robotics and Intelligent Systems (IRIS), ETHZ, Switzerland. He has worked as a postdoc in ISIR, UPMC, France for 2 years. Then, he joined to LPN-CNRS as CNRS tenured researcher. His research interest covers microfluidics, lab-on-a-chip, in-situ SEM characterization, nanomanipulation. Nicolas Andreff received the Engineer degree ( M.S. degree, 1994) in Computer Sciences and Applied Mathematics from ENSEEIHT, Toulouse, France and the Ph.D. (1999) in Computer Graphics, Computer Vision and Robotics from INPG, Grenoble, France. He is currently Professor at the AS2M department of Institut FEMTO-ST, Université de Franche-Comté, Besanc ¸on, France. He created the Biomedical Micro-/Nano-Robotics group at AS2M department in 2012. His research interests range from vision-based control and parallel kinematics to intracorporeal microrobotics. Stéphane Régnier received the Ph.D. degree in Mechanics and Robotics from the University of Pierre & Marie Curie, Paris, France in 1996. He is currently Professor at the Institute of Intelligent Systems and Robotics (ISIR), University of Pierre & Marie Curie, Paris, France. He has been head of the ISIR micromanipulation team since 2001. His research interests are focused on micro and nano manipulation, teleoperation and haptic feedback at the nanoscale, micromechatronics and biological cell characterization.


Morgan Bertin

François Hild

Stéphane Roux

Florent Mathieu

Hugo Leclerc
email: hugo.leclerc@lmt.ens-cachan.fr

Patrick Aimedieu
email: patrick.aimedieu@enpc.fr

Integrated Digital Image Correlation applied to elasto-plastic identification in a biaxial experiment

Keywords: Digital Image Correlation (DIC), Full field measurements, Identification, Constitutive law

The identification of the parameters of several constitutive laws is performed with the Integrated Digital Image Correlation (IDIC) technique in a biaxial experiment for a cruciform specimen made of stainless steel. The sought material parameters are assessed with the contribution of both reaction forces (from load sensors) and displacement fields (measured via digital image correlation). For each constitutive law a global residual quantifying the model error is assessed.

Introduction

The identification and validation of constitutive models are crucial issues for mechanical design. Challenges such as optimizing structures with respect to mass or using innovative materials, e.g., high performance steels or composite materials push toward complex and multiaxial constitutive models. However, standard tests [START_REF]Metallic materials -Tensile testing -Part 1: Method of test at room temperature[END_REF] related to the identification of material parameters require sample geometries for which the mechanical response must be homogeneous and uniaxial. Such limitations lead to numerous elementary tests even though more complex (i.e., multiaxial) tests individually provide much more data. To characterize the latter ones, inverse identification methods based on full field measurements [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF][START_REF] Grédiac | Full-Field Measurements and Identification in Solid Mechanics[END_REF] are developed.

Since Hill [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] proposed to account for anisotropic plasticity, several identification strategies have been proposed. After briefly reviewing some of them, the paper seeks to solve the problematic of the identification of elastoplastic laws with Integrated Digital Image Correlation (IDIC).

First the simplest identification is performed by tuning the constitutive parameters with uniaxial strain-stress curves and least-squares fit [START_REF]Metallic materials -Tensile testing -Part 1: Method of test at room temperature[END_REF]. However, this strategy does not provide internal checks of the accuracy of the resultant material parameters. Furthermore, it is assumed that the sample geometry behaves under statically determinate stress states [START_REF] Arcan | A method to produce uniform planestress states with applications to fiber-reinforced materials[END_REF]. This hypothesis is in opposition with the observed mechanical behaviors when complex and multiaxial experiments are performed.

Second approximately fifteen years after the appearance of the Finite Element Method, Kavanagh and Clough [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF] proposed the Finite Element Model Updating (FEMU) technique. They focused their work on the characterization of nonlinear elastic materials. FEMU consists of minimizing the sum of squared residuals, a residual being the difference between the numerical and experimental displacement fields (i.e., FEMU-U), strain data (i.e., FEMU-), load levels (i.e., FEMU-F), or combinations of the previous quantities (e.g., FEMU-UF) [START_REF] Grédiac | Full-Field Measurements and Identification in Solid Mechanics[END_REF].

This technique also provided new routes to tackle two main issues, namely, i)

the internal validation of the identified parameters, and ii) the use of statically indeterminate stress conditions.

Some authors proposed strategies aiming to assess elastoplastic properties of aluminium alloys [START_REF] Gelin | An inverse method for determining viscoplastic properties of aluminium alloys[END_REF] and metals [START_REF] Ghouati | Identification of material parameters directly from metal forming processes[END_REF] while combining experimental load measurements and FE simulations. However, full-field measurement techniques were not used and only reaction forces were considered to solve the identification problem.

The advent of DIC as an experimental tool [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF] has allowed elastic properties to be determined [10,[START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF][START_REF] Lecompte | Mixed numerial-experimental technique for orthotropic parameter identification using biaxial tensile tests on cruciform specimens[END_REF]. New strategies have been designed to tackle the increasing number of experimental data. Avril et al. [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF] summarize several techniques based on full-field measurements to identify isotropic linear elastic properties. Some of these identification strategies have been extended to elastoplasticity and nonlinear mechanical behaviors. Lecompte et al. [START_REF] Lecompte | Elastic and elasto-plastic material parameter identification by inverse modeling of static tests using digital image correlation[END_REF] identified Hill's parameters [START_REF] Hill | A theory of the yielding and plastic flow of anisotropic metals[END_REF] with the FEMU technique under biaxial experiments. Conversely, Grédiac and Pierron [START_REF] Grédiac | Applying the virtual fields method to the identification of elasto-plastic constitutive parameters[END_REF] have used the virtual fields method. The constitutive error gap was also extended to plasticity [START_REF] Latourte | Elastoplastic behavior identification for heterogeneous loadings and materials[END_REF]. Haddadi and Belhabib [START_REF] Haddadi | Improving the characterization of a hardening law using digital image correlation over an enhanced heterogeneous tensile test[END_REF] have investigated the characterization of a hardening law on a heterogenous tensile test using the weighted-FEMU technique with both reaction force and strain field. Réthoré et al. [START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3d kinematics[END_REF] considered quasi-3D measurements with stereo digital image correlation to perform the identification of elasto-plastic constitutive pa-rameters using the integrated digital image correlation technique. Schmaltz et al. [START_REF] Schmaltz | Comparison of different biaxial tests for the inverse identification of sheet steel material parameters[END_REF] made use of stereo-DIC measurements and FEMU-U to identify a plastic law for four different biaxial cross-shaped geometries.

All the previous methods use as input the measured displacement fields or derived strain fields. However, these kinematic fields may be seen as parameterized by the sought material properties. Such observation enables the mechanical identification to be integrated within the digital image correlation procedure.

This method is called Integrated-DIC [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties -a review[END_REF][START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF][START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated DIC[END_REF]. The latter is a global DIC technique that relies on a mechanically driven choice for the kinematic basis.

The measurement and identification of the sought parameters is thus performed in a single step. Since the identification is formulated at the pixel level, the process becomes mesh independent, provided the latter allows for a faithful description of the displacement fields as obtained from the mechanical modelling [START_REF] Lindner | On the evaluation of stress triaxiality fields in a notched titanium alloy sample via integrated digital image correlation[END_REF].

The same type of integrated strategy will be followed herein. In the approach followed herein, a commercial FE package will be used in a non-intrusive way to estimate the displacement and load sensitivities to the sought material parameters. Consequently, any built-in or user-defined material model can be probed. This type of implementation makes the present work very generic. Consequently different constitutive postulates will be assessed hereafter (e.g., linear and nonlinear kinematic hardening) to discuss model errors. The analysis will deal with a cruciform specimen in a biaxial experiment under load-controlled mode. The testing machine is a compact (2000 cm 3 ) biaxial device. The selected material is a precipitation-hardened stainless steel (i.e., 17-7 PH grade [START_REF] Ak-Steel | Mechanical characteristics 17-7 PH steel[END_REF]). It is used thanks to its excellent mechanical properties with respect to corrosion and fatigue [START_REF] Leybold-Herbert | Axial-load fatigue tests on 17-7 ph stainless steel under constant-amplitude loading[END_REF][START_REF] Tick-Hon | Cyclic plasticity of 17-7 precipitation-hardenable semiaustenic stainless steel[END_REF]. One key aspect of the present study is related to the small (i.e., sub-millimetric) thickness of the tested sheet since it will be used to make bellows.

The paper is divided into three main sections. In Section 2, the experimental tools and the theoretical background are presented. Integrated Digital Image Correlation (IDIC) and the chosen constitutive laws are also introduced.

Section 3 deals with the experimental procedures. The sample geometry, the prescribed loads and the sources of uncertainties are assessed. Section 4 shows the identification results and a discussion based upon the identification residuals is proposed.

Experimental protocol

The experimental procedure utilizes an electromechanical testing machine, Mini-ASTREE (Figure 1), and one digital camera with a telecentric lens. The machine prescribes arbitrary displacement and forces in two orthogonal directions, e 1 and e 2 , with opposite actuators that can be controlled in a symmetric fashion so that the specimen center is motionless [START_REF] Périé | Analysis of a multiaxial test on a c/c composite by using digital image correlation and a damage model[END_REF]. shown in Table 1. To avoid deflection during machining the sheet is maintained between two 3-mm thick aluminum sheets. The mechanical properties given by the manufacturer are the 0.2% yield stress (σ 0.2% 0 = 300 MPa) and Poisson's ratio (ν = 0.3). The Young's modulus is not mentioned and its value is assumed to be equal to E = 200 GPa, which is the value of the untreated steel. The loading history may influence the identification results. A so-called triangular loading path is prescribed. Since the assessment of elastoplastic laws with kinematic hardening postulates is one of the present objectives, a cyclic loading history is prescribed. This drawing is our property; it can't be reproduced or communicated without our written agreement. where each triangle consists of a first increase of F 1 at F 2 = 0, followed by an equal increase of F 2 at fixed F 1 , and finally and equibiaxial unloading down to the mesh is constructed from the reference image. Consequently, it follows very faithfully the boundaries of the sample. Further, a small element size was selected (i.e., 25 pixels) to capture the strain gradients as best as possible for regular T3-DIC. For the integrated approach, it could have been further refined [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF][START_REF] Lindner | On the evaluation of stress triaxiality fields in a notched titanium alloy sample via integrated digital image correlation[END_REF], but this option was not considered herein since both approaches are to be compared with the same underlying mesh. In the following the geometric coordinates shown in Figure 3(b) will be identical for any shown field.

F 1 = F 2 = 0

Identification strategy

The numerical tools and the identification strategy are summarized hereafter. First, global DIC and Integrated-DIC are formulated. Second, the constitutive laws and their respective material parameters are introduced. Last, the numerical implementation of the identification scheme is briefly recalled.

Digital Image Correlation (DIC)

Global DIC will be used for comparison purposes with Integrated-DIC.

Among various kinematic bases, finite element shape functions have initially been introduced to deal with regular meshes made of four-noded quadrilaterals [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to portevin-le châtelier bands[END_REF]. It was subsequently generalized to deal with unstructured meshes made of three-noded triangle (T3) elements [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]. DIC relies on the registration of an image f in the reference configuration and a series of pictures g in the deformed configurations. The problem consists of minimizing the global correlation resid-

uals χ 2 f = 1/N t Nt t=1 χ 2 f t (t)
, which is the sum of squared differences between the deformed image corrected by the measured displacement u(x, t) and the reference image (written for each time t independently) over the Region Of Interest

χ 2 f t (t) = 1 2γ 2 f N Ω Ω ((g(x + u(x, t), t) -f (x)) 2 (3.1)
with respect to the sought displacement fields u(x, t), where x is any considered pixel. In this expression Ω denotes the Region of Interest (ROI), N Ω its area expressed in terms of the number of pixels it contains, N t the number of time steps, and γ f the standard deviation (expressed in gray levels) of the white noise assumed to affect each image independently (including the reference one, which is responsible for the factor of 1/2 coming as a multiplicative term in this functional).

The displacement field is decomposed onto a basis of fields ψ n (x) that is selected at will

u(x) = n u n ψ n (x) (3.2)
where ψ n are (vector) T3 shape functions in the present case ((i.e., three noded triangles with a linear interpolation of displacements), and u n the unknown degrees of freedom. Ideally, the number of these fields should be kept as small as possible (i.e., to reduce the uncertainty1 on their amplitude) yet large enough to capture the anticipated variety or heterogeneity of the actual experimental field.

The minimization of χ 2 f is achieved by successive linearizations and corrections using Newton-Raphson's scheme

[M ]{δu} = {b} (3.3)
where [M ] is the DIC matrix, {δu} the vector gathering all increments of measured displacement amplitudes, and {b} the residual column vector. In the present case, the minimization is performed independently for each considered time, i.e., on χ 2 f t .

Integrated DIC

The integration of the mechanical identification at the measurement step consists of choosing as the generalized degrees of freedom the sought parameters associated with a chosen constitutive law. Consequently, the kinematic basis is chosen as the set of sensitivity fields gathered in a matrix [S U ] [28]

[S U ] = ∂{u} ∂{p} (3.4)
where {u} is the computed nodal displacement vector, {p} the vector gathering the sought parameters in the computation via forward finite differences in the present case. One of the interests offered by an integrated approach is that the mesh size does not prevent the convergence of computations since the sought material parameters are significantly fewer in comparison with the degrees of freedom associated with the underlying FE discretization [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF][START_REF] Lindner | On the evaluation of stress triaxiality fields in a notched titanium alloy sample via integrated digital image correlation[END_REF]. Consequently, errors coming from a coarse mesh can be avoided.

If the same mesh is used, and the DIC matrix has been computed as [M ],

and the residual vector {b}, the IDIC procedure simply consists of projecting the nodal displacement field onto the sensitivity fields [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated DIC[END_REF]. One main difference with the previous T3-DIC procedure is that only one global minimization is performed in a single spatiotemporal analysis [START_REF] Neggers | Timeresolved integrated digital image correlation[END_REF]. The identification of the sought parameters is achieved by solving iteratively until convergence linear systems

{δp} = 1 2γ 2 f [M IDIC ] -1 [S U ] t {b} (3.5)
where

[M IDIC ] = 1/(2γ 2 f )[S U ] t [M ][S U ]
is the weighted kinematic Hessian. If this unique quantity χ 2 f is minimized, it can be shown that it is equivalent to a weighted FEMU-U procedure provided the noise level is small [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated DIC[END_REF].

The applied load is also of importance for the identification because it provides an additional measured quantity, different in nature from the kinematics, and hence it contributes significantly to the identification and reduces the uncertainty of the determined material parameters. The reaction forces can also be measured and computed. Thus, the same approach is followed with the reaction forces for which χ 2 F is minimized

χ 2 F = 1 N F N t {F m -F c } t [C F ] -1 {F m -F c } (3.6)
where {F m } are the measured reaction forces and {F c } are the computed levels with respect to the chosen material parameter set, [C F ] the covariance matrix of the measured loads (in the present case it is assumed that the load measurements are uncorrelated so that [C F ] = γ 2 F [I], and N F the number of load data.

It is assumed that the load uncertainty is proportional to the magnitude of the load, namely,

γ 2 F = ρ 2 1 |F | 2 .
Moreover, it is easy at this stage to incorporate a minimum measurement uncertainty for the load cells by including an additional noise term whose variance ρ 2 0 is independent of the load level. This practically disqualifies all measurements of forces below F min = ρ 0 /ρ 1 . Thus, in the following it is assumed that

γ 2 F = ρ 2 1 |F | 2 + ρ 2 0 (3.7)
In the proposed procedure, the measured forces will be compared to the computed ones based on a numerical simulation that itself incorporates i) a constitutive law whose parameters will at convergence be adjusted to those of the material to be identified, and ii) Dirichlet (displacement) boundary conditions that are provided by the DIC measurements. The DIC measurement itself has an uncertainty that will induce an additional contribution, which is independent of the load magnitude as the DIC uncertainty is not dependent on the displacement amplitude, affecting computed forces rather than measured ones. However, as one will be interested in their differences, it is equivalent to transfer this additional force uncertainty on the measurement. The result is a term that cannot be distinguished from ρ 0 . In practice, load cells are designed to have a dynamic range adapted to the geometry and mechanical properties of the specimen, and hence, the minimum and maximum values of these sensors are not limiting. In contrast, the identification of elastic properties requires sufficiently small strains to ensure the relevance of linear elasticity. In this case, the contribution of DIC uncertainty to the boundary conditions, and hence to the load level, is generally expected to provide a level for ρ 0 that may be limiting.

The minimization of χ 2 F leads to the variation of the identified set of material parameters

{δp} = 1 γ 2 F [H F ] -1 [S F ] t {F m -F c } (3.8)
where The identification based upon both observables, i.e., displacement field and reaction force, is achieved by minimizing the global functional χ 2

[H F ] = γ -2 F [S F ] t [S F ]
I χ 2 I = N Ω N Ω + N F χ 2 f + N F N Ω + N F χ 2 F (3.9)
where the correlation and reaction force residuals have been introduced in Equations (3.1) and (3.6) respectively. The choice for the specific weight stems from a Bayesian foundation in the weighted quadratic difference including noise covariance. The minimization of the global residual (3.9) requires an iterative computation of the parameter increments

{δp} {δp} = [H IDIC ] -1 1 2γ 2 f [S U ] t {b} + 1 γ 2 F [S F ] t {F m -F c } (3.10)
where the global (i.e., kinematic and static) Hessian [H IDIC ] is the sum of kinematic [M IDIC ] and static [H F ] Hessians. The fact that images and load data are considered enables for the identification of the elastic parameters (i.e., the Young's modulus in particular) contrary to what was performed by Schmaltz et al. [START_REF] Schmaltz | Comparison of different biaxial tests for the inverse identification of sheet steel material parameters[END_REF].

The covariance matrix of the identified parameters reads

[C I p ] = {δp} ⊗ {δp} = [H IDIC ] -1 (3.11)
where • is the mean value of •. Another useful indicator is the correlation matrix (no index summation used)

(Cor p ) ij = (C p ) ij (C p ) ii (C p ) jj (3.12)
The diagonal terms of the correlation matrix are equal to 1 since one parameter is perfectly correlated with itself and off-diagonal terms vary between -1 and 1.

Constitutive laws

The objective of the present work is to assess several parameters for three von Mises flow rule. The total strain rate ˙ is written in terms of elastic and plastic strain rates ˙ = ˙ el + ˙ pl where ˙ el is the elastic strain rate tensor and ˙ pl the plastic strain rate tensor. The yield surface J 2 (σ -X) = σ 0 is defined such that J 2 is von Mises' stress, X the back-stress, and σ 0 the yield stress.

As a first approximation, a linear kinematic hardening model is chosen for the back-stress change [30]

Ẋ = 2 3 C ˙ pl (3.13)
where C is the hardening modulus. Under the assumption of exponential kinematic hardening, the back-stress becomes

Ẋ = 2 3 C ˙ pl -cX ṗ (3.14)
where C and c are material parameters, p the cumulative plastic strain [START_REF] Frederick | A mathematical representation of the multiaxial bauschinger effect[END_REF][START_REF] Lemaitre | Mechanics of Solid Materials[END_REF].

Numerical implementation

The numerical procedure has been implemented in a C ++ framework, which computes the sensitivity fields from the finite element computations performed with the commercial code Abaqus [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated DIC[END_REF]. The IDIC code provides to the latter the required input, namely, the finite element mesh, the current values of the material parameters, the chosen constitutive law and boundary conditions obtained from a DIC measurement. The parameterization is based on a log scale with a ratio between the current value of the parameter and its initial guess [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3d woven composite: Uncertainty analysis and regularization[END_REF]. where {•} ∞ denotes the infinite norm of the vector {•}. Finally, the numerical procedure uses a Levenberg-Marquardt regularization [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF], which is shown to be more robust than a pure Newton based scheme [START_REF] Schnur | An inverse method for determining elastic material properties and a material interface[END_REF][START_REF] Gelin | An inverse method for determining viscoplastic properties of aluminium alloys[END_REF][START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3d woven composite: Uncertainty analysis and regularization[END_REF] when seeking material parameters [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF].

As proposed by Gras et al. [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3d woven composite: Uncertainty analysis and regularization[END_REF], a regularization functional, namely χ R , is associated with the identification functional (χ I ) to enable for the evolution of the sought parameters only if their corresponding sensitivities are higher than a specific bound. This regularization prevents meaningless identification when the influences of material parameters on the observables are weak. χ R is a convex function reaching its minimum equal to zero and reads

χ 2 R = {q}[C R q ] -1 {q} (3.16)
where [C R q ] is the covariance matrix of the material parameters. Because χ I , and χ R are both dimensionless and normalized to one their summation is valid and the regularized functional becomes χ 2 = χ 2 I +χ 2 R . Finally, the regularization of the linear system reads

([H IDIC ] + λ * [I]) {δq} = {b IDIC } + λ * {q} (3.17)
where λ * is the regularization parameter whose choice is performed to cancel out the influence of noise on the change of material parameters [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3d woven composite: Uncertainty analysis and regularization[END_REF].

Most of the computational cost of such integrated procedures is associated with the computation of the various sensitivities, which each time requires an additional interrogation of the FE code. For the same number of iterations IDIC and FEMU procedures lead to very similar computation times [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated DIC[END_REF]. In the present case, an analysis consisting of 365 pictures required 6 hours on a PC with an 8-core Intel Xeon E5 processor for the elastoplastic nonlinear kinematic hardening law (i.e., five unknown parameters).

Identification results

This section presents the identification results and indicators of model errors.

First a resolution 2 analysis is carried out with an initial set of the sought parameters and the initial DIC solution. Second the assessment of the material parameters is performed. Last, a discussion is conducted on the kinematic, static and global residuals.

Resolution analysis

This part aims to estimate a priori the procedure ability to perform an identification. The DIC measurement is used to prescribe in the FEM analysis the experimentally measured displacement (Dirichlet) boundary conditions and get the sensitivity fields. Table 2 gathers the initial values of the sought material parameters. For the computation of Hessians, the parameters are ranked in the same order except for the eigenvalues ranked from smallest to largest. the reference configuration before loading the sample. Once the pictures have been registered, the variance γ 2 f is estimated as the root mean square of gray level differences at convergence (Table 3). The displacement uncertainty on the constrained boundaries is evaluated by measuring the displacement fields with DIC. The covariance matrix of the unknown degrees of freedom, namely [C u ] taken at the four boundaries provides this information. Furthermore, only the component in the displacement direction is kept. Therefore, the standard displacement uncertainty reads

σ 2 BC = 1 N BC ∆t t {V } t [C U ]{V } (4.1)
where {V } is the vector that cancels out the values of the covariance matrix when they are not related to a constrained node in the prescribed displacement direction. N BC is the number of degrees of freedom corresponding to the boundaries. The standard uncertainty is evaluated as σ BC = 1.3×10 -4 pixel. In these experimental conditions with a Young's modulus equal to E = 200 GPa, the related uncertainty provides a standard load uncertainty equal to ρ DIC 0 = 2.4 N.

Finally, the second parameter is evaluated while measuring the load with displacement control at several load levels F ∈ [0; 2000] N. The displacement is kept constant during 300 s. The evaluation of the covariance matrix (see Equation (3.11)) is the key quantity to study. The Hessians are evaluated with the linear elastic law over the first load cycle, which is assumed to be essentially in that regime

[H F ] =     10 1.4 1.4 0.4     × 10 3 [M IDIC ] =     0 0 0 8.2     × 10 5 [H IDIC ] =     0.
Cor([H IDIC ]) =    
1.00 0.02 0.02 1.00

    ; log 10 (λ I ) =        4.0 5.9        (4.3)
The off-diagonal terms denote a weak correlation between the two material parameters.

This first evaluation does not prove that the elastoplastic parameters can be assessed. The same approach is extended to the elastoplastic law with exponential kinematic hardening. The initial set of parameters is shown in Table 2.

[H IDIC ] is evaluated for the first loading cycle

[H IDIC ] =                4
4 × 10 -4 6 × 10 -5 -4 × 10 -4 -3 × 10 -5 3 × 10 -7                × 10 5 (4.4)
and the corresponding eigen values read

log 10 (λ I ) =                              -1.75 3
                             (4.5)
Since a plastic strain occurs a variation of the Young's modulus modifies the level of the plastic behavior (i.e., for a smaller E the elastic domain is greater since more strain is needed to reach the yield stress). Therefore both [H F ] and [M IDIC ] contribute for the identification of Young's modulus since both observables are influenced by the latter. Last the levels of the terms relative to C and c are very low denoting weak influences on the observables. The evolution of the corresponding material parameters will be avoided thanks to the regularized formulation. The corresponding correlation matrix reads

Cor([H IDIC ]) =                1.
               (4.6)
In the early stages of plasticity the differentiation between the elastic and plastic regimes is not an easy task. Therefore when the boundary between elastic and plastic regimes are unclear the value of the yield stress depends on the Young's modulus. As a consequence, the two parameters are anticorrelated (i.e., Cor([H IDIC ])(E, σ 0 ) ≈ -0.99). The resolution analysis over the first cycle shows that some parameters can be identified, especially elastic parameters.

However, a small incursion in the plastic regime is not sufficient to identify c and C independently.

The same analysis is now applied for the evaluation of [H IDIC ] over the entire loading history (Figure 2(a))

[H IDIC ] =                89
               × 10 5 ; log 10 (λ I ) =                              4
                             (4.7)
The levels of the diagonal terms show that the influences of all material parameters have increased (elasto-plastic parameters included). The five eigenvalues have all their levels greater than the uncertainty and therefore can lead to the identification of all the parameters. Last the correlation matrix is evaluated

from [H IDIC ] Cor([H IDIC ]) =                1.
               (4.8)
still showing an important anti-correlation between the parameters C and c.

The correlation between E and σ 0 has also been cut down.

Based on this resolution analysis two issues have been highlighted:

• Before performing the identification, it is known whether the test is discriminating to the sought parameters. In particular, it is shown that the five loading cycles are useful for the identification of the kinematic hardening parameters.

• The assessment of each parameter, its value and its uncertainty are related to the studied loading history.

In the following the results of the identifications over the first cycle and the entire loading history are analyzed for the three constitutive laws.

Parameter identification for the first loading step

The material parameters are assessed for the first cycle starting from the reference set of Table 2 and with the IDIC formulation (see Equation (3.9)).

For comparison purposes, a global DIC analysis is also run. In the latter no hypothesis is made on the underlying material behavior. The residuals χ f are close to the initial DIC residual (χ f = 2.06), themselves very close to the lower limit (i.e., 1) had the residuals only contained noise contributions. Figure 5 shows the change of IDIC residuals (color) and DIC residuals χ f (black) with time. After the first ten pictures captured before starting the experiment the residual levels increase. The level of DIC residuals can be explained by the mesh size equal to 25 pixels, which is a compromise between displacement resolution and spatial resolution [START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF]. The residual induced by the linear elastic law becomes slightly higher than the DIC residual for t ≥ 21 s, and may point out the onset of plasticity. This is confirmed by analyzing the load residuals, which are 15 times higher than the levels expected from noise alone. is also observed. The level of yield stress is found to be in accordance with known values [START_REF] Ak-Steel | Mechanical characteristics 17-7 PH steel[END_REF].

Figure 6 shows the change of the parameters with the iterations. The con-vergence is reached after few iterations, namely, 5, 7 and 7 for laws A, B and C respectively. As observed in the resolution analysis, the kinematic hardening parameters keep their initial values thanks to the Levenberg-Marquardt regularization. The corresponding standard uncertainty is significantly higher for the hardening parameters than those for the elastic parameters and the yield stress (Table 4) This observation validates the results shown in Table 4. The same analysis is now carried out over the entire loading history (Table 5).

Even if the IDIC residual of law C (χ = 6.3) is 2.5 times higher than the DIC residual (χ f = 2.5) it is the best of the three chosen postulates. The fact that the load residuals are still high is an indication that there remains a model error.

This is particularly true for the elastic law that induces very high residual errors.

The Young's modulus is found to be very small, and the Poisson's ratio tends to 0.5 (i.e., incompressibility as expected from dominant isochoric plasticity). All these trends clearly disqualify elasticity as a model able to describe the reported experiment. This is confirmed by Figure 8 that shows the change of IDIC and DIC residuals over time. Unlike the previous results (Figure 5), gaps between IDIC and DIC residuals are much more pronounced, especially when the level of plastic strain increases. Furthermore, the gain offered by each law with respect to raw DIC results is also clarified. (see Figure 9). Unlike Figure 6 all the parameters are modified and no regularization is needed in the present case. It is worth noting that the more freedom the identification code has (i.e., the number of unknowns is increased), the fewer iterations are needed. Further, the converged solutions are rather far from the initial guesses, which shows the robustness of the procedure.

Comparison with FEMU-F results

Table 6 shows the parameters assessed only by considering the reaction forces (Equation (3.6)). It corresponds to a FEMU-F procedure [START_REF] Mathieu | Estimation of elastoplastic parameters via weighted femu and integrated DIC[END_REF]. To compare both approaches, the gray level and global residuals are also computed. As expected, the static residuals are lowered for all laws since the residuals only depend on reaction forces. Conversely, the overall quality degrades in addition to that associated with gray level residuals. More importantly, the standard uncertainties of the identified parameters are significantly higher, which is caused by the limited number of data used per analyzed time.

Figure 11 shows the displacement field difference for component u 1 for the three constitutive laws with IDIC or FEMU-F approaches compared with DIC measurements at the end of the experiment. The choice of an elastoplastic A weighting that originates from a Bayesian foundation has been followed. Consequently, each pixel and load cell play the same role, when they respective noise level is accounted for.

An initial sensitivity analysis enables the impact of the loading history and the acquisition noise on the identification to be understood. It shows that the first loading cycle is not sufficient to identify plastic laws but it is sufficient to detect its onset with the initial set of the material parameters. The corresponding identification leads to an IDIC residual with similar level as raw DIC.

The identification is also performed over the entire loading history accounting for 5 loading/unloading cycles. Any of the chosen model does not match the available data with the same accuracy as observed for the first cycle. The three laws lead to three different residuals and the exponential kinematic hardening law provides the best results.

In the present study four different error indicators have been considered:

• When DIC and IDIC analyses are performed, the gap with respect to gray level conservation is computed. It estimates the registration quality when the pictures in the deformed configurations are corrected by the measured displacement field and subsequently compared with the picture of the reference configuration.

• DIC and IDIC can also be compared by computing the displacement differences measured by both approaches. When the mesh is identical, as in the present case, the comparison is straight forward.

• Load residuals are yet another way of characterizing the identification quality. In the present analyses, FEMU-F and IDIC could be compared since both included load residuals in their respective formulation.

• Last IDIC as developed herein has an overall quality indicator χ 2 IDIC that compares all sensor information (be it gray levels for each pixel or load level) in the same footing when normalized by the variance of corresponding noise.

This wealth of indicators allows the user to assess not only globally but also with each of the estimators the identification quality and the underlying model error (i.e., choice of constitutive law, finite element discretization).

The present study aims to pave the way to understand and design experimental protocols to not only identify but to (in)validate constitutive laws. In forthcoming works, the IDIC framework will be used to identify other laws for other materials such as polymers and composites. Other scales of observation may also be investigated by using other imaging systems such as confocal microscopy or scanning electron microscopy.

F 1

 1 will denote the load amplitude applied by the two coupled actuators along direction e 1 , and F 2 for the other two coupled actuators. Biaxial compressive and tensile tests can be run with a load range F 1,2 ∈ [-2000; 2000] N. Furthermore the stroke rate varies from 30 µm/min to 5 mm/min. 16-bit gray scale images are captured with a pco.edge camera. The effective magnification is 13.5 µm per pixel or 74 pixel/mm.

Figure 1 :

 1 Figure 1: Biaxial testing machine. A cross-shaped sample (center of picture)

Figure 2 (

 2 b) shows the F 2 versus F 1 cycles when the load amplitude is increased by 200 N between each of them. Furthermore, the load speed is constant and equal to Ḟ1 = Ḟ2 = 3000 N/min when different from 0.

Figure 2 :

 2 Figure 2: (a) specimen geometry (expressed in mm). (b) Prescribed loading

Figure 3 (

 3 Figure 3(a) shows the Region Of Interest (ROI) that contains the specimen

Figure 3 :

 3 Figure 3: (a) Reference image f with the Region Of Interest (ROI) in dark and

  is the static Hessian, and [S F ] = ∂{F c }/∂{p} the reaction force sensitivities defined in the same way as the displacement field sensitivities (see Equation (3.4)). If χ 2F is minimized alone, it corresponds to a load-based FEMU procedure, which is referred to as FEMU-F.

  constitutive laws (A) linear and isotropic elasticity, and (B, C) two elastoplastic laws. The latter ones correspond first to linear kinematic hardening with von Mises flow rule (B), and second to an exponential kinematic hardening (C) with

  It is defined such that the sought parameters in the new setting {q} are expressed from the initial basis as {q} = log{p/p 0 }, where {p} is the vector gathering the values of the parameters and {p 0 } their initial values. The procedure is considered to have converged when the change in the parameters reaches a chosen threshold {δq} ∞ ≤ 10 -4(3.15)

Figure 4 nt ifie d mode l γ 2 l = ρ 2 1 F 2 + ρ 2 0Figure 4 := 0 . 1 N

 42122401 Figure 4: Measurement of the load variance under static conditions as a function

Figure 5 :

 5 Figure 5: DIC and IDIC residuals for the three constitutive laws

Figure 6 :

 6 Figure 6: Change of material parameters with the iteration number for the three

Figure 7 (

 7 b) shows that law A induces very high differences in particular close to one of the connecting radii. This is an indication of plastic activity in the vicinity of stress concentration areas. When compared with Figures7(c-d) there are significant differences, again validating the identification with elastoplastic postulates that lead to equivalent differences with raw DIC. The fact that there remain displacement residuals shows that the plastic behavior has been only partially captured.

Figure 7 :

 7 Figure 7: Displacement field map u 1 measured by DIC and displacement differences between DIC and IDIC for the three investigated laws

Figure 8 :

 8 Figure 8: DIC and IDIC residual history for the three constitutive laws. Please

Figure 9 :

 9 Figure 9: Change of the material parameters for the three investigated laws

11 Figure 10

 1110 Figure 10: (a-b) Measured displacement fields with DIC at the last experimental

Figure 11 :

 11 Figure 11: Displacement differences u 1 between DIC and IDIC (first row) or FEMU-F (second raw) for the three investigated laws. The color bars correspond

29 9

 29 objective function (DIC) χ F finite element method updating objective function (FEMU-F ) χ I Integrated Digital Image Correlation objective function (IDIC) χ R Regularization objective function {p} vector gathering constitutive parameters in the initial basis {q} vector gathering constitutive parameters in the log basis x 2D or 3D coordinates in normal space e α unit vector along direction α = 1, 2 ψ(x) matrix gathering the shape functions of element e f (x), g(x) pictures in the reference and deformed configurations, respectively γ f standard deviation of gray levels γ F standard uncertainty of the load measurement Ω region of interest u(x) displacement field vector [M ] global DIC matrix DIC and IDIC residual history for the three constitutive laws. Please note that a semi-logarithmic scale is used . . . . . . . . . Change of the material parameters for the three investigated laws during the identification iterations. Each parameter is normalized by its initial value . . . . . . . . . . . . . . . . . . . . . . . . . . 30 10 (a-b) Measured displacement fields with DIC at the last experimental time step. (c) Corresponding strain field (normal component along direction 1) . . . . . . . . . . . . . . . . . . . . . . . . 31 11 Displacement differences u 1 between DIC and IDIC (first row) or FEMU-F (second raw) for the three investigated laws. The color bars correspond to the largest and lowest values of the corresponding residual map . . . . . . . . . . . . . . . . . . . . . . . 33

Table 1 :

 1 Chemical composition of 17-7 PH grade[START_REF] Ak-Steel | Mechanical characteristics 17-7 PH steel[END_REF]

	Composition	C	Mg	P	S	Si	Cr	Ni	Al
	(wt %)	0.09 1.00 0.040 0.030 1.00 16.00 -18.00 6.50 -7.75 0.75 -1.50

Table 2 :

 2 Initial value of the material parameters

	Parameter	E	ν	σ 0	C	c
	Value	200 GPa 0.3 300 MPa 10 GPa 10
	Load uncertainties are assessed for each iterations as discussed in Section 3.1
	and the displacement uncertainty is evaluated since 10 images are acquired in
	2 'Resolution: smallest change in a quantity being measured that causes a perceptible
	change in the corresponding indication. [...] The resolution can depend on, for example, noise

(internal or external) or friction.'

[START_REF]International Vocabulary of Metrology -Basic and General Concepts and Associated Terms, VIM, International Organization for Standardization[END_REF]

Table 3 :

 3 Standard gray level and load uncertainties

	Quantity	ρ 0	ρ 1	γ f
	Value	2.8 N 4×10 -4 302 gray levels

Table 4 :

 4 Table 4 gathers the corresponding values and the residuals for each law at convergence. The standard uncertainties γ p are obtained from the covariance matrices [H IDIC ]. Identified parameters and identification residuals for the three laws for the first cycle (γ p is expressed in % of p).

	law	χ	χ f	χ F	E	γ E	ν	γν	σ 0	γσ 0	C	γ C	c	γc
					GPa				MPa		GPa
	A	2.10 2.10	15	174	0.4	0.31	0.09	-	-	-	-	--
	B	2.09 2.09	8.9	195	0.21 0.297	0.1	306	0.25	10	-	--
	C	2.09 2.09	8.9	195	0.25 0.297	0.1	306	0.3	10	-	10 -

Table 5 :

 5 Identified parameters and identification residuals for the three laws for

	the whole loading history
	law	χ	χ f	χ F	E	γ E	ν	γν	σ 0	γσ 0	C	γ C	c	γc
				GPa				MPa		GPa
	A	15	14.9 1820	8.8	0.25 0.499 0.002	-	-	-	-	-	-
	B	7.1	7.1	113	157	0.15	0.31	0.04	480	0.15	6.7	0.15	-	-
	C	6.3	6.3	100	148	0.16	0.30	0.04	423	0.16	8.6	0.16 10.8 0.03

Table 6 :

 6 Identified parameters and identification residuals via FEMU-F for the three laws over the entire loading history cycle (no regularization prescribed)

	law	χ	χ f	χ F	E	γ E	ν	γν	σ 0	γσ 0	C	γ C	c	γc
					GPa				MPa		GPa
	A	19.1	19	1732	10	0.4 0.19	4	-	-	-	-	--
	B	8.1	8.0	11.9	188	0.7 0.23	4	330	0.4	3.7	1.2 --
	C	7.5	7.5	9.6	191	0.7 0.25	4	300	0.5	8.5	3.7 29 11

law instead of an elastic law decreases the displacement difference, thereby decreasing the model error. The residual level is lower when the identification is performed with IDIC rather than FEMU-F. An overall residual remains because i) the model does not correctly predict the plastic strains, ii) localization phenomena occur.

'The uncertainty [...] is a parameter, associated with the result of a measurement, that characterises the dispersion of the values [...][START_REF]International Vocabulary of Metrology -Basic and General Concepts and Associated Terms, VIM, International Organization for Standardization[END_REF].' In the present study, because the sole considered origin of uncertainty is that resulting from a Gaussian white noise, the uncertainty is quantified by the standard deviation of the concerned quantity.

[10] G. Geymonat, F. Hild, S. Pagano, Identification of elastic parameters

Acknowledgments

It is a pleasure to acknowledge the support of Région Ile de France (DICCIT and THERMOFLUIDE-RT projects).