I Mancisidor
email: imancisidor@ideko.es

M Zatarain

J Muñoa

Z Dombovari

Receptance Coupling for Tool Point Dynamics Prediction

In many applications, chatter free machining is limited by the flexibility of the tool. Estimation of that capacity requires obtaining the dynamic transfer function at the tool tip. Experimental calculation of that Frequency Response Function (FRF) is a time consuming process, because it must be done using an impact test for any combination of tool, toolholder and machine. The bibliography proposes the receptance coupling substructure analysis (RCSA) to reduce the number of experimental test. A new approach consisting of calculating the fixed boundary dynamic behavior of the tool is proposed in the paper. This way the number of modes that have to be considered is low, just one or two for each bending plane, and it supposes an important improvement in the application of the RCSA to the calculation of stability diagrams. The predictions of this new method have been verified experimentally.

INTRODUCTION

In recent years, high-speed machining (HSM) has been widely applied in aerospace industry, due to the good machinability provided by aluminum alloys and the improvements in machine and spindle designs. Under these conditions, the major limitation in the material removal rate are selfexcited vibrations. These vibrations, also known as chatter, cause a reduction in the surface quality and in the lifetime of mechanical elements and tools.

Since several years, many researchers have investigated methods to avoid the regenerative chatter [START_REF] Bediaga | Chatter avoidance method for milling process based on sinusoidal spindle speed variation method[END_REF][START_REF] Bilbao-Guillerna | Control laws for chatter suppression in milling using an inertial actuator[END_REF]. One of the most popular methods is to create a stability lobe diagram in order to determine the best cutting conditions [START_REF] Altintas | Analytical Prediction of Stability Lobes in Milling[END_REF]. The stability lobes separate the stable and unstable regions depending on the spindle speed and depths of cuts. These diagrams are built by a stability model [START_REF] Zatarain | Estudio comparativo de los modelos matemáticos de chatter en fresado: monofrecuencia, multifrecuencia y simulación en el tiempo[END_REF][START_REF] Munoa | Stability study of the milling process using an exponential force model in frequency domain[END_REF] fed with four inputs: tool geometry, cutting coefficients, cutting conditions and dynamic parameters.

To obtain good dynamic parameters, accurate measurements of FRF at the tool tip are required. The problem is that the measurements have to be repeated for each combination of tool/holder/spindle, which are time consuming and disturb the production.

In order to reduce the duration of the measurements and to increase the efficiency of the process, Schmitz et al. [START_REF] Schmitz | Predicting high-speed machining dynamics by substructure analysis[END_REF][START_REF] Schmitz | Tool Point Frequency Response Prediction for High-Speed Machining by RCSA[END_REF][START_REF] Schmitz | Three-Component Receptance Coupling Substructure Analysis for Tool Point Dynamics Prediction[END_REF] proposed a receptance coupling technique to predict the dynamic response at the tip of the tool. This method allows coupling of theoretical and/or experimental frequency responses of individual components to obtain the general response of the assembly.

In the simplest case, Schmitz et al. [START_REF] Schmitz | Predicting high-speed machining dynamics by substructure analysis[END_REF] considered only translational degree of freedom at the spindle portion of the substructure and introduced a interface flexibility between the tool and the spindletoolholder. The method was improved introducing the rotational degree of freedom with a torsional joint flexibility value [START_REF] Schmitz | Tool Point Frequency Response Prediction for High-Speed Machining by RCSA[END_REF]. The joint parameters were identified performing a correlation between theoretical and experimental values and updating the joint flexibility values to provide the best fitting. Park et al. [START_REF] Park | Receptance coupling for end mills[END_REF] proved that the rotational displacement of the tool at the joint cannot be neglected for accurate construction of FRF at the tool tip, and they proposed a new methodology to identify the receptance at the tool holder nose using blanks with different lengths.

The major advantage of RCSA is that it could save considerable time by means of the calculation of the FRF in the tool tip for different tools, combining the theoretical response of the different tools and the experimental results of the toolholder/machine assembly.

There are some ways to obtain the theoretical FRFs of the tools. Schmitz et al. [START_REF] Schmitz | Predicting high-speed machining dynamics by substructure analysis[END_REF][START_REF] Schmitz | Tool Point Frequency Response Prediction for High-Speed Machining by RCSA[END_REF] used analytical Euler-Bernoulli beam theory, which is simple and provides reasonable engineering approximations for many problems. However, this model tends to overestimate slightly the natural frequencies, especially for high frequency modes [START_REF] Han | Dynamics of Transversely Vibrating Beams using Four Engineering Theories[END_REF]. Also, the prediction is better for slender beams than non-slender beams.

Another classic beam theory was proposed by Timoshenko to overcome these inaccuracies, which adds the effect of shear as well as the effect of rotation to the Euler-Bernoulli beam. However, the Timoshenko beam model introduces a cut-off frequency that sometimes can be a low frequency. The results obtained above this frequency cannot be taken into account, and sometimes an insufficient number of modes is available.

Other way to obtain theoretical responses is by Finite Element Method (FEM). It applies a beam theory, normally Timoshenko's model, and the system is discretized in several elements. This method was used by Park et al. [START_REF] Park | Receptance coupling for end mills[END_REF] and Liljehren et al. [START_REF] Liljerehn | Dynamic substructuring in metal cutting machine chatter[END_REF]. The tools can be me meshed using 3D elements too.

All these methods are proposed in the bibliography with 'free-free' boundary conditions. Therefore, the modes are calculated considering the tool without any kind of of restraint or connection. This way, it requires the calculation of a high number of modes (natural frequencies) to obtain an accurate representation of the dynamics of the tool, because the mathematical description of strain distribution near the joint requires high frequency modes [START_REF] Schmitz | Predicting high-speed machining dynamics by substructure analysis[END_REF].

In this paper, this drawback is overcome proposing the fixed boundary approach. This model calculates a free-free beam response based on clamped-free boundary conditions beam results. Hence, the number of modes that have to be considered is very low and the Timoshenko beams can be used.

RECEPTANCE COUPLING

In RCSA, experimental or analytical FRFs of the individual components can be coupled to predict the final dynamic response of the assembly at any selected spatial coordinate. The objective of this paper is to couple an analytical model of the tool with an experimentally identified dynamic model of the toolholder-spindle assembly as illustrated in Figure 1. The RCSA will be explained briefly. According to Park et al. [START_REF] Park | Receptance coupling for end mills[END_REF], a rigid joint has been considered. Mathematically, the G ij matrix describes the frequency response function of the final assembly measured in point i and the excitation in point j (see equation ( 1)).

                   j j j i jy i iyj iyjy i i M f G G G G y      (1)
The Point FRF in the tool tip G 11 is one of the main inputs to feed the stability model [START_REF] Altintas | Analytical Prediction of Stability Lobes in Milling[END_REF]. By means of RCSA, this main receptance is estimated using receptances measured or calculated when the two substructures are disconnected (H ij ). Operating according to RCSA the next basic expressions can be obtained:

              21 1 33 22 12 11 11 H H H H H G       (2)               22 1 33 22 12 12 12 H H H H H G       (3)
where:

                          j i jy i iyj iyjy ij j i jy i iyj iyjy ij h h h h H ; g g g g G
