Jean-Yves Gauthier

Christian Lexcellent

Arnaud Hubert

Joël Abadie

Nicolas Chaillet

J.-Y Gauthier

Ni-Mn-Ga sigle

come

Introduction

The interest of Magnetic Shape Memory Alloys (MSMAs) compared with the classical Shape Memory Alloys (SMAs) is their possible activation not only by stress and temperature actions but also by magnetic field. A model of rearrangement process of martensite platelets in a non-stoichiometric Ni 2 MnGa MSMA single crystal under magnetic field and (or) stress actions has been proposed by the same authors in [1]. The aim of this following paper is to extend these works purpose to the anisothermal behavior when including the process of martensite platelets reorientation and phase transformation austenite ⇋ martensite.

Firstly, MSMA present the same properties as classical shape memory alloys but with the addition of a magnetic field sensibility. Several models are then devoted to the variant reorientation process. Some relatively old models are based on simple energy function [2], energy minimization [3,4] or using a magnetic stress to disconnect mechanical and magnetic behaviours [5]. One of the first thermodynamic approaches was built in [6] and the addition of internal variables in the thermodynamics of irreversible processes was proposed in [7,8] A very interesting experimental study concerning the shape memory and martensitic deformation response of Ni 2 MnGa single crystals is performed by Callaway et al. [9]. Few papers relate the global modeling of MSMA including temperature, stress and magnetic field effects in the same formalism. The thermodynamical approach proposed in this paper is a relevant way to model the complex behaviour of MSMA in a global and macroscopic form. a e-mail: christian.lexcellent@univ-fcomte.fr

Gibbs free energy expression associated with a magneto-thermo-mechanical loading

The Gibbs free energy G expression can be split into four parts : the chemical one G chem (generally associated with the latent heat of the phase transformation), the thermal one G therm (associated with the heat capacity), the mechanical one G mech and the magnetic one G mag .

G(Σ, T, H,z o , z 1 , . . . , z n , α, θ, α A ) = G chem (T, z o ) + G therm (T ) + G mech (Σ, z o , z 1 , . . . , z n ) + G mag (H, z o , . . . , z n , α, θ, α A ) (1)
where the state variables are:

-Σ = diag(0, σ, 0) the applied stress tensor, -H = H • x the magnetic field, -T the temperature.

The internal variables are:

-z o the austenite volume fraction, -z k the volume fraction of martensite variant Mk (k = 1 . . . n, i.e. the martensite presents n different variants), -α the Weiss domain proportion inside a martensite variant (see figure 1), -α A the Weiss domain proportion inside the austenite phase A, -θ the rotation angle of the magnetization vector M under the magnetic field H.

The coupling between mechanic and magnetism is made by the choice of internal variables, for instance z i , and not by the addition of an interaction term G mech-mag in the free energy. For a Ni-Mn-Ga single crystal presenting three variants of martensite and a mother phase (cubic austenite ⇋ tetragonal martensite), under thermal, magnetic loading and mechanical one, the Gibbs free energy expression can be:

M 1 M 2 (1-z) (z) x y α α - θ θ ¡ ¢
ρG(H, σ, T, z 0 , z 1 , z 2 , z 3 , α, θ, α A ) = u M o -T s M o + z o (∆U -T ∆S ) + C p (T -T o ) -T • log T T o - σ 2 (z 1 + z 3 ) β 2 a -1 + z 2 β 2 c -1 - 1 2 σ 2 E + Az o (1 -z o ) + K(z 1 z 2 + z 1 z 3 + z 2 z 3 ) -µ 0 m s (T ) z 1 (2α -1)H - m s (T ) 2χ a (2α -1) 2 + (z 2 + z 3 ) sin(θ)H - m s (T ) 2χ t (sin(θ)) 2 +z 0 (2α A -1)H - m s (T ) 2χ A (2α A -1) 2 (2)
14th International Conference on Experimental Mechanics

m 0 = tanh T c m s (T ) m O T
.

C p is the heat capacity, E the Young modulus, β a = a a 0 , β c = c a 0 . a 0 lattice parameter of A. a and c lattice parameters of M; A and K constant belonging to the mechanical part of free energy G. T c is the Curie temperature and m 0 the magnetisation at 0K. In order to simplify the model, m 0 and T c parameters are considered as the same for austenite and martensite. µ 0 , χ a and χ t are the vaccum magnetic permeabilty, the magnetic susceptibilities in the easy and transverse axis.

Clausius-Duhem inequality

In the classical frame of Thermodynamics of Irreversible Process, the total strain in y direction, the magnetization in x direction, and the entropy can be written as:

ε = - ∂ρG ∂σ µ o m = - ∂ρG ∂H s = - ∂ρG ∂T (3)
A magneto-thermal effect takes place in the entropy expression due to the temperature dependence of m s . This effect will be neglected in the present paper. The thermodynamical forces associated with the progression of the Weiss domains widths α, α A and rotation angle θ are:

∂ρG ∂α A = 0 ∂ρG ∂α = 0 ∂ρG ∂θ = 0 (4)
The choice of the free energy expression confirms that the pure magnetic behavior is considered as reversible. Actually, the magnetization curves of the two martensite variants have no hysteresis on the figure 2 taken from [10]. Finally, the thermodynamical forces associated with the z i martensite and austenite fractions are:

π f i = - ∂ρG ∂z i for i = 0, 1, 2, 3 (5)
The mechanical behavior is highly irreversible, hence the Clausius-Duhem inequality has to be written:

dD = -ρdG -µ o mdH -εdσ -sdT ≥ 0 ( 6
)
where dD is the dissipation increment. This expression can be reduced to:

dD = π f o dz o + π f 1 dz 1 + π f 2 dz 2 + π f 3 dz 3 ≥ 0 with 3 k=0 dz k = 0 (7)

kinetic equations (example with two variants)

In the previous inequality, three of the four variables are independent and the evolution depends on the configuration of the forces. For example, if the sample contains only two variants (z o = z 3 = 0) during the evolution at constant temperature, then z = z 1 = 1 -z 2 and the Clausius-Duhem inequality becomes: An external loop, e.g. a complete rearrangement from z = 0 to z = 1 (path a) and from z = 1 to z = 0 (path b), is reported on figure 3. Rearrangement begins when (π f 1 -π f 2 ) ≥ π cr (T ) for the path a and when (π f 1 -π f 2 ) ≤ -π cr (T ) for the path b. After the rearrangement start, the behavior is modeled according to the following kinetic equations:

dD = π f 1 dz 1 + π f 2 dz 2 ≥ 0 dD = (π f 1 -π f 2 )dz 1 ≥ 0 (8)
π f 1 -π f 2 = λż with: ż = ż1 = -ż 2 (9)
λ is considered as a constant value in the present paper. But, in [1], the value of λ was considered depending on the previous strain history. The concept of memorized point was introduced and a difference appears between partial loops and major loops. Moreover, π cr (T ) is a function of the temperature as it was shown in [11]. A linear dependence is used in this paper as A numerical simulation was done with the help of the Matlab/Simulink software. The powerful of the modeling approach presented here is that either mechanical (pseudo-elasticity and martensite reorientation), thermal, and magnetic effects can be taken into account in the same numerical simulation. Some different cases will be considered in this section.

π cr (T ) = π o cr + k cr (A o S -T ). Parameters are summarized in the table 1 with ∆S = 2∆U As+Ms , A = -∆S (As-Ms) 2 and λ A = -∆S (A f -Ms). A o S = 309.4 K M o S = 301.7 K A = 5.48.

Conclusion

The purpose of this paper was to propose a full single crystal modeling of Magnetic Shape Memory alloys including the magneto-thermo-mechanical coupling. In the frame of the thermodynamics of irreversible processes, a model was proposed using a Gibbs free energy expression. As it can be seen in ( 2), the model is quite complicated and its use may be difficult. Nevertheless in practical applications and for specific conditions, some reductions can be made (isothermal process on 2D motions) and in such a case, the usefulness of the model is shown. For example, in the paper [13], an extension of the quasi-static isothermal model to the dynamical case was made, including the magnetic circuit creating the magnetic field and a dynamical load applied to the MSMA sample. A thermodynamic approach with hamiltonian modeling was used by the authors to design actuators [13] and new control laws [START_REF] Gauthier | IEEE/RSJ International Conference on Intelligent Robots and Systems[END_REF]. Moreover, the thermodynamical approach proposed in this paper to model a specific MSMA can be easily extended to other materials, such as future MSMA monocristals and polycristals. Finally, this global model can be a basis to model MSMA in a Finite Element Analysis to design new mechanical structures for actuation and sensing applications.

Fig. 1 .

 1 Fig. 1. Representative Elementary Volume (REV) when the sample is composed by only two martensite variants M 1 and M 2 (z = z 1 and 1 -z = z 2 ).

  1, u M O and s M O the internal and entropy of martensite. m s (T )

Fig. 2 .

 2 Fig. 2. Magnetization curves measured in easy (free sample) and hard (sample constrained by stress) magnetization directions [10].

Fig. 3 .

 3 Fig. 3. Thermodynamical force [π f 1 -π f 2 ](σ, α, θ) as a function of the M 1 martensite fraction z ∈ [0, 1].

  10 5 J/m 3 a o = 5.84 Å a = 5.95 Å c = 5.60 Å E = 5.10 9 Pa λ M = 4.10 5 K = 0 χ a = 5 χ t = 1.05 χ A = 1.76 T c = 370 K m S 0 = 710 kA/m λ A = 1.26.10 6 J/m 3 π o cr = 12.10 3 J/m 3 k cr = 800 Pa/K

