# Approches multiéchelles en espace et en temps pour la simulation d'impacts de faible énergie sur les structures composites

Teddy Chantrait, Johann Rannou, Anthony Gravouil

► To cite this version:

Teddy Chantrait, Johann Rannou, Anthony Gravouil. Approches multiéchelles en espace et en temps pour la simulation d'impacts de faible énergie sur les structures composites. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717078

HAL Id: hal-01717078

https://hal.science/hal-01717078v1

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Approches multiéchelles en espace et en temps pour la simulation d'impacts de faible énergie sur les structures composites

Teddy Chantrait<fnref n="1" />, Johann Rannou<fnref n="1" />, Anthony Gravouil 2

**Résumé** — Les impacts de faible énergie peuvent être très néfastes pour les structures composites utilisées dans le secteur aéronautique notamment. En effet, ces derniers peuvent engendrer des endommagements importants (fissurations matricielles, ruptures de fibres, délaminages...) dans le cœur du composite ou sur la face opposée à l'impact alors que l'empreinte résiduelle qu'il laisse sur la face impactée est presque indétectable à l'œil nu. Les endommagements causés peuvent donc conduire à la ruine prématurée de la structure alors qu'ils peuvent passer inaperçus lors d'une inspection visuelle (notion de Barely Visibale Impact Damage). On comprend donc l'importance pour les constructeurs de maîtriser ce genre de situations. La simulation numérique de ce type d'impact est donc très utile mais elle reste pour le moment difficile à mettre en œuvre au niveau industriel car elles sont très délicates à gerer. En effet, d'un coté les sources de non-régularités introduites dans les modèles (contact, lois d'endommagement adoucissantes, modèles de zone cohésive...) rendent les solveurs implicite quasiment inopérant. Et de l'autre, le pas de temps nécessaire à la stabilité des algorithmes explicites sont si faible que les temps de calculs deviennent rédhibitoires pour des problèmes de grandes tailles. Nous présentons ici le développement d'une stratégie de couplage peu intrusive entre le code de calcul de dynamique implicite Zébulon (code co-développé par les Mines ParisTech, l'Onera et NW Numerics & Modeling) et le code de calcul de dynamique explicite Europlexus (code développé par le CEA et le Centre Commun de Recherche d'Ispra en Italie) basée sur la méthode GC [1] [2]. Ce couplage permet d'allier les performances des deux types d'algorithmes évoqués et il est appliqué sur l'exemple d'un panneau composite raidi impacté.

**Mots clés** — Couplage de codes, implicite, explicite, faible intrusivité, impact faible énergie, matériaux composites.

# 1 Introduction

