B Brogliato

P Orhant

Contact Stability Analysis of a One Degree-of-Freedom Robot

Keywords: free and constrained motion, compliant and rigid environment, quadratic Lyapunov stability, contact stability, impact Poincaré map

The aim of this note is to examine the conditions of stability of a simple robotic task: we consider a one degree-of-freedom (dof) robot that collides with a spring-like environment with stiffness k, the goal being to stabilize the system in contact with the environment. We study conditions on the feedback gains that guarantee quadratic Lyapunov stability of the task with a well-conditioned solution to the Lyapunov equation. It is shown that when the environment's stiffness k grows unbounded, those conditions yield unbounded values of the gains. Motivated by the stability analysis of the impact Poincaré map in the perfectly rigid case (k = +∞), we propose an analysis that is independent of k. It enables us to conclude on global asymptotic convergence of the system's state towards the equilibrium point. This work can also be seen as the study of stability of a contact (force control) phase, taking into account the unilateral feature of the constraint.

Introduction

The last fifteen years have witnessed important progresses in the theory of manipulator's control. It has been assumed that the robot evolves either in a free space (motion control), or that it remains in contact with a certain environment (hybrid force/position control). Both cases have been considered separately in the literature, and several solutions have been proposed in each case. However a real robotic task often implies phases of transition between free motion and constrained motion, and the so-called transition phase appears to be in most cases crucial for the system's stability. This is the case for hopping robots, walking machines, and manipulation with a robotic hand. Some interesting results for the transition phase control have been presented, see e.g. [START_REF] Volpe | A theoretical and experimental investigation of impact control for manipulators[END_REF] [25] [START_REF] Tornambé | Global regulation of planar robot arm striking a surface[END_REF] [START_REF] Tornambé | Modeling and controling two-degree-of freedom impacts[END_REF]. Also many studies have been devoted to investigate the so-called contact stability problem, due to the unilaterality of the robot's tip constraint [START_REF] Hogan | On the stability of manipulators performing contact tasks[END_REF] [24] [START_REF] Colgate | Robust control of dynamically interacting systems[END_REF] [1] [START_REF] Kazerooni | Contact instability of the direct drive robot when constrained by a rigid environment[END_REF]. It seems that Whitney [START_REF] Whitney | Force feedback control of manipulator fine motions[END_REF] was the first to point out and analyze such problems. However the first papers dealing with the transition phase containing a complete stability analysis for a n-degree of freedom robot controlled by a switching algorithm were recently proposed in [START_REF] Mills | Control of robotic manipulators during general task execution: A discontinuous control approach[END_REF] [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF].

The goal of this study is not to extend these works to more complex environments models or control algorithms, but rather to point out some problems related with stabilization of motion-controlled manipulators that come in contact with a compliant environment, in particular the sufficient conditions guaranteeing asymptotic convergence of the solutions towards the steady-state solution. Indeed we restrict ourselves to a simple continuous PD motion controller (in contrast with the sophisticated switching controllers studied in [START_REF] Brogliato | On the control of finite-dimensional mechanical systems with unilateral constraints[END_REF] [10] [START_REF] Mills | Control of robotic manipulators during general task execution: A discontinuous control approach[END_REF] [START_REF] Mills | Stability and control of robotic manipulators during contact/noncontact task transition[END_REF]) and to the case of a purely elastic environment. Basically our motivation is the following: a question a designer may ask himself when facing a real problem is: should the environment be considered as flexible or as rigid? In general one considers that rigid body collisions occur when the bodies show "sufficiently small" deformation so that they are geometrically rigid at a global observation scale [START_REF] Moreau | Some numerical methods in multibody dynamics: Application to granular materials[END_REF]. The answer is crucial for the choice of the control algorithm 1 and depending on it, the analysis of the whole robotic task may be quite different. Indeed rigid body dynamics involve models which completely differ from those used when compliant bodies are considered (i.e. respectively measure differential equations [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] [START_REF] Moreau | Some numerical methods in multibody dynamics: Application to granular materials[END_REF] instead of ordinary differential equations). In the rigid body impacting case, it is customary and convenient to study the so-called impact Poincaré map which is a discrete-time reduced-order system [START_REF] Wang | Dynamic modeling and stability analysis of mechanical systems with time-varying topologies[END_REF]. In the compliant body case (in general springdashpot contact models), one usually directly uses a Lyapunov-like stability analysis, the system being in this case a very simple example of a hybrid dynamical system [START_REF] Branicky | Stability of switched and hybrid systems[END_REF]. Note that the boundary between "flexible" and "rigid" is quite clear from a mathematical point of view, but not from a practical one: Besides clearly rigid environments made of hard materials (concrete, iron . . . ) and clearly flexible ones, some others might be considered to belong to one class or the other one depending on the task (masses of the bodies that collide, accuracy of the measurements, limits of the actuators . . . ). We believe the results in this paper may help in partially answering this question.

The note is organized as follows: in section 2, the system and the notations are introduced and we show with a particular Lyapunov function that the closed-loop system's fixed point is Lyapunov globally asymptotically stable (GAS) for a suitable choice of the feedback gains, and for k < +∞. In section 3 we analyze the problem related with quadratic stability of the task when the environment's stiffness grows unbounded. We show that in order to guarantee that the solution P to the system's Lyapunov equation satisfies λ min P ≥ δ > 0 for some δ, where λ min denotes the smallest eigenvalue, the feedback gains grow unbounded as k → +∞. Section 4 is devoted to analytically prove global asymptotic convergence of the position and velocity tracking errors towards zero, independently of the contact stiffness value. The analysis is shown to reduce to the stability analysis of the impact Poincaré map associated to the closed-loop system when k = +∞. Finally conclusions are drawn in section 5, and some technical results are developed and recalled in the appendices.

A Simple Example

The system consists of a simple mass moving horizontally without friction whose position is given by x(t), and a compliant environment at x = 0 whose model is a spring with stiffness k > 0 (see figure 1). The control law is given by u

= -λ 2 ẋ -λ 1 (x -x d ), x d ≥ 0, λ 1 > 0, λ 2 > 0.
We assume that contact is established at t = 0, with x = 0. Then the equations that govern our system are: Notice that convergence of the state (x, ẋ) towards the fixed point of the second equation in (1) may be investigated by considering the associated equivalent mechanical system to the closed-loop system in [START_REF] An | Dynamic stability issues in force control of manipulators[END_REF]. Notice however that the equivalent total mechanical energy for such a nonlinear discontinuous system is not simply the concatenation of both the (closedloop) energy functions for the free and contact phases. Indeed it corresponds to the so-called available storage function in dissipative systems theory. Calculations can be found in [START_REF] Lozano | Stability of interconnected passive systems[END_REF].

m ẍ + λ 2 ẋ + λ 1 x = λ 1 x d if x < 0 m ẍ + λ 2 ẋ + (λ 1 + k)x = λ 1 x d if x ≥ 0 (1)
In the sequel we focus on a particular stability property of that equilibrium point. The motivation for studying this type of stability is evident if one thinks of more complicated tasks as considered for instance in [START_REF] Mills | Control of robotic manipulators during general task execution: A discontinuous control approach[END_REF]. Also the equivalence with a mechanical system may no longer be possible in certain cases, e.g. when the feedback loop contains time-delays or for tracking control.

When the contact stiffness is finite, one can treat such a stabilization problem from different point of views, for instance: i) Study conditions that guarantee that after the first contact has occured, there is no rebound [START_REF] Marth | Stable transition phase for robot arm control[END_REF], ii) Relax the bounceless conditions by studying conditions that insure Lyapunov quadratic stability of the system, i.e. find a Lyapunov function V (x) such that along trajectories of the system V (x) = -x T Qx with Q > 0 (which does not a priori guarantee that the robot's tip will never take off the environment's surface), see e.g. [START_REF] Mills | Control of robotic manipulators during general task execution: A discontinuous control approach[END_REF]. Since these tools will generally provide sufficient conditions only, it is worth investigating whether these conditions are of any practical interest or not. In particular, if they yield lowerbounds on the feedback gains that are proportional to the environment's stiffness, it is clear that as soon as this stiffness becomes too large, the conditions become useless. It is then natural to seek a convergence proof that is independent of the stiffness as well as of the feedback gains values. For the sake of briefness of the paper, we shall study only approach ii). In fact it can be easily (and logically) concluded that bounceless conditions are impossible to obtain with finite force control, for nonzero contact velosity, as the stiffness k grows unbounded. It is note worthy that this work can also be seen as a study on conditions of stability of a force control scheme, taking into account the fact that the constraints are unilateral, i.e. the robot's tip may take off the surface and possibly start a sequence of rebounds.

Before going on with the stability analysis of system (1) when k < +∞, let us recall that as k → +∞, the solutions of (1) converge towards the solutions of the following dynamical system:

⎧ ⎪ ⎨ ⎪ ⎩ m ẍ + λ 2 ẋ + λ 1 x = λ 1 x d if x < 0 ẋ(t + k ) = -ẋ(t - k ) if x(t k ) = 0 ẍ = min(0, -λ 2 ẋ -λ 1 x + λ 1 x d ) if ẋ(t + k ) = 0 (2)
The t k 's generically denote the impact times. The proof of convergence can be found in [START_REF] Paoli | Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales: cas avec perte d'énergie[END_REF], together with a rigorous definition of the used notion of convergence. It is therefore legitimate to seek a stability analysis that encompasses both systems in (1) and (2), i.e. that works for all k ∈ [0, +∞]. We shall come back on the definition and on the stability analysis of the impact Poincaré map associated to (2) in section 4.

For the moment we shall analyze the stability of the task using a single Lyapunov function. To begin with, we show how the stability analysis of the closed-loop system in (1) can be led with a particular Lyapunov function candidate: Let us consider

