# Modélisation SPH 3D de l'impact basse vitesse sur plaque composite

Muhammad Ilyas, Jérôme Limido, Frederic Lachaud, Christine Espinosa,
Michel Salaün

► To cite this version:

Muhammad Ilyas, Jérôme Limido, Frederic Lachaud, Christine Espinosa, Michel Salaün. Modélisation SPH 3D de l'impact basse vitesse sur plaque composite. CFM 2009 - 19ème Congrès Français de Mécanique, Aug 2009, Marseille, France. pp.1-6. hal-01852247

HAL Id: hal-01852247

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

Submitted on 1 Aug 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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

This is an author-deposited version published in: http://oatao.univ-toulouse.fr/
Eprints ID: 3538

**To cite this document** : ILYAS, Muhammad, LACHAUD, Frédéric, ESPINOSA, Christine, MICHEL, Laurent, SALAÜN, Michel. Modélisation en dynamique rapide du délaminage des composites à l'aide d'éléments finis cohésifs. In: *19ème Congrès Français de Mécanique, Marseille, 24-28 Août 2009*

Any correspondence concerning this service should be sent to the repository
administrator: staff-oatao@inp-toulouse.fr
# Modélisation SPH 3D de l'impact basse vitesse sur plaque composite

M. ILYASª, J. LIMIDOᵇ, F. LACHAUDª, C. ESPINOSAª, M. SALAÜNª

a Université de Toulouse - ISAE, DMSM, 10 Avenue Edouard Belin, 31054 TOULOUSE Cedex 4, France

b IMPETUS Afea, Rue du Lanoux, 31330 Grenade Sur Adour, France

## Résumé:

L'objet de cette étude est d'évaluer l'utilisation d'une méthode Smoothed Particle Hydrodynamics (SPH) « stabilisée » pour modéliser le comportement macroscopique, ainsi que le comportement local mésoscopique, de plaques composites typiques de coupons aéronautiques, sollicitées par un impact basse vitesse. Les résultats présentés concernent la comparaison des courbes force-déplacement, et l'analyse couche par couche du comportement de la plaque.

## Abstract:

The present study aims to assess the ability of a stabilized Smoothed Particle Hydrodynamics method (SPH) to model the macroscopic and local mesoscopic behaviour of composites plates, typical of aeronautical parts, subjected to a low velocity impact. Compared results are the global force-displacement curve and the local ply by ply behaviour of the plate.

**Mots-clefs :** SPH, impact, plaque ; composite stratifié, comparaison essais-calculs

## 1 Introduction

L'emploi des matériaux composites dans les zones travaillantes des structures aéronautiques amène les concepteurs à devoir maîtriser les phénomènes d'initiation et de propagation des dégâts ayant une incidence sur leur tenue résiduelle, et donc sur leur certification. Des campagnes d'expérimentation essais/calculs sont réalisées dans ce but, qui utilisent principalement des modèles numériques basés sur les éléments finis (EF) pour représenter les modèles mécaniques de stratifiés composites. Les modèles de comportement intègrent les phénomènes multiples à l'origine des dégâts ou responsables de leur propagation. La méthode des éléments finis est parfois couplée avec des méthodes de liaisons cassables ou d'enrichissement pour représenter les dégâts : éléments cohésifs, contacts, liaisons nodales, XFEM, ...

Partant du constat que les délaminages et les fissures sont les ruptures significatives qui pilotent la résistance résiduelle des structures composites minces utilisées en aéronautique et que la méthode SPH a prouvé sa capacité à modéliser les ruptures [1], nous nous intéressons à l'apport de la méthode SPH dans la modélisation des stratifiés composites sous impact. Nous avons choisi dans un premier temps, d'améliorer la méthode SPH pour la rendre compatible avec : l'anisotropie des matériaux, la stratification, le comportement basse vitesse, un maillage EF autour du point d'impact. Les enjeux de méthode sont respectivement : l'introduction de l'anisotropie en général, la stabilisation de la méthode pour des basses déformations et vitesses de déformation, le couplage SPH/EF. Les travaux présentés sont le résultat d'essais menés au Département Mécanique des Structures et Matériaux de l'ISAE, de simulations EF et SPH menées au Département, de développements et de calculs SPH réalisés par la société IMPETUS Afea France pour le compte du LSTC propriétaire du code de calcul LS-DYNA. Une étude portant spécifiquement sur la comparaison essai calcul avec un modèle tout EF intégrant les dommages est présentée dans une autre publication [2], ainsi qu'une explication détaillée des développements réalisés pour stabiliser la méthode et la rendre compatible avec les éléments finis [3].