Les impacts de faible énergie peuvent être très néfastes pour les structures composites utilisées dans le secteur aéronautique notamment. En effet, ces derniers peuvent engendrer des endommagements importants (fissurations matricielles, ruptures de fibres, délaminages...) dans le cœur du composite ou sur la face opposée à l'impact alors que l'empreinte résiduelle qu'il laisse sur la face impactée est presque indétectable à l'œil nu. Les endommagements causés peuvent donc conduire à la ruine prématurée de la structure alors qu'ils peuvent passer inaperçus lors d'une inspection visuelle (notion de Barely Visibale Impact Damage). On comprend donc l'importance pour les constructeurs de maîtriser ce genre de situations. La simulation numérique de ce type de phénomène peut donc être une aide importante en particulier pour orienter les campagnes d'essais et pratiquer ainsi le "virtual testing". Divers travaux de recherche sont donc conduis dans la communauté scientifique et industrielle afin d'être en mesure de simuler ces impacts. L'utilisation de ces simulations reste cependant pour le moment difficile à mettre en œuvre au niveau industriel car elles sont très délicates à gerer. En effet, les sources de non-régularités introduites dans les modèles (contact, lois d'endommagement adoucissantes, modèles de zone cohésive...) rendent les algorithmes implicites difficiles à faire converger. Les simulations réalisées dans ce cadre sont donc majoritairement faites avec des schémas d'intégration en temps explicites [3] [4]. Ces derniers permettent en effet de mieux prendre en compte ces sources de non-régularités mais ils nécessitent d'utiliser un pas de temps dépendant notamment de la taille du plus petit élément du maillage afin de garantir leur stabilité. Ceci conjugué aux maillages utilisés ici, relativement fin à l'échelle de la struc-

<footnote n="1">¹Onera DMSM/MNU, Châtillon, France, {teddy.chantrait, johann.rannou} @onera.fr</footnote>

<footnote n="2">²Insa de Lyon LaMCoS, Villeurbanne, France, anthony.gravouil@insa-lyon.fr</footnote>
ture pour capter les non-régularités de la solution, conduit à réaliser un grand nombre d'incréments. Ces derniers sont alors réalisés sur la totalité du maillage ce qui devient rédhibitoire sur des problèmes de grande tailles. On observe par ailleurs que les non-régularités présentes dans cette problématique sont relativement localisées dans l'espace et dans le temps.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Schéma de principe d'une méthode multiéchelle en espace et en temps appliquée à la problématique des impacts de faible énergie sur composites.</figcaption></figure>

Il peut donc être intéressant de décomposer un maillage complet en plusieurs sous-domaines de manière à adapter la discrétisation en espace et en temps à la physique locale que rencontre le sous-domaine. La stratégie de cette décomposition est ainsi de focaliser l'effort de résolution numérique sur les zones qui rencontrent les phénomènes non-linéaires [5]. L'utilisation d'un maillage éléments finis fin ainsi qu'un petit pas de temps pour les parties de la structure où la solution est non-régulière comme la zone impactée est nécessaire et une formulation en temps explicite est quasiment indispensable, alors que pour les parties complémentaires plus régulières, un maillage plus grossier ainsi qu'un pas de temps plus grand sont suffisants. L'utilisation d'une formulation temporelle implicite peut donc s'avérer judicieuse. De cette manière, les incréments explicites ainsi que le temps passé à intégrer la loi de comportement peuvent être limités ce qui peut potentiellement faire gagner en temps de calcul. Les travaux réalisés ici sont basés sur la méthode GC [1] et ils visent à développer un couplage qui soit le moins intrusif possible [6] entre les codes de calcul implicite Zébulon (code co-développé par les Mines ParisTech, l'Onera et NW Numerics & Modeling) et explicite Europlexus (code développé par le CEA et le Centre Commun de Recherche d'Ispra en Italie) de manière à mettre en œuvre cette stratégie qui est schématisée sur la figure 1. On y voit notamment la tranche d'une plaque impactée ainsi qu'un exemple de maillage. Ce maillage est ainsi décomposé en deux zones. La partie impactée (au centre) qui est traitée par Europlexus et qui dispose d'une échelle de temps numérique fine et la partie complémentaire (autour) qui est traitée par Zébulon et qui dispose d'une échelle de temps numérique plus grande.

## 2 Pincipes de la méthode multiéchelles en espace et en temps adoptée

Différentes méthodes multiéchelles en espace et en temps [1] [7] [8] [9] qui peuvent être mises en place avec la famille des schémas d'intégration en temps de Newmark notamment existent dans la littérature. Ces dernières peuvent être vues comme des extensions des méthodes de décomposition de domaine duales sans recouvrement classiquement utilisées pour faire du calcul parallèle [10]. En effet, les méthodes de décomposition de domaine consistent à découper spatialement une structure en plusieurs sous-domaine et à rechercher la partie de la solution globale dans chacun des sous-domaines de la manière la plus indépendante possible (un exemple de décomposition en deux sous-domaines est montré dans l'encadré en pointillé de la figure 2). Le cœur de ces méthodes de décomposition de domaine
consiste à déterminer la condition aux limites particulière qui doit être ajoutée sur l'interface des sous-
domaines. Cette condition aux limites doit en effet assurer la continuité globale de différents champs
ainsi que l'équilibre au niveau de l'interface. Un problème additionnel qui a la taille de la discrétisation
spatiale de l'interface est alors résolu pour la déterminer. Pour formuler ce "problème d'interface" on
postule la continuité d'une quantité et on recherche la continuité de la quantité complémentaire assurant
l'équilibre interfacial. A ce stade on peut mentionner que du point de vu discret, il n'est pas équivalent
d'imposer la continuité du déplacement de la vitesse ou encore de l'accélération. Le problème d'interface
peut alors soit être résolu directement après avoir construit l'opérateur d'interface explicitement ou bien
être résolu de manière itérative qui ne nécessitent pas la construction explicite de l'opérateur d'interface
(méthode du gradient conjugué) comme c'est la cas dans les méthodes de type FETI.

