Miguel Aranda

Juan Antonio Corrales

Youcef Mezouar

Deformation-based shape control with a multirobot system

We present a novel method to control the relative positions of the members of a robotic team. The application scenario we consider is the cooperative manipulation of a deformable object in 2D space. A typical goal in this kind of scenario is to minimize the deformation of the object with respect to a desired state. Our contribution, then, is to use a global measure of deformation directly in the feedback loop. In particular, the robot motions are based on the descent along the gradient of a metric that expresses the difference between the team's current configuration and its desired shape. Crucially, the resulting multirobot controller has a simple expression and is inexpensive to compute, and the approach lends itself to analysis of both the transient and asymptotic dynamics of the system. This analysis reveals a number of properties that are interesting for a manipulation task: fundamental geometric parameters of the team (size, orientation, centroid, and distances between robots) can be suitably steered or bounded. We describe different policies within the proposed deformationbased control framework that produce useful team behaviors. We illustrate the methodology with computer simulations.

I. INTRODUCTION

Compared to a single robot, a multirobot team can manipulate larger and heavier objects and provide more refined and precise behaviors. The team's actions when manipulating a rigid or deformable body must be carefully coordinated, to avoid damaging it and encountering unpredictable motions [START_REF] Sieber | Iterative optimal feedback control design under relaxed rigidity constraints for multi-robot cooperative manipulation[END_REF]- [START_REF] Umlauft | Dynamic movement primitives for cooperative manipulation and synchronized motions[END_REF]. In particular, the robots typically need to ensure that the object remains close to a desired state. In this paper, we consider a scenario where a team of robots grasps rigidly a deformable object in 2D. We assume that the desired state of the object is encapsulated by a desired shape of the team. This assumption is reasonable in different practical cases, e.g.: 1) The polygon joining the grasping points represents the object's contour faithfully enough. 2) The object is highly deformable and thus its shape adapts to the team's shape. 3) The task does not require precise object shape control but, e.g., making the object fit in a given region of space. 4) The required object deformation is small. Therefore, we propose a method to control the shape of the robotic team. The approach minimizes a measure of the team's deformation relative to its desired shape. In this way, the controller aligns with the usual goal of keeping the manipulated object close to a desired state. Our measure of deformation is the function that one minimizes when solving

The authors are with Université Clermont Auvergne, CNRS, SIGMA Clermont, Institut Pascal, F-63000 Clermont-Ferrand, France. {miguel.aranda, juan-antonio.corrales-ramon, youcef.mezouar}@sigma-clermont.fr

This work was supported by project CoMManDIA (SOE2/P1/F0638) which is cofinanced by Interreg Sudoe Programme (European Regional Development Fund); and by the French government IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25) via project MaRoC. a certain Procrustes shape alignment problem [START_REF] Gower | Procrustes problems[END_REF]. In 2D space, the particular problem where rotation and scaling are optimized admits a closed-form solution. Our insight and contribution here is to appropriately formalize this solution and to define a gradient control strategy based on it. The controller is simple to compute, and amenable to analysis. Notably, we identify interesting properties for a manipulation task: the dynamics of the team's size, centroid and orientation can be controlled suitably both in short-term and asymptotic time scales. The robots' motions are tightly coordinatedas is typically desired in manipulation tasks-, balanced and efficient, and the distances between them can be bounded (to avoid overstretching, over-squeezing or collisions). We analyze these points and illustrate them via simulations.

A. Related work

Precise knowledge of object deformation is not a fundamental consideration for some multirobot manipulation approaches which deal with objects that are either essentially rigid [START_REF] Sieber | Iterative optimal feedback control design under relaxed rigidity constraints for multi-robot cooperative manipulation[END_REF], [START_REF] Umlauft | Dynamic movement primitives for cooperative manipulation and synchronized motions[END_REF], [START_REF] Bai | Cooperative load transport: A formation-control perspective[END_REF]- [START_REF] Petitti | Decentralized motion control for cooperative manipulation with a team of networked mobile manipulators[END_REF] or, on the other hand, highly deformable. In the latter case, it may be sufficient to simply respect upper and lower inter-robot distance bounds [START_REF] Tanner | Nonholonomic navigation and control of cooperating mobile manipulators[END_REF], [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF]. To ensure that the relative states of the robots maintain desired values, these approaches often exploit distributed techniques such as consensus or formation control [START_REF] Oh | A survey of multi-agent formation control[END_REF], which are very successful in many scenarios. Typically, they possess strong asymptotic convergence properties, but shortterm team dynamics are difficult to predict. Contrary to this literature, our method is based on a concept of global shape of the team. This creates closely coordinated motions, which is a singularly relevant property when the robots are manipulating a deformable object.

The robotic manipulation of deformable objects is a highly active area of research where a variety of methodologies, tasks, robotic platforms and object types have been considered [START_REF] Sanchez | Robotic manipulation and sensing of deformable objects in domestic and industrial applications: a survey[END_REF]. To keep the object's geometry under control, existing works typically exploit a certain degree of knowledge about how it deforms under the action of the robots. Some authors have used a prior model of the object's deformation [START_REF] Das | Autonomous shape control of a deformable object by multiple manipulators[END_REF], [START_REF] Jiménez | Survey on model-based manipulation planning of deformable objects[END_REF]. The models depend greatly on the type of object; linear objects, for instance, admit representations suitable for shape planning [START_REF] Bretl | Quasi-static manipulation of a Kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], [START_REF] Mukadam | Quasi-static manipulation of a planar elastic rod using multiple robotic grippers[END_REF]. Other authors favor modelfree approaches [START_REF] Berenson | Manipulation of deformable objects without modeling and simulating deformation[END_REF], [START_REF] Delgado | Tactile control based on Gaussian images and its application in bi-manual manipulation of deformable objects[END_REF], in some cases using a deformation Jacobian estimated online from sensory data [START_REF] Navarro-Alarcón | Modelfree visually servoed deformation control of elastic objects by robot manipulators[END_REF]- [START_REF] Hu | Three-dimensional deformable object manipulation using fast online Gaussian process regression[END_REF]. For feedback control, methods not reliant on prior models are considered more generalizable and computationally simpler. The methods cited in the paragraph have these general traits: 1) When multiple robots are used, the robots share a common global goal (e.g., object shape control), but they do not coordinate their motions with one another. In contrast, our method emphasizes the coordination/synchronization of the robots' actions. This can lead to higher efficiency in the team's motions; 2) The target configuration the robots pursue (e.g., corresponding to a desired object shape) is anchored, to the workspace or to a reference frame. In contrast, in our approach the target configuration is the one that minimizes, for the current configuration of the robots, the team's deformation measure. The motions are therefore not linked to external references and, being based purely upon deformation, they produce more efficient changes of shape. In addition, a free-floating target configuration -which our method provides-is adaptable and flexible for object transportation tasks. Furthermore, maintaining a minimum deformation can minimize damages to the object, and to the robots, because the contact forces will be reduced.