Nous avons choisi de présenter ici les résultats d'une étude de faisabilité qui montre la cohérence du comportement en flexion d'un modèle de plaque stratifiée composite tout SPH avec un modèle tout EF. Nous comparons ensuite des résultats d'essais d'impacts à 4J, avec ceux d'un modèle élastique tout éléments finis, et ceux d'un modèle élastique EF/SPH.
# 2 Etude de faisabilité

## 2.1 Formulation lagrangienne totale

La méthode SPH est une méthode numérique discrétisant le domaine par des particules qui ont des caractéristiques matérielles individuelles et se déplacent en interaction et interdépendance à l'intérieur de voisinages appelés sphères d'influence. L'état de la matière est décrit par la vitesse des particules, leur densité, les contraintes et flux d'efforts ainsi que la conservation de l'énergie dans la sphère d'influence. Cette méthode a évolué depuis sa création [4], notamment pour mieux s'appliquer aux cas de grandes déformations dynamiques [5]. Les méthodes SPH « standards » sont basées sur un noyau eulérien [6]. Le support d'intégration reste dans ce cas constant et les particules voisines dans la sphère d'influence sont déterminées à chaque cycle de calcul. Cette technique permet de modéliser les grandes déformations subies par la matière, mais est connue pour faire apparaître des instabilités en tension qui limitent souvent le champ d'application de la méthode aux grandes déformations [7]. Le noyau peut être choisi lagrangien. Dans ce cas, les particules voisines restent les mêmes au cours du calcul ce qui permet de gagner du temps de calcul, et le support du noyau évolue. Cette approche permet d'éliminer la majorité des instabilités en tension [8]. En contrepartie cette méthode n'est pas adaptée au traitement des grandes déformations [7]. La méthode que nous avons choisie ici est une formulation lagrangienne totale, cohérente avec les déformations constatées dans les essais réels d'impacts basse énergie sur plaques composites, et cohérente avec l'hypothèse simplificatrice du comportement élastique qui sera utilisée dans les modèles numériques. Pour plus de détails, le lecteur pourra consulter les références.

## 2.2 Prise en compte de l'anisotropie

Les études connues ayant utilisé la méthode SPH pour modéliser des composites concernent surtout des impacts énergétiques pour lesquels l'anisotropie locale n'est pas significative sur les résultats analysés, à savoir les limites d'écaillage et de perforation ou des seuils de perte totale de tenue [9]. Une étude de faisabilité réalisée il y a quelques années nous a permis de montrer la capacité de la méthode SPH à reproduire qualitativement les dégâts en fonction de la vitesse et de l'angle d'impact dans des plaques stratifiées 8 plis, en impact hypervéloce. Cette étude a permis de s'assurer de la stabilité de la méthode avec noyau eulérien dans les cas d'impacts anisotropes, et de trouver une méthode de modélisation des délaminages. La taille des trous était approximativement respectée, compte tenu des informations données par les auteurs des essais qui ont servi de base de comparaison, et des hypothèses simplificatrices choisies pour la modélisation 2D (figure 1) [10, 11]. Dans cette étude, nous sommes dans le cas où l'anisotropie des plis est prépondérante pendant toute la phase de charge par le projectile (initiation des dommages) puis lors de la déformation de la plaque suivant les modes de structure (propagation et extension des zones endommagées). Il est donc nécessaire de prendre en compte l'anisotropie des plis et de la stratification. L'anisotropie des plis est prise en compte via le modèle de comportement, ici un modèle 3D élastique orthotrope pour chaque pli. L'endommagement, possible avec ce modèle de matériaux du code LS-DYNA, n'a pas été utilisé dans nos modèles. Chaque pli est modélisé par une rangée de particules SPH. Les couches sont empilées en suivant la séquence d'empilement suivant le cas étudié. La repère d'anisotropie est suivi dans chaque particule au cours du temps. Les données matériaux sont résumées dans la table 1.

<figure><img src="image_3.png" /><figcaption>FIG. 1 - Répartition des contraintes axiales pendant l'impact, délaminage des couches externes pendant l'impact, comparaison des formes de trou dans l'épaisseur des plaques après impact</figcaption></figure>