Fig. 4 .

 4 Fig. 4. Simulation results of thermal action with or without magnetic field and mechanical stress.

Fig. 5 .

 5 Fig.5. Experimental results of strain and magnetization under magnetic action and constant stress (σ = -1 MPa) at different temperature levels (taken from[12]).

Fig. 6 .

 6 Fig.6. Simulation results of strain and magnetization under magnetic action and constant stress at different temperature levels (σ = -1 MPa, experimental results in F. 5).

Table 1 .

 1 Selected material parameters.

14th International Conference on Experimental Mechanics


Mohamed Kharboutly
email: mohamed.kharboutly@femto-st.fr

Membre Michaël Gauthier
email: michael.gauthier@femto-st.fr

Membre Nicolas Chaillet

Modeling the

I. INTRODUCTION

Manufactured products become always smaller and integrate more and more functionalities in small volumes. Several application fields are concerned such as bio-engineering, telecommunications or in a more general way Micro-Electro-Mechanical-Systems (MEMS). The assembly of these microproducts is a great challenge because of the microscopic size of the components. In fact, the major difficulties of micro-assembly come from the particularity of the microobject's behaviours which are more function of the surface forces than the volumic forces [START_REF] Lambert | Capillary Forces in Micro-assembly[END_REF], [START_REF] Gauthier | Forces analysis for micromanipulations in dry and liquid media[END_REF], [START_REF] Zhou | Ambient environment effects in micro/nano handling[END_REF]. The manipulation of a micro-object requires its handling, positioning, and release without disturbances of the surface forces such as electrostatic forces, van der Waals forces or capillary forces. The release is the most critical phase which is usually perturbed by adhesion phenomenon.

Several methods have been proposed in the last ten years to improve micromanipulations [START_REF] Sariola | Hybrid microhandling: a unified view of robotic handling and self-assembly[END_REF], [START_REF] Lambert | A case study of surface tension gripping: the watch bearing[END_REF]. The first approach deals with contact manipulation where the adhesion should be reduced or could be directly used for manipulation [START_REF] Hériban | Robotic micro-assembly of microparts using a piezogripper[END_REF], [START_REF] Lambert | Surface and contact forces models within the framework of microassembly[END_REF], [START_REF] Driesen | Micromanipulation by adhesion with two collaborating mobile micro robots[END_REF]. The release requires innovative methods overcome adhesion. Dielectrophoresis force, which is the force applied on a polarizable particle in a non uniform electric field, has nicolas.chaillet@femto-st.fr been recently used to induce a repulsive force on microobjects in order to release them [START_REF] Gauthier | Submerged robotic micromanipulation and dielectrophoretic micro-object release[END_REF]. The second approach consists in using non-contact manipulations as laser trapping [START_REF] Dholakia | Tutorial reviewoptical micromanipulation[END_REF] and dielectrophoresis [START_REF] Hunt | Dielectrophoresis tweezers for single cell manipulation[END_REF]. These principles are not disturbed by adhesion but the blocking force remains low.

In this paper we simulate the 3D behavior of micro particles under DEP force. This force is calculated in function of the electric potential applied on the electrodes. Using the vision capture, the control law of the system can be studied and a first step to a closed loop is proposed. The typical problem at this scale is the large difference between the high dynamics of the system (response time ≃ 1ms) and the low speed rate of the vision capture (≃ 1 image per 10 ms).

Thus, we propose to design a predictive control strategy based on the feedback of the vision sensor and a model of the DEP force. The introduction of a model will enable to control the object behavior despite the low speed rate of the vision capture. This paper is focused on the presentation of the model.

Because of the large variety of electrodes and their geometric complexity, it could be highly difficult to directly integrate analytic equations. Our method is based on merging preprocessed FEM simulations and analytic equations. Physical equations are computed in order to define the link between the electric potential and the DEP force.

Consequently, based on the database built by several preparatory FEM simulations, our simulator is able to provide the trajectory whatever the electric potentials applied on the electrodes (see Fig. 1). This approach enables to improve the simulation's time. Each iteration consists of the use of the equations linking the electric potential to the electric field and then to the DEP force. FEM software can simulate complex geometries but the simulation time remains high especially when the potential changes frequently. By using FEM simulator as a preprocessing simulation and integrating physical law in a specific simulator, the time of simulation can be highly reduced. In the next section, a numerical implementation of how to calculate the electric field in function of the electric potential applied on the electrodes are discussed. In the third section the principle of the dynamic model is presented and in the fourth section, simulations and experiment results are shown and discussed.

II. CONCEPT OF HYBRID SIMULATION

A. General principle

The general expression of the dielectrophoretic force applied to a micro object considered as a punctual point [START_REF] Chang | Multiple electrodes arrayeed dielectrophoretic chip with application on micro-bead manipulation[END_REF], [START_REF] Gascoyne | Dielectrophoresis-based sample handling in general-purpose programmable diagnostic instruments[END_REF] is:

---→

F DEP = 2πǫ 0 ǫ p r 3 Re[K(ω)] --→ ∇E 2 , (1)
where K(ω) is the Clausius -Mossotti factor:

K(ω) = ǫ * p -ǫ * m ǫ * p + 2ǫ * m , (2)
and

ǫ * = ǫ + σ jω , (3)
where ǫ are the permittivities, σ are the conductivities, index 0 refers to the vacuum, index m refers to the medium and index p refers to the micro object, r is the radius of the bead, ω is the angular frequency of the applied electric field, -→ ∇ is the gradient operator and E is the root mean square magnitude of the sinusoidal electric field. This paper focuses on the simulation of the trajectory of a spherical object in a non phased sinusoidal electric field where the DEP torque can be neglected. The DEP force is thus function of the magnitude of the electric field -→ E , and of the frequency ω of the electric field which modifies the value of the Clausisus -Mossotti factor [START_REF] Gauthier | Forces analysis for micromanipulations in dry and liquid media[END_REF]. The first challenge is to determine the spatial electric field -→ E according to the potentials applied on the electrodes (see Fig. 2). Let us consider the magnitude of the sinusoidal electric potentials applied on n+1 electrodes:

[ϑ 0 , ϑ 1 , • • • , ϑ n ] . (4)
The 3D calculation of the electric field -→ E with respect to the potentials ( 4) is quite complex. Analytic calculations is in the general case not possible because of its too high complexity. To get over this difficulty we use numeric simulations. The proposed idea consists in simulating the 3D behavior of micro beads under DEP force with the minimum of numeric simulations. Firstly, we consider n potential differences:

U = [U 1 = ϑ 1 -ϑ 0 , • • • , U n = ϑ n -ϑ 0 ] . (5)
The first goal is to find the relation between the electric field -→ E in a point M (x, y, z) and the potential difference U applied on these electrodes. This relation can be written as follows:

-

→ E (x, y, z) = -→ f E (U ). (6)

B. Numerical implementation

Different methods exist to numerically calculate the electric field. A tridimensional meshing of the space and a bidimensional meshing of the electrodes have been considered. The difference between this 2 methods is analyzed in [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF] and the relevance of electrode meshing is shown. Our study is consequently based on the bidimensional electrode meshing. In this case, the electrodes are meshed along two directions x and y using respectively e x and e y as the number of elementary steps through each direction (see Fig. 2). The electric field -→ E is the sum of contributions of the elementary charges Q r,s on the electrode surface at point P (r, s) (see Fig. 3). The electric field at point M (x, y, z) is thus:

-→ E (x, y, z) = 1 4πǫ m ex r=1 ey s=1 Q r,s . ------------→ P (r, s)M (x, y, z) ------------→ P (r, s)M (x, y, z) 3 ds, (7)
where ds is the elementary surface on the meshed electrode.

To calculate Q r,s , the superposition principle on the potential is applied:

Q r,s = n m=1 C r,s,m .U m , (8)
where C r,s,m represents the elementary capacitances which can be identified during preprocessing. Thus the relation between the electric field -→ E and the potential differences U applied on the electrodes is:

-→ E (x, y, z) = 1 4πǫ m ex r=1 ey s=1 ------------→ P (r, s)M (x, y, z) ------------→ P (r, s)M (x, y, z) 3 . n m=1 C r,s,m .U m ds. (9)
The C r,s,m matrix contains the elementary capacitances created by each potential difference. The electrodes are meshed into 2D (e x , e y ) elementary points on the surface of each electrodes and the elements of this matrix are calculated using 3D simulation. The equation ( 9) enables calculating the electric field vector -→ E in any point M (x, y, z) in the workspace (see Fig. 2).

C. Preprocessing

The capacitance matrix C defined in ( 8) is constant and is independent from the electric potential U and the charge density Q. To identify this matrix, n FEM simulations of the charge density on the electrodes in function of n algebraically independent vectors U are needed. Once the matrix C is calculated the relation between the electric field -→ E and the electric potential U is obtained through equation [START_REF] Gauthier | Submerged robotic micromanipulation and dielectrophoretic micro-object release[END_REF].

D. Static force model