Finally, Procrustes-based deformation metrics similar to the one we exploit -which appear frequently in shape alignment problems [START_REF] Kanatani | Analysis of 3-D rotation fitting[END_REF]-have been used previously to control mobile robot formations [START_REF] Aranda | Coordinate-free formation stabilization based on relative position measurements[END_REF], estimate their error [START_REF] Ze-Su | Formation control and obstacle avoidance for multiple robots subject to wheel-slip[END_REF], or determine optimal a priori destination formations [START_REF] Mosteo | Optimal role and position assignment in multi-robot freely reachable formations[END_REF]. To our knowledge, this is the first work to study the application of such metrics to cooperative manipulation tasks.

II. GENERAL PROBLEM STATEMENT

Consider a set, N = {1, 2, ..., N }, N > 1, of robots with the ability to move as kinematic point masses satisfying:

qi = u i , ∀i ∈ N (1)
(boldface font is used for multidimensional variables), where q i ∈ R 2 denotes robot i ′ s position in a given reference frame and u i ∈ R 2 its control input. We will in general not notate time dependence. {•} will stand for {• | i ∈ N }. The norm used in this paper is the Euclidean one. Our focus is on a scenario where the robots manipulate an object on a plane, each of them grasping it rigidly at a fixed contact point on the object's contour. Similar setups with free-floating rigidgrasping robotic actuators are assumed in, e.g., [START_REF] Berenson | Manipulation of deformable objects without modeling and simulating deformation[END_REF], [START_REF] Das | Autonomous shape control of a deformable object by multiple manipulators[END_REF], [START_REF] Navarro-Alarcón | Fourier-based shape servoing: A new feedback method to actively deform soft objects into desired 2-D image contours[END_REF]. We use the following notation for the relative position vectors: q ij = q iq j . We also work with a prescription of a desired geometry for the robotic team. This desired geometry is considered to be constant over time. It can be encapsulated by a layout of points c i ∈ R 2 in the workspace: {c i }. We assume that c i = c j for at least two distinct robots i, j ∈ N . We define the centroid of the current robot positions as q o , and the centroid of the desired geometry as c o .

We assume that the state of the object is defined by the configuration of the robotic team. Hence, we focus on controlling the team. We do not control the state of the object's full contour. This assumption is justified in the first paragraph of the paper. Still, one cannot ignore that there is an object being handled by the robots: this demands the team's motions to be closely coordinated and stable, and to satisfy size and distance bounds. Our deformation-based control approach allows to conform with these requirements.

We address the problem of shape control, defined as making {q i } satisfy q i = sRc i + t ∀i ∈ N for some scaling s ∈ R ≥0 , rotation R ∈ SO(2), and translation t ∈ R 2 . This condition implies that the team has the same shape as the desired geometry. Next, we describe the proposed solution.

III. CONTROL STRATEGY

Let us consider the following function defined from the two given point sets {q i } and {c i }:

γ g = 1 2 i∈N ||q i -H g c i + t g || 2 , (2)
where H g ∈ R 2×2 is a Euclidean similarity transformation -comprising rotation and uniform scaling-, and t g ∈ R 2 is a translation offset. γ g is a measure of the SSD (Sum of Squared Differences) error between the points in the sets when the points in {c i } are rotated and scaled by H g , and translated by t g . Finding the optimal H g and t g (i.e., the similarity and offset that minimize γ g for two given sets) is a type of Procrustes shape-fitting problem [START_REF] Gower | Procrustes problems[END_REF]. One can see that the optimal offset corresponds with making the centroid of the two point sets coincide. Therefore, we henceforth use q io and c io in (2), and choose t g = 0.

Let us define a similarity H ∈ R 2×2 having the form:

H = h 1 -h 2 h 2 h 1 (3)
h 1 = 1 c s i∈N q T io c io h 2 = 1 c s i∈N q T io c ⊥ io ,
where c s = i∈N ||c io || 2 , i.e., a strictly positive constant scalar, and the operator ⊥ on a 2D vector a is defined as a counterclockwise rotation of π/2 radians: a ⊥ = [(0, 1) T , (-1, 0) T ]a. It can be verified with simple manipulations that ∂γg ∂Hg(1,1) ⌋ Hg(1,1)=h1 = 0 and ∂γg ∂Hg(2,1) ⌋ Hg(2,1)=h2 = 0 and, via the second derivative test, that h 1 and h 2 are unique minimizers of γ g . Thus, the optimal choice in ( 2) is H g = H.

Taking the above into consideration, our control approach rests on the following cost function:

γ = 1 2 i∈N ||q io -Hc io || 2 . ( 4
)
Let us highlight that h 1 and h 2 satisfy ∂γ ∂h1 = 0 and ∂γ ∂h2 = 0. A visual interpretation of the cost function in the context of the shape control problem addressed is given in Fig. 1. We discuss important aspects about the chosen function next. 1) A crucial fact to notice is that γ can be interpreted directly as a measure of the deformation relative to the desired shape of the robotic team. H is such that {Hc io + q o } is the closest (in terms of the SSD in γ g ) possible team configuration to {q i } having the same shape as {c i }. The SSD for this optimal configuration is γ, which therefore quantifies purely how much the shape of {c i } is deformed in the current configuration {q i }. Then, our controller, based on the gradient of γ, directly minimizes the deformation with b) Current robot positions {q i } (hollow circles) and desired geometry in a) transformed by the optimal similarity H (3). H rotates and uniformly scales every c io . Therefore, the point set {Hc io } has the same shape as {c io }. Notice that γ = (1/2) i∈N ||d i || 2 . d i is robot i's motion vector, opposite to the gradient of γ with respect to q i , as explained in the text. respect to the desired shape. This is particularly appropriate for a task of manipulation of a deformable object. 2) γ is very suitable for our purpose: it is differentiable and simple computationally -which is important for feedback control-. It expresses unequivocally the achievement of the prescribed shape, and its gradients have closed-form expressions. Fundamental parameters (size, rotation and centroid of the team) can be steered and controlled, and exhibit dynamic behaviors that are interesting for manipulation tasks.