<table><caption>Table 1 : caractéristiques mécaniques choisies pour les plis unidirectionnels T700/M21.</caption><thead><tr><th>Masse volumique</th><th colspan="5">ρ 1550 kg/m<sup>3</sup></th></tr></thead><tbody><tr><td>Modules d'Young</td><td>E11</td><td>135 GPa</td><td>E22</td><td>8.5 GPa</td><td>E33</td><td>8.5 GPa</td></tr><tr><td>Modules de cisaillement</td><td>G12</td><td>4.8 GPa</td><td>G23</td><td>3.2 GPa</td><td>G31</td><td>4.8 GPa</td></tr><tr><td>Coefficients de Poisson</td><td>v21</td><td>0.019</td><td>v31</td><td>0.019</td><td>v32</td><td>0.4</td></tr></tbody></table>
## 2.3 Flexion d'un stratifié : comparaison SPH / EF

La stabilité de la méthode a été testée pour un chargement de flexion d'une poutre console en grand déplacement. La console est une plaque de dimensions 10mmx32mmx2mm encastrée à une extrémité et soumise à une force verticale à l'autre. La structure contient 8 plis T700/M21 à 0°. Elle est modélisée par huit couches d'éléments solides à 8 nœuds (3 degrés de libertés par nœud) et 1 point d'intégration dans le modèle EF et par huit couches de particules dans le modèle SPH. Les contraintes axiales sont comparées avec un modèle tout éléments finis. Les grands déplacements subis par la plaque SPH ne provoquent pas de rupture numérique. La méthode reste stable. Les amplitudes de contraintes axiales maximales sont identiques et leur répartition est similaire entre les deux modèles, dans le plan comme dans l'épaisseur. La plaque SPH fait apparaître un effet de Poisson légèrement plus marqué, ce qui est cohérent avec la rigidité supérieure des éléments finis utilisés pour le modèle EF en l'absence de rotations nodales. Les déformées sont quasiment identiques. Ce modèle valide à la fois la stabilité de la méthode SPH, et sa capacité à modéliser un comportement orthotrope, même en grand déplacement.

<figure><img src="image_4.png" /><figcaption>FIG. 2 - Répartition des contraintes axiales dans [-3GPa ; 5GPa] dans les modèles EF et SPH</figcaption></figure>

# 3 Application au cas d'impact basse vitesse

## 3.1 Les essais

Les essais d'impact sont réalisés au DMSM, grâce à une machine à poids tombant (figure 4).

<figure><img src="image_5.png" /><figcaption>FIG. 3 - Vues du système d'impact à poids tombant.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4 - micrographie de la coupe longitudinale et radiographie vue de dessus d'une plaque analysée</figcaption></figure>

Les cibles sont des plaques rectangulaires 100mmx150mmx2mm stratifiées [0/90]$_{2s}$ à base de plis unidirectionnels T700/M21, simplement posées sur un bâti ayant une fenêtre de 125mm par 75mm. Le projectile est un cylindre à embout hémisphérique de masse constante 1,369 kg, de rayon 8mm. Plusieurs vitesses d'impact ont été testées. L'essai que nous avons choisi d'illustrer est un essai à 4J. La fréquence
d'acquisition est de 1MHz pour les courbes historiques. Les plaques ont été analysées par des mesures destructives et non destructives. Les figures ci-dessous présentent ces deux analyses après un impact de 4J. Des délaminages entre les plis d'orientation différentes et des fissurations transverses dans les plis à 90° ont été relevés dans les plaques. Vue de dessus, la surface totale délaminée autour de l'axe d'impact est de forme ovale de dimensions 6mm par 15mm. Les dégâts sont peu importants et peu étendus. Les courbes d'efforts et de déplacement au cours du temps sont comparées ci-dessous aux résultats du modèle éléments finis.

## 3.2 Modèles EF et EF/SPH

Le modèle éléments finis a fait l'objet de calculs de mise au point et de recalages dans le but de réduire à minima l'utilisation de coefficients numériques purs dont le réglage peut notoirement modifier les résultats, en particulier force/déplacement. Ces travaux sont décrites dans une autre publication [2]. Le modèle éléments finis précédent sert de base au modèle qui contient les particules SPH. Le centre de la plaque est remplacé par des particules qui sont cinématiquement contraintes aux éléments finis à la frontière des deux discrétisations. La zone centrale fait une dimension de 32mmx32mm et contient environ 120000 mille particules. Une seule couche de particules est utilisée pour chaque pli composite. Les couches sont empilées en suivant la stratification de la plaque d'essai (figure 5).

<figure><img src="image_7.png" /><figcaption>FIG. 5 - Modèles de plaque composite stratifiée EF et EF/SPH</figcaption></figure>