La particularité des méthodes multiéchelles en espace et en temps consiste en plus, à doter chacun des sous-domaines de sa propre échelle de temps numérique. Un schéma de principe est montré sur la figure 2 où l'on a une échelle de temps fine associée au sous-domaine de gauche et une échelle de temps grossière associée au sous-domaine de droite. On note $\Delta t$ le pas de temps de l'échelle fine et $\Delta T$ celui de l'échelle grossière. On défini également le rapport entre ces deux pas de temps que l'on note m par la suite. La recherche de la condition aux limites particulière satisfaisant la continuité de différents champs ainsi que l'équilibre interfacial est donc plus délicate. En effet, on peut chercher à écrire le problème d'interface de manière à ce qu'il vérifie ces conditions à chaque piquet de temps de l'échelle la plus fine ou bien sur ceux de l'échelle de temps grossière. On peut également faire cela non pas à des instants précis (aux piquets de temps) mais plutôt en moyenne. Ces différentes possibilités conduisent à différentes stratégies multiéchelles en espace et en temps. Nous avons choisis pour notre problématique d'utiliser la méthode GC [1]. Dans cette méthode, les conditions d'équilibre et de continuités interfaciales sont recherchées à chaque instant de l'échelle la plus fine en temps. Cela permet notamment de mieux gérer la non-linéarité comparativement à d'autres méthodes qui nécessitent la reconstruction de l'opérateur d'interface ainsi qu'un processus itératif. La vitesse interfaciale du domaine le plus grossier en temps n'est pas connus sur l'échelle de temps fine, elle est donc interpolée [2]. Le problème d'interface à résoudre est donc légèrement modifié par rapport à une méthode de décomposition de domaine mono-échelle en temps mais il reste analogue. Nous choisissons pour notre problématique de construire explicitement ce problème d'interface. En effet, les opérateurs qui interviennent dans la construction de ce problème sont analogues à la masse pour la partie explicite. De ce fait, l'opérateur d'interface reste constant au cour du calcul tant que le pas de temps reste inchangé et tant que la partie implicite reste linéaire. Ce choix parait donc judicieux compte tenu du nombre important d'incréments à réaliser sous ces hypothèses.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Schéma de principe des méthodes multiéchelles en espace et en temps.</figcaption></figure>
## 3 Mise en œuvre de la méthode multiéchelles en espace et en temps GC

La particularité de l'approche proposée ici consiste à adopter une vision de type couplage de codes plutôt qu'une vision algébrique de la méthode multiéchelle en espace et en temps. Nous réalisons ainsi un couplage de codes de type structure/structure où chaque sous-structure est associée à un code de calcul et elle dispose ainsi de sa propre échelle de temps numérique. L'intérêt de ce type de couplage est de pouvoir étendre la méthode GC à des problèmes de grande taille tout en bénéficiant de la diversité des outils numériques déjà existant au sein de chacun des codes de calcul (algorithmes de contact, éléments de zone cohésive, modèles matériaux...). Nous pouvons ainsi voir le potentiel de la méthode en particulier sur la problématique des impacts de faible énergie où l'on met en jeux des sous-assemblages relativement représentatifs de l'industrie aéronautique ainsi que des problèmes de taille relativement importante (de l'ordre du million de degrés de liberté).

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Schéma de principe du couplage de code implicite/explicite développé sur la base de la méthode GC.</figcaption></figure>