3) Since we use a similarity (H), the desired physical size of the team is unprescribed, which may appear suboptimal or even risky. Various arguments support our choice: having the size parameter as a degree of freedom provides adaptability to changing task requirements or environment conditions, allowing to take advantage of the deformability of soft objects. Also, the similarity (3) is linear in the relative position coordinates, and this makes it easier to uncover its interesting dynamic properties. Finally, as we will show, the value of the team's size can still be controlled.

Let us define the scaling, s h , and rotation, R h , encapsulated by H. The control strategy does not use these parameters, but they provide us with a physical interpretation of the action of the matrix. We can write H = s h R h , with:

s h = h 1 2 + h 2 2 , (5)
R h = cos(α h ) -sin(α h ) sin(α h ) cos(α h ) , α h = atan2(h 2 , h 1 ), (6)
being by our definition α h = 0 when h 1 = h 2 = 0. Notice that there is no discontinuity in this alternative representation of H, because when h 1 = h 2 = 0, it holds that s h = 0. The size of the team is an essential property to take into account. We define it as the quadratic mean of the robots' distances to the centroid. The size of the configuration towards which the robots move at a given instant (which can be termed destination configuration) is thus defined as:

s d = 1 N i∈N ||Hc io || 2 = s h c s N . (7)
Observe that, as c s and N are constant, s d and s h are essentially equivalent ways of expressing this size. Similarly, we define the size of the current configuration as:

s q = 1 N i∈N ||q io || 2 . ( 8
)
We can already formulate a result that supports γ as being an appropriate cost function for our control purposes.

Lemma 1: {q i } is equal to {c i } up to translation, rotation and scaling if and only if γ = 0.

One can prove this lemma by considering the form of (3). We compute next the gradient of γ with respect to the position of one of the robots i ∈ N . Recall that the matrix H satisfies ∂γ ∂h1 = 0 and ∂γ ∂h2 = 0. Therefore, the gradient will be given directly by the partial derivative with respect to q i . In particular, it has the following expression:

∇ q i γ = ∂γ ∂q i = j∈N (q jo -Hc jo ) T ∂q jo ∂q i = (1 - 1 N )(q io -Hc io ) - 1 N j∈N j =i q jo -Hc jo = q io -Hc io , (9)
where we have separated in two different terms the case j = i and the case j = i, and used the fact that the centroids satisfy:

q o = i∈N q i /N , c o = i∈N c i /N .

A. Standard control law

For each robot, we define its standard control law in the direction opposite to the gradient of γ. The motion vector for robot i takes, from [START_REF] Petitti | Decentralized motion control for cooperative manipulation with a team of networked mobile manipulators[END_REF], the form that follows:

d i = -∇ q i γ = q oi -Hc oi . ( 10
)
This vector can be directly visually interpreted in Fig. 1. We propose the following control law for each robot i ∈ N :

u i = k i d i , (11)
where k i ≥ k min (k min being a positive constant) is a gain, constant or smoothly time-varying. The information needed to implement this control law is q oi ∀i ∈ N , which can be obtained from measurements of the relative positions of the team's robots. It can be seen (we omit the details for brevity) that this controller can be computed in a local reference frame. Note that H in ( 10) is a feedback matrix, recomputed at every instant according to [START_REF] Umlauft | Dynamic movement primitives for cooperative manipulation and synchronized motions[END_REF]. The controller has a closedform expression and, being based on minimizing a globally defined metric, it needs knowledge about all the robots.

IV. CONTROLLER ANALYSIS

Let us inspect how the multirobot system behaves under the controller. We first state a necessary assumption. Assumption 1: s q remains upper-bounded for all time. We now study the stability with respect to the desired shape.

Theorem 1: Under the action of the control law [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF], the robots converge exponentially to a configuration in which they form the desired shape and stay static.

Proof: As s q is upper-bounded under the controller (see Assumption 1), the relative positions of the robots remain finite [START_REF] Wang | Kinematic multi-robot manipulation with no communication using force feedback[END_REF]. We choose γ as candidate Lyapunov function, and analyze its dynamics when the controller [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF] is in operation:

γ = i∈N (∇ q i γ) T qi = - i∈N k i (t)||d i || 2 ≤ -k min i∈N ||d i || 2 = -2k min γ. (12)
From LaSalle's invariance principle, one can see that the system converges exponentially to a configuration in which γ = 0, i.e., where the robots are in the desired shape (Lemma 1). In addition, γ = 0 implies d i = 0 and thus, u i = 0 ∀i ∈ N , i.e., the robots' velocities converge to zero. Remark 1: The exponential decay of γ (12) implies that the robot speeds vanish exponentially. Since these speeds are thus integrable as time goes to infinity, the inter-robot distances will remain finite. We support Assumption 1 with this observation and with the analysis of the dynamics of s q provided in the following sections.

A. Properties of team size, rotation and centroid

It is essential to have some insight on how the team size is expected to evolve. Let us state a relevant result:

Lemma 2: The cost function γ (4) can be expressed as:

γ = N 2 (s q 2 -s d 2 ). (13)
Proof: The proof relies on algebraic manipulations of (4). We omit the details, for brevity.

There is an interesting fact to notice here: γ expressesas highlighted previously-the difference in shape; but, as [START_REF] Sanchez | Robotic manipulation and sensing of deformable objects in domestic and industrial applications: a survey[END_REF] shows, it additionally encapsulates the difference in size between the current and destination configurations. By minimizing γ, one is controlling size and shape simultaneously. We can state the following property regarding team sizes:

Proposition 1: It holds that s d ≤ s q for any current and desired configurations of the robots.

Proof: As γ ≥ 0, (13) proves the result. This means that the controller naturally moves the robots towards a configuration in which the size of the team is never larger than its current size. This suggests an interesting trend towards maintaining the energy of the system and avoiding divergent behaviors of the team. On a more practically specific note, this property can also be useful so as to prevent object overstretching. A result on the dynamics of H follows.

Proposition 2: If all robots apply the control law [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF] with the same gain k c , then H remains constant over time.

Proof: We examine the time-derivative of H:

ḣ1 = i∈N ∂h 1 ∂q i T qi . (14)
We can find that ∂h1 ∂q i = c io /c s , and substitute [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF] to get:

