I Iordanoff

D André

J L Charles

P Viot

2011, Presqu'île de Giens (Var) Méthode des éléments discrets : des problèmes multi-corps aux problèmes d'endommagement dynamique complexes

Keywords:

La méthode des éléments discrets est présentée comme une alternative aux approches de type mécanique des milieux continus pour aborder certains aspects liés aux problèmes dynamiques, notamment la multi fracturation de matériaux fragiles. Des exemples liés à l'usinage des composites et au surfaçage du verre sont présentés. Une extension de la méthode pour étudier le comportement des mousses est ensuite proposée et imagée par des premiers résultats.

Introduction

Il est important de commencer par rappeler une évidence : la méthode des éléments discrets est basée sur la description de l'état physique (vitesse, position, température, moment magnétique, potentiel électrique …) d'un grand nombre d'éléments discrets qui forment le milieu étudié. Cette méthode n'est pas basée sur une description continue de la matière. Elle est donc particulièrement bien adaptée pour décrire l'évolution d'un milieu subissant des phénomènes discontinus :

-comportement d'un matériau naturellement discret comme le sable, les poudres, les granulats …, -prise en compte des hétérogénéités et des anisotropies des matériaux, composites, mousses …, -problème de multi fracturation comme l'abrasion, les impacts de glace, l'usinage des composites, la rupture des matériaux fragiles, -génération de débris et prise en compte du rôle des débris dans les problèmes d'usure, d'avalanche de neige ou de roches … Il existe plusieurs types de méthodes par éléments discrets. On présente brièvement les différentes approches. Puis on détaillera particulièrement une approche basée sur des éléments discrets sphériques, des lois d'interactions régulières et une résolution explicite des équations de Newton. On montrera pourquoi cette approche est particulièrement bien adaptée pour aborder les problèmes dynamiques cités plus haut. Une extension de la méthode pour aborder le comportement dynamique des mousses sera ensuite présentée.

Méthodes des éléments discrets.

La méthode des éléments discrets consiste à simuler un matériau comme un assemblage d'un grand nombre d'éléments en interaction. La méthode la plus simple à mettre en oeuvre est la méthode du réseau de ressort ou 'lattice method' qui consiste à connecter un nuage de noeuds par des liaisons au comportement connu. On peut citer les travaux de Schlangenet al [START_REF] Schlangen | Fracture simulations of concrete using lattice models: computational aspects[END_REF], Van Mier et al [START_REF] Van Mier | Fracture mechanisms in particle composites: Statistical aspects in lattice type analysis[END_REF] ou Cusatis et al. [START_REF] Cusatis | Confinement-shear lattice model for concrete damage in tension and compression: II. Computation and validation[END_REF] sur la fracture des bétons. Généralement, les noeuds n'ont pas de masse et ces approches sont mal adaptées aux aspects dynamiques. D'autre part, la rupture ne génère pas de surface libre. Il n'est donc pas possible de gérer les éventuels contacts générés par la création des fissures. Une version plus complexe consiste à décrire la matière comme un ensemble de particules ayant des interactions de contact (liaisons unilatérales) et/ou des liaisons bilatérales. Les particules ont des formes sphériques [START_REF] Cundall | discrete numerical model for granular assemblies[END_REF][5], ou polygonales [START_REF] Kun | Fragmentation of colliding discs[END_REF] [START_REF] Delaplace | Discrete 3D model as complimentary numerical testing for anisotropic damage[END_REF]. La gestion des contacts peut se faire par des lois d'interaction régulières [START_REF] Cundall | discrete numerical model for granular assemblies[END_REF][5][6] [START_REF] Delaplace | Discrete 3D model as complimentary numerical testing for anisotropic damage[END_REF] ou non régulières [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. La méthode développée pour aborder les problèmes de rupture dynamique est basée sur une description par éléments discrets sphériques, des lois d'interaction régulières et un algorithme de résolution explicite des équations de Newton. Cette approche permet :

-de prendre en compte les effets d'inertie (par rapport aux méthodes lattice), -de gérer très facilement les contacts (simplicité de la recherche des contacts entre sphères), -de prendre en compte simplement la déformation des matériaux (lois d'interaction régulières), -de résoudre efficacement des problèmes aux petites échelles de temps (résolution explicite), Elle s'est montrée très efficace pour aborder des problèmes tribologiques complexes tel que l'usure ou le comportement thermo mécanique des interfaces. Cette méthode va être brièvement décrite.

Modèle mécanique

Les lois d'interaction définissent les propriétés mécaniques du matériau à étudier. Elles sont divisées en deux familles (Figure 1) : 1. Les joints cohésifs, qui agissent de façon bilatérale et qui représentent la cohésion d'un matériau solide non endommagé ; 2. Les lois de contact, qui agissent lorsque les deux particules en contact n'ont pas (plus) d'interaction de type joint cohésif.

Joint cohésif

Pour simuler le comportement mécanique d'un matériau, un joint cohésif entre les particules est créé. Le comportement de ce joint cohésif peur être plus ou moins complexe. Un joint cohésif de type micro poutre permet par exemple de retrouver de façon assez aisée le comportement mécanique statique et dynamique d'un matériau élastique. Ceci est détaillé par les travaux de D. André présentés dans ce même colloque [START_REF] André | Modèle par éléments discrets pour l'étude du comportement dynamique d'un matériau élastique[END_REF]. Les exemples détaillés dans cet article considèrent un joint cohésif plus simple, avec une composante normale F 1 (Figure 1b) et une composante tangentielle F 2 (Figure 1c). Le joint est créé initialement au point de contact des éléments discrets lorsque le matériau n'est pas sollicité. Le joint a une raideur tangentielle k et une raideur normale K. Dans cette étude, on considère une rupture fragile des joints cohésifs : un joint est rompu seulement lorsque la force équivalente F eq atteint une force critique F rt (en traction) ou F rc (en compression) avec

2 2 2 1 3F F F eq

Force de contact

La force de contact est divisée en trois contributions : répulsion, adhésion (les deux forces sont conservatives en énergie) et dissipation. La connaissance de l'interaction géométrique permet le calcul de la force de contact. La force de répulsion F r est simulée par un ressort élastique de raideur K, la force d'adhésion est caractérisée par une constante , l'énergie dissipée par le contact est représentée par une force visqueuse F d Avec le coefficient d'amortissement,  la vitesse d'impact et M ij la masse équivalente des deux éléments en contact.

* K F r a F  * . 2 ij d M K F

Intégration numérique

Un algorithme d'intégration explicite doit être utilisé pour aborder les problèmes de dynamique. L'algorithme choisi est un algorithme de Verlet vitesse qui permet une conservation de l'énergie.

Obtention de l'échantillon géométrique

Dans le cas de la fracturation des composites (par usinage), les fibres sont construites géométriquement grâce à des alignements d'éléments discrets sphériques. Lorsqu'on veut simuler un matériau isotrope, il est important de considérer un domaine d'éléments discrets le plus compact possible afin d'éviter les vides artificiels générés par un empilement de sphères. De plus, les éléments discrets sphériques doivent avoir une taille différente de façon à éviter les effets de structuration cristalline. Une répartition de taille de 15% autour d'une taille moyenne est la limite inférieure qui permet d'éviter toute organisation [START_REF] Poschel | Granular Gases[END_REF]. Il n'existe pas d'algorithme géométrique qui permettent d'obtenir directement un tel domaine. La méthode utilisée ici consiste à générer un domaine d'éléments discrets sans contacts. Puis de cisailler sous pression le domaine afin de le compacter à la compacité voulue et enfin de créer les joints cohésifs qui permettent de créer le matériau (Figure 2)

3 Applications.

Etude des dommages en sous surface créés par le procédé de doucissage de la silice.

Le matériau étudié est la silice utilisée pour réaliser les optiques des lasers de puissance. Les détails de cette étude sont présentés dans [START_REF] Iordanoff | A discrete element model to investigate sub-surface damage due to surface polishing[END_REF].

Usinage des composites à fibres longues

L'usinage des composites à fibres de carbone s'apparente plus à un processus de fracturation fragile qu'à la formation de copeaux sous l'effet de cisaillement plastique intense (cas classique des alliages ductiles). Pour modéliser les fibres de carbone avec le modèle simple présenté, on considère deux rangées d'éléments discrets par fibre ce qui permet de prendre en compte la flexion de la fibre. La matrice n'est représentée que par des joints cohésifs liants les éléments discrets représentant les fibres. L'outil est lui aussi modélisé par éléments discrets.

Les détails de cette étude sont donnés dans [START_REF] Iliescu | A discrete element method for the simulation of CFRP cutting[END_REF]. Les propriétés de l'outil, des fibres et de la matrice sont données dans le tableau 3. Le matériau composite étudié est monodirectionnel.