The response in the substructure A will be obtained theoretically with free-free boundary conditions. The direct FRF in the substructure B (h 3y3y ), can be obtained experimentally by measuring the spindle with a short blank, but the other two receptances (h 3y3θ , h 3θ3θ ) related to the rotational degree of freedom, are difficult to measure directly. However, following the methodology proposed by Park et al. [START_REF] Park | Receptance coupling for end mills[END_REF], this functions can be obtained using a long blank, measuring the direct response in the tool tip and finally subtracting the effect of the blank.

To simplify the formulation, the names of variables of equations 2 and 3 have been changed, obtaining the equations 4 and5

                         3 3 2 2 3 3 2 2 3 3 2 2 2 2 2 2 2 2 1 2 2 1 1 1 2 1 1 1     g b e c eg b e a u               2 (4)     g b e d eg b f b v               2 (5)
Two unknowns (β and δ) are going to be considered, and they can be solved by the utilization of the symbolic non-linear analytical toolbox:

   
                                2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2  (7)
After this extraction, h 3y3θ and h 3θ3θ receptances can be obtained:

        2 2 3 3 2 2 3 3 h h h h y y     (8)
Furthermore, the damping of the joint interface has been taken into account with this methodology.

FREE-FREE TOOL FRFS

In the receptance coupling, FRFs of substructure A (the tool) in free-free conditions are required. An analytical formulation is used to calculate this response, and predict the receptance of the complete system before the mounting of the tool. The first part of this section shows the development of this new model, while the second verifies theoretically the improvement introduced.

Mathematical development of the Fixed Boundaries Model

Before applying this approach, it is necessary to calculate the modes of the beam with clamped-free boundary conditions. A good strain distributions near the joint is provided, but it is not the required free-free beam response for RCSA. The next step is to add the effect of movements restrained in the clamped-free model (see Figure 2). Mathematically, the equation ( 9) represents the general equation of vibration, where a combination of rigid body and modal coordinates is going to be introduced [START_REF] Han | Dynamics of Transversely Vibrating Beams using Four Engineering Theories[END_REF]:

        f x M x K     (9)                                      n n n l l l k k k q q q q x   2 1 2 2 2 1 2 1 1 2 1           (10)
According with the last equation [START_REF] Han | Dynamics of Transversely Vibrating Beams using Four Engineering Theories[END_REF], the first term describes the movement of the joint and the second refers to clamped-free modes. n is the number of modes considered. This change drives to the next equation [START_REF] Liljerehn | Dynamic substructuring in metal cutting machine chatter[END_REF]:

        f q M q K T T T T T T                                                 (11)
If each part of this equation is analyzed a simple form of the stiffness matrix is obtained:

                                               2 2 2 2 2 1 diag 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 f n T T T T T T K K K K K                            (12)
However, the mass matrix presents more complexity in some of its terms:

                            M M M M M T T T T T T (13) Operating:                                                      1 0 0 0 1 0 0 0 1 0 0 2 0 1 0 0 2 0 1          nxn T L n L L L n L L T tot g tot g tot tot T I M dm x dm x dm x dm dm dm ' M M I x m x m m ' M M             
where x g is the distance from the gravity center of the beam and Φ are the mode shapes along x. Finally the equation ( 13) can be written as:

                nxn T T T I ' M ' M ' M M       (14)
If in the general equation ( 11) a sinusoidal force with a certain frequency () is considered, the next expression is obtained:

    f T ' M ' M ' M q T f T 1 2 2 2 2 2                           (15) where                     n n y y y n n y y y L T                     1 2 1 1 1 1 2 1 1 1 2 2 2 1 2 2 2 2 1 2 1 0 1 1 0 0 1    
Finally, the damping can be added introducing a loss factor, and obtaining the final receptance ready to introduce in the RCSA:

        T f T T i ' M ' M ' M T H 1 2 2 2 2 2 1                      (16)
Taking into account the numeration defined in Figure 2: This way, the free-free state response is obtained describing the correct strain distribution in the joint with a low number of modes.

                                1 1 1 1 2 1 2 1 1 1 1 1 2 1 2 1 1 2 1 2 2 2 2 2 1 2 1 2 2 2 2 2

h h h h h h h h h h h h

Mathematical model of the tool.

The previous solution has been combined with Timoshenko beam theory to calculate analytically the receptance of the tool.

The geometry of the tool should be modeled properly considering the complexity of the cross-section in his fluted part. Furthermore, the tool has been divided in two beams related to the fluted part and the pure cylindrical side. This way, different beam properties could be introduced along the different segments of the model.

The calculation of the inertia of the fluted part is an important issue to obtain accurate predictions. In this paper an approach proposed by Kivank et al. [13] has been used. They presented cross-section models of 2, 3 and 4 fluted mills. Each cross-section was divided in regions which were described by three variables: the radius r of the arc, the position of the center of the arc (a) and the diameter fd of the flute.

Theoretical verification

To verify theoretically the new fixed boundary model, a simple case has been selected to compare the result with the usual free boundary model. The results of both approaches to perform the RCSA have been correlated with the real theoretical response of the assembly. The Timoshenko beam formulation has been used in all the calculations of this comparison.

The substructure A which simulates the tool, is a simple 20mm diameter and 160mm length beam. The rigid substructure (B) is a 300mm diameter and 500mm length beam. The density and the modulus of two beams are 7820 kg/m3 and 207 Gpa, respectively.

In the substructure A, the cut-off frequency is reached analytically for the 12 th mode, therefore only 11 mode can be introduced using pure Timoshenko beams. This drawback can overcome using a FEM based model [START_REF] Munoa | Stability study of the milling process using an exponential force model in frequency domain[END_REF]. Figure 4 shows the effect of the number of considered modes in the free-free response for A substructure calculated by free boundaries formulation. Furthermore, in this case it is not possible to create an accurate result. However, when the new approach is used (Figure 5), two modes are enough to obtain an exact receptance.

Therefore, the great advantage introduced by the fixed boundaries model has been demonstrated.

EXPERIMENTAL RESULTS

In the first step, several experimental tests have been carried out to obtain the response of the assembly formed by the machine, the spindle-and the toolholder following the method proposed by Park et al. [START_REF] Park | Receptance coupling for end mills[END_REF]. This methodology has been chosen because it presents the possibility of characterize the damping of the interface between tool and toolholder in the experimental receptances of the assembly.

This way, the receptance h 3y3y is measured directly using a short blank cylinder, and a long blank cylinder is introduced to extract (equation 4-8) the other two responses of the machine-spindle-toolholder assembly (h 3y3θ , h 3θ3θ ). The properties of these blank cylinders are showed in Table 1.

On the other hand, three carbide tools with different properties have been tested. The material and geometrical properties, including the cross-section variables, of all of them are showed in the Table 1. The machine used in the different test is Soraluce SV6000 milling machine with a Kessler 34kW spindle.

Practical application

A thermal toolholder has been selected to perform the first part of the experimental verification (LAIP 1214121420 -HSK A63) due to their simple joint.

The next figures (Figures 678) shows the results of the receptance coupling compared with an experimental measurement of this response. In general, the results with three tools are good even though the signals are dynamically complex. The signals obtained by RCSA reflect the reality properly in frequency and amplitude, therefore it verifies that the estimation damping is correct. However in spite of these good results, there are some peaks that introduce an error. Thereby, the results with the more slender tool are better than with the other two tools. The reason of this is that when the tool is a slender tool, the mode shapes are clearer due to the dominant bending mode of the tool and in this way, the receptance coupling method obtains better results.

Hydraulic toolholder Collet chuck toolholder

Effect of toolholder