V = 1 2 m ẋ2 + 1 2 λ x2 + c x ẋ (3)
with λ = λ 1 + k + λ 2 c m , c > 0 is such that c 2 -λ 2 c -m(λ 1 + k) < 0 (since = λ 2
2 + 4m(λ 1 + k) > 0, and √ -λ 2 > 0, such a c can always be chosen arbitrarily small), and x = x -λ 1 x d λ 1 +k . λ and c guarantee that V is positive definite. Now we get along closed-loop trajectories:

• x < 0 (non-contact) V ≤ -λ 2 + c + 1 2 k 2 + 1 2 ẋ2 + - cλ 1 m + 1 x2 + 1 2 λ 1 k λ 1 + k x d 2 + 1 2 cλ 1 k (λ 1 + k)m x d 2 (4) or in compact form V = -a nc ẋ2 -b nc x2 + R (5) with R > 0. • x > 0 (contact) V = (-λ 2 + c) ẋ2 - λ 1 + k m c x2 = -a c ẋ2 -b c x2 (6)
Claim 1. For any stiffness 0 < k < +∞ there exist

P = P T > 0, Q = Q T > 0, λ 1 < +∞, λ 2 < +∞ such that λ 1 > λ 1 , λ 2 > λ 2 implies that for all t ≥ 0, V = z T Pz, V ≤ -z T Qz.
Thus the equilibrium point z = 0 is globally asymptotically stable (GAS).

The proof is given in appendix A.

Analysis of Quadratic Stability Conditions for Large Stiffness Values

We shall be content with these existence results on the feedback gains in the above analysis. However, let us note that if one takes the sufficient conditions for stability deduced from the above analysis, then the feedback gains λ 1 and λ 2 → +∞ as k → +∞. This suggests that in order to obtain quadratic Lyapunov stability of (1) one has to choose feedback gains proportional to the stiffness k as k becomes large. Let us rewrite (1) in state space form as

z ∈ (N C) = {x : x < 0} ż = A c z + 0 k m x z ∈ (C) = {x : x ≥ 0} ż = A c z (7)
where

z T = x - λ 1 x d λ 1 + k , ẋ (8)
A c = 0 1 -1 m (λ 1 + k) -λ 2 m (9)
Clearly the choice of the first component of z stems from the fact that we want to stabilize the robot in contact with the environment. Moreover from (1) one sees that the equilibrium point of the first equation belongs to (C), which means that the system in [START_REF] Hogan | On the stability of manipulators performing contact tasks[END_REF] possesses in fact only one equilibrium point, i.e. z T = (0, 0) (Note that the uniqueness holds for any value of x d ; when x d = 0 both equations in (1) have the same equilibrium point (x, ẋ) = (0, 0)). Stability of A c is independent of k since its eigenvalues are either real strictly negative or with real part equal to -λ 2 2m . Thus for any Q c = Q T c > 0 there always exists P = P T > 0 such that A T c P + P A c = -Q c . Since we want to stabilize the equilibrium point z = 0, we choose a Lyapunov function candidate as V = z T Pz. Along trajectories in (N C) we get

V = -z T Q c z + z T P 0 2k m x ( 10
)
For simplicity of the analysis, let us choose x d = 0. Then we can write

V = -z T Q c z + z T P K z = -z T Qc z (11)
with

K = 0 0 2k m 0 (12)
Simple calculations yield:

Q c = 2 λ 1 +k m p 12 λ 2 m p 12 + λ 1 +k m p 22 -p 11 λ 2 m p 12 + λ 1 +k m p 22 -p 11 2 λ 2 m p 22 -p 12 ( 13
)
Q nc = 2λ 1 m p 12 λ 2 m p 12 + λ 1 m p 22 -p 11 λ 2 m p 12 + λ 1 m p 22 -p 11 2 λ 2 m p 22 -p 12 ( 14
)
where Q nc is the symmetric part of the matrix Qc , that is independent of k. It is worth noting that only the skew-symmetric part of Qc depends on k. Thus a necessary and sufficient condition for Q c to be positive definite is that:

• λ 1 +k m p 12 > 0 • det(Q c ) = 4 λ 1 +k m p 12 λ 2 m p 22 -p 12 -λ 2 m p 12 + λ 1 +k m p 22 -p 11 2 > 0
For Q nc the conditions are the following :

• 2λ 1 m p 12 > 0 • det(Q nc ) = 4 λ 1 m p 12 λ 2 m p 22 -p 12 -λ 2 m p 12 + λ 1 m p 22 -p 11 2 > 0
Our aim in this section is to examine the conditions such that this simple task is Lyapunov quadratically stable, and in particular to find out which kind of conditions this implies on the feedback gains. It is clear that if one concludes for instance that λ 1 must be larger than k, then it means that this sort of stability analysis is completely meaningless and useless as soon as the environment's stiffness is too large; then one has to change the objectives (relax the stability conditions) or the model (consider that the environment's surface is rigid) to get satisfactory conditions on the feedback gains. As shown in appendix B, the following result is true:

Claim 2.
Consider the one-dof closed-loop equations in (1) with x d = 0. Then quadratic stability of the system implies conditions such that when the environment's stiffness k grows unbounded, then the feedback gains λ 1 and/or λ 2 have to be chosen of order ≥ k β , β ≥ 1 2 to guarantee that the solution P of the Lyapunov equation remains bounded away from singularities (i.e. λ min P ≥ δ > 0 for some δ) and that the matrices Q c and Q nc remain positive definite.

We reiterate that the only thing we have done is to study conditions on the feedback gains such that the Lyapunov equation A T c P + P A c = -Q c possesses a solution that is bounded-away from singularities and guarantees Q c > 0, Q nc > 0. The choice for such a stability analysis is quite natural: indeed it is the application of Lyapunov's direct method to a simple hybrid dynamical system [START_REF] Branicky | Stability of switched and hybrid systems[END_REF]. The result of claim 2 are consistent with those to be found in other studies, see for instance [START_REF] Mills | Control of robotic manipulators during general task execution: A discontinuous control approach[END_REF] and [START_REF] Tornambé | Modeling and controling two-degree-of freedom impacts[END_REF] [START_REF] Tornambé | Global regulation of planar robot arm striking a surface[END_REF], although the system we analyze is much simpler that those treated in these references.

A Stiffness Independent Convergence Analysis

Firstly let us consider the system in [START_REF] Branicky | Stability of switched and hybrid systems[END_REF]. Let us take the Poincaré section + = {(x, ẋ): x = 0, ẋ(t + k )}. Notice that if x(0) > 0, then the sequence of impact times {t k } is infinite (this can be easily shown by studying the vector field between the impacts, which forces the system to attain in finite time the constraint surface x = 0 whatever bounded initial conditions one may choose). The impact Poincaré map P : ẋ(t + k ) → ẋ(t + k+1 ) is thus well defined. However it is not explicitely calculable, despite of the simplicity of the dynamics. This is due to the nonzero dissipation during flight-times. Let us choose:

V (k) = 1 2 m ẋ2 (t + k ) (15)
We prove that P is Lyapunov stable with V as a Lyapunov function as follows. Consider the function

V = 1 2 m ẋ2 + 1 2 λ 1 (x -x d ) 2 (16)
Along free-motion trajectories of (2) one obtains

V = -λ 2 ẋ2 (17)
and at the impact times

σ V (t k ) = 1 2 ẋ2 (t + k ) -ẋ2 (t - k ) = 0 ( 18
)
where σ f (t k ) generically denotes the jump in the function f (•) at t k . Hence

V (t + k+1 ) -V (t + k ) = -λ 2 (t k ,t k+1 ) ẋ2 (τ )dτ ≤ 0 (19)
Now from the fact that

V (t + k+1 ) -V (t + k ) = V (k + 1) -V (k)
≤ 0, we conclude the proof. This stability result suggests 2 that one should be able to analyze the stability of the system in (1) for any k ≥ 0, without the drawbacks encountered in the previous section.

In the following, we propose a different convergence analysis to prove that the equilibrium point of ( 1) is asymptotically reached for any initial condition and any value of the feedback gains, independently of the value of k; the particular feature of the analysis is that it extends naturally to the rigid environment case (i.e. k = +∞), contrarily to the foregoing one. Roughly speaking, we consider a particular section of the phase-plane, = {(x, ẋ) : x = 0}. Then we analyze the mass velocity at the instants t i when the trajectories cross this section; we use the fact that these times define a sequence along which the kinetic energy is non-increasing; it follows that if {t i } is an infinite sequence, the velocity must converge to zero when i → +∞; if x d is strictly positive, this leads to a contradiction and there is a finite number of bounces, so that both x and ẋ converge to zero. To clarify the notations the instants t i and t k are depicted in figure 2.

We assume that the mass makes contact with the environment at t = t i , looses contact at t = t i+1 , i ∈ N, and that contact occurs at x = 0. Thus contact occurs on intervals [t 2i , t 2i+1 ], and free motion on intervals [t 2i+1 , t 2i+2 ]. Let us consider the positive definite functions and V nc as in [START_REF] Schmaedeke | Optimal control theory for nonlinear vector differential equations containing measures[END_REF]. On intervals [t 2i , t 2i+1 ],

V c = 1 2 m ẋ2 + 1 2 (λ 1 + k) x - λ 1 x d λ 1 + k 2 (20)
Vc = -λ 2 ẋ2 (21)
On intervals [t 2i+1 , t 2i+2 ],

V nc = -λ 2 ẋ2 (22)
Let T (t) denote the system's kinetic energy. From the fact that

V c (t 2i+1 ) -V c (t 2i ) = T (t 2i+1 ) -T (t 2i ) ( 23
)
and

V nc (t 2i+2 ) -V nc (t 2i+1 ) = T (t 2i+2 ) -T (t 2i+1 ) ( 24
)
we deduce that for all i, T (t i+1 ) -T (t i ) < 0, hence