The diagram presented in the Fig. 4 shows the procedure to calculate the DEP force in function of the potential difference U . During preprocessing, geometries, boundaries conditions and surroundings properties are specified in a FEM software and n simulations are launched. All data are arranged to create the capacitance matrix C. The link between the electric potential U and the electric field -→ E is then calculated. Thus the DEP force ---→ F DEP applied to the micro object can be obtained by using [START_REF] Lambert | Capillary Forces in Micro-assembly[END_REF]. The frequency ω of the electric field can also be modified during the trajectory simulation.

III. DYNAMIC TRAJECTORY MODELING

Let us consider a micro particle immersed in a liquid medium with a dynamic viscosity µ, in a non uniform electric field created by applying electric potential on the electrodes. The forces applied on this particle are ---→ F DEP , its own weight -→ P and the Stockes drag force ---→ F drag which verifies:

---→ F drag = -6πµr -→ V = -k µ -→ V (10)
where -→ V is the velocity of the particle. Using Newton's second law the particle's motion is defined by:

---→ where M is the mass of the particles and -→ a is the acceleration vector. In order to simulate this trajectory model, we have proved that the inertia M -→ a is a negligible volumic effect in the micro-world compared to the force [START_REF] Kharboutly | Modeling the trajectory of a micro particle in a dielectrophoresis device[END_REF]. Despite the fact that the weight is a volumic effect it can not be neglected. Indeed the equilibrium state is characterized by ---→ F DEP = -→ P , and without considering the weight, the equilibrium position is always toward z ∼ ∞. But, in fact, experimentally the micro particle stands on a finite and stable position.

F DEP + -→ P -k µ -→ V = M -→ a (11)

A. Dynamic model

As the inertia can be neglected, the equation [START_REF] Dholakia | Tutorial reviewoptical micromanipulation[END_REF] gives the velocity of the micro particle along the trajectory:

-→ V = ( ---→ F DEP + -→ P ) k µ . (12)
In the Fig. 5, a block diagram is presented to illustrate each part of the system acting between the electric potential U , the DEP force and the micro particle's position. An applied potential difference U on the electrodes creates the non uniform electric field -→ E which creates the DEP force used to manipulate the micro particle. The equation ( 12) manages the dynamical behavior of the micro particle under DEP force.

IV. SIMULATIONS AND EXPERIMENTS

Considering the geometry of the electrodes described in Fig. 6, 4 electrodes and 3 potential differences are used. After 3 simulations under the FEM software COMSOL, the matrix C is identified as described below, and the simulator is launched. For a silicon bead of radius equal to 50µm, submerged in the water, its motion under DEP is simulated from an initial position (x = 10µm, y = -10µm, z = 55µm) as it is proposed in Fig. 6. Applying the potential difference U (80V, 0V, 80V) from equation ( 5) between the electrodes using a frequency of 40KHz the Clausius -Mossotti factor k(ω) referred in (2) is equal to K(ω) = -0.4252, the DEP force is negative and the trajectory of the bead can be predicted.

A. Trajectory of a point-shaped object

If we consider the object as a point shaped object, there is no presence of the equilibrium state along the z axis (see Fig. 8). In the first phase of the motion, the bead goes Lateral view: The micro bead is in direct contact with the electrodes plane (Left view) when applying a tension of U (60V, 0V, 60V). The micro particle is in levitation up to the electrodes (Right view) when applying a tension of U (80V, 0V, 80V).

up when it moves toward the z axis. After reaching the equilibrium point in the XY plane (see Fig. 7), the bead starts the sedimentation phase and goes down. This motion is caused by the calculated DEP force, which is equal to zero on the micro particle's center. However, experiments performed on this geometry, using electrodes made from gold and deposit on a silicon wafer immersed in an ultra pure water and using a glass bead of radius equal to 50µm as a micro object, show that an equilibrium point exists in the center of the electrodes (see in Fig. 9) along the z axis (see in Fig. 10). In this figure, the micro bead and its reflexion on the wafer are shown on a lateral view. In the left view, the micro bead is on the wafer when an electric potential of U (60V, 0V, 60V) is applied. In this case, the z component of the DEP force is smaller than the weight of the micro bead. In the right view, the micro bead is in levitation above the electrodes plane when an electric potential of U (80V, 0V, 80V) is applied. This difference between the experimental result and the simulated result appears because the hypothesis used in (1) (micro particle is considered as a punctual point) is not available in the equilibrium point. This assemption induces a zero DEP force while the real force is not zero.

B. Trajectory of a volumic object

To solve this problem, the micro particle is divided into 4 identical parts around its vertical axes and then all forces applied on each part are calculated to find the general force by summing all forces together (see Fig. 11). Fig. 12 and Fig. 13 show the trajectory of the micro particle using the same initial conditions as above. Using this method, the simulated trajectory of the micro particle convinces more with the real trajectory and the simulation results show that the micro particle reaches an equilibrium point as observed in the experimental results.

Using our simulator, non linear trajectories can be simulated. Multi initials positions of a micro bead of radius equal to 10µm, under the same potential (U = 80V ), could be simulated and the projection on xy plane of the resulted trajectories are presented in figure 14.

C. Experimental validation

Using another geometry of electrodes described in 15, the static position of the micro particle of radius equal to 45µm, immersed in an ultra pure water, is calculated in function of the applied tension on the electrodes in both simulation and experiment. The applied tension is U (90V, 0V, U 3 ). When changing U 3 from 30V to 90V the static position of the x position (m) Experiment Simulation Fig. 16. static position of the micro particle in function of the applied tension on the electrodes micro particle changes. The Fig. 16 shows the static position of the simulated calculation and experimental position. In this figure, it is evident to find that when U 3 = 90V , the static final position is in the center of the electrodes. The difference between the experiment result ans the simulation is due to, firstly the deformed electrodes in the experiments and secondly the transformation from pixel to µm of the bead position in experiments.

V. DISCUSSION

Using a calculator with relatively high performance (Intel Xeon CPU 1.60 GHz (2 CPUs) with 3 GB of RAM) the measured time to simulate the electric field and calculate the DEP force applied on a particle using the electrodes described in Fig. 9 and its trajectory under the FEM software COMSOL 3.5 is 2 minutes.

Using the same calculator to simulate the charge density on the electrodes is 30 seconds. To create the base used in our simulator we need 3 FEM simulations witch means 1 minute 30 seconds.

Using this database, the time needed to calculate the trajectory is 10 seconds.

Thus comparing both methods, the time needed to simulate n trajectories of a micro object by applying n different tensions on the electrodes is:

1) For the FEM software: 2 × n minutes.

2) For our simulator: 1 minute 30 seconds + 11 × n seconds. Consequently, our approach, based on the use of FEM as a preprocessing, enables to clearly reduce the time of computation especially when different values of tension U have to be simulated on the same geometry. However, in this simulator, the calculated electric field is not limited in the space, it can be calculated in any point in the space contrary to the FEM software, where the space is limited to the meshed volume. Moreover this simulator does not take into account the DEP torque, but the general principle proposed here can be extended to angular positioning. This aspect will be explored in future works.

VI. CONCLUSION

We have proposed a new model and simulator using a hybrid method in which we combined both analytic and numeric calculation to simulate the behavior of a micro bead (spherical object) under DEP force. A pre-simulation under FEM software allows to obtain a database with respect to the electrode geometries. Based on the preprocessing, the simulator integrates the dynamical behavior equation to simulate the electric field under different tensions applied on the electrodes. An example of 4 electrodes have been discussed and compared to experiments. Current works are focused on the use of this model to study control laws able to move the micro particle along a specified trajectory.

VII. ACKNOWLEDGMENTS

Fig. 1 .

 1 Fig. 1. General principle of classic FEM simulation and hybrid simulation

E

  Fig. 2. Calculating → E for the 2D charge density

Fig. 3 .Fig. 4 .

 34 Fig. 3. Elementary charges in each electrode

Fig. 5 .

 5 Fig. 5. Block diagram representing the dynamic modeling of a bead under DEF force

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 6789 Fig. 6. Geometry of 4 electrodes studied and the initial position of the micro particle in simulation

Fig. 11 .Fig. 12

 1112 Fig. 11. Use of the punctual bead (left), use of the volumic particle(right). Only z component of the electric field is presented

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig. 13. z trajectory of the volumic particle under DEP force

This work is supported by the French research projects PRONOMIA N: ANR 05 -BLAN -0325 -01 and NANOROL N: ANR 07 -ROBO -0003 -01.


Catherine Aaron

Density support and intrinsic dimension estimation based on a hierarchical delaunay-type simplicial complex

Let X1, ..., XN , Xi ∈ R D be an uniform drawn on a compact d-dimensional manifold S with d ≤ D. Here is suggested a new way to estimate both S and d. The method is based on the computation of a set of simplicial complexes (one for each dimension d ≤ D) and on an inductive criterion to select the "good" one. Each computed complex is a subcomplex of Delaunay's complex computed using k-nearest neighbors restriction and local P CA. A proposition for the k value is given in the first part and the algorithm is detailed in the second part.

Introduction 1.Intrinsic dimension estimation