La plaque et le bâti sont modélisés par des éléments finis solides à 8 nœuds (3 ddls par nœud) et 1 point d'intégration. Le maillage a été raffiné dans la zone proche de l'impact afin de décrire plus finement les contraintes et les déplacements. Aux éléments de chacune des huit couches sont affectées les propriétés d'un pli unidirectionnel élastique orthotrope, avec un repérage de la direction des fibres. Les propriétés sont décrites en table 1. Un total de 30500 éléments solides composent la plaque et le bâti du modèle tout EF qui est supposé déformable élastique, avec les propriétés matériau d'un acier standard : $\rho$= 7850 kg/m3, $E$=210GPa, $v$=0,3. Dans les modèles présentés, le projectile est modélisé par une coque sphérique indéformable animée d'une vitesse verticale initiale de 2,42m/s, et de masse 1369g. Des algorithmes de contact en pénalité sans frottement sont introduits entre le projectile et la plaque, et entre la plaque et le bâti.

## 3.3 Comparaison essais / modèles

L'endommagement n'étant pas pris en compte dans le modèle, les résultats comparés sont les courbes temporelles du déplacement au centre de la plaque en face arrière, et d'effort entre le projectile et la plaque, et les répartitions des contraintes significatives.

Sur les courbes d'effort, le premier pic et le plateau sont correctement restitués par le modèle, ainsi que les fréquences d'oscillation jusqu'au pic. La durée d'application de l'effort de contact est également la même entre l'essai et la simulation. Les pics de déplacement face arrière et de force sont obtenus au même instant dans l'essai et les calculs, soit à environ 3 ms. Les valeurs maximales sont résumées dans le tableau 2. Les écarts entre les courbes en déplacement après le pic, et entre les courbes en force au pic, s'expliquent par le fait que les modèles ne sont pas capables de décrire l'endommagement, et sont donc plus raides. Mis à part la raideur élastique un peu supérieure à la réalité dans les deux calculs, les modèles sont relativement fidèles au comportement macroscopique de la plaque réelle faiblement endommagée.

<table><caption>Table 2 : valeurs de déplacements et forces maximales au pic (3 ms).</caption><thead><tr><td></td><td>Déplacement max (mm)</td><td>Force max (N)</td></tr></thead><tbody><tr><td>Essai</td><td>4,58 à 2,9 ms (réf.)</td><td>1900 (réf.)</td></tr><tr><td>Simulation EF seul</td><td>4,62 à 3 ms (+0.8%)</td><td>2350 (+23%)</td></tr><tr><td>Simulation EF / SPH</td><td>4,41 à 2,8 ms (-3.1%)</td><td>2250 (+18%)</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>FIG. 6 - courbes déplacement - temps et force - temps au centre de la plaque en face arrière (4J)</figcaption></figure>

<table><caption>Table 3 : contraintes σx, σy, σyz et σzx au pic (3 ms), dans les plis à 0° (externe face arrière) et 90°.</caption><thead><tr><th></th><th>SIGMA sens x<br>(MPa)</th><th>SIGMA sens y<br>(MPa)</th><th>SIGMA yz<br>(MPa)</th><th>SIGMA zx<br>(MPa)</th></tr></thead><tbody><tr><td>Dernier pli à 0° tout EF</td><td>-405 à 1900</td><td>-17 à 240</td><td>-23 à 23</td><td>-114 à 114</td></tr><tr><td>Dernier pli à 0° EF/SPH</td><td>-395 à 1950</td><td>-16 à 245</td><td>-16 à 16</td><td>-105 à 105</td></tr><tr><td>Dernier pli à 90° tout EF</td><td>-22 à 91</td><td>-230 à 2050</td><td>-160 à 160</td><td>-96 à 96</td></tr><tr><td>Dernier pli à 90° EF/SPH</td><td>-22 à 97</td><td>-220 à 2060</td><td>-173 à 173</td><td>-118 à 118</td></tr></tbody></table>

Les contraintes les plus élevées dans le pli à 0° face arrière et dans le dernier pli à 90° (le plus proche de la face arrière) sont résumées dans le tableau 3. Dans les deux modèles, l'étendue de chaque ellipse du délaminage observé en figure 6 entre les deux plis inférieurs (5mm de long) correspond à une zone du pli à 0° inférieur où se superposent des contraintes de traction sens x et de cisaillement transverse de 40 à 95 MPa, et des contraintes de traction sens y de l'ordre de 120MPa, et à une zone du pli à 90° où se superposent des contraintes sens x et y respectivement de l'ordre de 40 à 95 MPa et 115 MPa (figures 7et 8). L'enveloppe de niveaux des contraintes de cisaillement dans le pli à 0° est orientée dans le sens du délaminage observé entre les deux derniers plis dans les essais, et a la même forme ellipsoïdale. L'espacement entre les ellipses symétriques dans le calcul est également de l'ordre de grandeur de l'espace entre les ellipses délaminées dans l'essai (figure 4).