La mise en place de cette méthode a permis de dégager les points d'entrée minimaux que l'interface des codes de calcul doit fournir pour mettre en œuvre cette méthode. Avec cette interface minimale la méthode peut être facilement mise en place sans être intrusive. Seul le développement de cette interface peu rester légèrement intrusif. Il s'agit des points d'entrée suivants :

- Une méthode permettant d'imposer une condition limite en effort au début de chaque incrément sur une zone du maillage.

- Une méthode permettant d'extraire les quantités cinématiques nodales (au minimum la vitesse) d'une zone du maillage en fin d'incrément.

- Une méthode permettant de valider ou d'invalider un incrément.

Par ailleurs, un protocole doit être choisie pour les deux premières méthodes afin que les quantités manipulées soient identiquement ordonnées du point de vue externe. Pour faire cela une méthode supplémentaire réalisant ces opérations de trace peut être requise. Des méthodes permettant d'accéder aux paramètres des schémas de Newmark, au pas de temps et éventuellement à l'opérateur d'interface, peuvent être également utiles mais elle ne sont pas fondamentales.

Cette interface à ainsi été développée sur les deux codes de calcul évoqués (Zébulon et Europlexus) comme on peut le voir sur le schéma de principe de la figure 3. Il s'agit des encadrés situés autour de chacun des codes de calcul qui ont ici été développés dans le langage Python. Un script contenant la méthode GC ainsi que les outils pour réaliser la résolution du problème d'interface à également été construit. Ce script contient en particulier les méthodes permettant de réaliser les échanges entre les différents codes par le biais de la bibliothèque MPI. Il détermine également lors de l'exécution où réaliser la résolution du problème d'interface. En effet, comme ce dernier doit être résolu sur l'échelle de temps la plus fine, il est préférable de le faire au plus près du sous-domaine fin en temps de manière à localiser les données dans l'espace et dans le temps. Différents cas tests ont alors été réalisés sur des structures 2D et 3D. Certains proposés dans [7] ont notamment été reproduis afin de valider l'implémentation de la méthode. Il faut également noter que le script développé permet facilement d'interchanger des codes de calcul dès l'instant que l'interface minimale existe.
## 4 Applications de la méthode de couplage à l'impact d'un panneau composite raidi

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Exemple de décomposition d'un panneau composite raidis en deux zones (1 sous-domaine implicite et 9 sous-domaines explicite).</figcaption></figure>

Des simulations avec différents rapports de pas de temps ont été effectuées sur un panneau composite raidi relativement représentatif d'un sous-assemblage d'aéronef. La figure 4 illustre la géométrie et une décomposition de domaine adoptées pour une partie des simulations présentées. On retrouve sur cette figure la partie centrale du panneau impactée qui est traitée avec le code de calcul explicite Europlexus ainsi que la partie complémentaire qui comporte l'échelle de temps relativement grande et qui est traitée avec le code implicite Zébulon. Les empilements choisis sont $[90/45/0/-45]_s$ pour la peau et $[90/45/0/0/-45]_s$ pour les renforts. Du contact est également présent dans la partie impactée traitée par le code de calcul Europlexus apportant ainsi une source de non-linéarité. La taille des éléments dans le plan principal géométrique passe de 5mm à 0.25mm dans la zone d'impact en veillant à conserver un élément dans l'épaisseur des plis de toute la structure. La taille du problème est ainsi de l'ordre du million de DDLs avec 90% des noeuds environ dans la partie implicite. La taille de l'interface est quant à elle de l'ordre de 3000 DDLs. Un impact de 10J à été choisi pour les simulations présentées par la suite. Le pas de temps numérique de la partie impactée est globalement contraint par la condition de stabilité propre à l'algorithme explicite et nous le fixons à $2.10^{-8}$s pour cette étude. Nous faisons ainsi varier le pas de temps de la partie complémentaire traitée par Zébulon en fonction du rapport de pas de temps.