ḣ1 = k c c s i∈N c T io (q oi -Hc oi ) = - k c c s ( i∈N q T io c io - i∈N c T io Hc io ). (15)
Using that H is a similarity and substituting its expression (3), the last term of ( 15) can be readily seen to be:

i∈N c T io Hc io = i∈N ( ||c io || 2 c s j∈N q T jo c jo ) = i∈N q T io c io .
(16) Hence, ḣ1 = 0. An analogous analysis can be used to show that ḣ2 = 0 and, therefore, H remains constant. Thus, both s h and R h remain constant. This fact can actually facilitate the motion of the robots during the task, and help them attain efficiently the desired shape while avoiding unnecessary rotations/re-scalings of the object.

When all robots employ the same gain, the team's centroid is preserved under the action of the controller. One can verify this readily from [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF] and realizing that by definition the sums of all N vectors q oi or c oi are zero. Thus:

qo = 1 N i∈N qi = k c N ( i∈N q oi -H i∈N c oi ) = 0. ( 17
)
This fact is useful for navigation/transportation problems where we want to steer the centroid of the team: this centroid motion objective can be decoupled from the motions of the shape controller. Another important conclusion is that the changes in the manipulated object's shape are balanced, because they do not displace the center of the contact points.

B. Inter-robot distance bounds

Upper and lower inter-robot distance bounds are of great importance in practice because they can ensure that overstretching and over-squeezing/collisions, respectively, are avoided. To study this issue, let us first state a useful result.

Lemma 3: The following function:

γ p = 1 2 i∈N j∈N ||q ij -Hc ij || 2 (18)
satisfies γ p = 2N γ. Proof:

Defining v i = q io -Hc io , we have that γ = 1 2 i∈N ||v i || 2 and γ p = 1 2 i∈N j∈N ||v ij || 2 = 2N γ -v ij , where v ij = i∈N j∈N v T i v j = i∈N v T i j∈N v j .
Notice that by definition i∈N v i = 0. Thus, v ij = 0 and the stated result follows. Unlike γ, the function γ p is defined in terms of the pairwise inter-robot vectors. As both functions are equal up to a constant multiplicative factor, we can use γ to define bounds for the pairwise distances. Let us define the variable e(t) = 2 √ N γ(t). We can state the property that follows. Proposition 3: The distance at time t between any two robots i ∈ N , j ∈ N is bounded as follows:

s h (t)||c ij || -e(t) ≤ ||q ij (t)|| ≤ s h (t)||c ij || + e(

t). (19)

Proof: One can substitute (4), [START_REF] Delgado | Tactile control based on Gaussian images and its application in bi-manual manipulation of deformable objects[END_REF], and use the triangle inequality to obtain the result. Thus, by ensuring that γ and s h maintain suitable values, one has a means to keep the inter-robot distances bounded.

V. VARIANTS OF THE CONTROLLER

Modifications of the controller [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF] are presented next. These variants retain the underlying idea of maintaining the deformation under control, by minimizing γ. At the same time, they provide different dynamic behaviors, or allow to fix important geometric parameters, accommodating constraints that are relevant in practical applications.

A. Team motions in the rotation and scaling nullspaces

We next describe alternative control laws that preserve either the rotation or the scaling encapsulated in H. They are found by exploiting the nullspaces of these parameters.

1) Rotation-preserving motions: We define them as those for which Ṙh = 0. We can provide the following result.

Proposition 4: A motion strategy where the robots implement the control law that follows:

u i = k c (q oi -w i (t)Hc oi ), w i (t) ∈ R >0 , ∀i ∈ N , (20)
is rotation-preserving.

Proof: Let us express (20) for each robot as u i = u s i + u n i , where u s i is the standard control law [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF]. Then:

u n i = k c w * i (t)Hc oi , ∀i ∈ N , (21)
with w * i (t) = 1w i (t) a scalar. We will compute next the dynamics of h 1 and h 2 . From Proposition 2 we know that both terms are constant under u s i . Thus, we have, for h 1 :

ḣ1 = i∈N ∂h 1 ∂q i T u n i = k c c s i∈N c T io w * i (t)Hc oi = k c c s w(t)h 1 , (22) with w
(t) = -i∈N w * i (t)||c io || 2 . Analogously, for h 2 : ḣ2 = i∈N ∂h 2 ∂q i T u n i = k c c s i∈N c ⊥ io T w * i (t)Hc oi = k c c s w(t)h 2 .
(23) We assume the angle α h (6) is always well-defined, which it will be as long as the scaling satisfies s h > 0. We can express its dynamics as αh = (h 1 ḣ2 -ḣ1 h 2 )/s 2 h . Hence, substituting ( 22) and ( 23), it is concluded that αh = 0. This is interesting because if the object does not rotate, the manipulating robots' motions are more efficient and balanced. Also, a robot could use the degree of freedom given by w i (t) to satisfy other task criteria (i.e., avoid obstacles, or steer the normal of the object's contour).

2) Scale-preserving motions: These are team motions for which ṡh = 0. We find the result that follows.

Proposition 5: The following team motion policy:

u i = k c ((q oi -Hc oi ) + z i (t)Hc ⊥ oi ), z i (t) ∈ R, ∀i ∈ N , (24)
is scale-preserving.

Proof: Let us write [START_REF] Ze-Su | Formation control and obstacle avoidance for multiple robots subject to wheel-slip[END_REF] for each robot as u i = u s i +u m i , where u s i is the standard control law [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF]. Then:

u m i = k c z i (t)Hc ⊥ oi , ∀i ∈ N . (25)
We find the dynamics of h 1 and h 2 . The overall contribution of u s i is zero (Proposition 2), so we can write, for h 1 :

ḣ1 = i∈N ∂h 1 ∂q i T u m i = k c c s i∈N c T io z i (t)Hc ⊥ oi = - k c c s z(t)h 2 , (26) with z(t) = -i∈N z i (t)||c io || 2 . Analogously, for h 2 : ḣ2 = i∈N ∂h 2 ∂q i T u m i = k c c s i∈N c ⊥ io T z i (t)Hc ⊥ oi = k c c s z(t)h 1 .
(27) The dynamics of s 2 h (5) are d(s 2 h )/dt = 2(h 1 ḣ1 + h 2 ḣ2 ). Substituting ( 26) and ( 27), one finds that ṡh = 0, i.e., the scaling parameter remains constant. These motions result, as one could expect, in behaviors for which the team and the grasped object gyrate. This can be interesting if the object needs to be rotated rigidly in the workspace due to task demands (i.e., for a subsequent treatment, or to place it in a prescribed absolute orientation).