| ẋ(t i+1 )| < | ẋ(t i )| (25)
The same inequalities hold for V c and V nc . Now notice that there are two situations: either the sequence of instants t i is finite (the bounces stop after a finite time t 2N , N < +∞, and since x d > 0, x(t) > 0 for all t > t 2N ), or this sequence is infinite i.e. N = +∞.

• If N < +∞, then for t > t 2N the system is governed by the second equation in (1) (indeed each time the mass is "outside" the environment it necessarily collides again after a finite time) and we conclude that x → λ 1 x d λ 1 +k , ẋ → 0 asymptotically, globally and uniformly. • Assume that N = +∞. Since the kinetic energy is a positive definite function of the velocity that is non-increasing at times t i , T (t i ) converges as i → +∞, and so does ẋ(t i ). Suppose that | ẋ(t i )| → | ẋss | with | ẋss | ≥ δ > 0 (the ss subscript is for steady state value). Now δ > 0 and since sgn( ẋ(t i )) = -sgn( ẋ(t i+1 )) and x(t i ) = x(t i+1 ) = 0, one deduces that the length of the orbit between t i and t i+1 is strictly positive. Since the flow of both equations in (1) is exponential and bounded for bounded feedback gains and stiffness k, clearly μ[t i , t i+1 ] = μ i+1 > 0 and

T (t i ) -T (t i+1 ) = β i+1 = λ 2 t i+1 t i ẋ2 dt > 0 (26)
Note that for fixed and bounded coefficients in (1) μ i and β i depend only on δ (the other "initial" condition on the position needed to integrate the system on the interval [t i , t i+1 ] remaining fixed at the times t i ) so that in particular β i ≥ β(δ) > 0 for all i ≥ 0 and δ > 0. Since T (t i ) is non-increasing, its limit value is its minimum value and for all i ≥ 0,

| ẋ(t i )| ≥ | ẋss | > δ.
From the strictly positive variation of the kinetic energy we deduce that

ẋ2 (t i+1 ) = ẋ2 (t i ) - 2β i m ( 27
)
so that

ẋ2 i = ẋ2 0 - j=i-1 j=0 β j (28)
Therefore from the fact that the β i 's are strictly positive, we deduce that | ẋ(t i )| cannot converge towards a strictly positive | ẋss |. Since however T (t i ) and thus ẋ(t i ) converge, we deduce that the only possible limit value for the velocity is ẋss = 0. (Notice that if δ = 0, then both μ i and β i may asymptotically take arbitrarily small values and ẋ2

i = ẋ2 0 - j=i-1 j=0
β j no longer leads to a contradiction). Thus we have shown that if there is an infinite number of bounces, then the value of the velocity when contact is established or lost (x(t i ) = 0) is bounded and tends to zero.

Having proved that the velocity ẋ(t i ) tends to zero as i → +∞, we now show that the intervals i also converge to zero. Let us consider an arbitrarily large integer i such that | ẋ(t i )| is arbitrarily small, or in other words, for any ε > 0, there exists N (ε) > 0 such that i > N implies | ẋ(t i )| < ε. We shall denote i+1 = t i+1t i . First note that from any of the two dynamic equations in (1) we get i ≤ max < +∞ for some max since the "initial" velocities at times t i are bounded and tend towards zero. Now we use the fact that both vector fields in (1) are explicitely integrable; assume that we place ourselves at t 2i such that ẋ(t 2i ) = ε > 0, hence the system is in a contact phase for some time since ẍ(t i ) = λ 1 x d -λ 2 ε > 0 for some small enough ε > 0. We thus consider the second equation in [START_REF] An | Dynamic stability issues in force control of manipulators[END_REF]. If the negative roots r 1 and r 2 of the characteristic equation are real and separated, i.e. r 1 < r 2 , then the solution can be expressed as (recall that x(t i ) = 0 for all i):

x(t) = γ 1 e r 1 (t-t 2i ) + γ 2 e r 2 (t-t 2i ) + xd (29)
with xd = λ 1 x d λ 1 +k , and γ 1 = -γ 2 -xd , γ 1 r 1 = -γ 2 r 2 + ε. Since we assume a priori that the sequence {t i } is infinite, t 2i+1 exists and from (29) we get

ẋ(t 2i+1 ) = γ 1 r 1 (e r 1 2i+1 -(1 -ε)e r 2 2i+1 ) = γ 1 r 1 e r 1 2i+1 1 -1 - ε γ 1 r 1 e (r 2 -r 1 ) 2i+1 + 1 (30)
From the monotonicity of {| ẋ(t i )|} and its convergence, we deduce that 0 < | ẋ(t 2i+1 )| ≤ ε.

Assume now that the sequence { 2i+1 } does not converge towards zero, i.e. there exists > 0 such that 2i+1 ≥ for all i. Then we get for any ε > 0, ε small enough:

1 - ε γ 1 r 1 e (r 2 -r 1 ) 2i+1 -1 ≥ η(r 1 , r 2 , ) > 0 (31)
This inequality is true because when ε → 0 the left-hand-side of (31) tends to |e (r 2 -r 1 ) 2i+1 -1| which cannot be zero since r 1 < r 2 and 2i+1 ≥ > 0. Furthermore e r 1 2i+1 ≥ κ( max , r 1 ) > 0. Therefore from (30) we get:

|γ 1 r 1 κ( max , r 1 )η(r 1 , r 2 , )| < ε (32)
which cannot be true for ε small enough (Note that the roots as well as and max do not depend on ε). Since ε is strictly positive but arbitrarily small, we deduce that 2i+1 → 0 as i → +∞. A quite similar reasoning may be done for the case when r 1 = r 2 .

When the roots are complex conjugate r 1 = r + jω, r 2 = rjω, then the solution is given by:

x(t) = γ e r (t-t 2i ) cos(ω(t -t 2i ) + ϕ) + xd ( 33
)
with γ = -xd cos ϕ and tan ϕ = ε+ xd r xd ω . Now we get

ẋ(t 2i+1 ) = γ e r 2i+1 r 2 + ω 2 cos(ω 2i+1 + ϕ + ) ( 34
)
with tan = ω r . Using the same arguments as in the real roots case, one sees that for ẋ(t 2i+1 ) to be arbitrarily small, we must have cos(ω 2i+1 + ϕ + ) arbitrarily small, from which we deduce that ω 2i+1 + ϕ + is arbitrarily close to π 2 . Now for ε arbitrarily small, tan ϕ → r ω , and tan(ϕ + ) → +∞. But since 2i+1 is assumed to be bounded away from zero (and strictly positive by definition), tan( π 2 -ω 2i+1 ) is clearly bounded. Thus by contradiction we deduce that { 2i+1 } converges to zero. Now exactly the same reasoning may be done for the case of non-contact phases. It follows that if the velocities at times t i converge towards zero, so do the intervals i . Since again the sequence {t i } is infinite, if its limit is infinite also then (0, 0) is an equilibrium point of the system in [START_REF] An | Dynamic stability issues in force control of manipulators[END_REF]. Clearly this is not the case, except if x d = 0 (For the sake of briefness this case is not analyzed here; the analysis can be done using similar arguments).

In conclusion, we have proved that the sequence {t i } is either finite, or possesses a finite accumulation point. In both cases, we deduce that the equilibrium point of the system in (1) is asymptotically attained.

Relationship with the Case of a Rigid Environment

The only things that are modified in the rigid case are that since the intervals [t 2i , t 2i+1 ] → {t 2i }, the distinction between instants t 2i and t 2i+1 becomes worthless, and μ 2i+1 = β 2i+1 = 0 while μ 2i > 0 and β 2i > 0. One notes that in permanent contact (i.e. x = ẋ ≡ 0) and with k = +∞, then

V c = 1 2 m ẋ2 + 1 2 (λ 1 + k)x 2 -λ 1 x x d + 1 2 λ 2 1 x d λ 1 + k (35)
i.e.

V c = 1 2 kx 2 = 0 (36)
since the potential elastic energy vanishes as k → +∞ (indeed it can be verified that in this case, the roots r 1 and r 2 are necessarily complex conjugate, and that x(t) varies proportionally to

1 k ). Moreover σ V c (t k ) = σ V nc (t k ) = 0.
It follows that the contact phases reduce to the impact times t k . One retrieves the analysis done at the beginning of this section by studying the variation of V nc between the impacts (see [START_REF] Slotine | Adaptive manipulator control: A case study[END_REF]) and at the impact times (see [START_REF] Tornambé | Global regulation of planar robot arm striking a surface[END_REF]). Therefore the stability analysis for the compliant case k < +∞ naturally reduces to the stability analysis of the impact Poincaré map in the rigid case. In addition we have proved asymptotic stability.

We thus have proved the following: Claim 3. Consider the closed-loop equations in [START_REF] An | Dynamic stability issues in force control of manipulators[END_REF]. Then for any λ 1 > 0, λ 2 > 0, k ∈ [0, +∞], and for all initial conditions x(0), ẋ(0), x → λ 1 x d λ 1 +k and ẋ → 0 as t → +∞.

Remark 1.

A distinction has to be made between two different cases of analysis: We may consider i) Either an arbitrarily large but bounded k, ii) or a k that tends to infinity (that is implicitely a sequence of stiffnesses k n with unbounded limit together with the corresponding dynamics). Clearly claim 3 can be concluded from the analysis in sections 2.2 and 3 in case i), but not in case ii). The utility of the analysis proposed in this section is to enable us to draw conclusions in both the compliant and the rigid environment cases with a unique philosophy of stability analysis.

Remark 2. Since the solutions of (1) converge to those of (2), it would be interesting to reverse the reasoning done in section 4 as follows: if the fixed point of ( 2) is stable, then the fixed point of ( 1) is stable also for a large enough k. One could think of first studying the stability of the rigid body system impact Poincaré map, and then draw conclusions on the stability of an approximating compliant problem. This of course relies on the ability of proving the closeness of solutions of both systems for large enough k. Moreover it may be difficult to study the Poincaré map stability.