On étudie l'influence de l'orientation des fibres par rapport à la direction de coupe. Les résultats numériques sont comparés à une étude expérimentale menée sur une raboteuse instrumentée. Le procédé est filmé à l'aide d'une caméra rapide. Lorsque les fibres sont orientées à 90° (Figure 5 droite), la coupe s'effectue par flexion des fibres ce qui induit une profondeur de fissures entre les fibres très significative. Pour la coupe à 0°, le flambement des fibres conduit à des ruptures dans le sens de la coupe ce qui conduit à un bon état du matériau coupé. En revanche, les efforts sur l'outil en face de coupe sont très important (figure 5 gauche). La coupe à -45° s'effectue par flexion importante de paquets de fibres puis rupture bien en dessous de La coupe à 45° s'effectue par compression puis cisaillement des fibres (Figure 6 droite). Les débris formés sont très petits. Une légère flexion des fibres génère une usure par abrasion de l'outil en face de dépouille. la surface théorique de coupe (figure 6 gauche). Cette configuration est celle qui donne les plus mauvais états de surface liés à des défauts en profondeur. L'ensemble de ces résultats montrent la cohérence du modèle numérique avec les observations expérimentales. De nouveau, la complexité des phénomènes mis en jeu sont particulièrement bien captés par l'approche discrète.

Etude du comportement dynamique de matériaux alvéolaires.

L'endommagement de structures alvéolaires est un couplage complexe entre l'écrasement des cellules et la propagation de fissures dans le matériau. D'autre part, la structure alvéolaire inhomogène du matériau conduit à une hétérogénéité des contraintes qui participe aux phénomènes de localisation des déformations. Ces couplages dépendent bien sur des matériaux et de la géométrie de l'impacteur. Afin de prendre en compte l'écrasement des cellules, le modèle par éléments discrets à été étendu. Un élément discret représente une alvéole et doit avoir la propriété de s'écraser sous l'effet du chargement local.

Modèle 2D pour l'écrasement des cellules

Ce modèle présenté en détail dans [START_REF] Viot | Multiscale description of polymeric foam behavior: A new approach based on discrete element modeling[END_REF] est basé sur le calcul d'un tenseur de force calculé sur chaque élément discret :

z z L z F x L F z L x F x L x F D p 1 F
x et F z sont les composantes 2D de la force d'interaction entre deux éléments discrets, L x et L z les composantes du vecteur qui lie leurs centres, D le diamètre de l'élément discret considéré. A partir de ce tenseur, une force équivalente F eq , similaire à une contrainte équivalente de Tresca est calculée et comparée à une force limite d'écrasement F ec . Lorsque F eq atteint F ec , le diamètre de l'élément discret diminue jusqu'à obtenir l'équilibre F eq =F ec où lorsque D=D min , un diamètre minimum caractérisant la cellule complètement écrasée.

Influence de l'impacteur et du rapport force de rupture/force d'écrasement.

Dans cette étude qualitative, on travaille avec des données adimensionnelles afin de pouvoir comparer facilement les valeurs relatives des propriétés des matériaux. Les données physiques de référence sont le rayon moyen, la masse moyenne et la raideur des éléments discrets. Le domaine étudié a une hauteur de 130. L'impacteur est triangulaire avec 3 angles différents. La vitesse d'impact est de 0.001 et le déplacement maximum de 100. Les données des six simulations sont résumées dans le tableau 4. Lorsque la force à rupture des joints cohésifs est faible devant la force limite d'écrasement des cellules, le matériau s'endommage par fissuration (figure 7). L'apparition de la fissure est d'autant plus précoce (déplacement DL=18) que l'impacteur est plat. La force maximum à l'impact est aussi obtenue pour l'impacteur le plus plat (figure 9 gauche). Lorsque la force à rupture des joints cohésifs est égale à la force limite d'écrasement des cellules, le matériau s'endommage par écrasement des cellules puis par fissuration (figure 8). L'apparition de la fissure est d'autant plus précoce (déplacement DL=25) que l'impacteur est plat. La force maximum à l'impact est aussi obtenue pour l'impacteur le plus plat (figure 9 droite). Cependant, l'écrasement des cellules permet dans tous les cas de limiter la force maximale à l'impact et de retarder l'apparition de la fissuration.

Conclusion

Cette article présente les spécificités de la méthode des éléments discrets et développe plus particulièrement une méthode basée sur des éléments discrets sphériques. On montre comment la gestion de joints cohésifs et de contact permet d'étudier des phénomènes complexes liés à la multi fracturation. Enfin, une proposition est faite pour adapter ce modèle afin de modéliser des matériaux alvéolaires. Dans tous les cas étudiés, les résultats qualitatifs correspondent bien aux observations expérimentales. Cette étude montre donc la capacité des modélisations par éléments discrets à aborder des problèmes pour lesquels la mécanique continue est confrontée à des difficultés de modélisation. Un effort important doit maintenant être fait afin de caractériser correctement les paramètres micromécaniques du modèle qui permettront d'obtenir des résultats plus quantitatifs donc réellement prédictifs. C'est là un des enjeux majeur qui permettra à cette méthode d'apporter une réelle alternative aux méthodes continus pour la classe de problèmes présentés dans cet article.

Références

FIG. 1 -Fig. 2 -

 12 FIG. 1 -forces de contact et joints cohésifs

2 TAB. 2 -FIG. 4 -

 224 FIG.4-décroissance exponentielle des défauts dans la profondeur de la silice. Influence de l'abrasif.

TAB. 3 -

 3 FIG. 5 -coupe à 90° et coupe à 0°, comparaison numérique-expérience.

FIG. 7 -Fig. 9 -

 79 FIG.7-influence de la géométrie de l'impacteur pour un matériau à cellules ne s'écrasant pas

  La silice est modélisée par un ensemble d'éléments discrets liés par des joints cohésifs comme sur la figure 2.c. Les particules abrasives sont formées de quatre éléments discrets qui forment une plaquette. La surface simulée de l'outil est de 30 éléments discrets dans le sens de l'écoulement et de 10 éléments discrets dans le sens transverse. Le nombre et la taille des abrasifs testés sont résumés dans le tableau 2. Le calcul est mené jusqu'à ce qu'une demi-couche de silice (150 éléments discrets) soit usée par l'opération d'abrasion. Les joints cassés dans le volume de la silice sont alors comptés en fonction de la profondeur où ils se trouvent. La courbe figure4gauche montre que le nombre de joints cassés (de fissures) dans le cas d'une couche de petit grains d'abrasifs (CPA) est dix fois plus faible que dans le cas de l'abrasion menée avec peu de gros grains (PGA). L'influence de l'abrasif et la décroissance exponentielle des défauts dans la profondeur est qualitativement en adéquation avec les résultats expérimentaux de Sutatwala et al.[START_REF] Sutatwala | Sub-surface mechanical damage distribution during grinding of fused silica[END_REF] (figure 4 droite). Ce premier exemple montre la capacité de l'outil par éléments discrets à générer et étudier les réseaux de fissures.

	TAB. 1 -données matériau

L'outil est modélisé par une paroi rigide formée par une couche d'éléments discrets. Les données matériaux utilisées sont données tableau 1. Une pression normale de 1KPa est appliquée sur la surface supérieure de la silice et une vitesse de 1ms -1 est imposée sur l'outil. La taille des éléments discrets formant les particules abrasives est égale à la taille moyenne des éléments discrets formant la silice (petites particules) ou égale au double de cette taille (grosses particules).


A Micheletti
email: micheletti@ing.uniroma2.it-davide.cadoni@hotmail.it

D Cadoni

Presqu'île de Giens (Var) Design of Single-Layer Floating-Compression Tensegrities

Keywords:

Floating-compression tensegrities can conveniently be employed as deployable systems and as soft structures. Focusing on single-layer systems, we first identify the relevant parameters affecting their structural performances, by studying a simpler two-dimensional benchmark, then we derive new three-dimensional structures and present the design issues associated with them.

Introduction

Floating-compression tensegrities are prestressed systems composed by bars and cables such that bars are never connected to each other (see [START_REF] Motro | Tensegrity systems : The state of the art[END_REF][START_REF] Motro | Tensegrity : structural systems for the future[END_REF] for an introduction on tensegrities and [START_REF] Jáuregui | Controversial origins of tensegrity[END_REF] about the origin of these systems). As in conventional trusses, we can distinguish between single-layer and doublelayer systems : in the former case, by identifying a single surface with all the nodes of the system laying on it ; in the latter case, by identifying two closely spaced surfaces passing through all the nodes.

Since there are no bar-to-bar connections, these systems can be easily folded and it can be convenient to employ them as deployable structures [START_REF] Tibert | Deployable tensegrity structures for space applications[END_REF][START_REF] Motro | Tensegrity : structural systems for the future[END_REF]. For the same reason, their stiffness is often lower than that of conventional trusses, so that they can preferably be employed as soft/flexible structures. For example, they have been proposed for marine cages in aquaculture applications [START_REF] Jensen | Finite element analysis of tensegrity structures in offshore aquaculture installations[END_REF]. The tensegrities we derive here can also be designed to be employed as domes and roofs, given the architectural appeal of floatingcompression systems.