B. Fixing the team's size and orientation using a single robot

A single special robot can steer the matrix H for the full team. Without loss of generality, let us assume that robot 1 is the special robot and c 1o = 0. Consider the control:

u 1 = k c d 1 , with d 1 = q o1 -H 1 c o1 , (28)
where H 1 is a similarity. The following result holds: Proposition 6: If all robots use the same control gain k c , robots i = 2, ..., N implement the standard control law [START_REF] Alonso-Mora | Local motion planning for collaborative multi-robot manipulation of deformable objects[END_REF], and robot i = 1 uses (28), then H tracks H 1 , according to:

Ḣ = -k f 1 (H -H 1 ), with k f 1 = k c • ||c 1o || 2 /
= k c c s c T 1o H dif c o1 , ḣ2 = k c c s c ⊥ T 1o H dif c o1 . (29)
Note that H dif is a similarity. Using this constraint on its structure and developing the equations above, one can express the components of H back in matrix form and obtain:

Ḣ = -k f 1 H dif = -k f 1 (H -H 1 ). (30)
Thus, if H 1 is constant, then H converges to it exponentially. If H 1 is time-varying, then H tracks its variation. The special robot can in this way steer the value of H. It can also control the team's geometric parameters (size and orientation) individually. Two possible strategies based on the controllers in Section V-A are described next. S1: Drive the team size s h towards a desired value s h = s 1 , without changing R h : Robots use [START_REF] Navarro-Alarcón | Fourier-based shape servoing: A new feedback method to actively deform soft objects into desired 2-D image contours[END_REF], with w 1 (t) = s 1 /s h (s h > 0 is assumed) and w i (t) = 1 for i = 2, ..., N . S2: Rotate team towards a desired orientation R h = R(α 1 ), without changing s h : Robots use [START_REF] Ze-Su | Formation control and obstacle avoidance for multiple robots subject to wheel-slip[END_REF], with z 1 (t) = k r1 (α 1α h ), k r1 < 0, and z i (t) = 0 for i = 2, ..., N .

Notice that S1 and S2 are particular instances of the case studied in Prop. 6: robot 1 is enforcing a similarity H 1 , which may be time-varying, while the other robots' motions minimize the deformation metric γ. As H will track H 1 , robot 1 can suitably steer the team's parameters.

The special robot may be a robot with superior perceptual/decision making capabilities, e.g., global positioning and environmental awareness to manipulate/transport the object avoiding obstacles; or a privileged agent whom we want to lead the task (e.g., as in human-robot collaboration). The other robots, meanwhile, can have more modest capabilities and rely on local reference frames and simpler sensors.

VI. SIMULATION STUDY

We present illustrative results of simulations, performed using MATLAB . For our first example, we consider an inextensible elastic rod which rests in a planar workspace and is grasped rigidly at different points by four robots. Using gripping points along the rod -i.e., not just at its ends-makes it easier to steer this type of object around environmental obstacles and to control its shape precisely. We employ a quasi-static model of the rod's state based on minimal-energy configurations, as in previous works, e.g., [START_REF] Bretl | Quasi-static manipulation of a Kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], [START_REF] Mukadam | Quasi-static manipulation of a planar elastic rod using multiple robotic grippers[END_REF].

We control the orientation of the grippers/robots with a simple strategy aimed at avoiding sharp curvatures of the rod. We use for θ i ∈ (-π, π] -the orientation of robot i, equal to the tangent angle of the rod at that robot's position-:

θi = k ai (θ d i -θ i ), (31)
with k ai > 0 a control gain, θ d i = ∠(q i,i-1 + q i+1,i ) for i = 2, ..., N -1, θ d 1 = ∠(q 21 ), θ d N = ∠(q N,N-1 ). Here, ∠(•) denotes the angular polar coordinate of a vector and • the normalization to unit. Figure 2 illustrates the test's results.

Initially, the rod was fairly stretched out. The prescribed shape we chose was close to that of a square. The controllers used were those described in Section V-B. Robot 1 (leftmost at the start of execution) acted as the special robot. Three control phases were executed. First, robot 1 specified a team size (strategy S1). Then, it commanded a reduction of size with no rotation (S1). Finally, it made the team rotate without changing size (S2). The team parameters evolved as theoretically expected. The robot motions were efficient and the changes in the object's shape were smooth.

Our second simulation example, aimed at testing the robustness of the controller, included additive Gaussian noise both in the measurement and the actuation for all robots, and a systematic error added to the velocity of one robot. Five robots with a regular pentagonal desired shape were used. No object was considered in the simulation. Robots followed S1 so as to fix a desired team size s 1 = s d = 1. As is observable in Fig. 3, the controller was able to keep the team close to the desired shape despite the perturbations. The deformation metric γ was maintained in low values. A complementary video provides further illustration of the simulation results.

VII. CONCLUSION

We believe that the main appeal of the proposed method is that it applies multirobot formalisms to the control of metrics of deformation. This links it with object manipulation tasks and opens the door to multiple potentially interesting applications. However, several essential issues were not considered. We do not address the interaction with the object; for precise and safe operation, information of force and object state should be included in the control loop. Also, if the desired shape is very different from the initial one, it would be necessary to use additional planning. Finally, robot dynamics constraints, control in 3D space, and collision avoidance guarantees can be studied building on the presented results.

Fig. 1 .

 1 Fig. 1. Illustration of control strategy with four robots. Only the vectors for one robot (i = 1) are shown. The vectors for the other robots are analogous. The centroids are marked as squares. a) Desired geometry {c i }.b) Current robot positions {q i } (hollow circles) and desired geometry in a) transformed by the optimal similarity H (3). H rotates and uniformly scales every c io . Therefore, the point set {Hc io } has the same shape as {c io }. Notice that γ = (1/2) i∈N ||d i || 2 . d i is robot i's motion vector, opposite to the gradient of γ with respect to q i , as explained in the text.

Fig. 2 .

 2 Fig. 2. Simulation results for linear object. Top: snapshots of object and robots (circles) evolution from initial (left) to final (right) configuration. Bottom, left to right: desired geometry (robots as circles); time evolution of: team sizes and similarity angle, robot speeds, cost function, and inter-robot distances.

