M Collet

M Ouisse

M N Ichchou

Structural energy flow optimization through adaptive shunted piezoelectric metacomposites

Keywords: Adaptive Metamaterial, Distributed Shunted piezoelectric patches, Wave propagation optimization

In this paper, a numerical approach for modeling and optimizing 2D smart metacomposites is presented. The proposed methodology is based on the Floquet-Bloch theorem in the context of elastodynamics including distributed shunted piezoelectric patches. The dedicated numerical technique is able to cope with the multi-modal wave dispersions behaviour over the whole first Brillouin zone for periodically distributed 2D shunted piezo-mechanical systems. Some indicators allowing the optimization of the shunt impedance for specific performance objectives are presented, and applied for illustration purposes on the design of an adaptive metacomposite with specific functionalities. In order to validate the strategy, the designed metacomposite is integrated in a support structure, and a full 3D model is derived to illustrate the efficiency of the approach.

Introduction

The design of innovative structures with high degree of performances is a permanent challenge for both industrial and academic communities. The issue is often posed in many engineering sectors such as civil engineering and transportations (ground and aeronautics). In the transportation sector the question becomes drastic with regards to energy savings concerns and with respect to sustainable developments (see reports of Clean Sky and DREAM EU Projects, AIAA's Emerging Technologies Committee -ETC). This is the case, for instance, in the aeronautic sector where a fantastic effort is made in order to integrate lightweight structures in the design of new aircrafts. Advanced materials and structures such as composite materials are thus more and more used. Such composite materials make the structures lighter and stiffer. They are now intensively included in the design and manufacturing of new aircrafts and space shuttles. However, such composite structures bring new vibration and acoustic problems. Indeed, compared to metallic structures it is noticeable that composites structures evolves specific noise emissions and particular safety issues. The majority of the projects dedicated to smart materials have focused on adding smart functionalities to existing structures -piezoelectric patches to plates and shells for acoustic or structural health monitoring as an example. Nonetheless, these types of solutions tend to address efficiently only one particular problem at the time, and require increased complexity in terms of wiring, controllers and modifications of the baseline structure to tackle successfully structural, acoustic control and structural health monitoring. Hence, among the problems to be considered, the broadband nature (from the low frequency LF up to the high frequency HF and including the mid frequency MF range) of vibration and noise associated to composite structures and their assembly is still an open and challenging concern. This paper is a contribution to the challenges of designing, optimizing and manufacturing a new class of truly integrated and intrinsic smart active metamaterial, structures and systems, also called metacomposite with ground breaking engineering performance in terms of mechanical and vibroacoustics behaviour.

The concept of metacomposite has a specific meaning in this paper. This concept couples two different aspects in vibration control. The first concept is connected to periodic structures theories usually connected to metamaterial developments. In this case, it is well known that the dynamic behaviour is fully connected to periodicity ratios and existing pass bands and blocked bands can be of real use in vibration control. The second concept is associated to vibration control through piezoelectric and smart materials. Specifically, shunted piezoelectric smart materials are employed for the metacomposite achievement by integrating into the metamaterial electronics and numerical components allowing implementation of adaptive and controlled behaviour. We also tend to extend the notion of programmable matter within the meaning of work presented in to vibroacoustic programming. The paper main novelty is then the design through full numerical analysis of a smart structure with broad band control abilities. Wave based methods and numerical simulation tools are adapted to the proposed concept. The realized composite controlled dynamic is in the LF low frequency and mid frequency MF band as well. Development of a wave trap metacomposite can be the ultimate goal of the reported research and results.

Looking at the literature, there are intensive involvements in both the modelling and control communities. Indeed, researchers studied numerous methods to deal with undesirable vibrations and noise. These approaches can be generally classified as three catalogues which are passive control, active control and semi-active or semi-passive control. Several passive control ways were developed. Tuned mass damper (TMD) or harmonic absorber is one of the most popular passive control approach used since many years in many different engineering applications [START_REF] Ormondroyd | The theory of dynamic vibration absorber[END_REF][START_REF] Hartog | Mechanical vibrations[END_REF][START_REF] Spencer | State of the art of structural control[END_REF]. The corresponding device is normally mounted onto a host structure to absorb the structural mechanical energy when the structure vibrates at its resonant frequency. Passive damping technologies can also be considered as passive control strategies especially when optimization design is carried out [START_REF] Milman | Optimization methods for passive damper placement and tuning[END_REF][START_REF] Baz | Optimization of energy dissipation characteristics of active constrained layer damping[END_REF]). The main disadvantages of these passive strategies are also their narrow band of efficiency and the necessity to introduce massive elements (damping materials or massive absorber) inside the considered structure. Optimal efficiency in low and broad frequency bands can be reached by using active technologies for vibration or noise control [START_REF] Preumont | Vibration control of active structures: an introduction[END_REF][START_REF] Nelson | Active Control of Sound[END_REF][START_REF] Banks | Smart material structures: Modeling Estimation and Control[END_REF] with the drawback of introducing complex and massive numerical and electronics elements and external power supply systems.

Semi-active control or semi-passive control possesses partial advantages of both passive and active control. Piezoelectric shunt technologies are promising passive techniques for vibration control of flexible structures. Contrary to active control, the only external element to be used is a Passive Electrical Network (PEN) that is directly connected to the electrodes of the piezoelectric device. With this approach, the sensing element is not needed and the use of a passive network guarantees the stability of the coupled system. Furthermore, shunted piezoelectric systems offer a very attractive path for electro-mechanical integration and distribution. Numerous works have been published [START_REF] Park | A new approach to identify optimal of shunting elements for maximum damping of structural vibration using piezelectrical patches[END_REF][START_REF] Becker | Finite element-based analysis of shunted piezoelectric structures for vibration damping[END_REF][START_REF] Hurlebaus | Smart structure dynamics[END_REF][START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF][START_REF] Agnes | Active/passive piezoelectric vibration suppression[END_REF][START_REF] Dell'isola | Passive damping of beam vibrations through distributed electric networks and piezoelectric transducers: prototype design and experimental validation[END_REF][START_REF] Moheimani | Piezoelectric transducers for vibration control & damping[END_REF][START_REF] Cross | Shunted piezoelectrics for passive control of turbomachine blading flow-induced vibrations[END_REF], Beck, Cunefare & Ruzzene 2008a) that present analyses of the capability and efficiency of a single shunted piezoelectrical patch for structural stabilization and wave cancellation. Hagood and VonFlotow (Hagood & von Flotow 1991) provided the first analytical formulation for passive shunt networks. They demonstrated how a piezoelectric patch shunted through a single resistive-inductive (RL) circuit acts as a vibration absorber tuned at the resonance frequency of the circuit. Since then, more complex shunting circuits have been investigated to extend the effectiveness over broader frequency bands: multi-mode techniques (Wu 1998, Wu & Bicos 1997), state switching [START_REF] Corr | Comparison of lowfrequency piezoelectric switching shunt techniques for structural damping[END_REF], synthetic impedance [START_REF] Fleming | Synthetic impedance for implementation of piezoelectric shunt-damping circuits[END_REF] or negative impedance circuits [START_REF] Behrens | New method for multiple-mode shunt damping of structural vibration using a single piezoelectric transducer[END_REF][START_REF] Park | A new approach to identify optimal of shunting elements for maximum damping of structural vibration using piezelectrical patches[END_REF][START_REF] Park | Vibration control of beams with negative capacitive shunting of interdigital electrode piezoceramics[END_REF][START_REF] Fukada | Sound isolation by piezoelectric polymer films connected with negative capacitance circuits[END_REF][START_REF] Neubauer | Optimization of damping and absorbing performance of shunted piezo elements utilizing negative capacitance[END_REF]. Park and Palumbo decomposed the wave field on a beam to find a minimization of the reflected wave component [START_REF] Park | A new approach to identify optimal of shunting elements for maximum damping of structural vibration using piezelectrical patches[END_REF]. By noticing the link between the wave based approach of Park and Palumbo and the power input to the system, Cunefare developed a parameter selection technique based on the reactive input power [START_REF] Cunefare | Negative capacitance shunts for vibration suppression: wave based tuning and reactive input power[END_REF]). Periodic structures have been shown to improve the broadband performance of control systems on structures [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF], Thorp, Ruzzene & Baz 2001, Thorp, Ruzzene & Baz 2005). On all these works, piezoelectric shunt damping also appears as promising passive or semi-passive technique for vibration control of flexible structures (Beck et al. 2008a[START_REF] Cunefare | Negative capacitance shunts for vibration suppression: wave based tuning and reactive input power[END_REF]. Contrary to active control, the only external element to be used is a passive electrical network (PEN) that is directly connected to the electrodes of the piezoelectric device and can be integrated into the composite material [START_REF] Forward | Electromechanical transducer-coupled mechanical structure with negative capacitance compensation circuit[END_REF][START_REF] Fleming | Reducing the inductance requirements of piezoelectric shunt damping systems[END_REF][START_REF] Fleming | Improved Current and Charge Amplifiers for Driving Piezoelectric Loads, and Issues in Signal Processing Design for Synthesis of Shunt Damping Circuits[END_REF]. In this paper, thus shunted piezoelectric concept will be developed further for wave control purposes mainly. Unlike most of the existing papers, this contribution will concentrate on a deep multi-physical numerical modelling of composite structures coupled to shunted piezoelectric patches based on mathematical approach developed in (Collet, Ouisse, Ruzzene & Ichchou 2011). The numerical investigation is concerned with 2D wave treatment and optimization in such metacomposite. This is also an original aspect of the paper.