Remark 3. It is not realistic in general to assume that the interaction force F is measurable and can be compensated for, during the rebounds phase. Indeed if k is large, then F will be large also and its support (as a time function) very small (F approaches a Dirac distribution). This is a motivation to search for feedback control algorithms that are able to stabilize the transition phases without force feedback, see e.g. [START_REF] Tornambé | Modeling and controling two-degree-of freedom impacts[END_REF]. In fact it has been recognized [START_REF] Volpe | A theoretical and experimental investigation of impact control for manipulators[END_REF] that in general 3 distinct controllers have to be used for the control of a complete robotic task (i.e. a task involving free-motion as well as constrained motion phases, and transition phases). The algorithm analyzed in this paper could be used in this setting as a transition phase controller, as part of a more sophisticated controller (in the spirit of what has been proposed e.g. in [12] [3]). This could also be the starting point to a unified stability analysis framework encompassing both compliant and rigid body based models.

Conclusions

In this note, we have studied the stability of a simple robotic task that includes both contact and non-contact phases, when a PD motion controller is applied. It may also be considered as the stability analysis of a contact phase that takes into account the unilateral feature of the constraint surface. An interesting problem is to study conditions that guarantee that if the robot's tip happens to take off the surface (either because it has been initialized outside the environment or because of an oscillatory behaviour or due to a disturbance) then it will eventually attain its steady state value after a possible sequence of rebounds. We have studied what happens when the environment's stiffness grows very large and when one desires quadratic stability of the closed-loop system with well-conditioned solutions P to the Lyapunov equation, i.e. when the minimum eigenvalue of P is required to remain larger than some strictly positive constant. It appears that the velocity feedback gain must be chosen proportional to the environment's stiffness k, so that these sufficient conditions are clearly useless for practical purposes as soon as the contact stiffness is too large. The results indicate that the conditions that one may derive from an analysis done with continuous dynamic models may not be feasible for the environments which are too rigid; although we do not claim that such conclusions hold for any control algorithm for robot manipulators with any compliant environment, we note that this is quite consistent with the sufficient conditions found in the literature on the subject. Finally we have shown that one can conclude on global asymptotic convergence of the tracking errors towards zero independently of the values of the feedback gains and of the environment's stiffness. This is in accordance with the results obtained for the impact Poincaré map in the case of a rigid body model. We have proposed a convergence proof that encompasses both compliant and rigid cases, a property that is clearly not shared by some the other classical stability analysis.

A. Proof of Claim 1

λ 1 and λ 2 can be chosen such that a nc > 0 and b nc > 0. Thus we conclude that for all x:

V ≤ -α ẋ2 -β x2 + R ( 37
)
with α = min(a nc , a c ), β = min(b nc , b c ). Following the arguments in [START_REF] Corless | Continuous state feedback guaranteing unifrom ultimate boundedness for uncertain dynamical systems[END_REF], we deduce that the state ( x, ẋ) converges in finite time in a ball with radius r , with r → 0 as λ 1 and λ 2 tend to +∞. Therefore for all t ≥ t, t < +∞, we get | x| < r . Now notice (see ( 4)) that as λ 1 → +∞ then:

R → 1 2 k 2 x 2 d 1 + c 2 m 2 (38)
Since by taking λ 1 and λ 2 large enough r can be made arbitrarily small and since λ 1 x d λ 1 +k → x d as λ 1 → +∞, it follows that for λ 1 and λ 2 large enough, then | x(t)| > r for t ≥ t implies that x(t) > 0 for all t ≥ t. Then [START_REF] Corless | Continuous state feedback guaranteing unifrom ultimate boundedness for uncertain dynamical systems[END_REF] implies that both ẋ and x converge asymptotically to zero. Notice that outside some ball B R we have for some

Q = Q T > 0: V ≤ -z T Qz (39)
with z T = ( x, ẋ). This can be easily deduced by spliting α and β into α 1 and α 2 , β 1 and β 2 . Then one obtains:

V ≤ -α 1 ẋ2 -β 1 x2 -α 2 ẋ2 -β 2 x2 + R (40)
so that outside the ball B R with R = R min(α 2 ,β 2 ) we get:

V ≤ -α 1 ẋ2 -β 1 x2 (41)
Still λ 1 and λ 2 can be chosen large enough so that R is as small as desired. Thus we deduce that for 0 ≤ t ≤ t:

||z(t)|| ≤ V (0) λ min P exp - λ min Q λ max P t (42)
i.e. the ball B R is reached exponentially fast.

B. Proof of Claim 2

Starting from the Lyapunov equation, one may first fix Q c as a positive definite matrix and then try to calculate the unique corresponding positive definite P (see e.g. [START_REF] Vidyasagar | Nonlinear Systems Analysis[END_REF] lemma 42, chapter 5). A second way to attack the problem is to pick a P > 0 and study the properties of the resulting Q c ([20] p.198). Instead of choosing a Q c > 0 and solving the Lyapunov equation for P, we rather consider a matrix P and find conditions such that the corresponding Q c is positive definite, together with Q nc . Thus we prove that the only way for P not to tend towards a singular matrix while keeping Q c > 0 and Q nc > 0 when k increases is to take the gain λ 1 of the same order as k 2 .

The above determinants can be written in the following way:

det(Q nc ) =
• det(Q nc ) > 0 ⇔ 4λ 1 m|P| -Y 2 > 0 • det(Q c ) > 0 ⇔ Y 2 + 2kp 22 Y -4(λ 1 + k)m|P| + (kp 22 ) 2 < 0
We deduce that Y has to satisfy the following inequalities: k p 11 we can find P that satisfies (44) and that is positive-definite. From (43), Y satisfies the following inequalities:

-2 mλ 1 |P| < Y < 2 mλ 1 |P| -k p 22 -2 m(λ 1 + k)|P| < Y < -k p 22 + 2 m(λ 1 + k)|P| (43) Since -k p 22 -2 √ m(λ 1 + k)|P| < -2 √ mλ 1 |P|
-2 mλ 1 |P| < Y < min(-k p 22 + 2 m(λ 1 + k)|P|, 2 mλ 1 |P|) (45)
It is easy to prove that λ min (P) ≤ p 22 and λ max (P) ≥ p 11 from which we deduce that

λ max (P) λ min (P) ≥ p 11 p 22 (46)
From (44), we can write p 22 < 4 2 k p 11 , thus P has bounded entries when p 11 is bounded and the above conditions fulfilled. Then if p 11 is a finite real number the conditions of existence of Y imply that the coefficients p 12 and p 22 tend to zero when the stiffness of the environment becomes infinite, rendering the matrix P singular. Let us note that the stability analysis then becomes asymptotically (i.e. when k → +∞) meaningless since Q nc in ( 14) has bounded entries. The only way to avoid this problem is to increase the gain λ 1 such that the coefficient k does not tend towards zero when the stiffness increases, i.e. λ 1 has to be chosen of order ≥ k 2 . Assume that this is done so that P is well conditioned, and let us examine how λ 2 has to be chosen. λ 2 may be found by using (45):

mp 11 + λ 1 p 22 -Y max p 12 < λ 2 < mp 11 + λ 1 p 22 + 2 √ mλ 1 |P| p 12 ( 47
)
where

Y max = min(-k p 22 + 2 √ m(λ 1 + k)|P|, 2 √ mλ 1 |P|
). This implies that when λ 1 is of order k 2 and k grows unbounded, the gain λ 2 becomes infinite too.

Let us examine what happens if we allow p 11 to be proportional to k α , α > 1. Then p 22 may be chosen of order ≤ k α-1 from (44). Also p 12 will be of order ≤ k α-

α-1 ; then Y max is of order k α-1 2 so that λ 2 grows as k 1 2 . Now if Y max = A = -k p 22 + 2 √ m(λ 1 + k)
|P| that we assume > 0: then necessarily since p 22 ≥ 0, the second term in A is at least of the same order as kp 22 in k as k grows unbounded. Thus at most the order of Y max will be that of the second term 2 √ m(λ 1 + k)|P|, which is found to be at most k α . But if this is the case then this term dominates 2 √ mλ 1 |P| and asymptotically (in k) Y max will necessarily be equal to this last term, hence we are back to the previous case. Now if the order of A is k γ with γ < α then p 11 will asymptotically dominate Y max and the left-hand-side of (47) is asymptotically of order k 1 2 . Thus λ 1 may be chosen bounded but λ 2 will grow unbounded to guarantee λ min P ≥ δ > 0 for any arbitrarily small but fixed δ and Q c > 0, Q nc > 0.

C. Impulsive Dynamics

For the sake of briefness, we shall not recall here all the details about impulsive impact dynamics. Let us simply recall some basic facts: The interaction forces between two rigid bodies are mathematically modelled by Dirac distributions pδ t , whose magnitude p may be calculated from the velocities before the impact and using a physical law of percussion like Newton's restitution coefficient 0 ≤ e ≤ 1. e = 1 corresponds to the case when there is no loss of energy at the impact, so that in our case both velocities of the mass after and before the impact are equal in magnitude (and of course of opposite signs). In (2) we get p(t k ) = -2m ẋ(t - k ). Equations containing Dirac measures are called measure differential equations; the solutions of such equations are right-continuous functions of local bounded variation [START_REF] Schmaedeke | Optimal control theory for nonlinear vector differential equations containing measures[END_REF]. Stability in the sense of Lyapunov and Lyapunov functions can still be considered for measure differential equations [START_REF] Leela | Stability of measure differential equations[END_REF] 3 : One has however to consider in place of the traditional derivative V of a positive definite function V its right Dini derivative 4 in intervals between impacts (smooth dynamics), and its jumps σ V (t k ) = V (t + k ) -V (t - k ) at instants of discontinuities t k , i ≥ 0. The latter may also be seen as an application of a generalized chain rule for distributional derivatives to the function V (x(t), ẋ(t)) where ẋ(t) has local bounded variation and V (•, •) is continuously differentiable: Then at t = t k , DV = σ V (t k )δ t k (D usually denotes distributional derivatives of functions of bounded variation, see e.g. [START_REF] Schmaedeke | Optimal control theory for nonlinear vector differential equations containing measures[END_REF]). Rigorous convergence of sequences of continuous-dynamics problems (or "approximating problems") towards a nonsmooth dynamical problem when k → +∞ has been treated for instance in [START_REF] Paoli | Mouvement à un nombre fini de degrés de liberté avec contraintes unilatérales: cas avec perte d'énergie[END_REF].