The objective of this section is to validate this method for other types of toolholders with more complex interfaces. Therefore, two more toolholders have been tested: a hydraulic toolholder (Gühring GM 300 -HSK A63) and a collet chuck toolholder (ROMH CMBH 753216 -HSK A63).

Analyzing the results (see Figure 9), it is clear that the proper predictions of the RCSA are maintained for these toolholders too.

Comparing the results for the same tools in different tool holders, the results are not clear, In general, the amplitudes are the same for the first and third tools, but for the second tool the dynamic stiffness is lower with the thermal toolholder. In the other hand, it is clear that the collet chuck toolholder drives to lower frequencies due to its big mass. Anyway, considering that the toolholders present some nonlinearity a wider study is required to obtain general conclusions.

Stability prediction

The accuracy of the predictions obtained using the RCSA with the new approach has been measured using a stability model proposed by Altintas and Budak [START_REF] Altintas | Analytical Prediction of Stability Lobes in Milling[END_REF].

In this paper, the stability lobes have been calculated for slotting operation taking into account the cutting coefficients showed in Table 3. In the Figures 101112, the stability lobes diagrams obtained using the RCSA has been compared with the results provided by directly measured FRF. (see section chapter 4.1).

Cutting coefficients Value

The diagrams show that the approximation is good enough and the RCSA is able to predict the general shape of the diagram. The minimum stable depth of cut and the spindle speed related to the maximum stability point has been compared in the Table 4.

Deviation

CONCLUSION

A new approach to perform the receptance coupling (RCSA) has been proposed based on fixed boundary calculation. Applying this method, it is not necessary to calculate a high number of modes, and therefore they can be calculated analytically with Timoshenko beam theory. Furthermore, due to its good behavior with the receptance coupling method, this model is valid to predict dynamic parameters of different assemblies with different toolholders. Consequently, it is possible to obtain a fairly accurate stability diagram, especially reliable to determine the most stable spindle speeds.

On the other hand, this article shows the good results obtained by the methodology proposed by Park et al. [START_REF] Park | Receptance coupling for end mills[END_REF] to characterize the experimental data. In fact, this method provides a useful tool to consider the damping of the interface between tool and toolholder without any additional parameter.

Figure 1 :

 1 Figure 1: Tool and toolholder-spindle substructures

Figure 2 :

 2 Figure 2: Addition of fixed point movements

Figure 3 :

 3 Figure 3: a) Cross-sections (Kivank et al. [12]); b)Region considered for each cross-section.

Figure 4 :

 4 Figure 4: Results of RCSA with the usual free boundaries model.

Figure 5 :

 5 Figure 5: Results of RCSA with the new fixed boundaries model.

Figure 6 :

 6 Figure 6: Measured and predicted FRF of the tool 1 with thermal toolholder.

Figure 7 :Figure 8 :

 78 Figure 7: Measured and predicted FRF of the tool 2 with thermal toolholder

Figure 9 :

 9 Figure 9: Measured and predicted FRF of the tool 3 with thermal toolholder.

  Tangential direction, K t 700 N/mm 2 Radial direction, K r 245 N/mm 2

Figure 10 :

 10 Figure 10: Stability lobe diagram of the tool 1 with the thermal toolholder.

Figure 11 :

 11 Figure 11: Stability lobe diagram of the tool 2 with the thermal toolholder.

Figure 12 :

 12 Figure 12: Stability lobe diagram of the tool 3 with the thermal toolholder.

  :

	g	y	y	u	;	g	y	y	v	;	h	y	y	a	;	h	y	y		b	;		h	y	y	c	;
	h	y	y	d	;	h		y	e	;	h			f	;	h	y	y	h	y	y	g	;
	h		y	h	y			;	h			h				.

Table 1 :

 1 Properties of the blank cylinders and tolos.

		Short Blank	Long Blank	Tool 1	Tool 2	Tool 3
	Material	Carbide	Carbide	Carbide	Carbide	Carbide
	E (GPa)	580	580	580	580	510
	 (Kg/m 3 )	14500	14500	14350	14000	13800
	Diameter (mm)	20	20	20	20	20
	Total length (mm)	50	120	145	104	100
	Flute length (mm)	-	-	36	44	36
	Number of flutes	-	-	3	4	3
	r (mm)	-	-	8.52	7.30	9.07
	a (mm)	-	-	2.50	2.70	2.30
	fd (mm)	-	-	3.47	3.22	2.99

Table 2 :

 2 Cutting forcé coefficients.

Table 3 :

 3 Stability lobes results.

ACKNOWLEDGEMENTS