Considering numerical simulation of wave behaviours in elastodynamic media, two main strategies can be reported: the semi-analytical finite element (SAFE) and the wave finite element method (WFE). Both are of spectral nature but of different philosophies and practical use. Semi-analytical finite element method was investigated by many researchers in the past decades. This method was first proposed by [START_REF] Nelson | High frequency vibrations and waves in laminated orthotropic plates[END_REF] for plate waves to analyse natural vibrations of lam-inated orthotropic plates in which lamina is discretized by finite element. The (cross-sectionally) two-dimensional problem was first developed by [START_REF] Aalami | Waves in prismatic guides of arbitrary cross section[END_REF], which can be extended to the wave guides of arbitrary cross-sectional shape. Recently, Gavric calculated the dispersion relationship in a free rail by using triangular and quadrilateral elements, with those elements obtained from HamiltonŠs principle [START_REF] Gavric | Computation of propagatives waves in free rail using a finite element technique[END_REF]. [START_REF] Hayashi | Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example[END_REF] derived the SAFE formulation through virtual work principles and proposed a way to calculate the group velocity using the eigensolution at a given frequency. [START_REF] Damljanovic | Propagating and evanescent elastic waves in cylindrical waveguides of arbitrary cross section[END_REF] developed the linear triangular elements for SAFE method using Lagrange equations to investigate the elastic waves in wave guides of arbitrary cross-section. [START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF] examined the wave mode characteristics in damped waveguides by using SAFE method, where the elements were also derived from Hamilton's principle. This method has also been adopted to investigate the wave propagation characteristics for thin-walled structures [START_REF] Finnveden | Evaluation of modal density and group velocity by a finite element method[END_REF], where the polynomial interpolation is used in the propagation axis and thus leads to the polynomial eigenvalue problems. The method is extended to curved structures in [START_REF] Finnveden | Waveguide finite elements for curved structures[END_REF], where an isoparametric element is used.

Another spectral method has attracted many researchers, which is based on the technique dealing with the wave propagation problems of periodic structures. Mead [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF]) proposed a general theory in order to determine harmonic wave propagation characters, where both one-dimensional and two-dimensional periodic systems are considered. Recently, more studies focus on extending this idea to homogeneous structures rather than the generally periodic systems comprised of an arbitrary substructure (see [START_REF] Houillon | Wave motion in thin-walled structures[END_REF][START_REF] Bocquillet | Energetics of axisymmetric fluid-filled pipes up to high frequencies[END_REF]). This Wave Finite Element (WFE) method regards the homogeneous wave guide structure as a periodic system assembled by identical substructures. Thus the method for periodic systems can be seamlessly transplanted to the WFE method for eigenmode extraction. The identical substructures can be discretized with the aid of some commercial FE procedures, rather than the development of a relatively new FE code for specific elements. This allows the existing element libraries and powerful grid generation procedures to be used for many engineering structures. The WFE method is easy to apply due to its perfect connection with the standard FE method. This method has been used for structural vibration analyses [START_REF] Mace | Finite element prediction of wave motion in structural waveguides[END_REF], Zhou, Ichchou & Mencik 2009[START_REF] Hardy | A hybrid local energy formulation for plates mid-frequency flexural vibrations[END_REF][START_REF] Duhamel | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF][START_REF] Houillon | Wave motion in thin-walled structures[END_REF], wave propagation in elastic wave guides [START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF], Chen & Wilcox 2007), and calculations of the Green's functions for homogeneous acoustical and elastic media [START_REF] Duhamel | Finite element computation of green's functions[END_REF]. The method was also implemented for wave propagation and dynamic problems in homogeneous structures with internal fluid [START_REF] Bocquillet | Energetics of axisymmetric fluid-filled pipes up to high frequencies[END_REF], Mencik & Ichchou 2007[START_REF] Maess | Dispersion curves of fluid filled elastic pipes by standard fe models and eigenpath analysis[END_REF], where the studies were based on the WFE method for the one-dimensional wave propagation problem and concentrated on the lower frequency problems (far below the frequency of interest for NDT). One of the main problems all these approaches are facing is the difficulty to compute the damped wave numbers in complex media. Precisely, wave properties in the whole Brillouin domain for multiphysics problem incorporating shunted piezoelectric patches is not an easy task. This involvement is necessary for vibroacoustic behaviour optimization purposes of such periodic smart structures (Collet, Ouisse, Ruzzene & Ichchou 2011, Collet, Ouisse, Ichchou & Ruzzene 2011).

The objectives of the work reported in this document were threefold. First, the paper offers a suited numerical approach for computing wave dispersion in smart multiphysics metacomposite with shunted piezoelectric patches. Second, the achieved model is extensively used for electromechanical optimization studies of the metacomposite, and finally, the optimized configuration is used in a finite structure to check the applicability of the concept on a practical application. The paper is structured in four different sections. After recalling the Floquet-Bloch theorems in section 2, we introduce a numerical formulation to compute the multi-modal damped wave numbers dispersion in the whole first Brillouin domain. Based on this wave modelling, optimization of the electrical impedance of the shunted circuit is made in order to control group velocity of flexural waves or to increase damping induced by the electric circuit. The corresponding optimal impedances are then tested for controlling the mid frequency (MF) response of a semi-distributed finite system. The numerical achievements prove the effectiveness of the proposed metacomposite concept and draw future steps and perspectives.