Notes

1. For instance [START_REF] Youcef-Toumi | Impact and force control: Modeling and experiments[END_REF] show that an integral force feedback helps in stabilizing the impact phase when the environment is (sufficiently) rigid, whereas [START_REF] Volpe | A theoretical and experimental investigation of impact control for manipulators[END_REF] show that it is not suitable for a (sufficiently) flexible environment. 2. Notice that we have not proved the asymptotic stability of P .

3. Let us note however that the results in [START_REF] Leela | Stability of measure differential equations[END_REF] do not apply here since the author studies stability of systems of the form ẋ = f (t, x) + Du, u of local bounded variation, and Du is to be considered as a disturbance on ẋ = f (t, x) with globally asymptotically stable equilibrium point x = 0; In our case the impacts will drive the system to another point than the equilibrium point of the smooth dynamics (that may even possess no equilibrium point: The reader may think of the bouncing ball problem to illustrate clearly this). More details can be found in [START_REF] Brogliato | Nonsmooth Impact Mechanics. Models, Dynamics and Control[END_REF] chapter 7.

V (t) = lim h→0

+ sup 1 h [V (t + h) -V (t)].

Figure 1 .

 1 Figure 1. Controlled mass colliding an elastic wall.

Figure 2 .

 2 Figure 2. definition of the times t i (1.a) and t k (1.b).

  (λ 1 p 22 + mp 11 -λ 2 p 12 ) 2 -2 kp 22 m 2 (λ 1 p 22 + mp 11 -λ 2 p 12 ) Let us denote Y = λ 1 p 22 + mp 11 -λ 2 p 12 and |P| = p 11 p 22p 2 12 then:

	4	λ 1 m	( p 11 p 22 -p 2 12 ) -	λ 1 m	p 22 + p 11 -	m λ 2	p 12	2
	det(Q c ) = 4	λ 1 + k m	( p 11 p 22 -p 2 12 ) -	λ 1 + k m	p 22 + p 11 -	m λ 2	p 12	2
	1 m 2 + 4 = -λ 1 + k m	( p 11 p 22 -p 2 12 ) -	m kp 22	2

  1 2 from the second condition in (44). Now from (47) we have the following: If Y max ≤ 0 then obviously λ 2 is of order k 1 2 as k → +∞. If Y max > 0, let us analyze the case when Y max = 2 √ mλ 1 |P|: This value is maximum when p 12 is minimum, i.e. bounded, and when both p 11 and p 22 are maximum, i.e. respectively of orders k α and k


Xiaoli Guo
email: xiaoli.guo@ec-nantes.fr

Thomas Heuzé
email: thomas.heuze@ec-nantes.fr

Guillaume Racineux
email: guillaume.racineux@ec-nantes.fr

Ramzi Othman
email: ramzi.othman@ec-nantes.fr

Ramzi Othman Direct

very high strain rate and inverse identification of elastic-viscoplastic constitutive models

Keywords: Direct impact Kolsky bar, Inverse analysis, High strain rate, elastic-viscoplastic constitutive models

Introduction

The study of the material performance at high strain rate is an important topic in material sciences, since the materials usually exhibit strong viscous effects as the rate rises. In some material processes, the metallic materials may experience very high strain rates up to 10 6 s -1 , such as the laser shock peening [START_REF] Yu | Laser dynamic forming of functional materials laminated composites on patterned three-dimensional surfaces with applications on flexible microelectromechanical systems[END_REF]. Studying the dynamic behaviour of the materials requires to perform dynamic tests on a wide range of strain rate. The classical Split Hopkionson Pressure Bar (SHPB), also referred to as the Kolsky bar [START_REF] Davies | A critical study of the Hopkinson Pressure Bar[END_REF][START_REF] Kolsky | An investigation of the mechanical properties of materials at very high rates of loading[END_REF], is the most widely used apparatus to perform such dynamic tests. The classical SHPB device is capable to attain strain rates ranging from 10 2 to 10 4 s -1 [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF][START_REF] Ramesh | Part D : Chapter 33. High strain rate and impact experiments[END_REF]. Attempting to reach higher strain rate may yield the incident bar. Thus a direct-impact Hopkinson device has been developed [START_REF] Dharan | Determination of stress-Strain characteristics at very high strain rates[END_REF] by removing the incident bar. On this direct-impact device, a very high strain rate of the order of 10 5 s -1 is achieved by [START_REF] Dharan | Determination of stress-Strain characteristics at very high strain rates[END_REF] and [START_REF] Kamler | Measurement of the behaviour of high-purity copper at very high rates of strain[END_REF]. In this work, the Ti-6Al-4V alloy is tested on a direct-impact Hopkinson device at strain rates ranging from 3000 to 25000s -1 . Then an inverse identification is carried out to identify three elastic-viscoplastic constitutive models for the Ti-6Al-4V on this wide range of strain rates.

Direct-impact tests on the Ti-6Al-4V at very high strain rate

The dynamic tests are performed on a dedicated direct-impact device which design is solution of an optimization problem as detailed in [START_REF] Guo | Inverse identification at very high strain rate of the Johnson-Cook constitutive model on the Ti-6Al-4V alloy with a specially designed direct-impact Kolsky bar device[END_REF]. This device consists mainly of a transmitted bar, the projectile and the compulsory systems for the measurement. The transmitted bar is 1.2m length and 10mm diameter. To compress the specimen on a wide range of strain rate, four projectiles are manufactured, of diameter 15.8mm and lengths of 500mm, 125mm, 60mm and 30mm respectively. The bar and the projectiles are all made of a high strength steel, MARVAL X2NiCoMo18-8-5, with a yield stress of 1800MPa. Three Wheastone bridges of double gauges are mounted along the bar length to measure the elastic strain and monitor the dispersion. Two laser diodes are equipped to measure the impact velocity of the projectile. A high speed camera is also equipped to observe the compression of the specimen and the movement of the projectile during the impact.

The experimental plan comes down to determine the value of the specimen length l s , the projectile length l p and its impact velocity v p . However, the combination of these three parameters is not unique to attain the expected strain rate. Provided the expected strain rate εs and the allowable strain ε s in the specimen, their maximum values can be approximated through the following equations, c p being the sound speed in the projectile :

εsmax ≈ v p l s , ε smax ≈ εs 2l p c p (1)
The combination of the parameters for each objective strain rate is presented in table The measured data consist of the elastic strain ε b (t) in the transmitted bar and the impact velocity v p of the projectile. The force F s (t) applied by the bar on the specimen is computed from the recorded strain ε b (t) by

F s (t) = E b S b ε b (t),
where E b and S b are the Young's modulus and the cross-section of the bar. The force F s (t) computed for the tests listed in table 1 is plotted in the figure 1 for each projectile length. In these curves, some trays are observed during the unloading part. These arise from (i) the shorter length of the projectile than that of the bar and (ii) the mismatch of the generalized wave impedances [START_REF] Wang | Foundations of stress waves[END_REF] of the projectile (ρcS) p , the specimen (ρcS) s and the transmitted bar (ρcS) b . Moreover, when the projectile is long enough, the loading time is equivalent to the characteristic time, as shown in the figure 1(a) where the projectile is 500mm length. As the projectile length shortens, the actual period of the compression of the specimen lasts approximately twice the characteristic time or even longer as observed in the figures 1(b)-1(d), resulting from the different cross-sections of the projectile and the bar in these experiments. As a consequence, the specimen in these tests experiences a second-time compression before the unloading wave arrives, and thus obtained a greater residual strain than that of the test T1 as shown in the table 1, even though the greatest kinematic energy is achieved in this test.

Without any input bar to measure the input force, the strain rate in the specimen is assessed by the following equation [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF][START_REF] Guo | Inverse identification at very high strain rate of the Johnson-Cook constitutive model on the Ti-6Al-4V alloy with a specially designed direct-impact Kolsky bar device[END_REF] :

εs (t) = - v p + Sp+S b Sp c b ε b (t) l s (2)
where S p and c b are the cross-section of the projectile and the sound speed of the bar. The maximum compressive strain rate εmax assessed by this equation and the measured compressive strain ε max are presented in table 1 for each test.

According to this formula, the compressive strain rate gets its extremum equal to v p /l s at the beginning of the impact. However, this equation is only valid within the characteristic time [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF]. The stress-strain curves are thus assessed using the classical analysis of the Hopkinson tests [START_REF] Guo | On the direct impact Hopkinson system for dynamic tests at very high stain rates[END_REF] as shown in the figure 2. The curves exhibit a rising strain hardening as the strain rate increases. They are subject to scattering, as we 1: Force on the specimen in the direct-impact tests can observe through small oscillations and a smaller flow stress appear in the curves of the tests T8 and T9. Oscillations are also observed at higher strain rate in the figure 2(d). Consequently, the rate dependency is difficult to be assessed on these curves directly. An inverse analysis is carried out for the identification of the constitutive models.

Inverse identification of elastic-viscoplastic constitutive models