Fig. 3 .

 3 Fig.3. Simulation example with five robots. Left to right: initial (robots joined by dashed lines) and final (solid lines) shapes, with one same robot marked with a square in both shapes; evolution of similarity angle (thickest line) and team sizes (thinner lines, with sq > s d ); and of cost function.

  c s . Proof: We look at the time derivative of H. Let us define H dif = H -H 1 and u 1 = u s

	1 + u d 1 , where u s 1 obeys 1 = k c H dif c o1 . From Prop. 2, the dynamics of (11) and u d H is due only to u d 1 . Direct manipulations then lead to:
	ḣ1


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

NUMERICAL SIMULATION OF IMPACTS ON ELASTIC-VISCOPLASTIC SOLIDS WITH THE FLUX-DIFFERENCE SPLITTING FINITE VOLUME METHOD

Keywords: Flux difference splitting finite volume method, Elastic-viscoplastic solids, Dynamic impact, Curvilinear structured meshes

INTRODUCTION

The numerical simulation of hyperbolic initial boundary value problems including extreme loading conditions such as impacts requires the ability to accurately capture and track the wave front of shock waves induced in the medium. Indeed, this permits to correctly follow the path of waves and hence understand the mechanical phenomena occuring within that medium. For solid-type media, it allows also for an accurate assessment of the plastic strain field and hence that of residual stresses and distortions within the structure. High speed forming processes like electromagnetic material forming are some application examples of severe loading conditions in which the track of wave fronts is important both for understanding the development of plastic strain in the workpiece and optimizing its final shape. Hence, these problems require numerical schemes able to meet high orders of accuracy and a high resolution of discontinuity without any spurious oscillations.

The numerical simulation of impacts on dissipative solids has been and is again mainly performed with the classical finite element method coupled with centered differences or Newmark finite difference schemes in time [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. Though this approach has well-known advantages, classical time integrators introduce high frequency noise in the vicinity of discontinuities which is hard to remove with artificial viscosity without destroying the accuracy of the numerical solution. The finite volume method, initially developed for the simulation of gas dynamics [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], has gained recently more and more interest for problems involving impacts on solid media [START_REF] Barton | A high order eulerian godunov method for elastic-plastic flow in solids[END_REF][START_REF] Ortega | Numerical simulation of elastic-plastic solid mechanics using an eulerian stretch tensor approach and HLLD Riemann solver[END_REF][START_REF] Favrie | Mathematical and numerical model for nonlinear viscoplasticity[END_REF][START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: applications to dynamic fracture and fragmentation[END_REF]. This method shows some advantages to achieve an accurate tracking of wavefronts; among others (i) continuity of fields is not enforced on the mesh in its cell-centered version, that allows for capturing discontinuous solutions, (ii) the characteristic structure of hyperbolic equations can be introduced within the numerical solution, either through the explicit solution of a Riemann problem at cell interfaces, or in a implicit way through the construction of the numerical scheme, (iii) the amount of numerical viscosity introduced can be controlled locally as a function of the local regularity of the solution, so that to permit the elimination of spurious numerical oscillations while preserving a high order of accuracy in more regular zones.

We are here interested in elastic-viscoplastic systems, leading to a non-homogeneous system of partial differential equations, generating a system of weakly discontinuous waves beyond the viscoplastic yield, following a discontinuous (elastic) wave due to the transition between elastic and elastic-viscoplastic ranges. In this work, we use and extend a finite volume method applied to these dissipative media with bidimensional curvilinear structured meshes of quadrangles. The formulation is based on the flux-difference splitting method [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], has second order spatial accuracy through flux limiters, embeds the corner transport upwind method, and uses the second order accurate Strang splitting method to handle the right hand side of the system of balance laws. The approach is here derived with a Chaboche-type [START_REF] Lemaitre | Mechanics of solid materials[END_REF] elastic-viscoplastic solid within the small strain framework, and is illustrated on a problem of impact on a heterogeneous volume containing an inclusion. A comparison is performed with a finite element solution obtained with the finite element code Cast3M [START_REF]cast3M-cea[END_REF].

ELASTIC-VISCOPLASTIC INITIAL BOUNDARY VALUE PROBLEM

Let's consider an elastic-viscoplastic solid with an Odqvist creep law [START_REF] Lemaitre | Mechanics of solid materials[END_REF] (or Chaboche-type model) and a linear kinematic hardening, the associated constitutive equations written within the small strain framework are summarized below:

ε = εe + εp (1) σ = C : εe (2) f = σ eq (ξ) -σ y ; σ eq = 3 2 ξ : ξ ; ξ = s -X ; C = {(σ, X, εp )|f ≤ 0} (3) εp = 3 2 ṗ ξ σ eq (4) Ẋ = 2 3 D ˙ε p (5) ṗ = f K n ( 6
)
where equations ( 1), ( 2), ( 3), ( 4), ( 5) and ( 6) respectively refer to the additive partition of the strain rate ( ε) into elastic ( εe ) and viscoplastic ( εp ) parts, the elastic law, the yield criterion f (or viscous stress) associated to Mises' norm and the elastic convex C , the viscoplastic flow rule, the evolution law of the kinematic hardening and the creep law. One defines ξ = s -X, where s is the deviatoric part of Cauchy stresses σ, X the variable defining the center of the elastic convex, x = (x + |x|)/2 the positive part of x, C the elastic stiffness tensor, ṗ the rate of cumulated viscoplastic strain, and D the hardening parameter. To these constitutive equations, are added the equations of motion and the relations between the linearized strain rate tensor ε and the velocity field. One can show that we obtain the following non-homogeneous system of partial differential equations, written here for the bidimensional plane strain case:

∂u ∂t + ∂f x ∂x + ∂f y ∂y = S(u) ∀(x × t) ∈ (Ω×]0, t end ]) (7)
defined within a domain Ω and a time interval of duration t end , and has to be supplemented with appropriate initial and boundary conditions. The vector of unknowns u, fluxes f x and f y , and the right hand side S(u) are defined as follows:

u =       σ xx σ yy σ xy u v       ; f x =       -(λ + 2µ)u -λu -µv -σ xx /ρ -σ xy /ρ       ; f y =       -λv -(λ + 2µ)v -µu -σ xy /ρ -σ yy /ρ       ; S(u) = -3 µ σ eq f K n       ξ xx ξ yy ξ xy 0 0       (8)
where λ and µ denote the Lamé's parameters, u and v the velocity components, and ρ the mass density. The eigen-analysis of the jacobian matrix associated to the homogeneous part of ( 7)