2 Piezo-Elasto Dynamical Application of the Floquet-Bloch Theorem In this section the application of the celebrated Floquet-Bloch theorem is presented for piezoelastodynamic problems. Based on the well known results obtained by Floquet [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF] for monodimensional problems and later rediscovered by Bloch [START_REF] Bloch | Über die Quantenmechanik der Electron in Kristallgittern[END_REF] in multidimensional problems, an original application to bi-dimensional piezo-elastodynamical problem has been proposed recently [START_REF] Collet | Numerical tools for semi-active optimization of 2d wave's dispersion into mechanical system[END_REF]. This formulation leads to very general numerical implementation for computing waves dispersion for periodically smart distributed mechanical systems incorporating electronic components, damping effects or any frequency-dependent characteristics. The main ideas of the approach are recalled here, together with the specific points related to the inclusion of piezoelectric effects in the model.

The Bloch theorem

The Bloch theorem, in its original version, gives the form of homogeneous states of Schrödinger equation with periodic potential. This theorem can be considered as a multidimensional application of the Floquet theorem [START_REF] Joannopoulos | Photonic Crystals: Molding the Flow of Light[END_REF]. The periodic medium (or potential) properties satisfy M (x + R.m) = M (x), m ∈ Z 3 where R = [r 1 , r 2 , r 3 ] ∈ R 3×3 is a matrix grouping the three lattice's basis vectors (in 3D). The primitive cell is defined as a convex polyhedron of R 3 called Ω x . The reciprocal unit cell is denoted by Ω k and limited by the reciprocal lattice vector defined by the three vectors g j so that: r i .g j = 2πδ ij (δ ij being the Kronecker index). If Ω x is the irreductible primitive cell, Ω k corresponds to the first Brillouin zone of the lattice [START_REF] Kittel | Introduction to Solid State Physics[END_REF].

The Bloch Theorem stipulates that any function u

(x) ∈ L 2 (R 3 , C n ) can be expressed as u(x) = Ω k e ikx ũ(x, k)dk (1)
where the Bloch amplitude ũ(x, k) is Ω x -periodic and has the representations

ũ(x, k) = n∈Z 3 û(k + Gn)e iGn.x , = |Ω x | (2π) 3 n∈Z 3 u(x + Rn)e ik(x+Rn) (2)
where û(k) stands for the Fourier transform of u(x) and G = [g 1 , g 2 , g 3 ] is the reciprocal lattice matrix in the later. It can also be demonstrated that the mean value of the Bloch amplitude is the Fourier amplitude of u(x) for the corresponding wave vector: ũ(., k) Ωx = û(k). Using the Bloch theorem to represent the solutions of periodical partial derivative equations implies that all derivatives are shifted by k in the sense given by the spatial operator.

Based on that theorem, the expansion functions v m (x, k) can be defined. They are called the Bloch eigen modes, and can be used to represent the Bloch amplitudes of any solution of the corresponding partial derivative equation as

ũ(x, k) = m u m (k)v m (x, k) (3)
and at the same time diagonalize the partial derivative equations. The expansion coefficients u m (k) depend on the disturbance and on the induced wave vector (see [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF] for details). Let us consider a piezo-elastodynamic problem made of infinite periodic distribution of unitary cell described in figure 1. The harmonic homogeneous dynamical equilibrium of system is driven by the following partial derivative equations (t ∈ R):

Application to piezo-elastodynamic

ρ ẅ(x, t) -∇σ(x, t) = 0 ∀x ∈ Ω x -∇D(x, t) = 0 ∀x ∈ Ω x (4)
where w(x, t) ∈ R 3 (Ω x , R) is the displacement vector, σ represents the Cauchy stress tensor, = ∇ sym w = 1 2 (∇w T (x, t)+w(x, t)∇ T ) the Green strain tensor, D(x, t) the electric displacement. The linear constitutive material behavior relationships can be written as

σ = C E (x) -e T (x)E, (5)
D = e(x) + ε S (x)E, (6)
where E = -∇V is the electric field vector (V being the voltage), C E the elasticity tensor at constant electrical field, e T the piezoelectric coupling tensor and ε S the dielectric permittivity at constant strain. We add to this set of equilibrium equations an output expression

q o = - St D.ndS (7)
allowing the introduction of the charge measurement on the piezoelectric's top electrode and hence the dual counterpart of the imposed electrical Dirichlet boundary condition for applying the shunt impedance operator. The equations above are consistent for each kind of material to the extent that null piezoelectric and permittivity tensors can be used when passive materials are considered. All of these tensors also depend on the spatial location vector x and are Ω x periodic. By applying a Fourier transformation, the piezo-elastodynamic equilibrium can also be written as (ω ∈ R):

ρω 2 w(x, ω) + ∇C∇ sym (w(x, ω)) + ∇e T (x)∇V (x, ω) = 0 ∀x ∈ Ω x (8) -∇e(x)∇ sym (w(x, ω)) + ∇ε S (x, ω)∇V (x, ω) = 0 ∀x ∈ Ω x (9)
As the problem is 2D infinitely periodic, mechanical boundary conditions are included in the formulation, while electrostatic boundary conditions have to be considered on each cell (ω ∈ R):

   V (x, ω) = 0 ∀x ∈ S b V (x, ω) = V o (ω) ∀x ∈ S t D.n(x, ω) = 0 ∀x ∈ S l (10)
where S b is the grounded bottom electrode of the piezoelectric layer, S t is the top electrode connected to the external shunt and s l the lateral electrode. The top electrode applied feedback voltage V o depends on the shunt characteristic and on the collected charges q o given in eq. ( 7) and can be expressed in the Fourier space by:

V o (ω) = -Z(ω)q o (ω) (11)
In the following, the frequency dependency on ω of each functions in the Fourier space is removed for simplifying the notations. Considering a primitive cell of the periodic problem Ω x , the Bloch eigenmodes (3) and the dispersion functions can be computed by searching the eigen solutions of the homogeneous problem ( 8) and ( 9) with mechanical periodic boundary conditions and electric ones given by (10) as:

u(x) = w(x) V (x) = u n,k (x)e ik.x (12)
where

u n,k (x) = w n,k (x) V n,k (x)
is composed of Ω x periodic functions. By introducing expression (12) in the piezo-elastodynamic equations ( 8) and ( 9), w n,k (x), V n,k (x) and ω n (k) can be found by solving the generalized eigenvalues problem:

0 = ρω 2 n (k)w n,k (x) + ∇C∇ sym w n,k (x) + ∇e T ∇V n,k (x) +ik [(C∇ sym w n,k (x)).Φ + ∇(CΞ n,k (x))] +ik (∇e T V n,k (x)).Φ(e T ∇V n,k (x)).Φ -k 2 (CΞ n,k (x)).Φ + V n,k (x)(e T Φ).Φ ∀x ∈ Ω x , (13) 0 = -∇e∇ sym w n,k (x) + ∇ε S ∇V n,k (x) -ik [∇(eΞ n,k (x)) + (e∇ sym (w n,k (x))).Φ] +ik [(∇ε S V n,k (x)).Φ + (ε S ∇V n,k )(x).Φ] +k 2 [(eΞ n,k (x)).Φ -( S ΦV n,k (x)).Φ] ∀x ∈ Ω x , (14)
with the associated boundary conditions:

       w n,k (x -R.m) = w n,k (x) ∀x ∈ S r m ∈ Z 2 , V n,k (x) = 0 ∀x ∈ S b , V n,k (x) = -Z(iω)q o n,k ∀x ∈ S t , D.n = 0 ∀x ∈ S l .
In