Three elastic-viscoplastic constitutive models have been chosen for their potential abilities to fit correctly the experimental measurements and are identified by means of an inverse procedure involving a dynamic numerical analysis.

No thermal effect is addressed in these identifications. The Johnson-Cook model [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] in isothermal conditions reads :

σ(p, ṗ) = (A + Bp n ) 1 + C ln ṗ ṗ0 (3)
where p and ṗ are the cumulated plastic strain and strain rate respectively. The parameters A, B and n are identified on quasi-static stress-strain curves obtained at the strain rate of 10 -4 s -1 [10] : A=955MPa, B=770MPa and n=0.557. Thus only the parameter C is to be identified in the inverse process. An axisymmetric 2D Finite Element model in ABAQUS is created to proceed the dynamic simulation, consisting of the projectile, the specimen and the bar all meshed. The Zerilli-Armstrong model [START_REF] Zerilli | Armstrong. Constitutive relations for titanium and ti-6al-4v[END_REF] is calibrated on the same experimental database, using a unidimensional FE model in MATLAB developed by [START_REF] Andriamiseza | Rapport de stage de Master : Identification de modèles de comportement à haute vitesse de déformation sur des essais de Kolsky à impact direct[END_REF]. Two variants for bcc and hcp materials are considered, that read :

σ y = C 0 + C 1 exp(-C 3 T + C 4 T ln ṗ) + C 5 p n (bcc) (4) σ y = C 0 + C 1 exp(-C 3 T + C 4 T ln ṗ) + C 2 √ p exp(-C 3 T + C 4 T ln ṗ) (hcp) (5)
Assuming a constant T , the first formula (4) reduces to :

σ y = K 0 + K 1 p n + K 2 ṗK3 (bcc) (6)
where

K 0 = C 0 , K 1 = C 5 , K 2 = C 1 exp(-C 3
σ y = K 0 + K 1 ṗK2 + K 3 √ p ṗK4 (hcp) (7)
where

K 0 = C 0 , K 1 = C 1 exp(-C 3 T ), K 2 = C 4 T , K 3 = C 2 exp(-C 3 T ) and K 4 = C 4
T . K 0 is equivalent to the yield stress, hence the unknown vector consists of four parameters K 1 , K 2 , K 3 and K 4 during the inverse identification. The Nouailhas model [START_REF] Nouailhas | Unified modelling of cyclic viscoplasticity : application to austenitic stainless steels[END_REF] has been used with success for adhesively bonded assemblies to account for viscous effects [START_REF] Créac'hcadec | On modelling the non-linear behaviour of thin adhesive films in bonded assemblies with interface elements[END_REF]. The creep law of this model is expressed as :

ṗ = |σ| -R -σ y K n exp α |σ| -R -σ y K n+1 ( 8
)
where x denotes the positive part of x, that is x = x+|x| 2 . The parameter α enables to saturate the viscous effects. An isotropic hardening power law [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF] Ṙ = Q ṗm is adopted. The first three parameters σ y , Q, m are also calibrated on the quasi-static tests from [START_REF] Guo | On the direct impact Hopkinson system for dynamic tests at very high stain rates[END_REF] : σ y = 955MPa, Q = 770MPa, m = 0.557. The material parameters associated to the rate dependency are thus K, n and α.

The cost function f (x) is defined as the euclidean norm of the difference between the simulated strain ε sim (x, t) extracted at the location of the gauge and the recorded one ε exp (t). The identification is proceeded on each tests of the table 1 ; the calibrated strain ε sim of the tests T1, T4, T6, T8 and T11 are plotted and superposed on the experimental one ε exp in the figures 3(a)-3(e). An error bar is drawn associated to each mean value of the parameter C for the 2D identification, as shown in the figure 3(f), since at least two repetitions are carried out for each test. The simulated strains of different constitutive models are quite well superposed with the experimental one at strain rate up to 15000s -1 . A good concordance is observed at the second raising stage when the specimen is plastically deformed as shown in the figures 3(a)-3(c). Small discrepancy appears at this stage between the calibrated strain of the Johnson-Cook model and the measured one in the figure 3(d), whereas the Zerilli-Armstrong fits well the experimental data. As the strain rate goes higher, a greater strain is obtained from the test, while an equivalent amplitude of strain is simulated by the two constitutive models as presented in the figure 3(e). The figure 3(f) indicates that an approximately constant value of the parameter C of the Johnson-Cook model is identified at the strain rate ranging from 3000 to 18000s -1 . However a much greater value is obtained at very high strain rate in the tests T10 and T11. The discrepancies in these two tests may be due to (i) the very short projectile used in the test and/or (ii) the thermal effects not addressed in the identification. More repetitions of these two tests are to be carried out to enforce or not the results. The stress-strain curves of the identified constitutive models are plotted at the strain rate of 3000s -1 and 18000s -1 and superposed to those obtained from the quasi-static tests in the figure 4. The figures 4(a)-4(c) show a great raise of the flow stress from the quasi-static to the dynamic conditions. But a much slighter increase is observed when the strain rate goes from 3000 to 18000s -1 , especially for the Johnson-Cook model, where the two curves at these two strain rates are almost superposed in the figure 4(a). Indeed the rate dependency of this alloy has not been found greater than that predicted by the Johnson-Cook model at least in the studied range of strain rate, and so does for the other models.

Conclusion

In this work, the Ti-6Al-4V has been tested on a wide range of strain rate ranging from 3000 to 25000s -1 using the direct-impact Hopkinson device. An inverse analysis is carried out to identify three elastic-viscoplastic constitutive models on the dynamic experimental data. All three identified constitutive models fit well the experimental data at the strain rate ranging from 3000 up to 18000s -1 . However more experimental and numerical research at higher strain rate is required to refine these results.

F

F 2 :

 2 Stress-strain curves of the Ti-6Al-4V

  T ) and K 3 = C 4 T . The parameters K 0 , K 1 , n have the same value as the parameters A, B, n of the Johnson-Cook relation because of the equivalent form of strain hardening. Only two parameters K2 and K3 have been calibrated inversely. The second formula (5) reduces to :

F 3 :

 3 Error bars of parameter C Results of the identifications

F 4 :

 4 Stress-strain curves plotted within the range of 10 -4 to 18000s-1


Xiaoli Guo

Thomas Heuzé

Guillaume Racineux

Ramzi Othman

Direct-impact Hopkinson tests on the Ti-6Al-4V alloy at very high strain rate and inverse identification of the Johnson-Cook constitutive model

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

Direct-impact Hopkinson tests on the Ti-6Al-4V alloy at very high strain rate and inverse identification of the Johnson-Cook constitutive model X. GUO, T. Heuzé, G. Racineux. Research Institute in Civil and Mechanical Engineering (GeM, UMR 6183 CNRS), École Centrale de Nantes, 1 rue de la Noë, 44321 Nantes cedex 3, France Adresse(s) électronique(s) : {Xiaoli.Guo, Thomas.Heuze, Guillaume.Racineux}@ec-nantes.fr R. Othman. Mechanical engineering department, Faculty Engineering, King Abdulaziz University, P.O.BOX, 80248, Jeddah, 21589, Saudi Arabia Adresse(s) électronique(s) : Ramzi.Othman@ec-nantes.fr Mots clés : Direct impact Kolsky bar, Inverse analysis, High strain rate, Johnson-Cook constitutive model

INTRODUCTION

The study of the material performance at high strain rate is an important topic in material sciences, since the materials usually exhibit strong viscous effects as the rate rises. In some material processes, the metallic materials may experience very high strain rates up to 10 6 s -1 , such as the laser shock peening [START_REF] Yu | Laser dynamic forming of functional materials laminated composites on patterned three-dimensional surfaces with applications on flexible microelectromechanical systems[END_REF]. Studying the dynamic behaviour of the materials requires to perform dynamic tests on a wide range of strain rate. The classical Split Hopkionson Pressure Bar (SHPB), also referred to as the Kolsky bar [START_REF] Davies | A critical study of the Hopkinson Pressure Bar[END_REF][START_REF] Kolsky | An investigation of the mechanical properties of materials at very high rates of loading[END_REF], is a simple device that is widely adopted to perform such dynamic tests. The classical SHPB device is capable to attain strain rates ranging from 10 2 to 10 4 s -1 [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF][START_REF] Ramesh | Part D : Chapter 33. High strain rate and impact experiments[END_REF]. Attempting to reach higher strain rate may yield the incident bar. Thus a direct-impact Hopkinson device has been developed [START_REF] Dharan | Determination of stress-strain characteristics at very high strain rates[END_REF] by removing the incident bar. On this direct-impact device, a very high strain rate of the order of 10 5 s -1 is achieved by [START_REF] Dharan | Determination of stress-strain characteristics at very high strain rates[END_REF] and [START_REF] Kamler | Measurement of the behaviour of high-purity copper at very high rates of strain[END_REF]. In this work, the Ti-6Al-4V alloy is tested on a direct-impact Hopkinson device at strain rates ranging from 3000 to 25000s -1 . Then an inverse identification is carried out to identify the Johnson-Cook model [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] for the Ti-6Al-4V on this wide range of strain rates.

2 DIRECT-IMPACT TESTS ON THE TI-6AL-4V AT VERY HIGH STRAIN RATE 2.1 Design of the direct-impact configuration and experimental plan The dynamic tests are performed on a dedicated direct-impact device which design is solution of an optimization problem as detailed in [START_REF] Guo | Inverse identification at very high strain rate of the Johnson-Cook constitutive model on the Ti-6Al-4V alloy with a specially designed directimpact Kolsky bar device[END_REF]. This device consists mainly of a transmitted bar, the projectile and the compulsory systems for the measurement. The transmitted bar is 1.2m length and 10mm diameter. To compress the specimen on a wide range of strain rate, four projectiles are manufactured, of diameter 15.8mm and lengths of 500mm, 125mm, 60mm and 30mm respectively. The bar and the projectiles are all made of a high strength steel, MARVAL X2NiCoMo18-8-5, with a yield stress of 1800MPa. Three Wheastone bridges of double gauges are mounted along the bar length to measure the elastic strain and monitor the dispersion. Two laser diodes are equipped to measure the impact velocity of the projectile. A high speed camera is also equipped to observe the compression of the specimen and the movement of the projectile during the impact. The experimental plan comes down to determine the value of the specimen length l s , the projectile length l p and its impact velocity v p . However, the combination of these three parameters is not unique to get a given strain rate. Provided the expected strain rate εs and the allowable strain ε s in the specimen, their maximum values can be approximated through the following equations, c p being the sound speed in the projectile :