This research was partially supported by the OPENAER project (CENIT program of Technical Industrial Development Center (CDTI) of the Spanish government.


Teresa Cao-Rial

Peregrina Quintela

Carlos Moreno

NUMERICAL SOLUTION OF A TIME-DEPENDENT SIGNORINI CONTACT PROBLEM

The purpose of this work is to study the dynamic frictionless contact problem between an elastic body and a rigid foundation. In order to model the contact we consider Signorini conditions. A numerical algorithm is proposed to approximate the solution; the algorithm involves a contact multiplier, which is a fixed point of a nonlinear equation solved by using a generalized Newton method. We use one of the Newmark methods for time discretization and a finite element method for space discretization. The convergence of the method is numerically studied, and a simple test problem is used to validate the methodology.

1.

Introduction. Contact problems involving deformable bodies are present in many industrial processes, as well as in many aspects of everyday life. For this reason, in recent years they have been widely studied considering various constitutive laws and boundary conditions. Only some examples are the work of Kikuchi and Oden [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF], as well as Chau et al. [START_REF] Chau | Dynamic frictionless contact problems with normal compliance[END_REF], Laursen and Chawla [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF], Bécache et al. [START_REF] Bécache | A new family of mixed finite elements for the linear elastodynamic problem[END_REF] or Khenous et al. [START_REF] Khenous | Comparison of two approaches for the discretization of elastodynamic contact problems[END_REF] and the references therein. This paper studies a particular case: a dynamic frictionless contact problem between a linearly elastic deformable body and a rigid foundation. The contact is modelled by using Signorini conditions. The problem is analyzed theoretically and numerically. An existence result is presented and a numerical algorithm is proposed to obtain approximate solutions. The paper is organized as follows. In Section 2 we describe the problem, the contact conditions considered, and the mathematical formulation related to it. In Section 3 we present the functional framework considered and we establish an existence result. Section 4 is devoted to the numerical solution of the two-dimensional problem: its variational formulation as an inequality; its discretization in space by a finite element method; the derivation of a formulation as an equality, by means of the inclusion of a contact multiplier; and the time discretization. Also in this section, an algorithm in time and space is proposed to compute an approximate solution of the problem. Finally, in Section 5 the algorithm is tested on a simple problem, and numerical evidence of the convergence of the algorithm is presented. The body is assumed to be subjected to volume forces of density f 0 . The boundary of the body is partitioned into three mutually disjoint parts: Γ C , Γ N and Γ D , the latter with positive measure, mes(Γ D ) > 0. The body is under traction forces of density f 1 on Γ N , the displacements of the solid are prescribed on Γ D , and Γ C is the potential contact surface where we consider Signorini contact conditions. We denote by n the unitary outward normal vector. Then, the contact problem can be posed as follows:

Problem (P): Find (u, σ) satisfying: ρü -divσ(u) = f 0 in Ω × (0, T ), ( 1
)
σn = f 1 on Γ N × (0, T ), (2)
u = 0 on Γ D × (0, T ), ( 3
)
σ t = 0; σ n ≤ 0 on Γ C × (0, T ), ( 4
)
u n ≤ 0; σ n u n = 0 on Γ C × (0, T ), (5)
u(x, 0) = u 0 ; u(x, 0) = u 1 in Ω, (6)
where

f 1 ∈ W 2,∞ (0, T ; [L 2 (Γ N )] n ∩ [H -1 2 (Γ)] n ), f 0 ∈ W 2,∞ (0, T ; [L 2 (Ω)] n
) and σ(u) = Λ -1 ε(u), Λ being the elasticity tensor assumed to be time independent, symmetric and coercive. The usual notations have been used for the displacements, u, the strain tensor, ε(u), and the stress tensor, σ; the first and second time derivative are noted by u and ü respectively. The initial conditions u 0 and u 1 are assumed to belong to [H 1 (Ω)] n and satisfying

divΛ -1 ε(u 0 ) ∈ [L 2 (Ω)] n .
3. Functional framework and existence of a solution. Let V be the space defined by

V = {v ∈ [H 1 (Ω)] n ; v = 0 on Γ D }, and
V ad = {v ∈ V ; v n ≤ 0 on Γ C }, the closed and convex subset of admissible displacements. We consider the space of the stress fields

X = {τ = (τ αβ ) ∈ [L 2 (Ω)] n 2 ; τ αβ = τ βα }, ( 7
)
which endowed with the norm

τ X = Ω τ : τ dx 1/2 , ( 8
)
is a Hilbert space, τ : τ being the natural scalar product in X.

Let E be the subspace of X defined by

E = {τ ∈ X; div(τ ) ∈ [L 2 (Ω)] n }, (9)
which is also a Hilbert space with the norm

τ E = τ X + div(τ ) [L 2 (Ω)] n . ( 10
)
Given any function f ∈ [L 2 (Γ N )] n , let E ad (f ) be the set of admissible stresses

E ad (f ) = {τ ∈ E; τ t = 0 and τ n ≤ 0 on Γ C ; τ n = f on Γ N }.
Theorem 1 (Existence of a solution). Under the previous hypotheses for f 0 and f 1 , there exists a solution (u, σ(u)) of Problem (P) verifying:

• u ∈ L ∞ (0, T ; V ad ), u ∈ L ∞ (0, T ; [L 2 (Ω)] 2 ), and ü ∈ D (0, T ; [L 2 (Ω)] 2 ).
• The stress tensor σ(u) belongs to D (0, T ;

E ad (f 1 )) ∩ L ∞ (0, T ; [L 2 (Ω)] 4 ).
A sketch of the proof for this result is given in [START_REF] Cao | Existence of solutions for a dynamic Signorini's contact problem[END_REF].

4. Numerical solution of the problem. In the Section, we focus on the numerical solution of the problem in a two-dimensional domain. An algorithm for the solution of the discrete problem, the basic ideas of which come from the existence proof, is proposed.

4.1. Variational formulation. Following standard variational techniques we get a formulation of Problem (P) as a variational inequality given by:

Problem (VP): Find u : [0, T ] → V ad and σ : [0, T ] → E verifying a.e. t ∈ (0, T ), Ω ρü • (v -u)dx + Ω σ(u) : ε(v -u)dx ≥ F (v -u), ∀v ∈ V ad , ( 11
)
and the initial conditions [START_REF] Chau | Dynamic frictionless contact problems with normal compliance[END_REF]. The stress tensor is related to the deformation tensor by the constitutive law σ(u) = Λ -1 ε(u), and

F (v -u) = Ω f 0 • (v -u)dx + Γ N f 1 • (v -u)dγ. ( 12
)
4.2. Space discretization and treatment of the contact conditions. Let T h be a triangular mesh of the two-dimensional domain compatible with the boundary partition. Given one edge e on the mesh induced by T h on the boundary, we denote by n e the unit outward normal to Ω at the middle point of e.

Let S h = {e i , 1 ≤ i ≤ m h } denote the edges on Γ C .
We consider the discrete space

V h = {v h ∈ [C 0 ( Ω)] 2 ; v h|K ∈ [P 1 (K)] 2 , ∀K ∈ T h ; v h = 0 on Γ D }, ( 13
)
where P 1 (K) denotes the space of piecewise linear functions defined on the element K, and the subset of admissible discrete displacements

V ad h = {v h ∈ V h ; v h|e i • n e i ≤ 0, 1 ≤ i ≤ m h },
where v h|e i represents the value of v h at the middle point of e i . The discrete space considered for the stresses is given by:

E h = {σ h ; σ h|K ∈ [P 0 (K)] 4 , (σ h ) αβ = (σ h ) βα ∀K ∈ T h }, ( 14
)
P 0 (K) being the space of constant functions on K.

The variational formulation (11) can now be discretized as:

Problem (VHP): Find u h : [0, T ] → V ad h and σ h : [0, T ] → E verifying a.e. t ∈ (0, T ), Ω ρü h • (v h -u h )dx + Ω σ h (u h ) : ε h (v h -u h )dx ≥ F h (v h -u h ), ∀v h ∈ V ad h , (15) F h being defined by F h (v h -u h ) = Ω f 0h • (v h -u h )dx + Γ N f 1h • (v h -u h )dγ, ( 16
)
where f 0h and f 1h are the piecewise linear functions approximating f 0 and f 1 respectively. The discrete constitutive law is expressed as

σ h (u h ) = Λ -1 ε h (u h ) = λtr ε h (u h ) + 2μ ε h (u h ), (17)
where ε h (u h ) is a function constant per element calculated as the linearized strain tensor associated with u h|K , and λ and μ are the Lamé parameters of the material.

Treatment of the contact conditions.

In this subsection, we obtain a formulation of the problem as a variational equality involving a contact multiplier p h , which is a fixed point of a nonlinear equation. To approximate this multiplier, we consider the spaces:

P h = {q h ∈ L ∞ (Γ C ); q h|e i ∈ P 0 (e i ), 1 ≤ i ≤ m h },
and

Q h = {q h ∈ P h ; q h|e i ≤ 0, 1 ≤ i ≤ m h }.
We also define the operator B : V h -→ P h given by:

B(v h ) = p h with p h|e i = v h|e i • n e i , 1 ≤ i ≤ m h .
Then, following the papers of Bermúdez & Moreno [START_REF] Bermúdez | Duality methods for solving variational inequalities[END_REF] and Barral & Quintela [START_REF] Barral | A numerical method for simulation of thermal stresses during casting of aluminium slabs[END_REF], where maximal monotone operator techniques and subdifferential operators theory are used, we derive the following equivalent formulation of Problem (VHP):

Problem (DVP): Find u h : [0, T ] → V ad h and σ h : [0, T ] → E h verifying a.e. t ∈ (0, T ), for all v h ∈ V h , Ω ρü h • v h dx + Ω σ h : ε h (v h )dx = F h (v h ) - Γ C p h B(v h )dγ, (18)
p h = G λ c (B(u h ) + λ c p h ) , (19)
together with the discrete initial conditions u 0h , u 1h and constitutive law (17). In (19), λ c is an arbitrary positive real parameter and G λ c (s) is defined by

G λ c (s) = 1 λ c s -Π Q h (s) , (20)
Π Q h being the orthogonal projection operator over Q h . Note that

G λ c (s) = 0 if s ≤ 0, s λ c if s ≥ 0. (21)
Then,

p h = 0 i f B(u h ) + λ c p h ≤ 0, 1 λ c B(u h ) + p h if B(u h ) + λ c p h ≥ 0. (22)
4.4. Algorithm in space. To solve Problem (DVP) we propose an iterative algorithm based on a generalized Newton method. Given t ∈ [0, T ], and the initial values (u h0 , σ h0 , p h0 ), we compute successive approximations (u hk , σ hk , p hk ), k ≥ 1 of the dynamic solution (u h , σ h , p h ) at time t. To compute this approximation of the solution at iteration k, it is necessary to distinguish the effective contact region and its complement on Γ C . This is given by the following edge sets on the boundary:

Γ + C,k = {e i ∈ S h ; B(u hk ) + λ c p hk |e i > 0}, Γ - C,k = {e i ∈ S h ; B(u hk ) + λ c p hk |e i ≤ 0}. Now, from (22) we deduce that • B(u hk ) is approximated by zero if there exits contact, i.e., if e i ∈ Γ + C,k-1 . • p hk is approximated by zero on Γ - C,k-1 . 4.4.1.
Steps of the algorithm.

• Step 1: Given (u hk-1 , σ hk-1 , p hk-1 ), we compute (u hk , σ hk ) satisfying:

Ω ρü hk • v h dx + Ω σ hk : ε h (v h )dx= F h (v h ) (23) B(u kh )= 0 on Γ + C,k-1 , ( 24
)
for all v h in V h , such that B(v h ) = 0 on Γ + C,k-1 . In practice, the second equation is introduced in the first one as a penalization term, so (u hk , σ hk ) will be the solution of

Ω ρü hk • v h dx + Ω σ hk : ε h (v h )dx + 1 Γ + C,k-1 B(u hk )B(v h )dγ = F h (v h ), ( 25
)
for all v h in V h , being a small parameter.

• Step 2: Now, once computed (u hk , σ hk ) and since p hk = 0 on Γ - C,k-1 , from ( 19) and (25)we get that

p hk = 1 B(u hk ) on Γ + C,k-1 ,
and the sets Γ Then, we discretize equation (25) by using an implicit method of the Newmark family (see [START_REF] Hughes | The Finite Element Method[END_REF]), which consists of the following equations:

M ün+1 hk + (K + P n+1 k )u n+1 hk = F n+1 h , ( 26
)
where M represents the mass matrix, K is the stiffness matrix, P is the penalization matrix coming from the integral on Γ + C,k-1 in equation ( 25) and F is the force vector. The relations between displacements, velocities and accelerations are given by:

u n+1 h = u n h + Δt un h + Δt 2 2 [(1 -2β)ü n h + 2β ün+1 h ], (27)
un+1 h = un h + Δt[(1 -γ)ü n h + γ ün+1 h ]. ( 28
)
The parameters β and γ take the values β = 1/4, and γ = 1/2 for which the method is unconditionally stable. From here on we omit the subscript k to simplify the notation. There are several possible implementations for Newmark methods. We consider the a-form which computes the acceleration at time t n+1 from the data of the previous steps. This form of implementation consists of the following three steps:

• Defining predictors: having determined the displacement, velocity and acceleration fields at time n, we define

ũn+1 h = u n h + Δt un h + Δt 2 2 (1 -2β)ü n h , ( 29
) ũn+1 h = un h + (1 -γ)Δtü n h , ( 30
)
which will be used to compute the acceleration at time n + 1. • Computing accelerations: rewriting (26) in terms of ũn+1 h defined by (29), we get an equation only in terms of ün+1 h and the predictors, so we solve

(M + βΔt 2 (K + P n+1 ))ü n+1 h = F n+1 h -(K + P n+1 )ũ n+1 h . ( 31
)
• Computing velocities and displacements: from ( 27), ( 28), ( 29) and (30) we deduce that the displacement and velocity fields at time n + 1 are given by

u n+1 h = ũn+1 h + βΔt 2 ün+1 h , un+1 h = ũn+1 h + γΔtü n+1 h .
5. Numerical Results. In order to test the accuracy of the algorithm, we consider a simple test problem with a known solution and compute its numerical solution using the method described in the previous section.

We consider the two-dimensional domain Ω = [0, 0.5] × [0, 0.5] m 2 and the time interval [0, 2.e-06] s.

The material properties are the following:

• Young Modulus E = 7.41e+10 N/m 2 . • Poisson's coefficient: ν = 0.3302. • Density ρ = 2.7e+03 Kg/m 3 .
The contact boundary Γ C is considered to be the straight line [x 2 ≡ 0], Γ N = [x 2 ≡ 0.5] is the upper boundary, and

Γ D = Γ \ (Γ C ∪ Γ N )
. The problem to be solved is

ρü -divσ(u) = f 0 , in Ω × [0, 2.e-06], (32)
σn = f 1 , on Γ N , ( 33
) u = û, on Γ D , ( 34
)
σ t = h, σ n ≤ 0, on Γ C , ( 35
)
u n ≤ s(x 1 ), σ n (u n -s(x 1 )) = 0, on Γ C , ( 36
) u(x, 0) = u 0 (x), u(x, 0) = u 1 (x), in Ω, (37)
where s(x 1 ) is the initial gap between the deformable body and the rigid foundation given by

s(x 1 ) = ⎧ ⎨ ⎩ (0.3 -x 1 ) 3 if 0 ≤ x 1 ≤ 0.3, 0 i f 0 .3 ≤ x 1 ≤ 0.4, (x 1 -0.4) 3 if 0.4 ≤ x 1 ≤ 0.5.
Let {e 1 , e 2 } be the canonical basis for R 2 . Then, the volume and traction forces are given by As can be seen in Tables 1234, as the size of the mesh, h, decreases so does the relative error in the displacements, velocities, accelerations and stresses respectively. This convergence can also be perceived in Figure 3, where the graphic shows the t = 4.e-07 t = 8.e-07 t = 1.2e-07 t = 1.6e-07 t = 2.e-06 h = 0.0732 5.4410e-05 5.4443e-04 1.0896e-03 2.1819e-03 5.4752e-03 h = 0.0366 3.7908e-05 3.7955e-04 7.6015e-04 1.5244e-03 3.8410e-03 h = 0.0183 2.0618e-05 2.0647e-04 4.1355e-04 8.2928e-04 2.0844e-03 h = 0.0046 4.7353e-06 4.7385e-05 9.4748e-05 1.8881e-04 4.5545e-04 h = 0.0011 1.1837e-06 1.1753e-05 2.2956e-05 4.1675e-05 8.0119e-05 Table 2. Evolution of relative error in the velocities. convergence rate of the algorithm with respect to the mesh size. All these simulations have been made considering Δt = 2.e -08 and κ = 1.e + 3. The above algorithm behavior is reproduced when the acceleration increases, that is, for bigger values of κ 2 , as can be seen in Table 5, for κ = 1.e + 5. Similarly, in Table 6 we show that as the time step is reduced, so is the relative error in the displacements.

f 0 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 6A le 2 if 0 ≤ x 1 ≤ x 1 + l, -(λ + μ) l3 l 2 (3 l -4l)e 1 + 6x 2 A l2 l( l2 -4 ll + 2l 2 )e 2 if x 1 + l ≤ x 1 ≤ x 1 -l, 6Ale 2 if x 1 -l ≤ x 1 ≤ 0.5
Finally, Figure 4 illustrates the influence of the parameter κ on the errors for several Δt. The algorithm has been implemented in a software package named CRACKEW2D elaborated in MATLAB code and executed on a Pentium(R) 4, 3.00 GHz. Using a mesh of 32 elements (h = 0.0366) and a time step of 1.e-08 the CPU time needed was 0.4062 seconds and using a mesh of 512 elements (h = 0.0092) and the same time step, the CPU time needed was 1.2031 seconds.

2 .

 2 Statement of the problem. Let us consider an elastic solid with constant density ρ, initially occupying the bounded domain Ω ⊂ R n , n = 2, 3 of class C 1,1 .

5 .

 5 Time discretization. To compute a numerical solution of the problem we discretize the equations in time. We consider a regular partition of the time interval [0, T ] in I subintervals such that t 0 = 0, t i+1 = t i + Δt, Δt = T I , i = 0, ..., I -1.

Figure 1 .

 1 Figure 1. Stresses of the computed solution plotted on the deformed mesh at initial and final instant.

t = 4 Table 4 .Table 5 .

 445 .e-07 t = 8.e-07 t = 1.2e-07 t = 1.6e-07 t = 2.e-06 h = 0.0732 4.0802e-01 4.0777e-01 4.0748e-01 4.0688e-01 4.0490e-01 h = 0.0366 2.8429e-01 2.8449e-01 2.8470e-01 2.8508e-01 2.8593e-01 h = 0.0183 1.5463e-01 1.5478e-01 1.5490e-01 1.5500e-01 1.5412e-01 h = 0.0046 3.5512e-02 3.5478e-02 3.5335e-02 3.4718e-02 3.0496e-02 h = 0.0011 8.8759e-03 8.6621e-03 8.0236e-03 7.0398e-03 5.0352e-03 Table3. Evolution of relative error in the accelerations. t = 4.e-07 t = 8.e-07 t = 1.2e-07 t = 1.6e-07 t = 2.e-06 h = 0.0732 5.4699e-01 5.4721e-01 5.4745e-01 5.4793e-01 5.4938e-01 h = 0.0366 3.1676e-01 3.1693e-01 3.1711e-01 3.1748e-01 3.1860e-01 h = 0.0183 1.6807e-01 1.6817e-01 1.6827e-01 1.6849e-01 1.6913e-01 h = 0.0046 4.3723e-02 4.3749e-02 4.3777e-02 4.3836e-02 4.4013e-02 h = 0.0011 1.1034e-02 1.1040e-02 1.1048e-02 1.1063e-02 1.1109e-02 Evolution of relative error in the stresses with κ = 1.e + 3. t = 2.e-08 t = 2.e-07 t = 4.e-07 t = 1.2e-06 t = 2.e-06 h = 0.0732 5.4938e-01 5.7166e-01 5.9735e-01 5.6776e-01 4.7492e-01 h = 0.0366 3.1860e-01 3.3628-01 3.5826e-01 3.7799e-01 3.0112e-01 h = 0.0183 1.6912e-01 1.7935e-01 1.9225e-01 2.1811e-01 1.6490e-01 h = 0.0046 4.4009e-02 4.6800e-02 5.0348e-02 5.9031e-02 4.3563e-02 h = 0.0011 1.1107e-02 1.1818e-02 1.2724e-02 1.5014e-02 1.1026e-02 Evolution of relative error in the stresses with κ = 1.e + 5.

Figure 3 .Table 6 .

 36 Figure 3. Convergence rate in the displacements, velocities, accelerations and stresses at the first time step with κ = 1.e + 3, Δt = 2.e -08.

Figure 4 .

 4 Figure 4. Influence of Δt for different values of κ.

Table 1 .

 1 Evolution of relative error in the displacements.

	1,00E-04
	1,00E-05
	1,00E-06
	1,00E-07
	0,0011	0,0023	0,0046	0,0092	0,0183	0,0366	0,0732
				mesh size
	Dirichlet	Dirichlet and Neumann	Dirichlet, Neumann and Contact
	Figure 2. Relative error in the displacements at final instant for
	κ = 1.e + 3, Δt = 2.e -08 and different boundary conditions.
	t = 4.e-07 t = 8.e-07 t = 1.2e-07 t = 1.6e-07 t = 2.e-06
	h = 0.0732 5.4414e-09 5.4491e-07 2.1830e-06 8.7591e-06 5.5252e-05
	h = 0.0366 3.7910e-09 3.7980e-07 1.5223e-06 6.1139e-06 3.8674e-05
	h = 0.0183 2.0619e-09 2.0659e-07 8.2813e-07 3.3261e-06 2.1018e-05
	h = 0.0046 4.7356e-10 4.7429e-08 1.8995e-07 7.6045e-07 4.7067e-06
	h = 0.0011 1.1838e-10 1.1810e-08 4.6747e-08 1.7872e-07 9.3406e-07

Acknowledgements. This work has been supported by CICYT-FEDER (project DPI2004-01993) and the Xunta de Galicia (project PGIDIT05PXIC20701PN).

and

l3 l 2 (3 l -4l)e 1 + (λ + 2μ) l4 l 3 e 2 if x 1 + l ≤ x 1 ≤ x 1l x 2 = 0.5, -3μ l2 e 1 if x 1l ≤ x 1 ≤ 0.5, x 2 = 0.5, for all t ∈ [0, 2.e-6], μ and λ being the Lamé parameters of the material related to the Young modulus and Poisson's coefficient by

To simplify the notation, we have denoted A = (κ 2 ρμ), l = 0.3κtx 1 and l = x 1 + κt -0.4 and κ is a real constant.

The friction force considered is

and the initial conditions are given by

With these data, the solution of problem (32)-(37) for û(0, x 2 , t) = (0.3κt) 3 e 2 , x 2 ∈ [0, 0.5], t ∈ [0, 2.e-06], û(0.5, x 2 , t) = (κt -0.4) 3 e 2 , x 2 ∈ [0, 0.5], t ∈ [0, 2.e-06], is

The computed solution for κ = 1.e + 5 is presented in Figure 1 at initial and final instants respectively.

Numerical Convergence.

In this section we demonstrate the convergence of the proposed algorithm for the above test problem. Firstly, in order to show the accuracy of the contact algorithm, we present in Figure 2 In the following, we present several tables showing the evolution of the relative error for each mechanical variable as the size of the mesh is reduced and as the size of the time step decreases.


Martin Rosalie
email: martin.rosalie@uni.lu

Grégoire Danoy
email: gregoire.danoy@uni.lu

Pascal Bouvry
email: pascal.bouvry@uni.lu

Serge Chaumette
email: serge.chaumette@labri.fr

UAV Multilevel Swarms for Situation Management

Keywords: CCS Concepts, Computing methodologies → Multi-agent systems, Cooperation and coordination, cooperative UAV, multilevel swarm, mobility model, decisional autonomy

The development and usage of Unmanned Aerial Vehicles (UAVs) quickly increased in the last decades, mainly for military purposes. Nowadays, this type of technology is used in non-military contexts mainly for civil and environment protection: search & rescue teams, fire fighters, police officers, environmental scientific studies, etc. Although the technology for operating a single UAV is now mature, additional efforts are still necessary for using UAVs in fleets (or swarms).

This position paper presents the ASIMUT project (Aid to SItuation Management based on MUltimodal, MUltiUAVs, MUltilevel acquisition Techniques). The challenges of this project consist of handling several fleets of UAVs including communication, networking and positioning aspects. This motivates the development of novel multilevel cooperation algorithms which is an area that has not been widely explored, especially when autonomy is an additional challenge. Moreover, we will provide techniques to optimize communications for multilevel swarms. Finally, we will develop distributed and localized mobility management algorithms that will cope with conflicting objectives such as connectivity maintenance and geographical area coverage.

INTRODUCTION

The development and usage of Unmanned Aerial Vehicles (UAVs) quickly increased in the last decades, mainly for military purposes [START_REF] Valavanis | Handbook of Unmanned Aerial Vehicles[END_REF]. Nowadays, this type of technology is used in non-military contexts mainly for civil and environ-ment protection: search & rescue teams, fire fighters, police officers, environmental scientific studies, etc. Although the technology for operating a single UAV is now mature, efforts are still necessary for using UAVs in fleets (or swarms).

This position paper presents the ASIMUT (Aid to SItuation Management based on MUltimodal, MUlti-UAVs, MUltilevel acquisition Techniques) project that aims at developing surveillance services relying on multiple fleets of UAVs. It is planned to have collaborations among UAVs evolving at different altitudes: large (typically fixed-wing) UAVs covering the area from a high altitude with coarse grain acquisition capabilities, and low altitude multi-rotor UAVs in charge of precise close to the field measurements. Additionally, the considered scenario targets multimodal data collection, i.e. acquisition achieved with different types of sensors (e.g. electro-optical or infrared) and will thus allow retrieving meaningful information from the surveyed field. The ASIMUT consortium's vision entails a model proposing a balanced combination of different measures, framed by the types of threats and operator's reaction capabilities. The ASIMUT project aims at increasing detection, highlevel data fusion capabilities and UAV autonomy, by using innovative algorithms permitting to decrease the operator workflow by providing tools to enhance the capabilities of the system.

The remainder of the article is organized as follows. In section 2, a description of the ASIMUT project is given. Then in section 3 we provide a state-of-the-art on swarm decisional autonomy and on mobility management. Merging the latter fields, section 4 presents the innovations developed and implemented for the ASIMUT project. Finally, section 5 presents our conclusions.

ASIMUT PROJECT

The main goal of the ASIMUT project is to address scenarios where a mission is to be achieved based on data fusion coming from different sensors in a number of UAVs that constitute a swarm. The ASIMUT project aims at developing innovative algorithms based on learning techniques dedicated to fusion of data. From early warning process where data are provided by the UAV swarm payload, the target localization, monitoring and classification processing are performed using these techniques. These will improve the quality and significance of the pieces of information provided to an operator through detection and identification processes. It will be typically necessary to combine the information about multiple objects, events and their relations with con- textual information. To this end, the ASIMUT model concept integrates different approaches in three layers:

UAVs Collaboration: UAV swarm flight plan optimization according to external environment;

UAVs Detection: each UAV carries a specific payload to perform detection for air-ground surveillance. The autonomy of the UAV swarm permits to define which UAV(s) is(are) relevant to perform the analysis required by an operator;

Data Exploitation: automatic detection of targets in particular in a real-time FMV (Full Motion Video). The data fusion can be improved so as to decrease the workflow for an operator (Smart Information Management).

The ASIMUT system is composed of a Ground Control Station (GCS) and UAVs that constitute swarms (Figure 1). The GCS role is to communicate and transmit data between the swarms and the person(s) in charge of the situation management. Two types of swarms are distinguished: High Level Coordination Swarm (HLCS) and Low Level Swarm (LLS), composed of fixed-wing and multi-rotor UAVs respectively. The role of the HLCS is to manage the tasks assigned by the GCS or autonomously created by the HLCS itself in reaction to an event occurring in the environment. Finally the LLS role is to collect data and track vehicles in specific areas of this environment.

This distinction implies a multilevel cooperation between UAVs to achieve their missions. In this description the GCS is considered as global entity dealing with situation awareness, high-level data fusion and mission management of the ASIMUT project. The considered system targets multimodal acquisition, i.e. acquisition achieved with different types of sensors (e.g. electro-optical or infrared) and will thus allow retrieving meaningful information from the surveyed field. These data are transmitted to the GCS using a multilevel of collaboration between UAVs to ensure the communication performance. Figure 2 illustrates these relations between the GCS, the swarms (HLCS and LLSs) and the components of the system. This figure also underlines the autonomous capacity of the HLCS to take decisions facing critical situations identified from data collected by the LLSs.

Every considered aircraft will be equipped with a wireless transceiver to communicate with its neighbors, i.e. the aircrafts that are currently in the transmission range of the transceiver. Each UAV will thus be able to exchange information with its neighboring UAVs and for instance determine its flight plan solely based on the knowledge of its local environment and an estimate of the global situation. Indeed gathering a complete precise knowledge of the overall system would decrease the scalability of the approach, since in such networks this would require a continuous broadcasting process of all the aircrafts and no communication failure to maintain everyone updated. If feasible this would dramatically reduce the available bandwidth for other types of communications.

One of the addressed challenges in this project is therefore to handle several fleets, or swarms, of UAVs including their communications, networking and positioning. This thus motivates the development of multilevel cooperation algorithms which is an area that has not been widely explored, especially when autonomy is also a challenge. Moreover, we will provide techniques to optimize communications within a swarm and between swarms (including multilevel swarms). Finally, we will develop distributed and localized mobility management algorithms that will cope with conflicting objectives such as connectivity maintenance and geographical area coverage.

SWARM DECISIONAL AUTONOMY AND MOBILITY MANAGEMENT

This section first contains an analysis of the related work in UAV swarm decisional autonomy, followed by a survey on UAV swarm mobility management. In the first part we consider existing projects under the criteria of the swarm's decisional autonomy and on the fact that they rely, or not, on an existing communication infrastructure. We believe that the CARUS project [START_REF] Chaumette | an operational retasking application for a swarm of autonomous uavs: First return on experience[END_REF] went a step further than a stateof-the-art to deal with autonomous swarms of UAVs. In this project, five UAVs share the supervision of a grid of fifteen points of potential ground incidents. Although often using a motion capture system to localize the drone, the GRASP Laboratory generally supports a high level of autonomy in terms of decision as evidenced by their decentralized flight formation work [START_REF] Turpin | Decentralized formation control with variable shapes for aerial robots[END_REF], their achievement in collaboratively grasping different kinds of objects [START_REF] Mellinger | Cooperative grasping and transport using multiple quadrotors[END_REF]. The most remarkable achievement is a collaborative mapping experiment of an earthquake-damaged building of the Tohoku University by two rovers and one UAV carried by one of them [START_REF] Michael | Collaborative mapping of an earthquake-damaged building via ground and aerial robots[END_REF]. When the rovers have covered every part of the building they could reach, the UAV is autonomously launched at the right spots by the rover to map the specific parts that are only reachable by the UAV (which has less endurance than rovers).

The mobility management aspects of this project tackle issues encountered in both the robotic and the mobile networking domains, since for instance coverage and connectivity are shared objectives. Nowadays, these two scientific communities are working on very close problematics, sometimes without noticing it. In the mobile ad hoc networks domain, influencing the mobility of a device and its network equipment in order to enhance some networkingrelated quality metrics is currently being investigated due to advances in technology. In the robotic networking domain, the maintenance or the reinforcement of the connectivity has always been an important issue, as it greatly influences the fault tolerance of a system. When it comes to controlling the mobility of UAV swarms, many proposals exist in the literature. One can classify them as centralized/distributed and online/offline techniques. Historically, the first flight plan computation techniques are centralized and offline, meaning that the UAVs trajectory is computed by one single entity and distributed before the start of the mission [START_REF] Besada-Portas | Evolutionary trajectory planner for multiple uavs in realistic scenarios[END_REF]. However, these are not suitable due to the unpredictability of the potentially hostile environment in which the aircrafts will be deployed, e.g. harsh weather conditions or technical and mechanical problems.

Centralized and online mission planning techniques have therefore been investigated, allowing to recompute the swarm's flight plan at runtime, which is highly desirable when considering surveillance and tracking missions like in this work. However, updating the whole swarm flight plan via a unique access point, such as satellite connection, also has drawbacks in terms of scalability, robustness and communication costs. Indeed, what would happen if the central entity fails or is not able to communicate with the fleet frequently enough?

This work is thus interested in online and distributed mobility control techniques. In that case, UAVs are equipped with wireless communication interfaces that allow them to exchange information in a peer-to-peer fashion, i.e. no (costly) infrastructure is required. Such networks are a novel type of Mobile Ad Hoc Networks (MANETs) referred to as Flying Ad Hoc Networks (FANETs). That way, each UAV determines its flight plan solely based on aggregated knowledge from neighboring UAVs. Some promising approaches are inspired by nature, and more precisely by ant colonies. Ant Colony Optimization (ACO) is a set of probabilistic techniques for solving computational problems, which can be reduced to finding good paths through graphs. Initially proposed by Marco Dorigo in 1992 in his PhD thesis [START_REF] Dorigo | Optimization, Learning and Natural Algorithms[END_REF], ACO has found applications in many domains. The main idea is that ants use the environment to exchange information via pheromone deposit which concentration influences their behavior. This process is called stigmergy. UAVs deposit repulsive pheromones in visited geographical areas so as to prevent other UAVs to revisit the same places too rapidly [START_REF] Kuiper | Mobility models for uav group reconnaissance applications[END_REF]. However many of the proposed approaches are not realistic since they consider global pheromone maps while each UAV can only have a partial knowledge of the environment, i.e. a local pheromone map. Communicating this information therefore remains a challenge to ensure scalability, i.e. which information and when to broadcast it.

Additional objectives are to be considered such as network connectivity as recently proposed by Schleich et al. for surveillance scenarios [START_REF] Schleich | UAV fleet area coverage with network connectivity constraint[END_REF]. This is motivated by most safety-related regulations that currently require any UAV to be reachable at any moment in case of an emergency. That way maintaining a connected network would ensure that the furthest UAVs can still be controlled via multi-hop communications. However none of these existing works has considered multilevel swarms as proposed in the ASIMUT project. Novel mobility management techniques have thus to be developed to address the aforementioned objectives: coverage/tracking, communication, connectivity, and additionally energy-efficiency.

INNOVATION

The ASIMUT project innovates in several fields, including distributed load sharing, autonomous retasking, multilevel and global cooperation, and bio-inspired swarms mobility models. These are detailed in the following section.

In the ASIMUT project the system is fully decentralized. The GCS only sends missions requirements to the HLCS. The latter has the capability to first decide who will achieve this mission: the HLCS itself or in collaboration with one LLS. After that, the UAVs of the considered swarms will decide autonomously the actions they have to plan to complete this mission including communications. If the HLCS received too many missions, it has to decide to stop one or more missions or to wait until another mission ends before it can assign new ones. As a consequence: the swarms control the load-sharing. They are chosen depending on their capabilities, mainly based on the sensors they embed, and their relative position to the objective. This is done in order to reduce the number of communications and to optimize the system as a whole.

The Figure 3 details the global scheme of the autonomous retasking of swarms. In this diagram, the key point is that a UAV always receives information from the other UAVs and stores them. This information is regularly analyzed by the UAV before it takes a decision. This is the way UAVs collaborate at both levels because information can come from any swarms. The last steps are dedicated to the information of the other UAVs (position of an UAV, ongoing action, . . . ). We will also implement algorithms to select relevant information depending on their meaning and their duration impact. This participates to the communication optimization of the model.

In terms of communication optimization swarm members autonomously elect one of them to take care of the communication with other swarms. This election occurs regularly. Critical situations that the ASIMUT system faces might require collecting information far from the current position of the HLCS. Thus the HLCS can send several LLS to collect data while it remains over the critical situation area and help other LLS to handle it by directly post processing data to help them. This multilevel cooperation between swarm permits to face critical situations with massive amount of ressources involved while routine missions are still accomplished.

ASIMUT will focus on online and decentralized mobility management techniques. More precisely it will rely on bio-inspired techniques, to ensure the optimization of the autonomous UAV swarms movements. In order to comply with the targeted scenario, several possibly conflicting objectives and constraints should be respected by the UAVs, i.e. coverage, connectivity, randomness and energy efficiency. First of all, as the assigned mission is about surveillance, the considered area should be covered completely and on a sufficiently frequent basis. The exploration process will be based on Ant Colony Optimization (ACO). Its design should be light enough to represent the known environment accurately and be efficiently exchanged in a wireless ad hoc fashion between aircrafts. Extensions to this model will be proposed to optimize the mobility of the higher-level fleet, to maintain and optimize the connectivity between the two levels. Although this mechanism is antagonist to the exploration process, we will research how to limit its inhibition impact. In addition, the UAVs mobility should be as hard to predict as possible to prevent possible attacks. Two approaches are considered, either using stochastic or deterministic processes. For the surveillance movement, we plan to compare the usage of randomness to a behavior obtained from chaotic dynamical systems inspired by the mobile robots domain [START_REF] Fallahi | A cooperative mobile robot task assignment and coverage planning based on chaos synchronization[END_REF]. It should also be noted that in real conditions degraded modes must be considered. Loss of a UAV or loss of messages are to be taken into account not only in a military context but also in a civilian context such as disaster relief areas.

CONCLUSION

The ASIMUT project aims at developing surveillance services based on fleets of UAVs flying at different altitudes and equipped with various sensors (e.g. EO, IR, . . . ). The required processing of the collected information can be done at the level of each individual UAV, at the level of a higher altitude one or even on the ground. This thus requires intralevel and inter-level cooperation; the latter being without doubt an area that has not been widely explored, especially when autonomy is also a challenge. Therefore one of the innovations will be in the management and the optimization of the interactions between the swarms that participate in the mission and how their local decisions and work impact the behavior of the other swarms, i.e. of the global system.

Other innovations reside in the development of distributed and online mobility management algorithms. The latter deal with potentially conflicting objectives such as coverage, connectivity preservation and bandwidth efficiency. These innovations will be combined as a global framework for the ASIMUT project and evaluated via state-of-the-art simulations; a preliminary step before its real-world deployment.

Figure 1 :

 1 Figure 1: Multi-level swarms of UAVs.

Figure 2 :

 2 Figure 2: Interaction between system modules.

Figure 3 :

 3 Figure 3: Activity Diagram describing the communication process of a UAV.

Acknowledgement

The authors acknowledge the support of the ASIMUT project A-1341-RT-GP, which is coordinated by the European Defence Agency (EDA) and partially funded by 8 contributing Members (Austria, France, Germany, Italy, Luxembourg, The Netherlands, Poland and Sweden) in the framework of the Joint Investment Programme on Innovative Concepts and Emerging Technologies 2. The ASIMUT project consortium is composed of Thales, Fraunhofer IOSB, Fly-n-Sense, University of Bordeaux and University of Luxembourg.