(x) = 1 2 (w n,k (x).Φ T + Φ.w T n,k (x))
is the symmetric dyadic tensor or the dyadic product of the displacement w n,k (x) and direction vector Φ. S r are the interfaces of the cells continuum, and R is the matrix grouping the two lattice's basis vectors (in 2D in the considered problem). In the electrical boundary conditions, q o n,k is given by:

q o n,k = St [-e∇ sym w n,k (x) + ik∇eΞ n,k (x) +ε S (∇V n,k (x) + ikV n,k (x)Φ)].ndS ( 15
)
where n is the surface outpointing unitary normal vector. The proposed formulation is then based on the computation of the Floquet vectors given in equations ( 13) and ( 14)), instead of computing the Floquet propagators commonly used for elastodynamic applications. The full 2D waves dispersions functions can then be obtained, while damping and electrical impedance can clearly be introduced into the piezo-elastodynamic operator. The adopted methodology allows the computation of the complete complex map of the dispersion curves incorporating computation of evanescent waves and allowing the introduction of damping and shunt operator if any.

Weak formulation and computation of waves dispersion functions in periodical piezo-composite lattice

Let us consider the partial derivative equations ( 13), ( 14) on a unit cell Ω. It stands for a generalized eigenvalue problem leading to compute the dispersion functions ω n (k) and the corresponding Floquet eigenvectors u n,k (x). For computing the 2D dispersions curves, we need to introduce a suitable weak formulation.

If u n,k (x) is a solution of equations ( 13), ( 14), then:

∀ wn,k (x) ∈ H 1 (Ω, C 3 )/ wn,k (x -Rm) = wn,k (x) ∀x ∈ S r and ∀ Ṽn,k (x) ∈ {H 1 (Ω, C)/ Ṽn,k (x) = 0 ∀x ∈ S b and Ṽn,k (x) = V o ∀x ∈ S t } one has: 0 = Ω ρω 2 n (k) wn,k (x)w n,k (x) -(ε n,k (x) -ik Ξn,k (x))C(ε n,k (x) + ikΞ n,k (x)) +(ε n,k (x) -ik Ξn,k (x))e T (∇V n,k (x) + ikV n,k (x)Φ) -(∇ Ṽn,k (x) -ik Ṽn,k (x)Φ)e(ε n,k (x) + ikΞ n,k (x)) +(∇ Ṽn,k (x) -ik Ṽn,k (x)Φ)ε S (∇V n,k (x) + ikV n,k (x)Φ ] dΩ - Ṽ t n,k V t n,k Z(iω) (16)
This weak formulation is simply obtained by integrating equation ( 13), ( 14) projected onto any test function ũn,k (x). The boundary integral vanishes as the test functions are chosen so that wn,k (x -Rm) = wn,k (x) on S r . For a polyhedron cell, each boundary is generally a polyhedral plane sub-domain that can be associated with a parallel opposite one. The symmetry conditions called wn,k (x -Rm) = wn,k (x) explicitly link these associated surfaces.

Numerical computation of the Bloch's waves

The numerical implementation is obtained by using a standard finite elements method to discretize the weak formulation ( 16). The assembled matrix equation is given by:

0 = [K(Z(ω n )) + λL(φ, Z(ω n )) -λ 2 H(φ, Z(ω n )) -ω 2 n M u n (k, φ), (17)
where λ = ik, M and K(Z(ω n )) are respectively the standard symmetric semi-definite mass and stiffness matrices (the mass matrix is semi definite because elastostatic equation are condensed into the equation), L(φ, Z(ω n )) is a skew-symmetric matrix and H(φ, Z(ω n )) is a symmetric semi-definite positive matrix. u n is the generalized eigen vector defined on all degrees of freedom of the used finite element model.

When k and φ are fixed and Z does not depend on ω the system ( 17) is a linear eigen value problem allowing us to compute the dispersion functions ω 2 n (k, φ) and the associated Bloch eigenvector u n (k, φ).

This approach has been widely used for developing homogenization techniques and spectral asymptotic analysis like in the work of [START_REF] Allaire | Bloch waves homogenization and spectral asymptotic analysis[END_REF]. It can also be applied for computing wave's dispersion even if Floquet propagators is preferred for 1D or quasi 1D computation, as indicated in [START_REF] Ichchou | Guided waves group and energy velocities via finite elements[END_REF], Mencik & Ichchou 2005a). Nevertheless these approaches have been only developed for undamped or lightly damped mechanical systems. In these cases, most of the previously published works present techniques based on the mesh of a real k-space (i.e k or k and φ) following the boundary of the first Brillouin zone for obtaining the corresponding dispersion curves and the associated Floquet vectors. For undamped system, only propagative or evanescent waves exist corresponding to a family of eigen solutions purely real or imaginary. Discrimination between each class of waves is easy. If a highly damped system (K, L, H are complex frequency dependent) and a frequency dependent electrical shunt impedance are considered, the obtained eigenvalue problem is not quadratic and a complex specific numerical methodology has to be implemented. Furthermore, evanescent parts of propagating waves appear as the imaginary part of ω 2 n (k, Φ). It then becomes much more difficult to distinguish the propagative and evanescent waves as all solution appear complex.

Another much more suitable possibility for computing damped system, dedicated for time/space deconvolution and for computation of diffusion properties as defined by (Collet, Ouisse, Ruzzene & Ichchou 2011[START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF], Mencik & Ichchou 2005a), is to consider the following generalized eigen value problem:

0 = (K(Z(ω) -ω 2 M + λ n (ω, φ)L(φ, Z(ω)) -λ 2 n (ω, φ)H(φ, Z(ω)))u n (ω, φ). (18)
In this problem, the pulsation ω is a real parameter corresponding to the harmonic frequency. Wave's numbers and Floquet vectors are then computed. An inverse Fourier transformation in the k-space domain can lead us to evaluate the physical wave's displacements and energy diffusion operator when the periodic distribution is connected to another system as in [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF]. Another temporal inverse Fourier transformation can furnish a way to access spatio-temporal response for non-homogeneous initial conditions. As L is skew-symmetric, the obtained eigen values are quadruple (λ, λ, -λ, -λ) collapsing into real or imaginary pairs (or a single zero) when all matrices are real (i.e. for an undamped system). In this case a real pair of eigen values correspond to evanescent modes oriented in two opposite directions on the k-space and imaginary values to two traveling waves propagating in opposite direction as described in [START_REF] Garvic | Computation of propagative waves in free rail using a finite element technique[END_REF].

As previously mentioned, the real part of k = kΦ vector is restricted to stand inside the first Brillouin zone. In the quadratic eigen value problem (18) nothing restricts computation to only find eigen values satisfying this condition. For direction vector Φ orthogonal to the lattice facelets (i.e. for Φ p1 = [1, 0] T and Φ p2 = [0, 1] T in bi-dimensional rectangular cell), the periodical conditions expressed for one dimensional wave guide are still valid: if λ j (ω, Φ p ) is an eigen value associated to w j (ω, Φ p ) then ∀m ∈ Z 3 , λ+i.Φ T p (G.m) is also an eigen value associated to w j (ω, Φ p )e -i.Φ T p (G.m)x . Thus, for undamped systems, all obtained eigenvalues are periodically distributed in the k-space along its principal directions.

Computation of the evanescence and damped power flow criteria

One aim of this paper is to provide a numerical methodology for optimizing the piezoelectric shunt impedance Z(ω) for controlling energy flow into the periodically distributed piezo-composite structure. For doing this, we need to define suitable criteria evaluating the efficiency of the approach on the features of interest.

The first criterion which is considered for describing the capability of the metacomposite for transmitting structural energy is based on the computation of the wave group velocities. Indeed, they indicate how energy is transported into the considered system and allow to distinguish the 'propagative' and 'evanescent' waves. If a Bloch eigen solution (i.e u n (ω, φ), k n (ω)) is considered, the associated group velocity vector [START_REF] Maysenhölder | Körperschall-energie Grundlagen zur Berechnung von Energiedichten und Intensitäten[END_REF]) is given by:

C g n (ω, φ) = ∇ k ω = S e tot = I E tot ( 19
)
where : is the spatial and time average respectively on one cell and one period of time, S is the density of energy flow, I the mean intensity and e tot , E tot the total piezomechanical energy and its time average on a period (see [START_REF] Maysenhölder | Körperschall-energie Grundlagen zur Berechnung von Energiedichten und Intensitäten[END_REF]) for details). In this problem, we only consider mechanical energy transportation as the electrostatic coupling is decentralized on each cell and can not induce spatial energy transportation. The electromechanical coupling effect is also condensed into the mechanical behaviour as explicitly introduced in [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] and extended to finite element model in (Collet & Cunefare 2008).