J = ∂(F • n) ∂u = ∂f x ∂u n x + ∂f y ∂u n y = An x +Bn y = -       0 0 0 n x (λ + 2µ) λn y 0 0 0 n x λ n y (λ + 2µ) 0 0 0 n y µ n x µ nx ρ 0 ny ρ 0 0 0 ny ρ nx ρ 0 0       (9) gives the characteristic speeds λ 1 = -c L , λ 2 = c L , λ 3 = -c s , λ 4 = c s , λ 5 = 0, where c L =
(λ + 2µ)/ρ and c s = µ/ρ are the celerities associated to tensile/compression and shear waves respectively; and associated characteristic directions (eigenvectors)

K (p) , 1 ≤ p ≤ 5.
The finite volume method is based on subdividing the computational domain in elementary cells within which, for the cell-centered version, an approximation U i of the field u is defined in the cell i by integral averaging. This approximation is updated by means of balance laws [START_REF] Favrie | Mathematical and numerical model for nonlinear viscoplasticity[END_REF] written in integral form on each grid cell, coupled to an explicit time integration scheme.

Strang splitting for the right hand side

The solution procedure to account for the right hand side of the differential system (7) can be carried out by means of a splitting algorithm [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. For first order accuracy (Godunov splitting), a first step consists in solving the homogeneous part of ( 7) on the time step associated to the initial condition known at time t n (n referring to the time step number). The obtained solution ūn+1 is then used as an initial condition for the following initial value problem:

du dt = S(u) u(0) = ūn+1 (10)
This decomposition is exact for a linear advection-reaction equation. For the nonlinear system considered here, it is an approximation. Naming C et S operators associated to the solution of homogeneous and non-homogeneous parts, the update of the unknown field reads

u n+1 = S (∆t) C (∆t) u n (11)
In order to improve the accuracy of the computed viscoplastic strains, the second order accurate Strang splitting algorithm is used here:

u n+1 = S (∆t/2) C (∆t) S (∆t/2) u n (12)

Flux-difference splitting formulation

We use here a finite volume formulation based on flux-difference splitting [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], applied to curvilinear structured meshes of quadrangles. This version falls within the class of total variation non-increasing methods, that allow to meet both a high order of accuracy in zones where the solution field is regular and a high resolution of discontinuities without spurious numerical oscillations. Their strength relies on their ability to introduce a controlled amount of numerical viscosity locally, so that to adapt to the local regularity of the solution. This formulation is applied to the homogeneous part of the differential system (7) so that to be in accordance with a non-increasing total variation.

Let's consider the grid cell i shown in figure 1, of area |A i |, which edge s (1 ≤ s ≤ 4) joining points P s and P s+1 is of length L s , has an outward unit normal n s , and has left (L) and right (R) states known at time t n . The update of the cell state at time t n+1 reads:

U n+1 i = U n i - ∆t |A i | P k=1 L k A + k ∆U k + Q l=1 L l A - l ∆U l - ∆t |A i | P k=1 L k Fout k - Q l=1 L l Fin l (13)
The above formula consists of two groups of terms which are detailed below.

Godunov flux

The first group of terms in equation ( 13) shows fluctuations, denoted by operators A ± k ∆U k :

A ± k ∆U k = Mw p=1 λ ± p W (p) k = Mw p=1 λ ± p α (p) k K (p) k (14)
quantifying the effect of all waves travelling rightward or leftward respectively in the local frame of edge k. Theses fluctuations are defined with positive and negative parts of characteristic speeds λ ± p , and associated characteristic directions

K (p) k ≡ K (p) (n k ). Each wave is weighted with a coefficient α (p)
k determined by the projection of the jump of the unknown vector

∆U k = (U R -U L ) k across the edge k onto the characteristic basis K (p) k (1 ≤ p ≤ M w , M w = 5 here): ∆U k = Mw p=1 W (p) k = Mw p=1 α (p) k K (p) k = Kα k ( 15
)
The wavestrength is defined by

W (p) k = α (p) k K (p)
k . This projection amounts to solve a (linear) Riemann problem associated to edge k. The fluctuations associated to each edge are summed to compute the contribution of first order terms to the update of grid cell i; this summation is performed on negative fluctuations for edges having an outward normal, and on positive fluctuations for edges having an inward normal, so that P + Q = N , N being the number of edges of grid cell i. The formulation with fluctuations is actually equivalent [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] to applying the Godunov's method which interface fluxes are computed with the steady part (x/t = 0) of the Riemann problem solution.

Additional fluxes

The second group of terms in equation (13) involve additional fluxes defined at edges, either inward (in) or outward (out) depending on the edge normal, that consist of two types of contributions

Fin l = FHO l + Ftran l . (16)

High order fluxes

The first contribution allows to reach a higher order (HO) of accuracy (order two here), defined with wavestrength W (p) l that have been limited, hence denoted

W(p) l = α(p) l K (p) l : FHO l = 1 2 Mw p=1 |λ (p) l | 1 - ∆t ∆s l |λ (p) l | W(p) l (17)
where ∆s l refers to the distance between barycenters of grid cells sharing edge l. Waves are limited based on an upwind ratio θ (p) l defined for wave p on edge l as:

θ (p) l = W (p) J • W (p) l W (p) l 2 (18)
where J denotes the upwind edge, that is the opposed edge to l belonging to grid cell L located to the left in the local frame of edge l (see figure 1) if λ (p) l > 0, or the opposed edge belonging to grid cell R located to the right if λ (p) l < 0. However, for noncartesian quadrangles upwind and downwind edges do not necessarily have the same normal. Thus, the computation of the upwind ratio (18) has to be performed with wavestrengths computed in a same local reference frame, that of edge l. Wavestrengths W (p) J are hence computed with weighting coefficients α (p) J recomputed in the local frame of edge l:

α (p) J = K -1 (n l ) • ∆U J (19)
where ∆U J denotes the jump across edge J of the unknown vector. The upwind ratio (18) can be understood as a certain measure of the local regularity of the solution, the wavestrength of wave p associated to edge l is then limited by means of a limiting function φ(θ

(p) l ): α(p) l = φ(θ (p) l )α (p) l (20)
Many limiting functions exist and permit to obtain different known finite volume schemes [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]. Some of them enable the numerical scheme to satisfy a non-increasing total variation, so that the appearance of spurious numerical oscillations can be avoided in the vicinity of discontinuities. The limiter Superbee defined by φ(θ) = max(0, min(1, 2θ), min(2, θ)) falls in this family, and is used in this work. More generally, the limitation of wavestrength amounts to add locally some numerical viscosity, and to locally lower the order of accuracy to properly capture discontinuities. In zones where the solution field is more regular, the limitation is not active and an accuracy of order two can be reached.

Transverse fluxes

The second contribution to additional fluxes (16) enables to improve the stability of the numerical scheme, so that the Courant number can be set at one. These fluxes allows to account for informations travelling in bias with respect to the considered grid cell; this is the contribution to the grid cell to be updated of a cell only sharing a node (but not an edge) with it. This method is called the Corner Transport Upwind (CTU) method [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. These fluxes are of great importance to ensure numerical stability for elastic media; indeed elasticity couples strain components through Poisson's effect, so that a transverse information to the considered grid cell should be introduced in the numerical scheme.

Let's consider the patch of grid cells shown in figure 2. One focuses on the edge denoted (i) which local frame (n i , t i ) is shown. This edge gives rise to the computation of normal fluctuations A + i ∆U and A - i ∆U contributing to grid cells R et L respectively. These normal fluctuations leads to the computation of transverse fluctuations giving contribution to neighboring cells across edges (j) and (k) for cell L, and across edges (m) and (l) for cell R. These transverse fluctuations are computed by projecting normal fluctuations on the characteristic basis associated to the Riemann problem defined on the adjacent edge; it appears as a transverse Riemann solver. For instance the negative normal fluctuation is decomposed on the characteristic basis associated to edge (j) as

x y n i A + i ∆U A - i ∆U L R T (i) (j) (k) (l) (m)
B + j A - i ∆U B + k A - i ∆U B + m A + i ∆U B + l A + i ∆U n j n k n m n l
A - i ∆U = Mw p=1 β p K (p) j = K j β (21)
where K j accounts here for the normal n j of edge (j), but also of different material properties between grid cells L and T . Coefficients β p are determined analytically for bidimensional plane strain elasticity equations. For outward normals shown in figure 2, the transverse fluctuations are computed with the positive operator B + , that is only waves with positive characteristic speeds will contribute to this transverse fluctuation

B + j A - i ∆U = c L W (2) + c s W (4) (22)
An additional numerical flux defined at edges is hence built from these transverse fluctuations:

Ftran j = ∆t 2∆s j B + j A - i ∆U i (23)
which contributes to grid cell T . The flux Ftran l associated to edge l appearing in equation ( 16) thus denotes transverse contributions from adjacent grid cells to cell i. More generally, this finite volume scheme is linked to a Taylor expansion of the solution in the vicinity of a grid cell. In this framework, Godunov fluxes are first order terms, high order fluxes are second order terms, and transverse fluxes correspond to cross-derivatives.

SUDDEN VELOCITY LOADING AND UNLOADING OF A HETEROGENEOUS VOLUME

Let's consider a square heterogeneous volume, of side length 2a, with an inclusion of circular cross-section of radius R centered within the volume, which mesh is shown in figure 3. This volume in an initial natural state is suddenly loaded on its left side at time t = 0 with a constant first component of velocity ū. After time τ , the applied velocity is set to zero. Normal velocity components at the top and bottom sides of the volume are set to zero while transmittive boundary conditions have been set on the right side. We assume an arbitrary heterogeneous material, that consists of an inclusion in aluminium embedded in a matrix made of steel. The analysis is carried out in plane strain. Numerical data used for computations are summarized in table 1. The matrix is assumed to remain elastic, only the inclusion undergoes viscoplastic strains. For comparison purpose, a finite element solution is obtained with the finite element code Cast3M [START_REF]cast3M-cea[END_REF], computed with an implicit time discretization and bilinear quadrangular finite elements. This numerical solution is compared to the finite volume one in figures 4, 5, and 6 at times 1 × 10 -7 , 1.9 × 10 -7 et 3 × 10 -7 seconds respectively. The Courant number has been set at 0.9, the same time step being used for both numerical solutions. The comparison is performed on to the finite element solution. One can observe in figure 4 the slot of negative normal stress formed by the prescribed horizontal velocity on the left side. Its falling front has interacted with the front interface between the matrix and the inclusion, and has been both reflected and transmitted, generating an intermediate state of higher stress (lower in absolute value) and a higher first component of velocity since the impedance of the matrix is higher than that of the inclusion. The finite element stress field shows spurious numerical oscillations in the vicinity of discontinuities, especially close to the top and bottom sides where tensile spurious stress states appear. Cumulated viscoplastic strain are almost identical for both numerical solutions.

In figure 5, one can see that the reflected wave has been transformed into tensile wave after the interaction with the increasing stress front associated to velocity unloading; then has been reflected on the left side of the volume. The cumulated viscoplastic strains propagation is delayed with respect to the propagation of the negative stress slot, which is due to viscoplasticity. At last, one can see in figure 6 that this slot has interacted with the back side of the inclusion interface, generating intermediate states of normal compression stress, while cumulated viscoplastic strains grow more importantly close to this area. Generally speaking, the finite volume solution allows to obtain the same viscoplastic strains than these of the finite element solution without the spurious numerical oscillations on stresses obtained with the finite element solution. Figure 7 shows the comparison of the shear component of stress σ 12 at time t = 3 × 10 -7 s, on which these oscillations are visible.

CONCLUSIONS

It has been presented in this work the use and extension of a finite volume method for elasticviscoplastic solids with curvilinear structured meshes of quadrangles. The formulation is based on the flux-difference splitting method, and uses the second order accurate Strang splitting approach for the computation of the right hand side of the differential system of equations. Comparison is performed with a finite element solution obtained with the finite element code Cast3M [START_REF]cast3M-cea[END_REF] on a suddenly loaded then unloaded heterogeneous volume. The comparison shows identical cumulated viscoplastic strains computed. However, the finite volume solution enables to overcome spurious numerical oscillations on the stress and velocity fields thanks to flux limiters, applied on the homogeneous part of [START_REF] Favrie | Mathematical and numerical model for nonlinear viscoplasticity[END_REF], satisfying a non-increasing total variation.

Figure 1 :

 1 Figure 1: Quadrangular finite volume

Figure 2 :

 2 Figure 2: Normal and transverse fluctuations defined from edge i.

Figure 3

 3 Figure 3: Mesh

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Time t = 1 × 10 -7 s

Figure 7 :

 7 Figure 7: Shear stress component σ 12 at time t = 3 × 10 -7 s

Table 1 :

 1 Numerical data used for computations.