Un modèle matériau endommageable simplifié par rapport au modèle Onera Progressive Failure Model [11] a été utilisé. Il s'agit d'un modèle endommageable à effet retard identifié sur du T700/M21 dont le comportement est illustré sur la figure 5 pour un chargement de traction compression piloté en déformation. Un des résultats obtenu avec ce modèle pour un rapport de pas de temps de 500 est présenté sur la figure 6. On observe ainsi qualitativement un cône d'endommagement classiquement observé dans ce type d'impact situé au niveau de la zone d'impact.

Nous présentons sur la figure 7 l'évolution temporelle de la flèche du panneau composite pour deux rapports de pas temps (500 et 5000). On observe des allures globales relativement similaires. Néanmoins on remarque sur la courbe d'erreur quelques différences situées symétriquement approximativement autour des instants 0.0036 et 0.0052 s. Ces écarts peuvent s'expliquer par le grand rapport de pas de temps choisi. En effet, il conduit à réaliser uniquement 60 incréments sur la totalité du temps de simulation dans la partie implicite (Zébulon) ce qui est très grossier. Or à ces instants on observe des phénomènes de rebonds de l'impacteur sur la plaque qui sont donc mal captés pour m=5000. Malgré cela on voit sur la figure 7 que l'erreur reste relativement faible bien que l'on utilise des grands rapport de pas de temps. Cette erreur est définie par la valeur absolue de la différence entre la solution considérée comme référence (ici celle avec un rapport m de 500) et la solution courante normée par le déplacement maximum.
<figure><img src="image_6.png" /><figcaption>Fig. 5 – Comportements dans les sens fibre et matrice obtenus avec le modèle matériau élastique endomageable utilisé pour une charge décharge pilotée en déformation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Endommagement matriciel après impacte de 10J sur la panneau composite raidi avec l’algorithme de couplage mis en œuvre pour un rapport de pas de temps de 500.</figcaption></figure>

# 5 Conclusion et perspectives

Nous avons donc mis en œuvre la méthode GC au travers d'un couplage de codes de type structure/structure couplant ainsi les codes de calcul implicite Zébulon et explicite Europlexus. Ceci a permis de dégager l'interface minimale qui doit être publiée par le codes de calcul pour mettre en place cette méthode. Nous avons également mis en avant le potentiel de cette méthode sur un calcul de grande taille.
<figure><img src="image_8.png" /><figcaption>Fig. 7 – Evolution temporelle de la flèche du panneau composite raidi impacté à 10J pour différents rapports de pas de temps.</figcaption></figure>

On présente sur la figure 8 la répartition des temps de calcul obtenus sur ce panneau pour différents rapports de pas de temps (50, 500 et 5000) avec un modèle élastique et une décomposition de domaine légèrement différente de celle présentée.

<figure><img src="image_9.png" /><figcaption>Fig. 8 – Répartition des temps de calcul pour les différents rapports de pas de temps utilisés sur le panneau composite élastique linéaire impacté à 10J.</figcaption></figure>