The intensity vector I is expressed as:

I n = - ω 2 Re Ωx C(ε n (x) + ikΞ n (x)).(w * n (x)) dΩ V ol (20)
where . * is the complex conjugate, Re stands for real part and V ol the domain volume.

As the spatio-temporal average of the system Lagragian is null (see [START_REF] Maysenhölder | Körperschall-energie Grundlagen zur Berechnung von Energiedichten und Intensitäten[END_REF])), the total energy average is approximated by only computing the kinetic energy average:

E tot = 1 2V ol ( Ωx ρω 2 w n (x, ω, φ).w * n (x, ω, φ)dΩ). (21)
The group velocity vectors C g n (ω, φ) are computed for all wave numbers at each frequency. In order to focus our analysis on only flexural modes (S and SH ones) we introduce an indicator allowing to select them by computing the ratio of kinetic energy average on out of plane displacement as:

Ind(n, ω, φ) = 1 2V ol Ωx ρω 2 wz n (x)wz * n (x)dΩ E tot , (22)
wz n (x) being the (Oz) component of vector w n (x). Optimization of the shunt impedance Z(ω) is based on the minimization of the maximal group velocity collinear to the wave number vector ( 19) for waves having a ratio of transported flexural kinetic energy ( 22) greater than 0.8. The criterion can then be written as:

Crit 1 (Z(ω), φ) = max n/Ind(n,ω,φ)>0.8 (C g n (ω, φ).Φ). (23)
The second criterion for optimizing the structural energy absorption is based on the maximization of the damped electric power expressed as the active electrical power

P elec (n, ω, φ) = 1 2 Re(iωZ(ω)q o n,k q o * n,k ). (24)
In order to increase the damping effects inside the smart metacomposite, this term needs to be sufficiently large. In the second case, the criteria is then

Crit 2 (Z(ω), φ) = max n/Ind(n,ω,φ)>0.8 1 P elec (n, ω, φ) . ( 25
)
3 Optimization of the flexural energy flow inside the shunted periodic piezo-composite

The considered piezo-composite cell is presented in figure 1. The supporting plate material is standard aluminum with 0.1 % of hysteretic damping ratio and the piezoelectric material characteristics are given in table 1.

Optimization of the waves group velocities by using Crit 1

In this part, optimization of the transmission capability of the designed adaptive metacomposite is considered by using Crit 1 given in equation ( 23). The objective is to avoid any energy transportation when flexural waves (A 0 , A 1 , T H...) are excited into the periodically distributed shunted piezocomposite cells. The numerical optimization procedure is based on a multidimensional unconstrained nonlinear minimization algorithm (Nelder-Mead). The optimization is done by considering an active/reactive electronic circuit through an impedance Z(ω) ∈ C. The analysis is initialized with an arbitrary complex value of the shunt impedance. Optimization steps are then proceeded using criterion (23) by considering any frequency dependent complex impedance for describing the circuit behaviour. Figure 2 shows the dispersion curves using the propagative parts (here the real one) of the wave numbers kx n (ω) along (Ox) axis, for both initial and optimal impedance shunts.

Piezoelectric Material

Symbol

Value Property This figure underlines the ability of the procedure to catch all Bloch's solutions including evanescent ones that generally do not appear in literature because of their imaginary character. In the context of structural dynamics, including damping effects is mandatory for real-life applications and all wave numbers become complex, which constitutes a key point in the analysis as underlined in (Collet, Ouisse, Ruzzene & Ichchou 2011[START_REF] Garvic | Computation of propagative waves in free rail using a finite element technique[END_REF].

s E 11 = s E 22 = s E 33 11.6e -12 P a -1 s E 12 = s E 13 = s E 23 -3.33e
By using a suitable wave filter based on the group velocity, it is possible to retain for presentation only propagative solutions responsible for energy transport. For a sake of clarity, the following results will be presented by exhibiting only propagative waves. Above results will be filtered and presented as in figure 3, which shows the propagative parts of the wave number kx n (ω) along (Ox). It then becomes easy to observe branches similar to standard S 0 , A 0 , A 1 , T H waves.

The corresponding group velocities along (Ox) are presented in figure 4. It can immediately be observed that the optimization of the shunt impedance leads to a large decrease of the group velocity of the A 0 mode while the A 1 wave, which becomes propagative at 8.8 kHz, is not controlled by the optimal configuration (see figure 4). The bending waves also propagate energy with a very slow velocity and can be considered as evanescent. Flexural energy is, also, only transported by the A 1 mode after the cutting frequency.

The propagative part of the controlled wave is strongly modified, as illustrated in figure 3: from two complex conjugated solutions, the A 0 wave changes after control to 4 solutions (2 complex conjugates and their opposites). This situation is described in [START_REF] Garvic | Computation of propagative waves in free rail using a finite element technique[END_REF], it corresponds to high order evanescent waves solutions.

It can also be observed on the group velocity diagram 4 that a non dispersive symmetric wave is also affected by the control whereas it was not targeted by the criterion.

The optimization has been performed for several values of φ. In order to show the stability of the optimized configuration with this angle, the figure 5 shows the propagative parts of the wave number kx n (ω) along 6 directions forming angles with the (Ox) axis of nπ 20 , n being all integers from 0 to 5. It is clear on this figure that above comments are valid for all angles. In particular, the cluded). o: initial value of shunt; +: optimal value of shunt. Large marks correspond to flexural waves, small marks correspond to other waves high order of evanescence of controlled waves in the optimized configurations can be observed in all situations, while some non dispersive symmetric waves are systematically affected by the control whereas they were not targeted by the criterion.

The real and imaginary parts of the optimal impedance are plotted in figure 6 for all angles φ. The optimal impedance values almost correspond to a constant negative capacitance in all directions. The corresponding average value is -233.66 pC.V -1 . Equivalent resistances corresponding to the active part of the shunt impedance (figure 6) are negative which indicates that the optimization leads to provide energy to the system for controlling mechanical damping effects introduced by hysteretic damping in the model. The optimized configuration tends to converge towards a fully conservative system. The obtained mean value of the resistance is -1.5319 Ω. This is confirmed in figure 7 where the electrical dissipated energy appears negative when the optimal shunt is connected to the patch.