Dimension estimation is a challenging problem that has many statistical applications in data analysis, the most obvious applications are dimension reduction methods (for instance isomap [START_REF] Silva | A global geometric framework for nonlinear dimensionality reduction[END_REF], LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], HLLE [START_REF] Donoho | Hessian eigenmaps : new locally linear embedding techniques for high-dimensional data[END_REF] or SOM algorithms [START_REF] Kohonen | Self-Organizing Maps[END_REF] all need the initial choice of a dimension), but dimension knowledge is also helpful for modeling problems as in time series [START_REF] Verleysen | Engineering Applications of Bio-Inspired Artificial Neural Networks, chapter Forecasting financial time series through intrinsic dimension estimation and non-linear data projection[END_REF] or regression [START_REF] Cadre | Dimension reduction for regression estimation with nearest neighbor method[END_REF]. Dimension estimation is also related to other mathematical fields as neural network, signal processing [START_REF] Hero | On local intrinsic dimension estimation and its applications[END_REF] and physics [START_REF] Almoga | A vectorized algorithm for correlation dimension estimation[END_REF] so giving an exhaustive bibliography is vain (for a complete review see Cutler's chapter of [START_REF] Ray | DIMENSION ESTIMATION AND MODELS[END_REF]). We are only going to deal here with the main ideas of intrinsic dimension estimation.

Mostly we can find two types of dimension estimation methods :

• Topology based dimension estimation : They are mainly based on the Hausdorf dimension simplification based on Grasberg and Proccacia's work [START_REF] Procaccia | Measuring the strangeness of strange attractors[END_REF]. Let us define :

C(r) = lim N →∞ 2 N (N -1) 1≤i≤j≤N
1 ||Xi-Xj ||<r d = lim r→0 lnC(r) ln(r) Since 1992, the limitation of such a method has been studied in [START_REF] Eckmann | Fundamental limitations for estimating dimensions and lyapounov exponents in dynamical systems[END_REF] (in this paper the needed N according to the dimension and diameter of the set is studied) but there exists two kinds of limitations with such a method : 1 -It is strictly based on uniformly randomized data (and there is little hope to remove this hypothesis)

-The first C(r) can't be computed (because it is only possible to observe a finite number of points) so if the d formula is applied, the estimated dimension will be 0 (the true dimension of a discrete set).

A good value for a small r value must be found or the method has to be improved (as in [START_REF]Intrinsic Dimension Estimation Using Packing Numbers[END_REF] or [START_REF] Bozhangc | Intrinsic dimensionestimationofmanifoldsbyincisingballs[END_REF]) for instance)

• Local PCA : A radically (and more statistical) approach is the local P CA's approach [START_REF] Olsen | An algorithm for finding intrinsic dimensionality of data[END_REF]. The idea is to compute tangent T space via local P CA and to observe the decreasing of the variance of the projection on T ⊥ . The main problem here is to define well the neighborhood of each point. The present work may be useful for this via the theoretical result in section 2.

Density support estimation and its topological properties

Density support estimation has a lot of application fields, for instance in medical diagnosis, machine condition monitoring, marketing and econometrics as noticed Biau and Pelletier in [START_REF] Cadre | Exact rates in density support estimation[END_REF]. The support density estimation using union of small balls centered on observation has been studied in [START_REF] Rodrguez-Casal | A nonparametric approach to the estimation of lengths and surface areas a nonparametric approach to the estimation of lengths and surface areas[END_REF], [START_REF] Cadre | Exact rates in density support estimation[END_REF], [START_REF] David | Asymptotic normality in density support estimation[END_REF]. Such a method has great asymptotical properties but a huge inconvenient : the topological properties of the estimated support may not be those of the "true" support : for instance the dimension of the estimated support is D the dimension of the embedding R D in which the support is and not the dimension of the support. There can be holes in the estimated support that do not exist in the true supports so the estimated homology groups might differ from the "true" one. The importance of the computation of the homology groups for application (and the methods for it) can be seen in [START_REF] Carlsson | Topology and data[END_REF] and asymptotical properties of an estimation of Betti numbers can be found in [START_REF] Meckes | Limit theorems for betti numbers of random simplicial complexes[END_REF]. In [START_REF] De Silva | On the local behavior of spaces of natural images[END_REF] it has been applied and a Klein bottle shape has been observed in a real data base.

The proposed method

This paper's aim is to manage to build a complex on the data to estimate the density support. This complex is expected to have the same dimension and the same topological invariant as the unknown true support of the density. Contrary to [START_REF] Zomorodian | Computing persistent homology[END_REF] the final result is a complex on which the homology can be computed (and not a set of complexes which implies to find "the persistent one").

Section 2 is dedicated to the theoretical search of a k to restrict Delaunay's complex by a k-nearest neighbors and not to create an undesirable hole (which could skew the computation of Betti Numbers). Section 3 presents an algorithm that computes a complex for each supposed dimension d using Delaunay's complex and local P CA.

Section 4 gives an indicator to choose a dimension and so a complex.

2

Finally section 5 presents some results on simultated data.

2 Majoration of the probability to suppress an inside edge by restricting Delaunay's complex by k-nearest neighbor 2.1 Introduction

let X 1 , ...X N be a sample on R D in a d dimensional submanifold S ⊂ R D (D ≥ d).
One of the goal of this paper is to build a "good" complex that links the points of the sample. S will thus be estimated as the union of all the simplexes.

We choose to use initially Delaunay's complex. Let us denote Delaunay's complex T . It satisfies :

• The D-dimensional simplexes of T are the (X i1 , ..., X iD+1 ) such as the hypersphere circumscribed to the points does not contain any X j

• The k-dimensional simplexes of T are the k-sub-simplexes of the D-dimensional simplexes

• The 1-dimensional simplexes of T will be called edges of the complex

• Property : Delaunay's complex gives an estimation of the convex hull of the sample

• Corollary : Since S is not convex it is necessary to remove some simplexes to get a good estimation of S using the complex (see figure 1)

This section focuses on the case D = d and the way to remove simplexes from T to estimate correctly S. We choose the well-known rule of constraining D by k-nearest neighbors : the removed simplexes will be those that contain at least one edge [A, B] such as B is not a k-nearest neigbor of A and A is not a k-nearest neighbor of B.

The choice of a "good" value for k is fundamental to get good results. k must be small enough to respect the local non-convexity of S but large enough to avoid unexpected "holes" (which may be a problem if using results to compute the homology group, for instance).

The following part of this section is dedicated to proving that the probability to remove an inside edge when D = d and the drawn is uniform is majored by :

P (N, k) ≤ a k-1 d (N -2) 3/2 π/2 + a N -2 d with : a d = 2 d -1 2 d
Unfortunately our definition of an inside edge is a little more restrictive than the creation of an undesirable hole but it gives an approximation (see appendix for a discussion about this definition). The constrained graph is plain and black ; Delaunay's complex is dashed and red. (c) : Delaunay's complexe contrained by 14-nearest neighbors. The constrained graph is plain and black ; Delaunay's complex is dashed and red. 9 undesirable holes can be observed

Theoretical Study

In all this part we will assume that X 1 , ..X N is a uniform drawn in R d on S a d-dimensional bounded manifold. Without loss of generality, we also will assume that V (S) the volume of the manifold is 1.

We will first start by two elementary lemmas : Lemma 1. Let us denote :

• B 1 = B(O, r) • X a point of the boundary of B 1 • O ′ such as d(O, O ′ ) = d(O ′ , X) = r ′ • B 2 = B(O ′ , r ′ ) • B ′ 1 = B 1 \ (B 1 ∩ B 2 ) then : V (B ′ 1 ) V (B 1 ∪ B 2 ) ≤ 1 - 1 2 d = a d Proof. Let first us define O ′′ ∈ [O, O ′ ] with d(O, O ′′ ) = r/2 and B 3 = B(O ′′ , r/2) (see figure 2). Let us denote c d = π d/2 /Γ(1 + d/2) (the volume of the unit d-ball). B 3 ⊂ B 1 ∩ B 2 ⇒ V (B ′ 1 ) ≤ c d r d -c d (r/2) d B 1 ⊂ B 1 ∪ B 2 ⇒ V (B 1 ∪ B 2 ) ≥ c d r d Lemma 2. S n = 1 n n-1 k 1 (1/n)(1 -k/n) < π O O ′ X B 1 B 2 B 3 O"
(x) = 1/ x(1 -x) defined on ]0, 1[.
We define g n stepwise by : 3) finishes to prove it. Hypothesis, definitions and notations