On observe premièrement que le choix de résoudre le problème d'interface sur l'échelle de temps fine n'est pas trop pénalisant du point de vu du temps de calcul pour cette taille d'interface (environ 3000 DDLs). On constate également que le temps de calcul global diminue lorsque l'on augmente le rapport de pas de temps. Cela montre ainsi que dans cette situation nous avons réussis à faire focaliser l'effort de résolution dans la partie explicite. En effet, pour un rapport de pas de temps de 500 et de 5000, les temps passés dans la partie implicite deviennent très faibles (mois de 3h pour m=500 et de 1h pour m=5000).
Cependant, ce gain de temps tend vers une limite qui correspond au temps passé dans la zone explicite. Des travaux sont donc en cours pour améliorer l'équilibrage des charges de manière à réduire le temps de calcul passé dans la partie impactée. Nous avons en particulier adapté le couplage développé pour l'étendre à une décomposition mettant en œuvre plus de deux sous-domaines. Les calculs réalisés sur le panneau composite ont été réalisés avec ces développements. Nous avons décomposé la partie explicite en plusieurs sous-domaines illustrés sur la figure 4 par la présence ou non du maillage (la partie explicite comporte ainsi 9 sous-domaines en forme d'anneau autour du point d'impact). Ils ont cependant tous le même pas de temps (m=1) ce qui, dans ce cas, reviens à réaliser une parallélisation mono-échelle en temps classique du code explicite Europlexus. Cette parallélisation est cependant réalisée de l'ex-térieur par l'intermédiaire de l'algorithme développé. Nous continuons également à enrichir le modèle de la partie impactée (éléments de zone cohésive, modèle OPFM) dans le but de réaliser des compa- raisons expérimentales/numériques dans le cadre du Projet de Recherche Fédérateur intitulé "transition statique/dynamique pour les matériaux composites" en cours à l'Onera.

**Remerciements** : Les auteurs de cet article tiennent à remercier Vincent Faucher du CEA ainsi que Steve Belon et Roland Ortiz de l'Onera (CL/DADS/CRD) pour leur aide relative à Europlexus.

Références

[1] A. Combescure, A. Gravouil. *A time-space multi-scale algorithm for transient structural nonlinear problems*, Mécanique & Industries, Elsevier, 2(1), p. 43-55, 2001.

[2] A. Gravouil, A. Combescure. *Multi-time-step explicit-implicit method for non-linear structural dynamics*, International Journal for Numerical Methods in Engineering, Wiley Online Library, 38(1), p. 199-225, 2001.

[3] C. Lopes. *Damage and Failure of non-conventional composite laminates*, PhD thesis, TU Delft, Faculteit Luchtvaart-en Ruimtevaarttechniek, 2009.

[4] C. Bouvet. *Étude de l'endommagement dans les structures composites*, habilitation à diriger des recherches, Université de Toulouse III, Paul Sabatier, 2009.

[5] C. Dupleix-couderc. *Sur une stratégie multi-échelle d'analyse des grands délaminages en dynamique transitoire*, Thèse, LMT-Cachan, 2011.

[6] L. Gendre, O. Allix, P. Gosselet, F. Comte. *Non-intrusive and exact global/local techniques for structural problems with local plasticity*, Computational Mechanics, Springer, 44(2), p. 233-245, 2009.

[7] A.Prakash, K.D. Hjelmstad, *A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics*, International journal for numerical methods in engineering, Wiley Online Library,61(13), p. 2183-2204, 2004.

[8] N. Mahjoubi, A. Gravouil, A. Combescure, N. Greffet. *A monolithic energy conserving method to couple heterogeneous time integrators with incompatible time steps in structural dynamics*, Computer Methods in Applied Mechanics and Engineering, Elsevier, 2010.

[9] A. Ghanem, M. Torkhani, N. Mahjoubi, T.N. Baranger, A. Combescure. *Arlequin framework for multi-model, multi-time scale and heterogeneous time integrators for structural transient dynamics*, Computer Methods in Applied Mechanics and Engineering, Elsevier, 2012.

[10] C. Farhat, L. Crivelli, F.X. Roux. *A transient FETI methodology for large-scale parallel implicit computations in structural mechanics*. International Journal for Numerical Methods in Engineering, Wiley Online Library, 37(11), p. 1945-1975, 1994.

[11] F. Laurin, N. Carrere, J.-F. Maire. *A multiscale progressive failure approach for composite laminates based on thermodynamical viscoelastic and damage models*, Composites Part A, Applied Science and Manufacturing, Elsevier, 38(1), p. 198-209, 2007.