In order to check this point, an optimization run has been performed on a highly damped system, using an hysteretic coefficient of 5 %. Figure 8 show the corresponding dispersion curves for φ = 0. The optimized configurations leads to a behaviour which is very similar to the one illustrated in figure 3 for lower structural damping value. The same observation can be done for other angle values. This similar behaviour is only related to propagative parts of the wave numbers, while the dissipative effects are modified as illustrated in figure 9 which shows the new optimal impedances. It can be underlined that the electric circuit is fully active with larger resistive parts of average value -74.88 Ω. The observed proportionality between the optimal average resistances (-1.5319 and -74.88 Ω) and the chosen mechanical damping terms (0.1 % and 5 %) indicates that the optimization leads to annihilate mechanical damping in order to reach a total cancellation of the group velocity of A 0 mode. The final controlled system is then almost fully reactive.

This last point is of first importance in both physical and technological points of view. Indeed, this results tends to show, on this example, that a pure band gap effect is only reachable if the Large marks correspond to flexural waves, small marks correspond to other waves system is completely reactive, which means that no damping occurs in the system. The technological implementation of such interface allowing cancellation of structural transmissibility is then unrealistic since it leads to the implementation of only 'Lyapunov stable' system: the stability margins have to be drastically minimized by using negative resistance for obtaining transmission control.

Optimization of damped power flow inside the electric shunts by using Crit 2

Another strategy for optimizing the adaptive metacomposite consists in focusing on the damped power flow inside the electric shunts by using the criteria Crit 2 given in equation ( 25). In this case the objective is to improve the absorption capability of the smart structure.

The propagative parts of the wave number kx n (iω) along (Ox) axis appear as in figure 10. The corresponding group velocities along (Ox) are presented in figure 11, while the real and imaginary parts of the optimal impedance are plotted in figure 13 and the electrically damped power in figure 12.

The optimal impedances of the electric shunts are plotted in figure 13 for various values of φ.

The first observation is that the optimization of the shunt impedance for improving the absorption characteristics of the system induces modifications of the group velocities of the controlled waves (see figure 11), while the propagative part of the wave numbers remain unchanged (figure 10). This can be explained by a large improvement of the ratio between the real and imaginary parts of the waves numbers, which physically corresponds to the forcing of propagating effects to increase damping effects: energy can propagate inside the periodically distributed set of active cells for allowing electrical energy conversion. The figure 12 shows that a the dissipated power is largely increased when optimal shunt is connected to the patch.

The optimal impedance values correspond to an almost constant negative capacitance in all directions (see figure 13). The corresponding average value is -237.43 pC.V -1 . Equivalent resis- with n being all integers from 0 to 5. o: initial value of shunt; +: optimal value of shunt. Large marks correspond to flexural waves, small marks correspond to other waves tances corresponding to the active part of the shunt impedance (figure 13) are positive, which is in accordance with the fact that a damping effect is awaited.

Validation on a periodically semi-distributed set of adaptive cells

Up to that point, the optimizations of the smart metacomposite which have been proposed were performed by considering infinite periodic structures. The objective of this section is to illustrate the applicability of the designs for integration into real-life finite structures. The optimal impedance is then applied on a finite set of shunted piezo-composite cells distributed on a part of a totally free plate system submitted to a point force in one corner, as described in figure 14. Material damping is introduced throughout a constant imaginary part of the elasticity tensor of 0.1%. The harmonic response of this system is then computed on the [500, 10000] Hz frequency band when optimal impedances, obtained by criteria Crit 1 or Crit 2 , are connected or not to each patches. The figure 15 presents the mean kinetic energy on the considered system (i.e on the whole domain From the observation of figure 15, the main particularities of the synthesized optimal impedances can be emphasized. The first configuration, obtained by using Crit 1 , lead the adaptive interface to avoid any flexural energy transfer. The mechanical response of the system is strongly modified, and the resonance peaks are moved far form the open circuit response, which indicates a modification of the internal dynamics. The second configuration, based on Crit 2 optimal shunt, forces energy dissipation into the resistive part of the electric shunt. In that case, the conservation of the modal architecture into the energy spectrum can be observed, even if some damping effects limit the peak's amplitudes.

Ω = Ω s ∪ Ω i ∪ Ω t )
To go further, one can take a closer look at the energy repartition into the considered system in figure 16, which shows the ratio of the kinetic energy finally transmitted to the sub-domain Ω t . The minimization of the flexural wave group velocity (Crit 1 ) leads to decrease the part of kinetic energy transmitted to the non-excited part of the plate located behind the adaptive interface while the maximization of the damped energy (Crit 2 ) only slightly modify this distribution. The modification of the energy distribution into the system is very important after 6kHz and increase with the frequency when the energy is globally transported by A 0 mode. This is not the case in the lower frequency band when the interface works in the near field domain of the applied point force.

The figure 17 shows the corresponding results in terms of velocity root mean square. These numerical results clearly show a strong influence of the modifications of the wave dispersion characteristics on the standing wave responses plotted in figure 17 especially when Crit 1 is used. In particular, at 2550Hz, a large attenuation of the dynamic response of the system can be observed when optimal shunts are connected. In this case, when the shunts are optimized for cancelling flexural energy flow (Crit 1 ), the location of the nodal line are also modified. With an optimization aiming at damping energy in the electric resistance (Crit 2 ), only to a standard attenuation can be observed. At 8000 Hz when Crit 1 optimal shunts are connected, a large part of the system energy remains in its left part where the applied force is located. The smart metacomposite interface also changes the system admittance and filter wave diffusion by increasing is reflexivity properties. The energy diffusion is also clearly condensed into the left part of the system with a largely decreased

Conclusions

This article presents a numerical procedure able to compute the damped wave's dispersion functions in the whole first Brillouin domain of multi dimensional piezo-elastodynamical wave guides. The method was applied for determining the optimal shunt impedance allowing the minimization of the group velocities of the flexural waves and the damping capability of the distributed system. The smart metacomposite can then be included in a finite structure to confer on it new functionalities. Some numerical tests on a finite dimension system incorporating a semi-distributed set of shunted piezo-composite cells has been performed. A strong influence of the designed shunt circuits has been underlined in the dynamical response of the system and energy distribution into each structural sub domains. Even if the link between the obtained wave properties are not clearly established especially in the low frequency range where the adaptive interface is located in the near field domain of the input force, it has been demonstrated that the numerical procedures can be used for optimizing the energy diffusion operator of such adaptive mechanical interface. To do so, additional work has to be done for optimizing the complete interface scattering and for controlling the evanescent waves playing an important role in the finite system dynamical response in the low frequency band.

The proposed methodology can also be used for studying particular dissipation phenomenon such as those induced by complex shunted piezoelectric patches as proposed by [START_REF] Beck | Broadband vibration suppression assessment of negative impedance shunts[END_REF]) and [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF]), or even foams or complex polymers behaviours. The proposed method furnishes an efficient tool for future optimization of distributed smart cells as proposed in the case of 1D wave guide by [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF]).

-

Figure 1 :

 1 Figure 1: Generic 3D piezocomposite periodic cell

  these equations k = k where φ represents the direction angles into the reciprocal 2D lattice domain and Ξ n,k

Figure 2 :

 2 Figure 2: Propagative parts of the wave number kx n (iω) along (Ox) (all waves included). o: initial