<figure><img src="image_9.png" /><figcaption>FIG. 7 - contraintes σx EF et EF/SPH [-400;1950] MPa (haut) et σy EF et EF/SPH [-16 ;240] MPa (bas).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8 - contraintes σzx et σyz [-120MPa ; 120MPa] dans le repère matériau du modèle EF, σzx dans le repère global en SPH: lieu et étendue des délaminages et des fissures transverses</figcaption></figure>

## 4. Conclusion

Les travaux présentés ont démontré la capacité de la méthode SPH à modéliser le comportement élastique orthotrope de plaques composites stratifiées ; la méthode est stable et le couplage EFS/PH est continu. L'étude de faisabilité a permis de comparer un modèle SPH de flexion et un modèle éléments finis. Le modèle éléments finis présente une raideur plus importante tandis que le modèle SPH permet de prendre en compte les effets de Poisson. Les valeurs de déplacements et contraintes sont quasi identiques. Dans le cas d'application d'un impact à 4J, nous avons comparé les résultats de mesures d'essais et les résultats de calculs tout FE et couplé EF/SPH. Les pics en forces et déplacements sont obtenus au même instant dans les essais et les modèles. Les déplacements maximaux en face arrière sont quasi identiques. Les valeurs maximales de forces restent trop élevées avec les deux modèles élastiques par rapport à l'essai. Les valeurs maximales de contraintes dans les plis sont tout à fait similaires entre les éléments finis et le modèle couplé EF/SPH, et l'analyse des contraintes permet de retrouver qualitativement l'origine des fissurations transverses observées dans les plis à 90°, et les délaminages. Une comparaison quantitative a été menée entre les valeurs des contraintes significatives et l'étendue de la zone délaminée entre les plis de la face arrière. Les résultats sont prometteurs. L'analyse des contraintes dans les différentes couches permet d'envisager d'étendre le modèle SPH à la prévision du délaminage et de l'endommagement dans les plis, grâce à des modèles de comportement de composite intégrant l'endommagement en particulier.

## Remerciements

Nous remercions de leur appui, Michel Labarrère du DMSM et Jean Luc Lacome de la société LSTC.

## Références bibliographiques

[1] J. Limido, Etude de l'effet de l'usinage grande vitesse sur la tenue en fatigue de pièces aéronautiques, Thèse de l'Université Paul Sabatier, Toulouse, 2008

[2] M. Ilyas, F. Lachaud, Ch; Espinosa, M. Salaün, Simulation of dynamic delamination and mode I energy dissipation, 7th European Conference LS-DYNA, Salzburg, 2009

[3] J.-L. Lacome, J. Limido, Ch; Espinosa, SPH Formulation with Lagrangian Eulerian adaptive kernel, 4th SPHERIC Workshop, Nantes, 2009

[4] L.B. Lucy, A numerical approach to testing of the fission hypothesis, Astron. J. 82, 1977

[5] J.J. Monaghan, SPH without a Tensile Instability, Journal of Computational Physics 159, 290-311, 2002

[6] G.R. Liu, M.B. Liu. Smoothed Particle Hydrodynamics – A Meshfree Particle Method. World Scientific Publishing Co., ISBN 981-238-456-1, 2003

[7] S. P. Xiao, T. Belytschko. Material stability analysis of particle methods. Adv. Comput. Math., 171-190, 2005

[8] A. Combescure, M. Brunet, B. Maurel, G. Nefussi, F. Abed Merail. Les simulations numériques de la mise en forme sont elles fiables? Congrès Français Mecamat, Aussois, 2008

[9] R. Salome, V. albouys, Ch. Le Floch, D. Sornette, J.-P. Vila, High pressure composite tank behaviour under an hypervelocity impact, 3rd European Conference on Space Debris, Darmstadt, 2001

[10] V.V. Silvestrov, A.V. Plastinin, N.N. Gorshkov, Hypervelocity impact on laminate composite panels, International Journal of Impact Engineering, Vol. 17, 751-762, 1995

[11] E. Jacquet, Etude numérique de l'endommagement induit par l'impact de debris spatiaux sur des structures composites, DEA de l'ENSICA, 2004