g n (x) = min{f (t), t ∈ [(k -1)/n, k/n]} when x ∈ [(k -1)/n, k/n] Obviously : g n (x) ≤ f (x) so 1 0 g n ≥ 1 0 f = π. S n + f (0.5)/n = 1 0 g n (see figure
• X = {X 1 , .., X N } is a uniform sample on a d-dimensional manifold S with V (S) = 1.
For all that follows Delaunay's complex is Delaunay's complex computed on X.

• X j is the k(X i , X j ) the neighbor of X i . • k * (X i , X j ) = min( k(X i , X j ), k(X j , X i )). • If t = [X i , X j ] is included in Delaunay's complex then k * (t) = k * (X i , X j ). • an edge t = [X i , X j ] of Delaunay's complexe is an inside edge if B(X i , d(X i , X j ))∪ B(X j , d(X i , X j )) ⊂ S (see appendix for a discussion of this notion). X 1 X 2 B 1 B 2 k -1 points 0 points N -1 -k points S
P (k * (t) = k) ≤ a k-1 d √ N -2 √ 2π (k -1)(N -k -1) exp 1 12(N -2)
.

Proof. Let us first assume that :

• t = [X 1 , X 2 ], • k(X 1 , X 2 ) = k.
As t is in Delaunay's complexe there exists an empty ball B 2 such as X 1 and X 2 are on the boundary of

B 2 . As k(X 1 , X 2 ) = k then B 1 = B(X 1 , d(X 1 , X 2 ))
contains k -1 other points (see figure 4).

Knowing B 2 , the probability to get such a configuration is :

P B2 (k * (t) = k) = C k-1 N -2 (V (B 1 \ (B 1 ∩ B 2 ))) k-1 (1 -V ((B 1 ∪ B 2 ) ∩ S)) N -1-k .
Let us denote x = V (B 1 ∪ B 2 ) and apply lemma 2 :

P B2 (k * (t) = k) ≤ C k-1 N -2 a k-1 d x k-1 (1 -x) N -1-k .
A maximisation of the expresion leads to :

P (k * (t) = k) ≤ a k-1 d C k-1 N -2 (k -1) k-1 (N -k -1) N -k-1 (N -2) N -2 .
Finally, we use the Stirling unequality (for k ≥ 2 and k ≤ N -2) to get :

P (k * (t) = k) ≤ a k-1 d √ N -2 √ 2π (k -1)(N -k -1) exp 1 12(N -2)
.

Theorem 1. Let X = (X 1 , ..., X N ) be a uniform sample on a d-dimensional manifold S with V (S) = 1, the probabilty P (N, k) that an inside edge of Delaunay's complex is suppressed when restricting by the k-nearest neighbors graph satisfies :

P (N, k) ≤ a k-1 d (N -2) 3/2 π/2exp 1 12(N -2) + a N -2 d .
Proof. An inside edge t is suppressed when restricting by the k-nearest neighbors graph if k * (t) ≥ k so :

P (N, k) = N -1 j=k P (k * (t) = j)
Lemma 3 gives

P (N, k) ≤ N -2 j=k a j-1 d √ N -2 √ 2π (j -1)(N -j -1) exp 1 12(N -2) + a N -2 d . P (N, k) ≤ a k-1 d N -2 j=k √ N -2 √ 2π (j -1)(N -j -1) exp 1 12(N -2) + a N -2 d
, and lemma 2 leads to the conclusion.

Corollary 1. The restriction of the delaunay complex by the k-nearest neighbor graph with :

k ≥ 1 + ln(ε) -3 2 ln(N -2) -1 2 ln(π/2) ln(a d ) = k 0 (N, d, ε)
creates an inside edge with a probability :

P ≤ εe 1 12N + a N -2 d ∼ ε 2.3 Numerical results

Graphical results

We prensent here 4 graphes, illustrating Delaunay's restriction to k 0 (N, 2, ε = 0.1) for the same holed squares as in figure 1 (N varying in {100, 200, 500, 1000}). First, it can be observed that there is no undesirable hole creation in the obtained complexes. In the first case N = 100 doesn't allow to obtain a complex that respects the topology of the initial set S : there is no hole at all. We think that this is due to the fact that the size of the hole is too small in regard to the density around the hole. The convergence of the computed complex to the density support seems here to occur (the convergence of the algorithm will be studied in a further paper).

Simulation on Spheres

Let K d,N (X 1 , ...X N ) be the maximum number of neighbors observed on inside edges of Delaunay's complex when the density support is d-dimensional and the sample is uniformly distributed on S. Our majoration implies that :

k 0 (N, d, ε) ≥ Q 1-ε (K d,N ) (with Q 1-ε (K d,N ) the (1 -ε) percentile of K d,N
). It is expected that this majoration is not too big. To compute K d,N (X 1 , ...X N ) without boundary effect and so avoiding the problem of looking for the inside edges, the X i has been uniformly randomized on d-dimensional ball S d which lies in R d+1 . The d-dimensional Delaunay's complex is here the boundary of the (d + 1)-complex. Results of our computation are presented in figure 6 for dimension 1 (1000 draws), 2 (1000 draws) and 3 (500 drawns) and each time N ∈ {100, 200, 300, 500, 1000, 2000}. Plain lines represent the "theoretical" k 0 for ε ∈ {0.1, 0.05, 0.01, 0.005} and dashed lines the simultated percentiles for same ε. It can be observed that :

• The majoration is verified and is not too big • The growth speed seems to be the right one (for the dimension 3 the Q 0.99 and Q 0.995 simulated values might not be good because they are only computed on 500 samples 3 Computation of a simplicial complexe on the data

Let X = (X 1 , ..., X N ), X i ∈ R D be a uniform draw on a compact d dimension manifold S with d ≤ D.
Let us denote T the D-dimensional Delaunay's complex of X. We will denote : T = p∈{1,..,D} i t p i with :

• the t D i are the D-dimensional Delaunay-simplexes : i.e. such as the D-ball circumscribed of t i is empty (doesn't contain any points of X),

• t k i is the face of a t k+1 j .

The goal is now to extract of T , T * a "good" subset of simplexes. By "good" we expect that :

• T * is a d-dimensional simplex with d the true dimension of S,
• T * gives a correct estimation for S.

For that we will compute D subcomplexes of T (one for each supposed dimension from 1 to D) and the choice of the final complex will be done afterwards.

Computation of the D-dimensional complex

According to section 2 the D-dimensional simplex is the restriction of T by the k 0 (N, D, ε)-nearest neighbor graph.

Computation of the p-dimensional complexes for p < D

The idea that leads us to the following algorithm is very simple : let us assume that S is a p-dimensional smooth manifold. Then locally the manifold and its tangent hyperplan are close and we are going to consider some simplexes of Delaunay's compexe of the local projection on tangent hyperplan. More precisally :

for each point i :

• 1-search V i = {X j1(i) , ..., X j k(i) (i) } a neighborhood of X i
-Practically : according to section 2, we propose to use the k 0 (N, d, ε)-nearest neighbors of X i as V i .

• 2-H (p) i the hyperplan tangent to S at the point X i according to the hypothesis that it is p-dimensional:

-Practically : using local P CA (P CA on V i ).

• 3-W i is the set of all the projections of the V i points on H (p) i .

• 4-compute Delaunay's complex of W i : T i .

• 5-keep the p-dimensional simplexes of T i that satisfy the following properties :

i is in the simplex, -the simplex is in the set of Delaunay's simplexes T .

The final p-dimensional complex will be the union of all the p-dimensional simplexes kept in the algorithm.

Choice of the intrinsic dimension and complex

It is now needed to choose a complex in the set of all computed complexes (one for each supposed dimension). We suggest here two ways to choose. The first one is only classical local P CA method and the second one is based on the geodesic distances. The local P CA method is not linked to the simplicial approach and can be used beforehand to reduce the set of tested dimension (which can be usueful to reduce the computational time). The second method depends on the simplicial approaches and can be used afterwards to confirm the choice.

Local PCA

The local PCA method for estimation of the dimension is well known [START_REF] Olsen | An algorithm for finding intrinsic dimensionality of data[END_REF]. We here only apply it using the neighborhood according to section 2. Each tested dimension d is associated to a k 0 (N, d, ε)-neighborhood and for each point of the sample the eigenvalues of a P CA on its neighborhood can be computed. Boxplotting these eigenvalues can help to choose possible dimensions

Geodesic distance recognition

Let us suppose that the true dimension is d, then the k 0 (N, d, ε)-nearest neighbor graph and the d-dimensional complex's graph obtained with our method might both be used to compute approximation of the geodesic distance correctly and so geodesic distances computed using the two graphs might be close.

For each dimension, the plot of the geodesic distances computed with nearestneighbor method and simplicial method will be plotted and a dimension has to be chosen between dimensions that lead to a plot near the bisector.

5 Some examples

Examples for dimension D = 2

We computed here two examples for D = 2 (N = 500 and ε = 0.01) : the holed square and the circle. For the holed square, let us first look at the local P CA results : they both indicate a dimension 2. Looking at the geodesic distances scatterplots also indicates a dimension 2. The complex associated to the dimension 2 respects the topology of the (known) density support.

For the circle example, the local P CA method and the simplicial method also agree to decide for a dimension 1 (if the hole observed for the dimension 2 had not existed, only the PCA method would have been discriminant to conclude).

Examples for dimension D = 3

We computed here three examples for D = 3 (N = 500 and ε = 0.01) : the cylinder, the sphere and the spiral.

The sample on the cylinder has been realized as follows : θ U[0, 2π] r U[0.7, 1] and X 3 U[0, 1]. X 1 = cos(θ), X 2 = sin(θ) and X 3 . Local P CA hesitates between dimension 2 and 3 (the third eigenvalues is realy small compared to the two first ones, we can wonder if it is due to a "thin" 3-dimensional For the sphere example, the geodesic criterion hesitates between dimension 2 and 3 ; the local P CA helps to conclude to dimension 2.

Conclusions : For all the presented examples (dimension D = 2 and D = 3) the topological properties of the chosen complexes are the good ones. The choice of the elected complex has to be done according to both local P CA and geodesic approach.

Conclusions and persepectives

Our theoretical value for k is coherent with the known good properties that are expected for the choice of k, the k-nearest neighborhood (generally expected to satisfy k → ∞ and k/N → 0 ( [START_REF] Crou | On the rate of convergence of the bagged nearest neighbor estimate[END_REF]). Fortunately, our majoration in section 2 that can be considered as a quite strong majoration is not so far from the simulated value for k. Applications using this value gives quite good results for not Two main furter theoritical axes can be envisaged :

• On k value to avoid undesirable holes : can we improve the definition of an inside edge to be closer to the hole creation ? Can we prove that 1.5ln(N )/ln(a d ) speed is really the good one ? Can our result be extended to non-uniform samples ?

• On the S estimation : does our ŜN converge ? do Betty number's estimation converges ? what is the speed ? Graphical results of section 5 are encouraging. The intuition is that it is required that N > c(S)k 0 with c(S) a constant reflecting the complexity of S.

For the applied part : can the complex building be improved ? can it be adapted to sets S where the local intrinsic dimension is not constant [START_REF] Hero | On local intrinsic dimension estimation and its applications[END_REF]?

Figure 1 :

 1 Figure 1: 3 complexes computed on an holed squared sample (500 points U([0, 3] 2 \ [1, 2] 2 )). (a) : Delaunay's complexe. (b) : Delaunay's complexe contrained by 42-nearest neighbors.The constrained graph is plain and black ; Delaunay's complex is dashed and red. (c) : Delaunay's complexe contrained by 14-nearest neighbors. The constrained graph is plain and black ; Delaunay's complex is dashed and red. 9 undesirable holes can be observed

Figure 2 :

 2 Figure 2: Illustration for lemma 1. Remark for d > 2 this graphic corresponds to the projection in the plane containing O, O ′ and X

Figure 3 :

 3 Figure 3: Illustration for lemnm 2 with n even and n odd. The red line represents the function, the step function g n and the grey area's surface S n

Figure 4 :Lemma 3 .

 43 Figure 4: Illustration for lemna 3

Figure 5 :

 5 Figure 5: ε = 0.1, N = 100, N = 200, N = 500, N = 1000 : black edges are those of the restricted Delaunay's complex, dashed red ones are those of the Delaunay's complex. It can be observed that the removed simplex are on the boundary and that the topological properties of the estimated support are the same as the true support since N ≥ 200

2. 5 . 1

 51 Some values for k 0 and ε = 10 -2 Some numerical values for k 0 are presented in figure7: the empty cells correspond to k 0 (N ) < N . Cells in Italics correspond to k 0 < 2N and bold k 0 > 10N

Figure 6 :

 6 Figure 6: Theoretical (plain) and simulated (dashed) values for k and different values for ε for 1, 2 and 3.

Figure 7 :

 7 Figure 7: Some Values for k 0 (N, d, 10 -2 ).

Figure 9 :

 9 Figure 8: example of the holed square

  Figure 10: example of the cylindre

Figure 11 :

 11 Figure 11: example of the ball Finally for the spiral example it is only the local P CA approach that helps to conclude.

Figure 12 :

 12 Figure 12: example of the spiral

Figure 13 :

 13 Figure 13: Illustration for inside, semi-inside and not-inside edges

A Discussion of the inside edge notion

The aim of this appendix is to show that not removing an inside edge (that allows to get the theorem 1) is not so far from not creating an undesirable hole (that is the really interesting notion). The discussion here is not (at all) a proof but helps to understand why, at least for dimension 2, the two notions are close.

Let us remember the definition :

And Let us add two other definitions :

Illustrations for these 3 possible cases can be seen in figure 13. Our proof in section 2 can easily be adapted for semi-inside edges because the volume majoration in lemma 1 can be adapted (see figure 13) and because there is at least k point in each presented ball (but the writing is a little more difficult).

Let us now focus on the not-inside edge case and the dimension 2. Let us assume that t = [X 1 , X 2 ] is a not-inside, Delaunay's edge, that satisfies k * (t) ≥ k. and that removing t create an undesirable hole. As there is a creation of an undesirable hole t is not on the boundary of the complex. So there are two Delaunay triangles [X 1 , X 2 , X 3 ], and [X 1 , X 2 , X ′ 3 ]. Let us assume that [X 1 , X 2 , X 3 ] is the one "closest to the boundary" (see figure 14). [X 1 , X 2 , X 3 ] is removed by our algorithm so it is not on the boundary of the complex (otherwise we will not create a hole). Iterating such reasoning leads to the fact that there exists X 4 and X 5 as on figure 14 which satisfies [X 4 , X 5 ] is an edge of Delaunay's complex which is not removed by restrinction to k-nearest neighbors. As the construction implies that d(X 4 , X 5 ) > d(X 1 , X 2 ) the fact that k * ([X 4 , X 5 ]) < k * ([X 1 , X 2 ]) may be small.


Roba El Khoury

LIST, Moussa Cea

Mathieu Grossard

Mehdi Boukallel

Nicolas Chaillet

Arnaud Hubert

Observation-Oriented Design of a Monolithic Piezoelectric Microactuator with Optimally Integrated Sensor

This paper presents an evolution in an optimization method, called FlexIn, developed for the optimal design of piezoelectric compliant smart structures. FlexIn is based on a flexible building block method that uses a genetic algorithm approach, to optimize a truss-like planar structure made of piezoelectric passive, active and, with the work reported in this paper, sensitive building blocks. The model of these blocks is established by means of a finite-element electromechanical formulation. The main contribution of this paper is to present a new observation-oriented criterion, along with a static electromechanical one, considered in the optimization procedure for the optimal placement of piezoactuating and piezosensing parts in the compliant micro-structure. In order to underline the interests of such a criterion, performances of three pseudo-optimal piezoelectric smart structures are drawn. Their analysis and comparison illustrate the role of the optimization method and the observation-oriented criterion, in the design of smart structures and in simplifying their control afterwards.

Introduction

In many robotic applications, the use of compliant smart mechanisms is being widely propagated. In particular, when considering small-scale systems, e.g. for microrobotic operations like micro or nano-positioning, microgripping, etc., there are many advantages of compliant smart mechanisms, such as simplified manufacturing, reduced assembly costs, reduced kinematic noise, high precision , and miniaturization by actuation and sensing functions integration. One type of smart material typically used in compliant structures design is the piezoceramic PZT. An interesting advantage of such material is the reversibility of its electromechanical coupling effect, which explains its potential use in microrobotic applications as actuator and/or sensor [START_REF] Lin | Diagnostic lamb in an integrated piezoelectric sensor/actuator plate: analytical and experimental studies[END_REF], [START_REF] Campolo | Development of piezoelectric bending actuators with embedded piezoelectric sensors for micromechanical flapping mechanisms[END_REF]. Few studies consider the topology optimization (shape) of monolithic PZT active structures [START_REF] Silva | Design of piezoelectric transducers using topology optimization[END_REF] with exteroceptive sensors [START_REF] Grossard | Mechanical and control-oriented design of a monolithic microgripper using a new topological optimization method[END_REF]. Furthermore, a number of papers only address the problem of optimally designing PZT active structures with optimally integrated sensors [START_REF] Leleu | Piezoelectric Actuators and Sensors Location for Active Control of Flexible Structures[END_REF], [START_REF] Collet | Active modal control with piezoelectric layers optimisation[END_REF], [START_REF] Han | Optimal placement of piezoelectric sensors and actuators for vibration control of composite plate using genetic algorithms[END_REF]. In each one of these studies aiming for vibration control, the optimal placements of piezoelectric actuator and sensor on a flexible plate, based on controllability and observability criteria respectively, are achieved separately. The resulting modes' observability and controllability could disagree, thus eliminating the possibility of integrating both smart functions in the structure. Therefore, we developed two criteria that can help the designer to optimally place the piezoelectric parts simultaneously in the compliant structure. The first one considers purely static electromechanical aspect, i.e. the maximization of the amount of electric charges induced by the piezoelectric sensor's deformation, while the second one is based on modal observability aspects. This last criterion is a useful tool to optimize dynamic operating flexible actuators with integrated piezoelectric sensors, and to ensure their efficient control afterwards. It is based on the structure's balanced-gramian and observability gramian. From these two indicative characteristics, the criterion aims to optimize the observability of the 'k' first modes specified by the designer, under a condition that the control authority of those same modes is highly dominant comparing to the residual ones. Both criteria are integrated in a more global systematic design approach, based on topology optimization of the structure, as well as that of its frequency responses, to design compliant integrated smart micromechanisms. This method is based on the flexible building block method, called FlexIn ("Flexible Innovation") [START_REF] Bernardoni | A new compliant mechanism design methodology based on flexible building blocks[END_REF] and [START_REF] Grossard | Mechanical and control-oriented design of a monolithic microgripper using a new topological optimization method[END_REF], which uses an evolutionary approach, to optimize a trusslike planar structure including passive, active and now sensitive building blocks made of piezoelectric material. This paper is organized as follows: we will first briefly review the underlying idea of the FlexIn methodology for the design of smart compliant mechanisms. The electromechanical approach, based on a finite element (FE) formu-lation, is established for the model of the piezoactuating and piezosensing building blocks. From this computation, the piezosensing model leads to the first static criterion described earlier. In a second part, the state model approach used in FlexIn is presented, where the mechanism is discretized on its modal components. This particular representation simplifies the computation of a gramianbased criterion, taking into account simultaneously, controllability aspects for actuators and observability aspects for sensors in flexible structures. Hence, the two criteria ensure static and dynamic performances, required to optimally synthesize a compliant monolithic microactuator with integrated piezoelectric sensor. Optimization results are presented in the last part. Performance comparisons between three of the obtained pseudo-optimal structures exhibit the interests of the proposed optimization method and obervability-oriented criterion, for the design of smart structures.

2 FlexIn: A Compliant Smart Mechanisms Stochastic Design Methodology

Principles of the method

FlexIn is an optimization software for the design of planar compliant micromechanisms. It uses a multiobjective genetic algorithm approach, which consists of searching for an optimal distribution of allowed building blocks, chosen in three specific libraries, composed of 36 passive blocks, 19 piezoactuating blocks, and 19 piezosensing blocks, made of beams assembly (Fig. 1). In addition to topological specifications, the optimization problem appoints an optimal set of boundary conditions (fixed frame location, contacts, actuators, sensors, end-effectors, etc.), dimensions and materials, based on the optimization criteria selected by the designer. More detailed descriptions of the method can be found in [START_REF] Bernardoni | A new compliant mechanism design methodology based on flexible building blocks[END_REF] and [START_REF] Grossard | Mechanical and control-oriented design of a monolithic microgripper using a new topological optimization method[END_REF].

Electromechanical FE model of the piezoelectric structures

In the optimization procedure, the computation of different criteria requires the FE model of each block of the libraries. To obtain the FE formulation of the piezoelectric blocks, a model of a piezoelectric beam is needed first, exploiting both direct and inverse piezoelectricity effects separately for sensing and actuation purposes respectively (Fig. 2). To do so, the compliant mechanism is assumed to undergo structural planar deformation, leading to considering Navier-Bernoulli beams type FEs. Figure 2: Thickness-polarized piezoelectric beam transducer with electroded surfaces, and orientation in the material reference frame (e 1 , e 2 , e 3 ). (ϕ 1 , ϕ 2 ) and (q 1 , q 2 ) denote respectively the electric potential for actuation case, and the electric charges for sensing case, of the two electrodes.

FE formulation of the piezoelectric beam

From modified Hamilton's principle for a general electromechanical system [START_REF] Preumont | Mechatronics: Dynamics of Electromechanical and Piezoelectric Systems (Solid Me-chanics and Its Applications)[END_REF], piezoelectric beam model is established through two equations representing respectively the actuator and sensor functioning modes, as follow:

M b ηb + K b η b = G b Φ b + Fr b q b = G t b η b + C b Φ b (1)
where

η b = (u A , v A , ω A , u B , v B , ω B ) t
Rp is the nodal displacement vector in the beam coordinate system R p = (A, x p , y p , z p ) (see Fig. 3). M b , K b , G b and C b are respectively, the mass, stiffness, electromechanical coupling and electric capacity beam matrices. Φ b = [ϕ 1 ϕ 2 ] t and q b = [q 1 q 2 ] t are the vectors representing the electric potentials and the electric charges respectively on the upper and lower faces of the piezoelectric beam (see Fig. 2). The forces vector Fr

b = R x A , R y A , H z A , R x B , R y B , H z B t
Rp is due to the external mechanical work (see Fig. 3). Corresponding matrices are expressed as follows:

M b = ρA          L 3 0 0 L 6 0 0 13L 35 11L 2 210 0 9L 70 -13L 2 420 L 3 105 0 13L 2 420 -L 3 140 L 3 0 0 (sym.) 13L 35 -11L 2 210 L 3 105          Rp (2) K b = Y         A L 0 0 -A L 0 0 12I L 3 6I L 2 0 -12I L 3 6I L 2 4I L 0 -6I L 2 2I L A L 0 0 (sym.) 12I L 3 -6I L 2 4I L         Rp (3) G b = Y hd 31 1 0 0 -1 0 0 -1 0 0 1 0 0 t Rp (4)
C b = hLε 33 e p 1 -1 -1 1 Rp ( 5
)
where

I = eph 3
12 is the inertia moment of the beam cross section A = e p h, ρ the density of the beam, Y = its Young's modulus, and ε33 its modified electric permittivity, function of the piezoelectric material properties. The second equation in [START_REF] Lin | Diagnostic lamb in an integrated piezoelectric sensor/actuator plate: analytical and experimental studies[END_REF] represents the electric charges generated on the sensing beam electrodes via the direct effect of piezoelectricity. To measure these charges experimentally, charge-voltage converters are usually required using operational amplifier (see Fig. 4). In ideal mode, upper and lower electrodes of the piezoelectric beam are considered to be shorted-circuit, so that zero-voltage Φ b is applied in [START_REF] Lin | Diagnostic lamb in an integrated piezoelectric sensor/actuator plate: analytical and experimental studies[END_REF]. This assumption ensures that sensing and actuation functions of piezoelectric beams are decoupled.

FE model of piezoelectric structures

The piezoelectric blocks' matrices are calculated by the association of corresponding beam matrices in the global coordinate system of the structure. Each flexible structure is then defined as a finite-dimensional linear system modelled by:

M g ηg + K g η g = E g u y co = δ = F g η g y ob = q g = L g η g (6)
The foregoing second-order differential matrix equations represent the undamped dynamic behaviour of such a system, where M g and K g are the structure's mass and stiffness matrices respectively, obtained by the assembly of the matrices of all the blocks constituting the structure. Considering the integers p, s, and r, as respectively the numbers of DOFs of the structure, number of inputs (i.e. actuators), and number of tip displacement outputs, η g is then the p × 1 nodal displacement vector and u is the s × 1 input vector. The p × s input matrix E g reflects the location of the actuated DOFs, while y co is the r × 1 controlled output vector representing the output tip displacement δ through r × p output displacement matrix F g . The third equation expresses the electric charges (y ob = q g ) obtained by the integrated sensing function from the direct piezoelectric effect. It will serve as a static criterion maximizing the amount of sensing electric charges (see criterion J 2 in section 4.1). Note that L g is the 1 × p single output matrix indicating the placement of piezoelectric sensor in the structure. Furthermore, it is important to note that controlled output variable (y co ) is not the observed output one (y ob ), as in usual other microrobotic systems, where the controlled tip of the piezocantilever is observed through external sensor.

A New Criterion for Matching Observability and High Control Authority Optimization

In order to successfully achieve suitable dynamic openloop performances for further closed-loop control (see section 3.3), an optimal topology design strategy is derived taking into account control-observability criterion. The latter based on modal balanced gramians and observability gramians interpretations will be defined to help optimally integrating actuators and sensors in the microstructure. To do so, the physical coordinate base representation ( 6) is firstly transformed in the modal base to display the flexible modes. Then, we design an observer that will estimate the flexible modes included in the state vector through the measured electric charge vector q g on the sensing blocks of our optimal device. Finally we propose to control the output tip displacement δ through a state feedback corrector.

Modal representation of flexible structures

The harmonic solutions of the first equation in [START_REF] Collet | Active modal control with piezoelectric layers optimisation[END_REF] give the eigenvectors matrix Ψ and natural frequencies ω i of the system. Details of the modal representation computation is given in [START_REF] Grossard | Mechanical and control-oriented design of a monolithic microgripper using a new topological optimization method[END_REF]. We then obtain:

z + diag(2ξ i ω i )ż + diag(ω 2 i )z = Ψ t E g u y co = δ = F g Ψz y ob = q g = L g Ψz (7)
where z is the p × 1 modal displacement vector, and ξ i is the i th modal damping ratio introduced using Basil's hypothesis. One interesting 2p × 1 state vector x typically used for flexible structures, and whose advantages are revealed in [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF], consists of modal velocities and frequency-weighted modal displacements:

x = ż1 ω 1 z 1 . . . żp ω p z p t (8)
Since controlled and observed output vectors are not the same, the modal state-space representation can be written as follows:

ẋ = Ax + Bu, y co = C co x, y ob = C ob x. (9)
which leads to two matrices triplets (A, B, C co ) and (A, B, C ob ) designating the control and observation state space models respectively. The matrices take the forms A = diag(A 1 , . . . A p ), B = (B t 1 , . . . B t p ) t , C co = (C co1 , . . . C cop ), and C ob = (C ob1 , . . . C obp ), with, for i = 1, ..., p,

A i = -2ζ i ω i -ω i ω i 0 (10)
B i = b i 0 (11)
C co i = 0 cco i ωi (12)
C ob i = 0 c ob i ωi [START_REF]PI Piezo Ceramic Technology[END_REF] where b i , c co i , and c ob i are the i th components of Ψ t E g , F g Ψ and L g Ψ respectively. Matrix A revolves around the structural parameters (eigen frequencies and damping ratio), whereas matrix B depends on the location of actuated DOF, matrix C co on the location of desired displacement output, and matrix C ob on the location of integrated piezoelectric sensors. Fig. 5 shows the control-observation diagram principle to be applied at the end on our system.

Computation of the observability and balanced gramians

Observability gramian (W ob ) between state x and measured output q g is found to be convenient to characterize the modes' observability by the mean of electric charge q g . Its energetic and geometric interpretations are demonstrated in [START_REF] Leleu | Piezoelectric Actuators and Sensors Location for Active Control of Flexible Structures[END_REF] and [START_REF] Jonckheere | Principal components analysis of flexible systems, Open-loop case[END_REF]. For stable A, W ob is obtained from algebraic solution of following Lyapunov equation:

A t W ob + W ob A + C t ob C ob = 0 (14)
Assuming that the damping ratios are infinitely small and the natural frequencies well spaced, which is widely accepted for flexible structures, the block diagonal forms of (A, C ob ) couple can be exploited to give closed-form analytical solution for the expression of the modal observability gramian [START_REF] Gawronski | Balanced Control of Flexible Structures[END_REF]. It is diagonal and equal to:

W ob = diag W ob11 , . . . , W obpp (15)
with, for i = 1, ..., p,

W ob ii = γ qii 4ξ i ω 3 i I 2 = α i I 2 ( 16
)
where γ qii = c t ob i c ob i , and I 2 is the 2 × 2 identity matrix. For a given mode (ξ i , ω i ), γ qii scalars represent the way the i th mode is seen through the piezoelectric sensor blocks. On the other hand, the observer should be able to reconstruct the dominant modes of δ/u transfer by measuring q g (see Fig. 6). These dominant modes are symbolized by high Hankel singular values (HSVs) defining the balanced gramian W eδ of (A, B, C co ) system as follows:

W cδ = W oδ = W eδ = diag (σ i ) ( 17
)
where σ i are the HSVs of the (A, B, C co ) system. Note that balanced gramian is a useful tool for quantifying the joint controllability and observability of the system. It is shown that when the damping ratios decrease to zero, the physical modal state coordinates are approximately balanced in this asymptotic situation, and the approximate i th Hankel singular value for flexible structure is given by [START_REF] Gawronski | Balanced Control of Flexible Structures[END_REF]

σ i = c t co i c co i b i b t i 4ξ i ω 2 i ( 18
)
HSV describes the degree of the corresponding modal state's input-output energy flow through the system.

A gramian-based criterion for optimal placement of piezoelectric sensor within a piezoelectric microactuator

Considering the statements made in Section 3.2, this kind of active structures which are to be finely controlled are confronted with two main issues: 1) A reduced model of the structure must be developed, which includes the few dominant low frequency modes, without destabilizing the system by rejecting the residual modes (i.e. high roll-off after the dominant modes).

2) If the dominant modes are not all observed, the reconstruction of δ will not be guarenteed in an optimal way by the observer. As mentioned earlier, each piezostructure suggested in the optimization procedure is evaluated according to the specified criteria. Hence, we developed a new criterion to help overcome these two difficulties: the first problem is overcome by forcing the optimal structure to have k first dominant modes of the δ/u transfer in order to reduce the model without the residual modes affecting it. Then, to surmount the second problem, we guarantee a high-level of observability of these dominant modes by means of electric charge q g . It is thus presented by a procedure formulated as follows:

If σ min i=1→k σ max j=k+1→p ≥ thv, ( 19
) then J k,n,thv = σ min i=1→k σ max j=k+1→p . k i=1 α i σ i σ max j=1→k n ( 20
)
where k is the number of the first dominant modes specified by the designer, σ min i=1→k and σ max j=k+1→p are respectively the level of the least dominant mode within the first k and the level of the most dominant mode within the residual ones. thv is a threshold value specified by the designer. Hence, the condition (19) represents the domination of the k first modes by at least thv times compared to the residual modes (see Fig. 6). Furthermore, in the numerical expression (20) of the criterion itself, α i = γ q ii 4ξ i ω 3 i corresponds to the coefficient of the i th observability modal gramian

W ob ii . σ i σ M ax j n is a weighing ratio ∈ [0, 1].
When the latter's value is close to 1, the i th mode in question is a dominant mode within the first k ones. Note that the exponant n helps emphasizing the most dominant modes by accelerating the convergence towards σ M ax i . Thus, the corresponding α i is privileged compared to other modes' observability. In other words, by maximizing this criterion, we privilege the modes where good observability of δ coincides with its dominant modes (Fig. 6). In this way, it appears that the procedure given by ( 19) and ( 20) is able to solve the problems mentioned above.

Optimal Synthesis of Monolithic

Piezoelectric Mechanism with Integrated Actuator and Sensor

Optimization problem specifications

We consider here the optimal synthesis of a monolithic micromechanism with integrated actuator and sensor, made of the piezoelectric material PIC151 from PI Piezo Ceramic Technology [START_REF]PI Piezo Ceramic Technology[END_REF].

The structure is considered to have a maximal size of 15mm × 18mm, and a constant thickness of 200µm. It is defined to be made of either, passive, active or sensitive blocks inside a 2 × 3 mesh (see Fig. 7). For the optimal synthesis run, the number of active (resp. sensitive) blocks in the structure, chosen among blocks given in Fig. 1, will be allowed to vary between 1 and 4. The size ratio of the blocks can vary as b max /b min ∈ [ [1; 2]] and a max /a min ∈ [ [1; 2]] (see Fig. 7 for details about a and b parameters definitions). When external voltages are applied to the actuating blocks' electrodes, the chosen output node of the structure has to move along the x-axis (see Fig. 7).

For evaluation of static mechanical criteria, the potential difference between upper and lower face is taken equal to 200V . The number of blocked nodes is comprised between 1 and 3 among the locations permitted, which are reported in Fig. 7.

Finally, three numerical criteria to be maximized simultaneously with FlexIn are:

• J 1 : free mechanical displacement δ x at the output node in x-direction,

• J 2 : amplitude of the sensing electric charges induced on the piezoelectric blocks, numerically computed by means of the third equation in ( 6),

• J k,n,thv : modal observability of the mechanism output δ by the observed sensitive blocks charges q g . In this example, we chose k = 2, n = 2 and thv = 3.

Results and performances comparison

After setting the optimization to run, if during 200 subsequent generations, the genetic algorithm does not find new pseudo-optimal solutions, the optimization procedure is automatically stopped, resulting into Pareto fronts such as in Fig. 8, showing the best compromises kept. From these fronts, the designer can chose the structure that is most appropriate to his study. In our case, performances of three selected pseudo-optimal solutions are compared in Fig. 9. Structure C exhibits the best output displacement (highest J 1 value), whereas structure A has the best dynamic characteristics (highest J k,n,thv value). As for J 2 criterion, it is found to be essential in increasing the range of electric charges induced on the piezoelectric sensor parts, for quasi-static measurements. In fact, in the latter mode, electric charges less than 10 -9 C require special electronic circuits, because they desert the sensor surface as soon as they are generated. Hence, even though the structures selected do not have the best J 2 values among all the pseudo optimums, structures A and C present measurable charges quantities, via circuits based on charge-voltage converters as mentioned in section 2.2.1.

To illustrate the dynamic performances, bode diagrams of control and observation transfers of the three structures are presented in Fig. 10. The control authority of the first k = 2 selected modes in the three structures is properly insured by the criterion condition (19), as shown in Fig. 10.a. As for J k,n,thv value itself presented by (20), it symbolizes significantly the same modes' observability. In fact, as expected, structure C has the poorest observability, demonstrated by a small peak in the amplitude of the observation transfer q g /u of the second mode (see Fig. 10.b). For a greater value of J k,n,thv , structure B shows better observability of the first two modes. One problem remains: the existence of the first anti-resonance complicates the computation of modal representation and thus the control of such structure via a state observer. Finally, structure A having the highest J k,n,thv value, it exhibits ideal observability of the first two modes, as reported by the high amplitude peaks of q g /u transfer (see Fig. 10.b).

Conclusion

A contribution to an existing optimal compliant micromechanisms design method has been presented in this paper. In this method, structures monolithically machined, composed of an assembly of building blocks made of PZT, optimally integrate actuating and sensing blocks by means of two new criteria: first, a static one, maximizing the amount of electric charges induced on the sensor surfaces via the direct piezoelectric effect; the second one is a control-observation oriented gramian-based criterion, considering the open-loop dynamic control and observation transfers' performances of the structure. In fact, for optimally integrating both actuation and sensing functions in the structure, the second criterion implies matching high control authority modes from the balanced gramian matrix, with high observation ability from the observation gramian matrix.

An optimization problem is specified to illustrate the mentioned methodology. The analysis of three resulting pseudo-optimal solutions brings to satisfying conclusion: the structures have good static characteristics due to J 1 and J 2 criteria, but more importantly, a high value of J k,n,thv insures good observability of the selected dominant modes. This successfully sought property will simplify the control of such structures via state observer. A prototype of the compliant pseudo-optimal micromechanism, structure A, will be developed for experimental validations and control.

Figure 1 :

 1 Figure 1: Passive (blue), piezoactuating (red) and piezosensing (green) libraries of compliant building blocks, for planar compliant mechanisms synthesis using FlexIn.

Figure 3 :

 3 Figure 3: Curvilinear coordinates of the piezoelectric beam A -B, and its orientation in the global coordinate system R ′ = (0, x, y, z). R and H represent the in-plane nodal force and moment at the beam extremities.

Figure 4 :

 4 Figure 4: Basic electronic circuit for electric charges measurements based on charge-voltage converter.

Figure 5 :

 5 Figure 5: Control-observation diagram for compliant mechanism with integrated piezoactuator and piezosensor where F co and K ob are the control and observator gains respectively.

Figure 6 :

 6 Figure 6: Example of desired form of Control and Observation transfers.

Figure 7 :

 7 Figure 7: Mesh of the PZT compliant micromechanism with imposed and permitted boundary conditions. a and b optimization parameters define the relative height and width of the blocks.

Figure 8 :

 8 Figure 8: Pareto fronts of compliant mechanisms synthesized using FlexIn.

Figure 9 :

 9 Figure 9: FlexIn representation, 3D commercial simulation software representation and performances of the A, B and C compliant structures.

Figure 10 :

 10 Figure 10: Frequency responses of A, B and C structures of a) control transfer δ/u and b) observation transfer q g /u.