Figure 3 :

 3 Figure 3: Propagative parts of the wave number kx n (ω) along (Ox) (only propagating waves in-

Figure 4 :

 4 Figure 4: Group velocities along (Ox) direction. o: initial value of shunt; +: optimal value of shunt.

Figure 5 :

 5 Figure 5: Propagative parts of the wave number kx n (iω) along 6 directions forming an angle with (Ox) axis of nπ 20 with n being all integers from 0 to 5. o: initial value of shunt; +: optimal value of

Figure 6 :

 6 Figure 6: Optimal electric impedance, represented as equivalent resistance and capacitance, obtained

Figure 7 :

 7 Figure 7: Damped power inside the electric circuit. o: initial value of shunt; +: optimal value of

Figure 8 :

 8 Figure 8: Propagative parts of the wave number kx n (iω) along (Ox), for high hysteretic damping of

Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 9101112 Figure 9: Optimal electric impedance for high hysteretic damping of 5 %, represented as equivalent resistances and capacitances, obtained along 6 direction forming an angle with (Ox) axis of nπ 20 withn being all integers from 0 to 5.

Figure 13 :FrequencyFigure 15 :Figure 16 :

 131516 Figure 13: Optimal electric impedance, represented as equivalent resistances and capacitances, ob-

Figure 17 :

 17 Figure 17: Velocity root mean square level of the system submitted to an external point force at

Table 1 :

 1 Piezoelectric patch characteristics

		-12 P a -1	Compliance matrix
	s E 44 = s E 55 = s E 66	45.0e -11 P a -1
	η	0.1 %	Hysteretic Damping ratio
	d 31 = d 32	-6e -11 C/N
	d 33	15.2e -11 C/N	Piezoelectric matrix
	d 24 = d 15	730e -12 C/N
	ρ	7600kg/m 3	Density
	ε T 11 = ε T 22	504.1 ε o C/V /m	Dielectric Permittivity
	ε T 33	270 ε o C/V /m

  with open circuit and Crit 1 , Crit 2 optimal shunts connected to each piezoelectric patch.

			-1
			-2
	R opt [Ω]	-3 -4									0 π/20 π/10
			-5									3π/20 π/5
												π/4
		-6 -2.32	x 10 -11
		-2.33
	C opt [F]	-2.34
		-2.35
		0 -2.36	1000	2000	3000	4000	5000	6000	7000	8000	9000	10000
								Frequency [Hz]

This work is a collaborative effort supported by the French Research Agency under grant number N T 09 -617542.


E Cervera

P Martinet

Combining Pixel and Depth Information in Image-Based Visual Servoing

Image-based visual servoing has been found generally satisfactory and robust in the presence o f c amera and hand-eye calibration errors. However, in some cases, singularities and local minima may arise. We propose a modi cation of the feature v e ctor which alleviates this problem without need o f i n t r oducing further information, i.e., only pixel coordinates and depth estimations are used. Using the task-function approach, we demonstrate the relationship between the velocity screw of the camera and the current and desired p oses of the object in the camera f r ame. Experimental results on a real robotic platform illustrate the presented approaches.

Introduction

This work aims at improving the behavior of imagebased visual servoing. As pointed out in 1, 5], in some cases, convergence and stability problems may o c c u r . We propose a new feature vector based on the same information used by the classic approach, i.e. image point coordinates and depth estimation.

Using the task-function framework 8, 3 ], the velocity screw of the camera can be obtained. We w i l l demonstrate that, for small changes of orientation, the computed screw depends only on the current and desired poses of the object in the camera frame, and consequently there are no singularities nor local minima in the trajectory of the camera.

Later, we study the case of directly using the 3D coordinates of object points, which turns to be a particular case of the previous one. Nonetheless, stronger results for the velocity screw are obtained for some particular geometric con gurations of the target object. Though 3D visual features have been studied before (see e.g. 7]), we propose a new approach a n d present for the rst time analytical results for the velocity screw based on a set of feature points.

The remainder of this section brie y reviews the classic image-based visual servoing approach a n d presents some de nitions. In Section 2 the new feature vector is presented, together with the development o f the interaction matrix and the velocity screw of the camera. Section 3 discusses the case of using the 3D coordinates of the points of the object. Experimental results are described in Section 4. Finally, Section 5 draws some conclusions and open issues.

Image-based visual servoing

In image-based visual servoing, a feature vector s has to reach a desired value s . Usually, s is composed of the image coordinates of several points of the object.

The key issue in visual servoing is to nd the relationship between the derivative of the feature vector and the velocity screw of the camera v = ( t ! t ) t :

_ s = L s v (1)
where L s is the interaction matrix or Jacobian matrix.

If the feature vector s i is composed of the image coordinates (u i v i ) of a single point p i , then, as shown in 6], the interaction matrix is the presented in Equation 2 (shown at the top of next page) where Z i is the z-coordinate of p i , A is the matrix of the camera intrinsic parameters (see e.g. 4] for more details)

A = 2 4 f k u ;fk u cot u 0 0 f k v sin v 0 0 0 1 3 5 = 2 4 u uv u 0 0 v v 0 0 0 1 3 5
(3) and where (x i y i ) t are obtained from the pixel coordinates by means of the relationship

x i y i = 1 u ; uv u v 0 1 v u i ; u 0 v i ; v 0 (4) L si (s i Z i A) = u uv 0 v ; 1 Zi 0 xi Zi x i y i ;(1 + x 2 i ) y i 0 ; 1 Zi yi Zi (1 + y 2 i ) ;x i y i ;x i (2)
If several image points are used, the interaction matrix is obtained by simply stacking the matrices for each point L s = (5)

De nitions

Unless otherwise stated, the frame of reference is the camera frame. All the coordinates are expressed with respect to this frame of reference. Thus, let p i = ( X i Y i Z i ) t be the coordinates of the i point of the object, expressed in the camera frame, at the current camera pose. Let p i be its coordinates at the desired camera pose.

Without loss of generality, let us de ne a frame F b attached to the object, whose origin p is the center of gravity of the object. Then, assuming that the object is composed of n points, p = 1 n n X i=1 p i [START_REF] Hutchinson | A tutorial on visual servo c o n trol[END_REF] Let the rotation matrix R be the rotation component o f F b , i.e., the rotation between the camera and object frames. Identically, let p and R be respectively the position and orientation of the object in the desired camera frame.

The coordinates of the points can be readily expressed with respect to the object frame. Let b p i be such a coordinates vector, which is related to the coordinates in the camera frame by

p i = R b p i + p (7) p i = R b p i + p (8)
From the previous de nitions, it follows that The main advantage of image-based visual servoing is the direct utilization of image features, without the need to compute the three-dimensional pose of the target. However, in the original approach, depth information corresponding to each image point m ust be introduced in the Jacobian matrix, as well as the camera calibration parameters. In fact, if these informations are known, the three-dimensional coordinates of the points can be readily obtained, and 3D image features can be used.

The feature vector

However, the explicit calculus of the coordinates is not necessary. Assuming a camera model without distortion, the relationship between the pixel coordinates (u i v i ) t and the image coordinates (x i y i ) t is linear:

u i v i = u uv 0 v x i y i + u 0 v 0 ( 12
)
Since the image coordinates are obtained from the perspective projection of the coordinates of each p o i n t p i in the camera frame,

x i y i = 1 Z i X i Y i (13)
the following relationship is readily obtained:

2 4 u i Z i v i Z i Z i 3 5 = 2 4 u uv u 0 0 v v 0 0 0 1 3 5 p i = Ap i (14) Let s i = ( u i Z i v i Z i Z i
) t be the feature v e ctor which corresponds to point i. Then:

s i = Ap i (15) and its derivative i s _ s i = A _ p i (16)

The interaction matrix

The velocity o f p o i n t i of the object, expressed in terms of the velocity screw of the camera v