We present a study on single-layer system, such as two-dimensional or three-dimensional tensegrity balloons : systems where the characterizing surface is closed, with cables laying on the boundary of the convex hull of the structure, and bars inside. The term reflects the well-known balloon analogy by Fuller [START_REF] Buckminster | Tensegrity[END_REF], which fully applies here. Cables and bars are analogous to the membrane in tension and the air under pressure in common balloons. In two dimensions, these systems have been referred to as cablestrut polygons by Connelly [START_REF] Connelly | Second-order rigidity and prestress stability for tensegrity frameworks[END_REF]. A couple of examples is shown in Figure 1 (left and center). Externally FIG. 1 -Tensegrity polygons (left and center). Externally constrained system obtained by sectioning a cable-strut polygon (right). constrained structures can be easily obtained from tensegrity balloons by sectioning them in various ways, a two-dimensional example is shown in Figure 1 (right).

After discussing the stability of these systems, we show that there are two main parameters affecting their properties, namely the number of bars per unit area and the bars' overlap. In particular, buckling is more easily induced by external loads when the number of bars increases and when the overlap decreases. In order to prevent this effect, we evaluate the benefits of having additional internal bracing cables.

The study is continued by performing the design of two new types of tensegrity structures. The first type is derived from classical spherical tensegrities by modifying the typical layout of the set of cables, arranging them in a triangulated net. The resulting structures possess a higher stiffness if compared to classical ones. The second type is entirely original and it is obtained by a parallel juxtaposition of cablestrut polygons in three dimensions.

Modeling and computation

By modeling tensegrities as simple pin-connected bar frameworks, in the linear theory, the properties of solutions of the equilibrium and compatibility equations are determined from the fundamental subspaces of the equilibrium operator [START_REF] Pellegrino | Matrix analysis of statically and kinematically indeterminate frameworks[END_REF]. The system is characterized in terms of self-stress states, internal stresses balanced by null external loads, and mechanisms, nodal velocities causing no change in length of any member. We recall here the extended Maxwell's rule for a free-standing system, respectively, in two-and three-dimensions [START_REF] Pellegrino | Matrix analysis of statically and kinematically indeterminate frameworks[END_REF] :

2n n -3 -n e = n m -n s (2D) ; 3n n -6 -n e = n m -n s (3D) . (1)
This rule relates the number of nodes, n n , and members, n e , with the number of independent mechanisms, n m , and self-stress states, n s . Roth, Connelly and Whiteley were the firsts to develop the theory of tensegrity frameworks [START_REF] Roth | Tensegrity frameworks[END_REF][START_REF] Connelly | Rigidity and energy[END_REF], which are composed of rigid bars and inextensible, tension-only, cables. Tensegrity frameworks can be classified in terms of the numbers n m and n s [START_REF] Micheletti | A marching procedure for form-finding for tensegrity structures[END_REF]. If n s = 0 then the system cannot be stable without the action of external forces causing tension in all cables. When n s > 0 and n m = 0, then the system is stable under no external forces if there exists a self-stress state for which all cables are in tension. When n s > 0 and n m > 0, then the system is stable under no external forces if the following condition holds for all mechanisms and self-stress states :

∑ members t i j l i j (v i -v j ) 2 > 0 . (2)
In this condition, which is referred to as prestress-stability, t i j and l i j are the axial force and the length of the member connecting nodes i and j, while v i and v j are the velocities of these nodes.

When the members of a tensegrity system are considered to be linear elastic, then the above stability conditions are not enough to guarantee the actual physical stability of the system. A configuration will be stable if the total potential energy is at a minimum value. In other words, the tangent stiffness operator needs to be positive definite. We recall here that this operator can be decomposed additively into its material and geometric parts,

K T = K M + K G , (3)
the former depending on members' material properties, the latter depending on members' stresses. When a system under no external loads is stable regardless of its material properties and level of self-stress it is called super-stable [START_REF] Connelly | Tensegrity structures : Why are they stable ?[END_REF] (cf [START_REF] Zhang | Stability conditions for tensegrity structures[END_REF]).

Finding the stable configuration(s) of a tensegrity system constitutes the most important task when dealing with kind of structures. In the literature, this is always referred to as form-finding problem, which can be formulated in many non-equivalent ways (see [START_REF] Ruscica | Tensegrity : from didactical tasks to innovative marketing proposals[END_REF]) and it has been solved with a large variety of methods [START_REF] Tibert | Review of form-finding methods for tensegrity structures[END_REF][START_REF] Hernandez Juan | Tensegrity frameworks : Static analysis review[END_REF]. Here, we solved the form-finding problem for linear elastic tensegrity balloons by shortening cables in a large-displacement elastic analysis, and we employ a Newton-Raphson procedure. A similar conceptual strategy, the cocoon method, has been proposed and demonstrated in [START_REF] Sakantamis | The cocoon method : a physical modelling tool for designing tensegrity systems[END_REF] by making use of physical models. In that method, a deformable envelop is stretched against elongating bars, to find stable configuration according to the balloon analogy.

We employ the same Newton-Raphson procedure to perform static analyses. We remark that there are many other possibility for the choice of form-finding and static analysis methods, which can be more effective in many diverse situations. However, due to the smooth behavior of tensegrity balloons, this simple procedure was sufficient for our purposes.

A two-dimensional parametric benchmark

In this section we study cable-strut polygons, whose mechanical properties can be useful to understand more complicate three-dimensional cases. An important theoretical result have been established by Connelly about these and other similar polygonal systems [START_REF] Connelly | Rigidity and energy[END_REF] (cf [START_REF] Connelly | Tensegrity structures : Why are they stable ?[END_REF]), proving that any self-stressed convex polygon, having a positive stress in external members and negative (or zero) otherwise, is superstable (Polygon theorem [START_REF] Connelly | Tensegrity structures : Why are they stable ?[END_REF]).