εsmax ≈ v p l s , ε smax ≈ εs 2l p c p (1)
The combination of the parameters for each objective strain rate is presented in table 1.

Measurements from the direct-impact tests

The measured data consist of the elastic strain ε b (t) in the transmitted bar and the impact velocity v p of the projectile. The force F s (t) applied by the bar on the specimen is computed from the recorded strain ε b (t) :

F s (t) = E b S b ε b (t) (2)
where E b and S b are the Young's modulus and the cross-section of the bar. The force F s (t) computed for the tests listed in table 1 is plotted in the figure 1 for each projectile length. In these curves, some trays are observed during the unloading part. These arise from (i) the shorter length of the projectile than that of the bar and (ii) the mismatch of the generalized wave impedances [START_REF] Wang | Foundations of stress waves[END_REF] of the projectile (ρcS) p , the specimen (ρcS) s and the transmitted bar (ρcS) b . Moreover, when the projectile is long enough, the loading time is equivalent to the characteristic time, as shown in the figure 1(a) where the projectile is 500mm length. As the projectile length shortens, the actual period of the compression of the specimen lasts approximately twice the characteristic time or even longer as observed in the figures 1(b)-1(d). This results from the different cross-sections of the projectile and the bar in these experiments. FIGURE 1 -Force on the specimen in the direct-impact tests Without any input bar to measure the input force, the strain rate in the specimen is assessed by the following equation [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF][START_REF] Guo | Inverse identification at very high strain rate of the Johnson-Cook constitutive model on the Ti-6Al-4V alloy with a specially designed directimpact Kolsky bar device[END_REF] :

εs (t) = - v p + Sp+S b Sp c b ε b (t) l s (3)
where S p and c b are the cross-section of the projectile and the sound speed of the bar. The maximum compressive strain rate εmax assessed by this equation and the measured compressive strain ε max are presented in table 1 for each test. According to this formula, the compressive strain rate gets its extremum equal to v p /l s at the beginning of the impact. However, this equation is only valid within the characteristic time [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF].

Since the force equilibrium of the specimen can not be checked, and the actual loading period is much longer than the characteristic time, an inverse analysis is carried out to identify the dynamic behaviour of the Ti-6Al-4V alloy.

INVERSE IDENTIFICATION OF THE JOHNSON-COOK CONSTITUTIVE MODEL

A constitutive model is first postulated to describe the behaviour of the tested material. The parameters of this constitutive model are then identified so that some given quantities extracted from the numerical simulation fit experimental data. The Johnson-Cook model is widely adopted to describe the dynamic behaviour of metallic materials, and is here used for the Ti-6Al-4V alloy. The temperature effects are not addressed in this work, the Johnson-Cook equation thus reads :

σ(ε p eq , ε p eq ) = A + B(ε p eq ) n 1 + C ln εp eq ε0 (4)
The parameters A, B and n are identified on quasi-static stress-strain curves obtained at the strain rate of 10 -4 s -1 : A=955MPa, B=770MPa and n=0.557. A reasonable range of search of 0.005 ≤ C ≤ 0.05 is set for the parameter C. Then a 2D axisymmetric finite element model is established in ABAQUS to run the inverse analysis in order to identify the parameter C. The cost function f (x) is defined as the euclidean norm of the difference between the simulated strain ε sim (x, t) extracted at the location of the gauge and the recorded one ε exp (t).

f (x) = ε sim (x, t) -ε exp (t) 2 (5)
The figure 2 presents the superposed plots of the recorded and calibrated strain curves of some tests, and also shows the values of C identified on each test with respect to the experimental strain rates assessed by equation 3. In the figures 2(a)-2(c), the calibrated strain curves fit well the recorded ones in the loading stage. However some discrepancies can also be observed during the unloading period. In the figure 2(d), slight discrepancies appear in the loading stage, when the strain rates exceed 10 4 s -1 . This kind of discrepancies increases as the strain rate goes beyond 20000s -1 as shown in the figure 2(e). In the figure 2(f), the parameter C is almost constant in the range of 3000 to 15000s -1 . As the strain rate increases, the results of the identification are subject to more scattering. A much greater value of C of about 0.03 is obtained when the strain rate is over 20000s -1 , whereas a smaller one is identified on the tests T7 and T9. For many ductile materials, it has been observed that the stress increases more rapidly as the strain rates excess the level of 10 3 s -1 [START_REF] Rule | A revised form for the Johnson-Cook strength model[END_REF]. However here, more tests and identifications are required to come to a reliable conclusion for the highest strain rates.

CONCLUSION

In this work, the Ti-6Al-4V has been tested on the direct-impact Hopkinson device on a wide range of strain rate ranging from 3000 to 25000s -1 . Due to the mismatch of the material impedances and crosssections, a longer loading period has been observed in the force curves. Thus an inverse identification procedure has been adopted to identify the Johnson-Cook model for the Ti-6Al-4V on the experimental data. It has been shown that this model works well at least up to 20000s -1 , even though more tests are required for strain rates greater than this value.

FIGURE 2 -

 2 FIGURE 2 -Recorded and calibrated strain curves

TABLE 1 -

 1 Experimental plan and results

	Test Obj. ε n o (s -1 ) (m/s) (mm) (mm) (s -1 ) v p l p l s εmax ε max Test Obj. ε n o (s -1 ) (m/s) (mm) (mm) (s -1 ) v p l p l s εmax ε max
	T1	3000	15	500	5	3064 11%	T7 18000	36	60	2 17970 33%
	T2	5000	25	125	5	4642 16%	T8 18000	30	60	1.5 18350 16%
	T3	7000	28	125	4	6925 17%	T9 20000	30	60	1.5 21222 23%
	T4 10000	30	60	3 10740 17% T10 20000	40	30	2 19659 24%
	T5 12000	36	60	3 12040 25% T11 25000	37.5	30	1.5 25050 25%
	T6 15000	30	60	2 15010 11%


Xiaoli Guo

Cheikh-Tidiane Sow
email: cheikh-tidiane.sow@ec-nantes.fr

Chady Khalil
email: chady.khalil@ec-nantes.fr

Thomas Heuzé

Guillaume Racineux

Material constitutive behavior identification at high strain rates using a direct-impact Hopkinson device

Keywords: Direct impact Kolsky bar, inverse analysis, high strain rate, elastic-viscoplastic, constitutive models

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

To further develop electromagnetic or electrohydraulic forming, crimping or welding processes, accurate simulation of theses process is needed. This goes by the study of the material performance at high strain rates and to define a precise dynamic behavior by using constitutive models. The strain rates involved can be up to 10 5 𝑠 -1 . It is well known that the behavior of many metallic materials usually depends on the strain rate; therefore, a good representation of this behavior is required by the mean of dynamic tests on a wide range of strain rates. The most popular experiment is the classical Split Hopkinson Pressure Bar (SHPB) [START_REF] Andriamiseza | Rapport de stage de Master: Identification de modèles de comportement à haute vitesse de déformation sur des essais de Kolsky à impact direct[END_REF][START_REF] Kolsky | An investigation of the mechanical properties of materials at very high rates of loading[END_REF] capable to attain strain rates ranging up to 10 4 s -1 [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF][START_REF] Ramesh | Part D: Chapter 33. High strain rate and impact experiments[END_REF]. For higher strain rates, another system of the SHPB is available: a direct-impact Hopkinson device [START_REF] Dharan | Determination of stress-Strain characteristics at very high strain rates[END_REF]. In this device, the incident bar is removed and very high strain rates of the order of 10 5 𝑠 -1 can be achieved [START_REF] Dharan | Determination of stress-Strain characteristics at very high strain rates[END_REF][START_REF] Kamler | Measurement of the behaviour of high-purity copper at very high rates of strain[END_REF]. In the present work, the Ti-6Al-4V alloy is tested at strain rates varying from 3000 to 25000 𝑠 -1 using a direct-impact Hopkinson device and an elastic-viscoplastic constitutive model is identified by using an inverse identification.

Direct-impact tests on the Ti-6Al-4V at very high strain rate

This device is composed of a transmitted bar, a projectile and three Wheatstone bridges of double gauges mounted on the bar for the measurement of deformations and dispersion monitoring (Figure 1). The transmitted bar is 1.2m length and 10 mm diameter. Four projectiles of diameter 15.8 mm and lengths of 500 mm, 125 mm, 60 mm and 30 mm respectively are used.

Figure 1: direct impact device

The higher the objective strain rate is, the lower the length of the projectile must be to avoid crushing the specimen. The bar and the projectiles are made of MARVAL X2NiCoMo18-8-5, with a yield stress of 1800 MPa. The impact velocity of the projectile is measured by two laser diodes. For each objective strain rate, the value of the specimen length 𝑙 𝑠 , the projectile length 𝑙 𝑝 and its impact velocity 𝑣 𝑝 should be determined. The maximum values of the expected strain rate 𝜀 𝑠 ̇ and the allowable strain 𝜀 𝑠 in the specimen can be approximated by the following equation:

𝜀̇𝑠 𝑚𝑎𝑥 ≈ 𝑣 𝑝 𝑙 𝑠 𝜀 𝑠𝑚𝑎𝑥 ≈ 𝜀 𝑠 ̇2𝑙 𝑝 𝑐 𝑝 (1)
where 𝑐 𝑝 denotes the sound speed in the projectile. 𝑐 𝑝 is defined by the equation / E  . E and  are respectively the Young's modulus and the density of the projectile material.

The parameters for each objective strain rate (Obj. S  ) are listed in Table 1. The elastic strain 𝜀 𝑏 (𝑡) in the transmitted bar and the impact velocity 𝑣 𝑝 of the projectile are measured. From the recorded strain 𝜀 𝑏 (𝑡), the force 𝐹 𝑠 (𝑡) applied by the projectile on the specimen is computed by:

𝐹 𝑠 (𝑡) = 𝐸 𝑏 𝑆 𝑏 𝜀 𝑏 (𝑡) (2)
where 𝐸 𝑏 is the Young's modulus and 𝑆 𝑏 is the cross-section of the bar. The results are plotted in Figure 2. The trays observed during the unloading part of these curves arise from the shorter length of the projectile than that of the bar and the mismatch of the generalized wave impedances between the projectile 𝑍 𝑝 = (𝜌𝑐𝑆) 𝑝 , the specimen 𝑍 𝑠 = (𝜌𝑐𝑆) 𝑠 and the bar 𝑍 𝑏 = (𝜌𝑐𝑆) 𝑏 [START_REF] Wang | Foundations of stress waves[END_REF]. When the projectile is long enough, the loading time and the characteristic time are equivalents, as we can see in the Figure 2(a) where the projectile is 500 mm length. As the projectile shorten the actual period of loading of the specimen lasts approximately twice the characteristic time or even longer as is observed in Figure 2(b) and Figure 2(d), these results from the different cross-sections of the projectile and the bar. Following this, the specimen undergoes a second compression before the unloading starts. Thus, it is subjected to a higher strain than that of the test T1 as shown in the Table 1, even though the projectile of the test T1 has a greater kinetic energy before impact.

For the direct impact device, the strain rate in the specimen is assessed by the following equation [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF][START_REF] Guo | Inverse identification at very high strain rate of the Johnson--Cook constitutive model on the Ti-6Al-4V alloy with a specially designed direct-impact Kolsky bar device[END_REF]:

𝜀̇𝑠(𝑡) = - 𝑣 𝑝 + 𝑆 𝑝 +𝑆 𝑏 𝑆 𝑝 𝐶 𝑏 𝜀 𝑏 (𝑡) 𝑙 𝑠 (3)
where 𝑆 𝑝 is the cross-section of the projectile and 𝐶 𝑏 the sound speed of the bar. The maximum compressive strain rate 𝜀̇𝑠 𝑚𝑎𝑥 computed by using this equation and the measured compressive strain 𝜀 𝑠𝑚𝑎𝑥 for each test are listed in Table 1. This equation states that the compressive strain rate gets its maximum equal to 𝑣 𝑝 /𝑙𝑠 at the beginning of the impact and is however only valid during the characteristic time [START_REF] Gorham | An improved method for compressive stress-strain measurements at very high strain rates[END_REF]). The stress-strain curves are computed using the classical analysis of the Hopkinson tests [START_REF] Guo | On the directimpact Hopkinson system for dynamic tests at very high strain rates[END_REF]. The results are shown in the Figure 3. As the strain rate increases, the curves exhibit a rising strain hardening. However, they are subjected to scattering, as we can observe through small oscillations. Consequently, the rate dependency is difficult to be assessed on these curves directly. An inverse analysis is carried out for the identification of the constitutive models.

Inverse identification of elastic-viscoplastic constitutive models

The elastic-viscoplastic constitutive models used to represent the behavior of the material must be able to fit correctly the experimental measurements. They are identified by an inverse procedure that uses a dynamic numerical analysis. No thermal effect is addressed in these identifications. Three models that fit well to this requirement have been chosen:

a) The Johnson-Cook model [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF])

𝜎(𝑝, 𝑝) = (𝐴 + 𝐵𝑝 𝑛 )(1 + 𝐶𝑙𝑛 𝑝ṗ ̇0) (4)
where p and 𝑝̇ are the cumulated plastic strain and strain rate respectively. The parameters A, B and n are obtained by quasi-static classical tension test [10]: A=955 MPa, B=770 MPa and n=0.557. The parameter C is the only one to be identified by the inverse process. An axisymmetric 2D Finite Element model in ABAQUS is used to compute the dynamic simulation. The model consists of the projectile, the specimen and the bar.

b) The Zerilli-Armstrong model [START_REF] Zerilli | Constitutive relations for titanium and Ti-6Al-4V[END_REF] This model is calibrated using a unidimensional FE model in MATLAB developed in [START_REF] Andriamiseza | Rapport de stage de Master: Identification de modèles de comportement à haute vitesse de déformation sur des essais de Kolsky à impact direct[END_REF]. Assuming a constant temperature T, two variants for bcc and hcp materials are considered in Eq. 5 and Eq. 6 respectively:

𝜎 𝑦 = 𝐾 0 + 𝐾 1 𝑝 𝑛 + 𝐾 2 𝑝̇𝐾 3 (5)
where 𝐾 0 = 𝐶 0 , 𝐾 1 = 𝐶 5 , 𝐾 2 = 𝐶 1 exp(-𝐶 3 𝑇) and 𝐾 3 = 𝐶 4 𝑇. The parameters, 𝐾 0 , 𝐾 1 and 𝑛 have the same value as the parameters A, B and n of the Johnson-Cook relation because of the equivalent form of strain hardening. The parameters 𝐾 2 and 𝐾 3 have been calibrated inversely.

𝜎 𝑦 = 𝐾 0 + 𝐾 1 𝑝̇𝐾 2 + 𝐾 3 √ 𝑝𝑝̇𝐾 4 (ℎ𝑐𝑝) (6)
where 𝐾 0 = 𝐶 0 , 𝐾 1 = 𝐶 1 exp(-𝐶 3 𝑇) , 𝐾 2 = 𝐶 4 𝑇, 𝐾 3 = 𝐶 2 exp(-𝐶 3 ′ 𝑇) and 𝐾 4 = 𝐶 4 ′ 𝑇. 𝐾 0 is equivalent to the yield stress. Four parameters 𝐾 1 , 𝐾 2 , 𝐾 3 and 𝐾 4 are computed during the inverse identification.

The cost function 𝑓(𝑥) is defined as the Euclidean norm of the difference between the simulated strain 𝜀 𝑠𝑖𝑚 (𝑥, 𝑡) and the recorded strain 𝜀 exp (𝑡). The identification is done on each tests of the Table 1. The calibrated strain 𝜀 𝑠𝑖𝑚 of the tests T1, T4, T6, T8 and T11 are plotted and superposed to the measured strain 𝜀 𝑒𝑥𝑝 in the Figure 4(a)-Figure 4(e). An error bar is plotted associated to each mean value of the parameter C for the 2D identification, as shown in the Figure 4(f). At least two repetitions are carried out for each test. The simulated strains of the three different constitutive models fit well to the experimental data at strain rate up to 15000 𝑠 -1 . We observe at the second raising stage when the specimen is plastically deformed that the results are concordant as shown in the Figure 4(a)-Figure 4(c). A small difference appears at this stage between the calibrated strain of the Johnson-Cook model and the measured one in the Figure 4(d). The Zerilli-Armstrong model fits well the experimental data. The Figure 4(f) indicates that an approximately constant value of the parameter C of the Johnson-Cook model is identified at the strain rate ranging from 3000 to 18000 𝑠 -1 . At very high strain rate as observed in the tests T10 and T11, a much greater value of C is obtained. The differences in these two tests may be due to (i) the very short projectile used in the test and/or (ii) the thermal effects not addressed in the identification. More repetitions of these tests must be carried out to clarify the results.

Conclusion

In this work, the Ti-6Al-4V has been tested on a range of strain rate ranging from 3000 to 25000 𝑠 -1 using the direct-impact Hopkinson device. To identify three elastic-viscoplastic constitutive models, an inverse analysis has been carried out on the experimental data. The three identified constitutive models fit well the experimental data at the strain rate ranging from 3000 up to 18000 𝑠 -1 . However more experimental and numerical research at higher strain rate is required to refine these results.

Figure 2 :Figure 3 :

 23 Figure 2: Force on specimen in the direct impact tests

  Error bars of parameters

Figure 4 :

 4 Figure 4: Results of identifications

Figure 5 :

 5 Figure 5: stress-strain curves plotted within the range of 10 -4 to 18000 s -1

Table 1 :

 1 experimental plan and results

	Test	Obj. 𝜀 𝑠 ̇	𝑣 𝑝	𝑙 𝑝	𝑙 𝑠	𝜀̇𝑠 𝑚𝑎𝑥	𝜀 𝑠𝑚𝑎𝑥	Test	Obj. 𝜀 𝑠 ̇	𝑣 𝑝	𝑙 𝑝	𝑙 𝑠	𝜀̇𝑠 𝑚𝑎𝑥	𝜀 𝑠𝑚𝑎𝑥
	n°	(s -1 )	(m/s)	(mm)	(mm)	(s -1 )		n°	(s -1 )	(m/s)	(mm)	(mm)	(s -1 )
	T1	3000	15	500	5	3064	0.11	T7	18000	36	60	2 17970	0.33
	T2	5000	25	125	5	4612	0.16	T8	18000	30	60	1.5 18350	0.16
	T3	7000	28	125	4	6925	0.17	T9	20000	30	60	1.5 21222	0.23
	T4	10000	30	60	3 10740	0.17	T10	20000	40	30	2 19659	0.24
	T5	12000	36	60	3 12040	0.25	T11	25000	37.5	30	1.5 25050	0.25
	T6	15000	30	60	2 15010	0.11