= t ! t ] t , is _ p i = ; + p i ] ! (17)
Thus, the derivative of the feature vector is

_ s i = ;A A p i ] v (18)
For a feature vector composed of a set of points, i.e., s = ( : : : s t i : : : ) t , the interaction matrix L s is:

L s = 2 6 6 4 . . . ;A A p i ]
. . . where each p i can be obtained from the corresponding s i since, from Equation 15, p i = A ;1 s i (20

)
where

A ;1 = 2 4 1 u ; uv u v ; ; v u0+ uvv0 u v 0 1 v ; v0 v 0 0 1 3 5 (21)
which always exists, since neither u nor v are null.

Then, the interaction matrix can be written in terms of the feature vector only as L s = 2 6 6 4 . . . ;A A A ;1 s i . . . Usually, the pseudo-inverse of the interacion matrix is computed numerically at each iteration of the control algorithm. However, it may be derived symbolically, a s w e show in next section. For that purpose, the expression of the matrix in Equation 19 is utilized.

We do not intend to use the symbolic result in the control loop, but it will be useful to calculate the velocity screw of the camera.

We will demonstrate that, for small rotations of the camera, such v elocity screw can be approximated by a compact expression which only depends on the current and desired p oses of the object in the camera frame.

Computation of the pseudo-inverse matrix

In general, matrix L s is not square, thus it is not invertible. However, a pseudo-inverse matrix L + s can be calculated as, L + s = ; L t s L s ;1 L t s (23) so that L + s L s = I. Let us calculate the product L t s L s of the interaction matrix by its transpose: ;nA t A p] n p] A t A ;n p] A t A p] ; S (24) where,

;A t ; p i ] A t
S = RMR t (25) M = n X i=1 b p i R t A t AR b p i (26)
The inverse matrix of the above product L t s L s is

1 n (A t A) ;1 + p] S ;1 p] ; p] S ;1 S ;1 p] ;S ;1 (27)
which exists as long as M has an inverse, since S ;1 = RM ;1 R t . Finally, the pseudo-inverse matrix is

L + s = ; 1 n (A t A) ;1 A t + p] Ti Ti ( 28
)
where

T i = RM ;1 b p i R t A t (29)

Computation of the velocity s c r e w

Prior to the calculation of the velocity screw, we must de ne the error vector e between the current feature vector s and the desired one s :

W = n X i=1 M ;1 b p i R t A t A(R ; R ) b p i (33)
According to the properties of rotation matrices, the above di erence is R ; R = ;R u] 2 (1 ; cos ) ; R u] sin (34) where u are the axis and the angle which correspond to the rotation matrix product R t R .

If is small, we c a n t a k e the approximations sin and cos 1. Then, (38) Thus, the center of gravity of the object moves along a straight line trajectory from its initial to its nal position in the camera frame. Consequently, t h e object is most likely to remain in the camera eld of view during the whole task.

W n X i=1 M ;1 b p i R t A t A ; ;R u] b p i M ;1 n X i=1 b p i R t A t AR b p i u (

The use of 3D point coordinates

The presented control law utilizes the pixel coordinates and depth estimation of each point, which a r e the same requirements as those of classic image-based visual servoing. However, additional interesting results can be obtained if the feature vector is composed of the coordinates of the points of the object, i.e., s = ( : : : p t i : : : ).

One should note that any additional information is needed for the computation of such coordinates, which are obtained from pixels and depth estimates as shown in Equation 20.

The interaction matrix

Instead of repeating the development, all of the previous theoretical results can be used by means of simply replacing A by the identity matrix. Thus, the new interaction matrix is L s = (39) and its pseudo-inverse is

L + s = ; 1 n I + p] Ti Ti (40)
where

T i = RM ;1 b p i R t (41) M = n X i=1 b p i 2 (42)

The velocity screw

The same approximation of the velocity screw of the camera can be obtained. Nonetheless, an exact computation of the velocity screw for any rotation can be obtained if M is a diagonal matrix (some objects for M is a diagonal matrix are, e.g., the tetrahedron, the square and the cube).

From the properties of skew-symmetric matrices,

M = n X i=1 b p i b p i t ; n X i=1 b p i t b p i I (43)
Since both M and the second sum are diagonal, then the rst sum is diagonal too, i.e., where is a positive scalar. Then, from Equations 33, 34 and 44,

W = ;M ;1 n X i=1 b p i u] 2 b p i (1 ; cos ) + M ;1 n X i=1 b p i 2 usin = ;M ;1 n X i=1 b p i uu tb p i (1 ; cos ) + usin = M ;1 u] n X i=1 b p i b p i t u(1 ; cos ) + usin = M ;1 u] u(1 ; cos ) + usin = usin (45)
The velocity screw of the camera is, v = ; (p ; p) + p] Rusin Rusin (46) which i s v alid for all u .

As shown in Equation 38, the object will follow a straight line path in the camera frame. However, care must be taken if 2 < , since sin decreases. Consequently, the camera accelerates during the rst phase of the motion until = 2 then decelerates until convergence. This may cause unexpected behavior if the initial velocity of the camera is near to the maximum limit in order to improve c o n vergence times.

Finally, i f = then the camera does not rotate at all, which will produce inconsistent results if the translational motion does not change the orientation of the object with respect to the camera, e . g . if the object is positioned along the straight line which joins the current and desired camera positions.

Experimental results

The control law has been implemented on a robotic platform with 6 degrees of freedom and an eye-in-hand con guration. The target object is composed of four points which de ne a tetrahedron. Three-dimensional coordinates are obtained from the pose of the object, which is extracted from the images and an internal model with Dementhon's algorithm 2].

Initial and desired poses of the camera are shown in Table 1. As explained in 6], this is a di cult task since the displacement is important and the object moves towards the border of the image. Figures 1, 2 and 3 depict the trajectories of the object in the image plane. Since the involved rotation is high, the points of the object do not follow straight lines from the initial to the nal positions. Though trajectories are more curved in the proposed approaches than in image-based visual servoing, the object reemains in the eld of view of the camera during the complete servoing task.

In addition, the trajectory of the camera in the proposed approaches is not as elongated as in the imagebased scheme, thus preventing the manipulator from reaching its joint limits. Camera trajectories in an absolute frame are depicted in Figure 4.

Conclusion

The use of new visual features in image-based visual servoing has been proposed, without the need of additional information, i.e., only the pixel coordinates and the estimated depth of each point is used.

Such new features result in a linear control law, from which the computed screw can be obtained, resulting in a motion of the object along a straight path in the camera frame thus, the object is most likely to remain in the eld of view of the camera during the visual servoing task.

The use of 3D coordinates of the points of the object, obtained from the image features, the estimated depth and the camera parameters, is proposed too. We demonstrate that it is a particular case of the rst proposed scheme. Nevertheless, an exact value of the velocity screw is obtained, provided that an additional geometric property of the target object is met.

Both proposed approaches exhibit a better behavior than the classic image-based one, since the trajectories of the camera in the absolute frame are less elongated when the change of orientation is high. Thus, the risk of the robot going out of its joint limits is lowered.

Future work includes a comparative study of the e ciency and robustness of the proposed schemes and the classic image-based and position-based visual servoing approaches, particularly in face of camera calibration errors.

2

  Combining pixels and depth information

PoseFigure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: Object trajectory in the image when using image point features.

 7 7 5 (30)

Acknowledgements

Support from Fundaci o Caixa Castell o -Bancaixa for a temporal stay of E. Cervera in Blaise Pascal University is gratefully acknowledged.