A n-bars symmetric cable-strut polygon, with n > 2, possesses a n-fold axis of cyclic-symmetry, which is orthogonal to the plane of the polygon, and n axes of reflection symmetry, which belong to that plane (Figure 2 circle circumscribed to the polygon, the geometry of this system can be uniquely identified by the angle at center α depicted in Figure 2 (left). This angle measures the overlap between bars, a crucial structural parameter.

Symmetric cable-strut polygons always possess a single states of self-stress, n s = 1. Since there are 3n members and 2n nodes, by the extended Maxwell's rule there are n m = n-3 independent mechanisms. The self-stress state it is easily computed from the graphic representation of axial forces. We have

t b = -t c1 , t c2 = 2 cos π n t c1 , (4)
where t b is the force in bars, t c1 and t c2 are the forces in cables (Figure 1). Notice that the self-stress does not depend on α.

Although cable-strut polygons are super-stable, it can be seen that they can be affected by buckling under an external loading. To test their resistance against buckling, we performed numerical computations by loading these systems radially, as shown in Figure 2 (right). Results can be summarized as follows :

-when α (the overlap) increases, the buckling load increases ; -when n increases, then the buckling load decreases ; -when the self-stress level increases, the buckling load increases. Intuitively, the first two effects can be explained geometrically since the distance between a bar and the parallel cable close to it decreases, resulting in a 'thinner' structure. We also verified that any 'weak' cable-strut polygon can be significantly strengthened with additional bracing cables, as shown in Figure 3. Notice that these cables need only to be slightly in tension to be effective. Notice also that a high tension in these bracing cables would affect and lower the tension in the other external cables. We can obtain spherical systems from a net of cables with the layout similar to that of a geodesic sphere, and then by adding bars to it following the arrangement found in classical tensegrities by Fuller (Figure 4). As introduced earlier, by taking linear elastic elements, we progressively decreased the rest length of cables until a self-stressed system is obtained. To improve structural performances, it is beneficial to avoid cables with low self-stress, since they can be slackened by the action of external loads. For this reason we iteratively shortened the rest lengths of low stressed cables until a more uniform self-stress state is obtained. Figure 5 shows the result of this procedure. FIG. 4 -Layout of cables and bars for an improved spherical system.

The structure obtained in this way can be further modified by changing the layout of cables. In Figure 6 we replaced a 'star' of six cables, and the common node connected to them, with a set of three cables, with no significant geometrical modification of the rest of the structure.

At this point, we tried obtained domes from these structure by sectioning a part of it and by adding external constraints. This operation can be done in different ways. For example, we can section the members with a plane and pin the intersection points to ground, as shown in Figure 7 (left), or we can 'cut' along the cables, section the bars, and add external constraints in a similar way, as shown in Figure 7 (right).

Both domes may look too complicate with respect to their foundations, as in the first we have a large number of anchorages and in the second the foundations are not horizontal. We also tested a system where the starting layout of the cable-net is that of a geodesic dome, then we added bars to it and performed the form-finding procedure. However, this way we obtained poor results, since a suitable uniform self-stress state was much more difficult to achieve. In the cases we presented, the high symmetry of the system allowed us to easily obtain structures whose properties are intrinsically closer to the optimal ones.

Figure 8 shows the system obtained from a finer geodesic subdivision of the sphere. The structural performances of spherical systems are quite good. The interlocked layout of bars make these systems quite strong against overall buckling and the triangulated net of cables provides a higher stiffness with FIG. 5 -Form-finding results : self-stressed system obtained from the layout in Figure 4. A sphere has been plotted inside the structure for a better visualization of the members. FIG. 6 -Form-finding results : system obtained from a slightly modified arrangement of cables. A sphere has been plotted inside the structure for a better visualization of the members. FIG. 7 -Two different ways of obtaining a dome from the system on the left in Figure 5.

FIG. 8 -System obtained from a finer geodesic subdivision of the sphere, with the same kind of modification performed in the system in Figure 6 (right). A sphere has been plotted inside the structure for a better visualization of the members. respect to classical tensegrities. Moreover, these structures are highly redundant, i.e. they have a large number n s of self-stress states, so that breaking a cable would not cause the collapse of the entire system. On the other side, during construction, to precisely establish the desired self-stress state it is necessary to monitor the tension of at least n s members.

A new system

The idea for a new single-layer system originated from the seed structure in Figure 9 (left). By looking at this system in projection (Figure 9, right) we recognized a cable-strut polygon. Figure 10 shows a system where such a polygon is visible in a clearer way. Bars lay alternatively on different planes, and these planes are orthogonal to a central longer bar. By extending this concept further it is possible to obtain a whole family of structures, like the one in Figure 11. Figure 12 shows a dome obtained with the operation previously described. The stability properties of these systems are not as good as those of the spherical ones, hence, we needed to add a set of bracing cable in order to strengthen them. FIG. 12 -A dome obtained from the system in Figure 11.

Discussion and conclusions

From the analysis performed in this study, the design of single-layer floating-compression tensegrities results from the balance of different factors. It is clear that it is not possible to increase indefinitely the number n of bars in a cable-strut polygon, since this would increase the chance of buckling under external load. On the other side the number n cannot be too small to avoid an excessive bar length and to gain internal space. Due to buckling as well, the overlap parameter α cannot be too small and, at the same time, to avoid an excessively short length of cables, cannot be too large. Regarding the self-stress state, it should be designed to be at a level guaranteeing both the tension in cables under the envisaged load conditions and an adequate buckling strength, without causing excessive compression in bars. Lastly, the tension in bracing cables, when they are necessary, should be low to not affect the rest of the cables.

In this study we devised a methodology for the form-finding, the analysis and the design of singlelayer tensegrities, by identifying the relevant parameters which are typical of these systems and by pre-senting two new types of systems. This methodology can be easily applied to discover and realize new systems, both regular and irregular. Several variation are possible, involving connectivity operations and numerical procedures. Future work can regard the automatic generation and optimization of these systems, with the aid of the criteria outlined here. Another important issue is the establishment of stability results in the three-dimensional case.

  FIG.2-Overlap angle and self-stress of a cable-strut polygon (left). Buckling under a radial load (right).

FIG. 3 -

 3 FIG. 3 -Possible arrangements for bracing cables.

FIG. 9 -

 9 FIG. 9 -The seed structure. Perspective view (left), top view (right).

FIG. 10 -

 10 FIG.10 -Asmall system derived from the seed structure. Perspective view (left), top view (right).


N Bel Hadj

L Rhode-Barbarigos
email: giosef.rhode-barbarigos@epfl.ch

I F C Smith
email: ian.smith@epfl.ch

9-13 mai 2011, Presqu'île de Giens (Var) Active tensegrity structures with sliding cables -static and dynamic behaviour

Keywords: Tensegrity structures, sliding cables, active and deployable structures, dynamic relaxation

Tensegrities are spatial, reticulated and lightweight structures that are increasingly investigated as structural solutions for active and deployable structures. Tensegrity systems are composed only of axially loaded elements providing opportunities for actuation and deployment through changing element lengths. In cable-based actuation strategies, the deficiency of having to control too many cable elements can be overcome by connecting several cables. However, replacing sequences of cable elements by continuous cables sliding over joints significantly changes the mechanics of classical tensegrity structures. Furthermore, challenges emerge for structural analysis, control and actuation. In this paper a modified dynamic relaxation (mDR) algorithm is employed to investigate actuated tensegrity structures with continuous cables. The dynamic behavior of such structures is also investigated.

Introduction

Tensegrities are spatial, reticulated and lightweight structures that are composed of struts and cables. Stability is provided by a self-stress state established in tensioned and compressed elements. The tensegrity concept has applications in fields such as sculpture, architecture, aerospace engineering, civil engineering, marine engineering and biology [START_REF] Skelton | Tensegrity systems[END_REF]. Tensegrity structures have a high strength-to-mass ratio leading to strong but lightweight structural systems [START_REF] Wang | Free-Standing Tension Structures From Tensegrity Systems to Cable-Strut Systems Spon Research[END_REF][START_REF] Hadj | Design optimization and dynamic analysis of a tensegrity-based footbridge[END_REF]. Furthermore, tensegrities are flexible and easily controllable with low energy requirements [START_REF] Skelton | An introduction to the mechanics of tensegrity structures[END_REF]. Therefore, tensegrity structures are particularly attractive for active and deployable structures.

As a special type of prestressed pin-jointed framework, tensegrity structures are composed of axially loaded elements providing opportunities for actuation and deployment through changes element lengths. Length changes can be made to struts and/or cables through various actuation strategies. Strut-based actuation, employing telescopic members, has already been used in active tensegrity control applications. Fest et al. [START_REF] Fest | Active Tensegrity Structure[END_REF] experimentally explored shape control of a five-module large-scale active tensegrity structure. The actuation strategy was based on controlling the self-stress state of the structure through small movements of ten telescopic struts. This actuation was also used for self-diagnosis, self-repair and vibration control [START_REF] Adam | Active tensegrity: A control framework for an adaptive civil-engineering structure[END_REF][START_REF] Hadj | Dynamic behavior and vibration control of a tensegrity structure[END_REF]. Tibert and Pellegrino [START_REF] Tibert | Deployable Tensegrity Reflectors for Small Satellites[END_REF] numerically and experimentally investigated use of telescopic struts for the deployment of tensegrity reflectors. Generally, strut-based actuation becomes difficult to implement under conditions where internal forces are substantial, and required changes in shape are large. Furthermore, when strut-actuation is used for deployment, the structure may have no stiffness until it is fully deployed.

Cable-based actuation has been investigated in many research projects involving active and deployable structures. Bouderbala and Motro [START_REF] Bouderbala | Folding tensegrity systems[END_REF] studied folding of expanded octahedron assemblies and showed that cable-mode folding was less complex than strut-mode, although the latter produced a more compact package. Djouadi et al. [START_REF] Djouadi | Active Control of Tensegrity Systems[END_REF] developed a cable-control strategy for vibration damping of a tensegrity structure. Sultan and Skelton [START_REF] Sultan | Deployment of tensegrity structures[END_REF] proposed a tendon-control deployment strategy for tensegrity structures. Actuation is conducted in such way that the structure goes throughout successive equilibrium configurations. Smaili and Motro [START_REF] Smaili | Folding/unfolding of tensegrity systems by removal of self-stress[END_REF] investigated folding of tensegrity systems by activating finite mechanisms. A cable-control strategy was applied to a double-layer tensegrity grid. The proposed strategy was extended to the folding of curved tensegrity grids [START_REF] Smaili | Foldable/unfordable curved tensegrity systems by finite mechanism activation[END_REF].

Most research studies of deployment of tensegrity structures showed that cable-actuation strategy directs tensegrity structures to maintain stiffness as they move from one equilibrium position to another. There are, however, a few disadvantages with this approach. Tibert and Pellegrino [START_REF] Tibert | Deployable tensegrity masts[END_REF] argued that controlling cables is complicated, because of all the additional mechanical devices that are necessary. The deficiency of having to control too many cable elements can be overcome by connecting several cables together and using only one motor to control them [START_REF] Sultan | Deployment of tensegrity structures[END_REF]. This suggests that groups of individual active cable elements could be combined into continuous active cables. A single continuous cable can slide over multiple nodes through frictionless pulleys. This strategy has the advantage that fewer actuators are necessary for control. However, using continuous cables significantly changes the mechanics of classical tensegrity structures. Specifically the number of selfstress states may decrease and the number of internal mechanisms may increase eventually resulting in unstable configurations [START_REF] Motro | Tensegrity Systems[END_REF]. This leads to significant challenges for structural analysis, control and actuation.

Finite-element formulations for sliding cable-elements have been developed for modeling of suspension systems [START_REF] Aufaure | A finite element of cable passing through a pulley[END_REF][START_REF] Zhou | Finite element formulation for modeling sliding cable elements[END_REF][START_REF] Pauletti | Modelling sliding cables and geodesic lines through dynamic relaxation[END_REF][START_REF] Hincz | Nonlinear Analysis of Cable Net Structures Suspended From Arches with Block and Tackle Suspension System, Taking into Account the Friction of the Pulleys[END_REF] and fabric structures [START_REF] Pargana | Fully integrated design and analysis of Tensioned Fabric Structures: Finite elements and case studies[END_REF]. Chen et al [START_REF] Chen | Formulation and application of multi-node sliding cable element for the analysis of Suspen-Dome structures[END_REF] presented a formulation of multi-node sliding cable-element for the analysis of Suspen-Dome structures. Genovese [START_REF] Genovese | Strutture tensegrity -Metodi di Analisi e Ricerca di Forma[END_REF] investigated form-finding and analysis of tensegrity structures with sliding cables. The complete formulation of such systems was provided by Moored and Bart Smith [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF]. Moored and Bart-Smith [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF] formulated the potential energy, equilibrium equations and stiffness matrix for tensegrity structures with continuous (clustered) cables. The equilibrium equations of a tensegrity structure are non linear. Analysis can thus be carried out in an iterative manner through use of the transient stiffness method. Matrix methods generally require iterative assembling and inversion of large stiffness matrices. As a vector-based method, the dynamic relaxation method (DR) does not require such complexity. This method introduced by Otter [START_REF] Otter | Computations for prestressed concrete reactor pressure vessels using dynamic relaxation[END_REF] and Day [START_REF] Day | An introduction to dynamic relaxation[END_REF] in the mid-1960s is particularly attractive for modeling nonlinear structural behaviour. DR is an explicit iterative method for the static solution of structural-mechanics problems [START_REF] Underwood | Dynamic Relaxation[END_REF]. When the DR method is used, the static problem is transformed into a pseudo-dynamic one by introducing fictitious inertia and damping terms in the equation of motion. DR traces the motion of each node of the structure until it comes to rest in static equilibrium due to artificial damping. One of the advantages of this method is that the formulation of the structure stiffness matrix is not needed and hence the method is particularly suitable for problems with material and geometrical nonlinearities. DR has been used by many researchers to solve a wide variety of engineering problems. Furthermore, Barnes [START_REF] Barnes | Form Finding and Analysis of Tension Structures by Dynamic Relaxation[END_REF] showed that DR is particularly efficient for form finding and analysis of tension structures. Hundreds, perhaps thousands of structures such as cable-stayed bridges and large tent structures have been designed and then analyzed using DR. For a tensegrity structure with continuous cables, the uncoupled nature of the DR process makes it particularly straightforward to implement [START_REF] Wakefield | Engineering analysis of tension structures: theory and practice[END_REF].

In this paper, a modified dynamic relaxation (mDR) algorithm applicable to the analysis of tensegrity structures with continuous cables is proposed. The algorithm is validated by simulating load response and actuation of an active tensegrity beam. Results are compared with those obtained employing a stiffness-based algorithm to show effectiveness of the proposed methodology. The effects of introducing continuous cables sliding over frictionless pulleys into tensegrity structures are also investigated.

Modified dynamic relaxation method

The DR method follows from augmenting static equilibrium equations (Eq. ( 1)) by including inertial and damping terms (Eq. ( 2)):

In Eq. ( 1) and ( 2), u and v are the vectors of nodal displacements and velocities, M and C are the mass and damping matrix, F int is the vector of internal forces and F ext is the vector of external forces. Introducing the residual force vector R as the difference between external and internal forces at any time t, Eq. ( 2) becomes In DR, Eq. ( 3) is iteratively resolved. Due to damping, nodal velocities and accelerations decay to zero as the solution is approached. The transient response is attenuated leaving the steady state solution for the applied load. The static equilibrium is thus attained and the out-of-balance or residual forces come to zero.

To obtain the DR basic equations, the following central difference approximations are used for temporal derivatives: Using these approximations, Eq. ( 3) can be re-arranged to give the recurrence equations for nodal velocities where subscript i,x refers to the i th node and direction x (respectively, for directions y and z):

The velocities are then used to predict displacements at time (t + Δt):

The iterative process of DR method consists of a repetitive use of Eq. ( 5) and [START_REF] Adam | Active tensegrity: A control framework for an adaptive civil-engineering structure[END_REF]. The process continues until the residual forces are close to zero. The values of masses M and damping C have to be chosen to ensure that the recurrence scheme converge to the static equilibrium [START_REF] Underwood | Dynamic Relaxation[END_REF]. Generally, a diagonal mass matrix is used along with a mass proportional damping matrix. This strategy involves the determination of a critical viscous damping coefficient. An alternative damping approach is the use of kinetic damping. This approach is adopted in this study [START_REF] Hadj | Analysis of clustered tensegrity structures using a modified dynamic relaxation algorithm[END_REF].

The modification introduced to DR in order to adapt it to the analysis of tensegrity structures with continuous cables is concerned with the calculation of the residual forces. With the assumption that all sliding cables run over frictionless pulleys, each continuous cable in the structure can be assumed as a string of cable sub-elements all carrying the same tensile force. At each time step t, Eq. ( 6) is used to determine current node coordinates of the structure. The new member-length vector is thus easily determined. For continuous members, current element lengths can be calculated using Eq. [START_REF] Hadj | Dynamic behavior and vibration control of a tensegrity structure[END_REF].

int ext F (u) = F (1)
{ } int ext + + M v C v F (u) = F (2) t t t M v + C v = R (3) /2 -/2 2 t t t t t +Δ Δ = v +v v , /2 -/2 t t t t t t +Δ Δ = Δ v -v v (4) , , /2 /2 , , , , , , , 2 1 ( )
( ) 2 2 i x i x t t t t t i x i x i x i x i x i x i x M t C v v R M t C M t C +Δ -Δ Δ - = + Δ + Δ + (5)
/2 , , ,

.

t t t t t i x i x i x u u t v +Δ +Δ = +Δ (6)
= l Sl [START_REF] Hadj | Dynamic behavior and vibration control of a tensegrity structure[END_REF] where S is a transformation matrix. Moored and Bart-Smith [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF] called this transformation matrix the clustering matrix and showed that this matrix can be defined as follow: Subsequently, current internal forces in the m th -member of the tensegrity may be determined as follows:

where 0,m l and t t m l +Δ are rest and current length of continuous member m. m E , m A and 0 m t are Young modulus, cross-section area and initial prestress of continuous member m. Once the vector of internal forces in the tensegrity elements ( t ) is determined, the vector of internal forces in the tensegrity sub-members ( t ) can be computed employing Eq. ( 10).

The residual forces can thus be calculated. For any node i, the residual force in x-direction R i,x is calculated as the sum of the external force f ext,i,x and the x-component of the resultant force induced by the contributions of the N members meeting at node i (Eq. ( 11)). Similar equations may be written for the y and z coordinate directions.

A stiffness based mass matrix and a kinetic damping are adopted in this work [START_REF] Hadj | Analysis of clustered tensegrity structures using a modified dynamic relaxation algorithm[END_REF].The iterative process continues until static equilibrium is attained when the norm of the vector of residual forces goes below a fixed precision value.

Case study: a tensegrity beam

The performance of the proposed modified DR algorithm is demonstrated using the topology of a tensegrity structure studied by Moored and Bart-Smith [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF]. The structure is a tensegrity beam composed of an assembly of three prismatic modules commonly known as quadruplex modules. A quadruplex unit comprises four struts held together in space by 12 cables. When only rigid-body movements are blocked, this tensegrity unit has a unique state of self-stress and three infinitesimal mechanisms. The quadruplex unit can be reinforced by adding four reinforcing cables which remove mechanism modes from the structure. Three reinforced modules are connected together with no bar-to-bar connections, forming a class I tensegrity structure. A perspective view of the tensegrity beam is given in Fig. 1 where grayed lines denote bars and thin lines denote cables. Three of the end nodes of the structure (nodes 1, 3 and 6) are pinned forming a cantilever beam. Ten cable elements of the top surface and ten cable elements of the bottom surface of the structure are grouped into four continuous cables. In Fig. 2

( ) 0 0, 0, t t m m m m m m m E A t l l t l = - + (9)
T t = S t The tensegrity beam used in this study has a length of 212 cm, a width of 80 cm and a height of 30 cm. Struts are made of aluminum hollow tubes with a length of 85cm. Saddle, vertical and reinforcing cables have a length of 60, 48 and 40 cm, respectively. All cable members are made by stainless-steel. Detailed characteristics of used members are summarized in Table 2 Actuated bending deformation of the tensegrity beam is first studied. Bending deformation can be obtained through antagonist actuation of the top and bottom sliding cables. Actuation is performed by changing the effective rest length of actuated cables. For example, a prescribed actuation stroke of 20% is defined as a change in the rest length of 20%. Prior to actuation, top continuous cables are contracted by 2% in order to introduce self-stress in the structure and counteract deflection induced by self-weight. The tensegrity beam is then actuated through modifying lengths of the four actuated cables. Top sliding cables are actuated with 10% contraction while the bottom sliding cables are expanded by 10%. Contraction and elongation of actuated cables is conducted progressively in steps of 1%. Note that actuation is deliberately performed through small and slow steps so that inertia effects can be neglected when the structure is in motion. The actuation response obtained by the proposed DR method and transient stiffness method based on Moored and Bart-Smith formulation [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF] are compared in Fig. 3. Displacements at the top node 18 of the beam are displayed with respect to the actuation ratio in active elements (Fig. 3). It is observed that the result of the present study are in agreement with those obtained by Moored and Bart-Smith [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF]. When the top continuous cables have been actuated with 10% contraction while the bottom clusters have been expanded by 10%, this results in a 55cm tip deflection in the positive z-direction (Fig. 6). Load-displacement curves of the tensegrity beam obtained employing modified and unmodified dynamic relaxation are compared in Fig. 6. This comparison reveals that predicted displacement are under estimated if the effect of sliding cables is not considered in the analysis. Results indicate that the tensegrity beam with sliding cables is about 55% more flexible than an equivalent configuration having discontinuous cables. An eigenfrequency analysis showed also that including continuous cables induced a decrease in the tensegrity natural frequencies approaching 50% (Table 3). Furthermore, the nonlinear behaviour induced by large displacements is more accentuated in tensegrity structures with continuous cables.

(10) ( ) , ,, , , 1
t N t t t m i x ext i x j m i m t m m t R f x x l = = + - ∑ (11)

Conclusion

Dynamic relaxation is an attractive static analysis method for tensile and tensegrity structures. The classic DR method is extended here to accommodate tensegrity structures with continuous cables. In cable-based actuation of tensegrity structures, the deficiency of having to control too many cable elements can be overcome by connecting several cables. Continuous cables are thus used instead of discontinuous members. Sliding cables are assumed to run without friction through structural nodes. The concept of cable grouping is a scalable solution that can be employed for active structures that incorporate many actuated elements in order to reduce the number of actuators needed for active control and deployment. However, grouping cables significantly changes the mechanics of classical tensegrity structures and this leads to new challenges for structural analysis, control and actuation. This study shows that the uncoupled nature of the DR process makes it particularly attractive to apply to structures with sliding cable elements.

6 FIG. 2 -

 62 FIG. 1 -A perspective view of the tensegrity beam

FIG. 3 -

 3 FIG. 3 -Actuation response of the tensegrity beam

23 FIG. 4 -

 234 FIG. 4 -Deformed shape of the tensegrity beam due to 10% actuation of active cables

FIG. 6 -

 6 FIG. 5 -Load-displacement curves of the tensegrity beam

  , sliding cables are shown in dashed lines in a top view of the tensegrity beam. Each sliding cable is attached to two end nodes and runs frictionnlessly through four intermediate nodes. Details of the four continuous cables are given in Table1.

	S	ij	=	1, 0, ⎧ ⎪ ⎨ ⎩ ⎪	j if sub element e is part of i continuous element e otherwise -	(8)

  .

		TAB. 2 -Material characteristics for the tensegrity beam
	Member	Material	Cross-section area	Young modulus	Specific weight
			(cm 2 )	(kN/cm 2 )	(kN/cm 3 )
	Struts	Aluminum	2.55	7000	2.7 10 -5
	Cables	Stainless-steel	0.5026	11500	7.85 10 -5


G Anciaux

S B Ramisetti

J-F Molinari

Presqu'île de Giens (Var) Coupling molecular dynamics with finite elements Projective methods issues at finite temperature

Keywords: finite elements, molecular dynamics, finite temperature, coupling

In the context of coupling methods that match molecular dynamics with finite elements, we focus on the finite temperature and heat rates treatment done by Bridging Domain approaches. The energy dissipation occurring due to the projection of the fine scale onto the coarse scale is demonstrated during a sliding simulation of contacting surfaces. We describe the problems brought by the dissipation and propose a method that does not alter the temperature of the atoms in the coupling region. A coupling of the thermal fluxes with a continuum temperature field is also introduced and tested on simple examples.

Introduction

The main motivation for multiscale methods in the context of condensed matter simulations is to reduce MD simulation times. In many cases, such as realistic indentation simulation, a full atomistic model is out of reach (0.1µm 3 ≃ 10 9 atoms for crystalline copper). Thus a multiscale approach is mandatory as well as Parallel High Performance Computing efforts. A large body of work exists [START_REF] Curtin | Atomistic/continuum coupling in computational materials science[END_REF] to circumvent the difference in nature between molecular dynamics and continuum models, naturally inducing pathological wave reflections at the coupling interface. For instance, the Bridging Domain method consists in having an overlap region where the degrees of freedom have a mixed dynamic evolution. These methods have the merit to deal with wave reflections that occur because of the change of scales. In the overlap zone, the formulation of the Bridging Domain method damps (with an acceptable control) the kinetic energy not forwarded to the continuum. It has been proven [START_REF] Anciaux | Ghost force reduction and spectral analysis of the 1d bridging method[END_REF] that some parameters, such as element sizes, allow to damp problematic waves in order to ameliorate the reflective behavior of coupling interfaces. Thus, the coupling region can be seen as a parameterizable filter.

Nevertheless, a critical issue of coupled multiscale methods for solid mechanics appears in finite temperature configurations. Because temperature plays an important role in material response, efforts have been made to extend to finite temperature the Quasi-continuum method [START_REF] Dupuy | Finite-temperature quasicontinuum : Molecular dynamics without all the atoms[END_REF], CGMD [START_REF] Broughton | Coarse-grained molecular dynamics : Nonlinear finite elements and finite temperature[END_REF], and the Bridging scale [START_REF] Wagner | Coupling of atomistic and continuum simulations using a bridging scale decomposition[END_REF]. But these techniques only provide a correct accounting of the local temperature with no fluctuations between the two models. To the best of our knowledge, only few methods allow a two-way coupling of thermal fluxes with thermal energy provided by MD [START_REF] Wagner | An atomistic-to-continuum coupling method for heat transfer in solids[END_REF][START_REF] Kulkarni | A variational approach to coarse graining of equilibrium and nonequilibrium atomistic description at finite temperature[END_REF][START_REF] Mathew | Concurrent coupling of atomistic and continuum models at finite temperature[END_REF].

In this presentation, we will focus on the Bridging domain method and more generally on the direct coupling that uses an overlapping with a projection of the fine solution onto the coarse finite element mesh. There are two problems concerning this approach. First it is hardly acceptable to loose energy when the molecular dynamics ensemble (NVE) should be conservative. Secondly, the damping has a dramatic effect on the temperature field. Indeed, the Bridging Domain strategies are very efficient in dissipating the waves that are not representable by the coarse mesh. However it always brings the atoms in the coupling area to a temperature very close to absolute zero which creates very steep temperature gradients. We will demonstrate this behavior on a sliding contact simulation [START_REF] Anciaux | Contact mechanics at the nanoscale, a 3d multiscale approach[END_REF] at low temperature. Because heat is generated at an important rate during the sliding motion of two contacting surfaces, the damping effect of the enforced Bridging Domain coupling is particularly evidenced. Thus thermal fluxes were not playing the important role we are emphasizing now.

We will present a much simpler coupling strategy with no projection of the atomic part onto the smooth mesh field. The simple idea is to separate spatially the overlapping zone in two as plotted on figure 1. In Ω 1 , the atoms are dictating the nodal fields in least square sense whereas in Ω 2 the atoms have displacement and velocities interpolated from the embedding finite elements. The idea is that no dissipation of the lattice vibrations should be imposed. Such a method shall not introduce any spurious temperature gradient or ghost forces. We demonstrate the effectiveness of the approach by using the simple coupled situation. Boundary conditions are periodic on the lateral sizes whereas fixed Dirichlet conditions are imposed on the top and bottom boundaries. Atoms are positioned with taking into account the thermal expansion coefficients ensuring a thermodynamic equilibrium with no pressure acting on the walls. When no load is applied, it is clear that equilibrium should remain stable during a dynamic simulation. In order to test these points the LibMultiScale (http://libmultiscale.gforge.inria.fr/) parallel framework which is an investigation tool for 3D parallel concurrent coupling approaches was used. With a Bridging Domain method, a strong gradient is created in the entire domain while our approach allows to maintain the thermodynamic equilibrium. The absence of projection of fine to coarse scale leads to this improvement. The crucial conclusion is that damping strategies, as efficient as they can be in limiting wave reflections at low temperature, are inaccurate in the context of finite temperature simulations.

Still, wave reflections need a specific treatment so that the energy flux due to the highest frequency modes is not trapped in the atomistic region. Thus, in order to recover a proper energy balance, a coupling of the thermal energy flux with a continuum temperature field is introduced. As a pure phonon heat transport is assumed (dispersion heat equation) within the body, a set of thermostat can be used to drive the energy exchange between molecular dynamics and finite elements. An energy balance equation allows to derive the scaling parameters of such thermostats. Finally we will validate the methodology by the analysis of a coupled domain with heat flows and mechanical waves passing through the overlapping region.

FIGURE 1 -

 1 FIGURE 1 -Schematic of the overlap regions when no projection is introduced.


Dr M Collet

Dr M Ouisse

Presquîle de Giens (Var) New methodology for vibroacoustic energy diffusion optimization by mean of distributed electromechanical devices

Keywords: Distributed control, 2D Waves Dispersion, Bloch Theorem, Shunted piezoelectric, Mid-Frequency Optimization

The paper proposes a synthesis of new methodologies for developing a distributed, integrated shunted piezo composite for beams and plates applications able to modify the structural vibro acoustical impedance of the passive supporting structure so as to absorb or reflect incidental power flow. This design implements tailored structural responses, through integrated passive and active features, and offers the potential for higher levels of vibration isolation as compared to current designs. Optimization tools for designing shunt circuits is based on using dedicated numerical approach able to compute the multi-modal wave dispersions curves into the whole first Brillouin zone for periodically distributed 2D shunted piezo-mechanical systems. By using a specific indicator evaluating the evanescent part Bloch's waves, we optimize, in a second time, the piezoelectric shunting electrical impedance for controlling energy diffusion into the proposed semi-active distributed set of cells. A 3D modeling of semi-distributed distribution of the optimal smart metamaterial is used for validating the obtained cell design.

INTRODUCTION

Tailoring the dynamical behavior of wave-guide structures can provide an efficient and physically elegant approach for optimizing mechanical components with regards to vibration and acoustic criteria, among others. However, achieving this objective may lead to different outcomes depending on the context of the optimization. In the preliminary stages of a product's development, one mainly needs optimization tools capable of rapidly providing global design directions. Such optimization will also depend on the frequency range of interest. One usually discriminates between the low frequency (LF) range and the medium frequency (MF) range, especially if vibration and noise are considered. However, it should be noted that LF optimization of vibration is more common in the literature than MF optimization. For example, piezoelectric materials and other adaptive and smart systems are employed to improve the vibroacoustic quality of structural components, especially in the LF range [START_REF] Preumont | Vibration control of structures : An introduction[END_REF][START_REF] Elliott | Active Control of Sound[END_REF][START_REF] Banks | Smart material structures Modeling Estimation and Control[END_REF] even if ditributed transducers are used [START_REF] Batra | Multimode vibration suppression with passive twoterminal distributed network incorporating piezoceramic transducers[END_REF][START_REF] Tzou | A study of segmentation of distributed piezoelectric sensors and actuators, part i : Theorical analysis[END_REF]. Recently, much effort has been spent on developing new multi-functional structures integrating electro-mechanical systems in order to optimize their vibroacoustic behavior over a larger frequency band of interest [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF][START_REF] Thorp | Attenuation and localization of wave propagation in rods with periodic shunted piezoelectric patches[END_REF][START_REF] Beck | Broadband vibration suppression assessment of negative impedance shunts[END_REF][START_REF] Beck | Experimental analysis of a cantilever beam with a shunted piezoelectric periodic array[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF][START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF]. However, there is still a lack of studies in the literature for MF optimization of structural vibration. To that end, the focus of this study is to provide a suitable numerical tool for computing wave dispersion in 2D periodic systems incorporating controlling electronics devices. The main final aim is to allow their optimization in order to optimize vibroacoustic diffusion in 2D wave guides. Two numerical approaches can be distinguished for computing that dispersion : the semi-analytical finite element method (SAFE) and the wave finite element (WFE) method. The main disadvantage of the SAFE method is that FE used are not standard so they must be specifically defined for each application. Nevertheless, a large amount of FE has been developed since 1975 to compute dispersion curves of rails [START_REF] Akrout | Comportement dynamique déterministe et large bande des structures guidées[END_REF], laminated composite plates [START_REF] Manconi | The Wave Finite Element Method for 2-dimensional Structures[END_REF][START_REF] Berthaut | Multi-mode wave propagation in ribbed plates : Part i k-space characteristics[END_REF] and viscoelastic laminated composite plates [START_REF] Floquet | Sur les équations différentielles linéaires à coefficients périodiques[END_REF]. To avoid development of specific FE, the WFE method considers the structures as periodic in order to model, with standard FE, a period of the structure. By using the periodic structure theory (PST) introduced by Mead [START_REF] Akrout | Comportement dynamique déterministe et large bande des structures guidées[END_REF], an eigenvalue problem can be formulated from the stiffness and mass matrices of the FE model to find wave numbers of all the propagating waves. The WFE method has been successfully used to deal with wave propagation in two dimensional structures [START_REF] Bloch | Über die Quantenmechanik der Electron in Kristallgittern[END_REF][START_REF] Collet | Numerical tools for semi-active optimization of 2d waves dispersion into mechanical system[END_REF]. One of the main problem of all these approaches is the difficulty to compute the damped wave numbers in the whole Brillouin domain necessary for optimizing vibroacoustic behavior of smart periodic structures. After recalling the Floquet-Bloch theorems, we introduce a new numerical formulation for computing the multi-modal damped wave numbers dispersion in the whole first Brillouin domain of a periodical smart structure made of periodically distributed shunted piezoelectric patches. Based on this wave modeling, optimization of the electrical impedance of the shunted circuit is made in order to decrease group velocity of flexural waves. The obtained optimal impedance is also tested in controlling the HF response of a semi-distributed system.

PIEZO-ELASTO-DYNAMICAL APPLICATION OF THE FLOQUET-BLOCH THEOREM

In this section the application of the celebrated Floquet-Bloch theorem is presented for piezo-elastodynamic problems. Based on the well known results obtained by Floquet [START_REF] Joannopoulos | Photonic Crystals : Molding the Flow of Light[END_REF] in one-dimensional and later rediscovered by Bloch [START_REF] Kittel | Introduction to Solid State Physics[END_REF] in multidimensional problems, we propose an original application to bi-dimensional piezo-elastodynamical problem leading to very general numerical implementation for computing waves dispersion for periodically smart distributed mechanical systems incorporating electronic components and damping effects [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF].

The Bloch Theorem

The Bloch theorem gives the form of homogeneous states of Schrödinger equation with periodic potential. This theorem can be considered as a multidimensionnal application of the Floquet theorem [START_REF] Allaire | Bloch waves homogenization and spectral asymptotic analysis[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). We can also define the primitive cell as a convex polyhedron of R 3 called Ω x . The reciprocal unit cell is denoted by Ω k limited by the reciprocal lattice vector defined by the three vectors g j so that : r i .g j = 2πδ i j (δ i, j the Kronecker index). We note G = [g 1 , g 2 , g 3 ] the reciprocal lattice matrix in the later. If Ω x is the irreductible primitive cell, Ω k corresponds to the first Brillouin zone of the lattice. One can see [START_REF] Ichchou | Guided waves group and energy velocities via finite elements[END_REF] for details.

The Bloch Theorem stipulates that any functions 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 , u(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). One can also demonstrate 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 used spatial operator.

Based on that theorem one can define the expansion functions v m (x, k), called the Bloch eigen modes, such that they 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. One notes that the expansion coefficients u m (k) depend on the applied disturbance and also on the induced wave vector (see [START_REF] Houillon | Wave motion in thin-walled structures[END_REF] for details).

Application to Piezo-Elastodynamic

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 equation :

ρ ẅ(x) -∇σ(x) = 0 ∀x ∈ Ω x -∇D(x) = 0 ∀x ∈ Ω x (4)
where w(x) ∈ R 3 (Ω x ) is the displacement vector, σ represents the Cauchy stress tensor, ε = ∇ sym w = 1 2 (∇w T (x) + w(x)∇ T ) the Green strain tensor, D(x) 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 the electric field vector (V 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 = - S t 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. The piezo-elastodynamic equilibrium can also be written as :

ρω 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 ∈ Ω (9)
As the problem is 2D infinitely periodic, only electrostatic boundary conditions have to be considered on each cell : [START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF] 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 less boundary. The top electrode applied feedback voltage V o depends on the shunt characteristic and on the collected charges q o (7) and can be expressed in the Fourier space by :

   V = 0 ∀x ∈ S b V = V o ∀x ∈ S t D.n = 0 ∀x ∈ S l
V o (iω) = -Z(iω)q o (iω) (11)
By considering a primitive cell of the periodic problem Ω x and by using the Bloch theorem, we can compute the associated Bloch eigenmodes (3) and the dispersion functions by searching the eigen solutions of the homogeneous problem ( 8) and ( 9) as :

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

w n,k (x) V n,k (x)
, Ω x periodic functions. By using the same approach as proposed in [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF],

one can demonstrate that the finite element discretization of u n,k (x) are solutions of the generalized algebraic eigenvalues problem :

(K(Z(iω n (λ, φ)) + λL(φ, Z(iω n (λ, φ))) -λ 2 H(φ, Z(iω n (λ, φ))) -ω 2 n (λ, φ)M)u n,k (φ) = 0, (13)
where

k = k   cos(φ) sin(φ) 0   = kΦ, λ = ik, M and K(Z(iω 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(iω n (λ, φ))) is a skew-symmetric matrix and H(φ, Z(iω n (λ, φ))) is a symmetric semi-definite positive matrix. These matrices are linked to specific weak term of the problem's weak formulation as described in [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF]. When k and φ are fixed and Z does not depend on ω the system ( 13) 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] Mencik | Multi-mode propagation and diffusion in structures through finite elements[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] Maysenhölder | Körperschall-energie Grundlagen zur Berechnung von Energiedichten und Intensitäten[END_REF][START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF][START_REF] Collet | Modal Synthesis and Dynamical Condensation Methods for Accurate Piezoelectric Systems Impedance Computation[END_REF]. Nevertheless these approaches have been only developed for undamped mechanical systems that is to say represented by a set of real matrices. In this case, most of the previously published works present techniques based on the mesh of a real k-space (i.e k or λ and φ) inside the first Brillouin zone for obtaining the corresponding frequency dispersion and the associated Floquet vectors. For undamped system only propagative or evanescent waves exist corresponding to a of eigen solutions purely real or imaginary. Discrimination between each class of waves is easy. If a damped system is considered (K, L, H are complex frequency dependent) or frequency dependence of the electrical shunt impedance is considered, the obtained eigenvalue problem is not quadratic and a complex specific numerical methodology has to be implemented. Furthermore, evanescent part of propagating waves appear as the imaginary part of ω 2 n (λ, Φ). It then becomes very difficult to distinguish the propagative and evanescent waves but also to compute the corresponding physical wave's movements by applying spatial deconvolution.

Another much more suitable possibility for computing damped system, dedicated for time/space deconvolution and for computation of diffusion properties as defined by [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF][START_REF] Collet | Modal Synthesis and Dynamical Condensation Methods for Accurate Piezoelectric Systems Impedance Computation[END_REF], is to consider the following generalized eigen value problem :

(K(Z(ω) -ω 2 M + λ n (ω, φ)L(φ, Z(ω)) -λ 2 n (ω, φ)H(φ, Z(ω)))u n (ω, φ) = 0. (14)
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.

Computation of the group velocity and evanescence criterion

The main 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 a suitable criterion. The waves group velocities indicate how energy is transported into the considered system and allow to distinguish the 'propagative' and 'evanescent' waves. If one Bloch eigen solution (i.e u n (ω, φ), k n (ω)) is considered, the associated group velocity vector [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF] is given by :

C g n (ω, φ) = ∇ k ω = S e tot = I E tot ( 15
)
where : is the spatial and time average respectivelly on one cell and one period, S is the density of energy flux defined in [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF], I the mean intensity and e tot , E tot the total piezomechanical energy and its time average on a period (see [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF] for details). In this problem, we only consider mechanical energy transportation as the electrostatic coupling is decentralized and can be condensed as a mechanical interface as proved in [32] and generally computed in [33]. So we also compute the intensity vector I by :

I n = - ω 2 Re( Ω x C(ε n (x, ω, φ) + ikΞ n (x, ω, φ)).(w * n (x, ω, φ))dΩ/V ol (16)
where . * is the complex conjugate and V ol the domain volume. As the spatio-temporal average of the system Lagragian is null (see [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF]), the total energy average is approximated by only computing the kinetic energy average :

E tot = 1 2 ( Ω x ρω 2 w n (x, ω, φ).w * n (x, ω, φ)dΩ/V ol (17)
The group velocity vectors C g n (ω, φ) is 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 2 ( Ω x ρω 2 wz n (x, ω, φ)wz * n (x, ω, φ)dΩ)/V ol E tot (18)
with wz n (x, ω, φ) being the (Oz) component of vector w n (x, ω, φ).

Optimisation 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 is P1-91 PZT. The used methodology for optimizing the shunt impedance Z(iω) is based on the minimization of the maximal group velocity collinear to the wave number vector [START_REF] Akrout | Comportement dynamique déterministe et large bande des structures guidées[END_REF] for waves having a ratio of transported flexural kinetic energy (18) greater than 0.8. The used criterion can also be written as :

Crit(Z(iω), φ) = max n/Ind(n,ω,φ)>0.8 (C g n (ω, φ).Φ) (19)
The used numerical optimization of the criteria is based on a multidimensional unconstrained nonlinear minimization (Nelder-Mead).

We also optimize the criterion by considering any frequency dependent complex impedance. We present in figure 3(a) and 3(b) the obtained real parts and imaginary parts of the wave number kx n (iω) along (Ox) axis. The red circles mark the dispersions curves for Z = 0 and the blue crossed the optimal dispersion. The corresponding group velocities along (Ox) are presented in figure 4 while the real and imaginary parts of the optimal impedance are plotted in figure 5.

We immediately observe that the optimization of the shunt impedance leads to greatly modify the group velocity of the A o mode in the first part of the spectrum (i.e before the first bending band gap FIGURE 5 -Real and imaginary parts of the optimal impedance between 22 and 27 kHz). The bending waves also propagate energy with a lower velocity and can even been very low as evanescent waves before 6 kHz. The corresponding real and imaginary parts of the dispersion curves are slightly modified. We notice an increase of the imaginary values indicating an increase of the spatial decay rates. The optimal impedance values are almost real, and correspond to those obtained if a constant negative capacitance is used. The corresponding average value is -150.05 pC.V -1 . Some imaginary parts of the optimal impedance are negative which indicate that the optimization leads to provide energy to the system for controlling mechanical damping effect introduced with hysteretic damping ratios into the model, and, also, obtain a fully conservative system.

Validation on a periodically semi-distributed set of adaptive cells

In order to study the translation of the wave properties obtained by optimizing the previously described criterion on the dynamical response of a finite dimension system, we apply the optimal impedance on a finite set of shunted piezo-composite cells semi-distributed onto a plate system as described in figure 2. The harmonic response of this system is also computed at different frequency when optimal impedance is connected or not to each patches. We plot in figure 6 the obtained results in term of kinetic energy density.

These numerical results clearly show up a strong influence of the modifications of the wave dispersion (i.e on the group velocities) in the standing wave responses plotted in figure 6. We observe at 5 kHz an increase in the dynamic response of the system when optimal shunt is connected. A large part of the system energy remains in its left part where the applied forces is located. The semi-distributed interface also change the system admittance and filter wave diffusion by increasing is reflexibility property. At 22 kHz, the energy diffusion is clearly condensed into the left part of the system with a largely decreased amplitude compared to this obtained with open circuit. The structural dynamical admittance has been decreased by connecting the shunt circuits. Finally, at 60 kHz, the energy is concentrated on the right part of the system between the adaptive interface and the clamped boundary condition. We observe something similar to a wave trap effect. The average value is, one more time, largely decrease by using the shunt circuits.

Conclusions

This paper presents a numerical procedure able to compute the damped wave's dispersion functions in the whole first Brillouin domain of multi dimensionnal piezo-elastodynamical wave guides. The method was applied for determining the optimal impedance allowing to minimize the group velocities of the flexural waves. Based on this approach, some numerical test on a finite dimension system incorporating a semi-distributed set of shunted piezo-composite cells has been performed. We underline a strong influence of the designed shunt circuits in the dynamical response of the system. Even if the link between the obtained wave properties are not clearly established, we also demonstrated that our developed 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. 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 [34], or even foams or complex polymers behaviors. 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].

acknowledgments

This work was carried out with a grant of French agency ANR number NT 09-617542. We gratefully acknowledge the French ANR and CNRS for supporting this program.

FIGURE 2 -

 2 FIGURE 1 -Generic 3D piezocomposite periodic cell

FIGURE 3 -

 3 FIGURE 3 -Real parts (a) and imaginary parts (b) of the wave number kx n (iω) along (Ox), the red circles mark values obtained for Z = 0 and the blue crossed the optimal ones

FIGURE 4 -4

 4 FIGURE 4 -Group velocities along (Ox) direction ; the red circles mark values obtained for Z = 0 and the blue crossed the optimal ones

FIGURE 6 -

 6 FIGURE 6 -Kinetic Energy density for different source frequency ((a) 5 kHz, (b) 22 kHz, (c) 60 kHz). On the left with optimal impedance connected to the piezo-patches and on the right with open circuit

Références