# Développement en dynamique d'un élément fini multicouche avec interfaces imparfaites

van Anh Duong, Jean-François Caron, Gilles Forêt, Dang Truc Nguyen

► To cite this version:

van Anh Duong, Jean-François Caron, Gilles Forêt, Dang Truc Nguyen. Développement en dynamique d'un élément fini multicouche avec interfaces imparfaites. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422162

HAL Id: hal-01422162

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

Submitted on 23 Dec 2016

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
# Développement en dynamique d'un élément fini multicouche avec interfaces imparfaites

V.A. Duong, J.F. Caron, G. Foret, D.T.Nguyen

Paris Est, Institut Navier ENPC-LCPC

Ecole Nationale des Ponts et Chaussées (ENPC) - 77455 MARNE LA VALLEE CEDEX 2

caron@enpc.fr

**Résumé** — Un élément fini de type layer-wise, une cinématique par couche, est proposé ici. Ce modèle comporte 5n champs cinématiques pour un multicouche à n couches. Cet article a pour but de présenter les récents développement de MPFEAP qui propose désormais, un module dynamique qui permet de calculer des modes propres et d'aborder des problèmes d'impact. Il permet également désormais la prise en compte d'interface non linéaires ou imparfaites.

**Mots clés** — mot éléments finis, composites, impact, interface.

## 1 Introduction

L'optimisation fine des structures composites nécessite de mettre en place des outils de modélisation du comportement mécanique de plus en plus sophistiqués et prenant en compte les spécificités de ces matériaux-structure. Les points clés à maîtriser sont essentiellement les phénomènes d'interface, tels que délaminages, glissements ou autres endommagements. De nombreux travaux s'intéressent à ces aspects, introduction de la mécanique de l'endommagement dans les interfaces (ou dans les couches) [1] [2], calcul raffiné des champs de contraintes singuliers [3], comportements non linéaires des interfaces par exemple [4]. Ici un élément fini à une cinématique par couche, est proposé afin de prévoir la mécanique du multicouche ayant des interfaces "imparfaites". Il est basé sur les travaux de l'équipe multicouche de l'institut Navier [5][6] [8] qui s'inspirent du modèle de Pagano [7]. Ce modèle comporte 5n champs cinématiques pour un multicouche à n couches et porte le nom "M4-5N". Il approche chaque couche par une plaque de Reissner et intègre des efforts d'interface généralisés. Il rejoint ainsi la famille des modèles layerwise ([9],[10]...). La première version du code d'élément fini appelé MPFEAP (MultiParticle Finite Element Analysis Program) a été implémentée [11]. Cet article a pour but de présenter les récents développement de MPFEAP qui propose désormais, un module dynamique et une prise en compte d'interface non linéaires ou imparfaites.

## 2 Modes propres, influence du degré d'anisotropie.

Une méthode classique de sous-espace est utilisée pour le calcul des modes propres. Pour démontrer l'efficacité du modèle et du code de calcul, des exemples qui comprennent des plaques isotropes, des plaques sandwichs sur appuis simples ont été calculés (en particulier $(0/90)_n$, $(45/ - 45/45/ - 45)$.... Les résultats ont été comparés à des solutions analytiques et à d'autres codes d'éléments finis. Les résultats sont très satisfaisants. On ne cite ici que les résultats concernant les plaques multicouches orthotropes $(0/90)_n$. Le degré d'orthotropie est définit par le rapport des
modules d'Young suivant deux directions du matériau $E_1/E_2$. Le nombre de couches varie de 2 à 6. La première fréquence propre non-dimensionnelle obtenue par MPFEAP est validée par comparaison avec celle de la théorie 3D [12] et le code élément fini donné par Khare-Kant [13] et Putcha, Reddy [14], les solutions de la théorie 3D sont utilisées comme références. Les effets du facteur $E_1/E_2$ sont illustrés dans la fig. 1. Le pourcentage d'erreur de chaque modèle se trouve dans un groupe indépendant et sur un axe logarithmique. A travers cet exemple, on voit ainsi l'efficacité et la stabilité du modèle $M4-5n$ dans le calcul des valeurs propres des plaques dont le degré d'orthotropie est important. L'erreur qui oscille autour de 1% est acceptable, et surtout non influencée par le degré d'anisotropie.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Pourcentage d'erreur de la première fréquence propre non-dimensionnelle des différents modèles (de gauche à droite, MPFEAP - [13]- [13]raffiné- [14], pour une plaque asymétrique (0/90)$_n$ ($n = 2$ ici)</figcaption></figure>

### 3 Impacts de plaques composites.

MPFEAP s'est également vu enrichi par la prise en compte possible d'impacts gérés par la loi de Hertz ou celle de Yang et Sun (calcul implicite type Newmark).

Pour valider le code de calcul, une plaque isotrope impactée (loi de Hertz) est comparée à la solution analytique de Karas [15] au niveau de la force de contact, du déplacement de la plaque et de l'impacteur et de la vitesse de la plaque et de l'impacteur. Les résultats comparés sont très proches. Pour la même plaque, l'étude pour différents maillages et différents pas de temps montre une convergence très rapide. Des résultats concernant des plaques multicouches avec 3 conditions aux limites différentes et 3 vitesses différentes d'impact (loi de Yang et Sun) sont très proches des solutions numériques de Choi [16]. L'élément fini proposé semble donc une solution assez robuste et peu couteuse (2D) pour l'étude des impacts ou chocs sur les composites. Son originalité provient du fait que compte tenu des grandeurs (déplacements, déformations, cisaillements et arrachements) d'interface présentes dans le modèle, il permet l'étude des interfaces lors d'impact, ce qui est compte tenu de la nature des endomagements des multicouches est potentiellement très intéressant.

### 4 Interfaces imparfaites.

Le modèle $M4-5n$ est donc spécifiquement dévolu à l'étude des phénomènes d'interfaces. L'introduction des efforts et déformations d'interface en fait un outil puissant pour l'étude des délaminages ou autres glissements. Une extension du modèle original (considérant l'interface comme parfaite) a déjà été réalisée par Diaz [17][18]. Une discontinuité a été introduite dans la
formulation qui permet une non linéarité au niveau de l'interface, jeu, glissement élastique, plastique parfaite ou autre. Dans le développement présenté, on propose l'introduction de ces aspects dans MPFEAP en se concentrant sur deux problématiques particulières, celle des jeux élastiques introduits par des connecteurs dans une poutre mixte par exemple et celle du collage avec une interface particulière élastoplastique.

## 4.1 Modèle de connecteur élastique

L'interface frontière est ici représenté par trois ressorts élastiques indépendants entre eux dans les trois directions principales de l'interface. La rigidité de chaque ressort représente la relation linéaire entre le glissement de l'interface et la contrainte suivant la direction correspondante.Il n'existe pas de couplage entre les directions. Ce modèle est applicable dans les calculs des poutres mixtes (en bois-béton ou en béton-acier) où différents types de connecteurs sont utilisés qui présentent tous des souplesses à prendre en compte (des recommandations sont faites en ce sens dans les eurocodes).

Les efforts à l'interface se composent de trois composantes : $\sigma_{13}$, $\sigma_{23}$, $\sigma_{33}$ correspondant respectivement à deux contraintes de cisaillement et la contrainte normale. La relation entre les efforts et les glissements d'interface élastique $\delta_1$, $\delta_2$, $\delta_3$ s'écrit sous la forme :

$$ \delta_{\alpha} = S_{\alpha\beta}\sigma_{\beta\gamma} \qquad \delta_{\beta} = S_{\beta\gamma}\sigma_{\gamma\alpha} \qquad \text{avec } \alpha, \beta = 1, 2 \qquad (1) $$

Ce modèle est simple à introduire dans les codes d'éléments finis (généralement sous la dénomination de "modèles cohésifs") et peut également fournir des informations importantes sur l'interface. En fonction de son besoin, chaque auteur choisit une forme de S à implémenter dans son modèle. Quelques auteurs ([19]...), acceptent une rigidité de l'interface calculée à travers les propriétés des couches, d'autres [18] calculent la rigidité de l'interface à partir des propriétés de l'interface, d'autres enfin [20] utilisent des coefficients expérimentaux. Ce sera le cas ici puisqu'il s'agit de caractériser des connecteurs. Dans tous les modèles, l'interface n'a pas d'épaisseur.

On considère une poutre mixte bois-béton sur appuis simples (fig. 2). La couche supérieure est en béton et celle inférieure en bois. Une charge uniforme est imposée à la surface supérieure de la poutre.

<figure><img src="image_3.png" /><figcaption>Figure 2 – La poutre mixte en bois-béton considérée et son maillage 2D</figcaption></figure>

On choisit une interface assez faible pour bien mettre en valeur le glissement. La rigidité choisie de l'interface correspond à celle identifiée par le modèle analytique pour un connecteur métallique discret classiquement utilisé.

Les contraintes de cisaillement $\tau_x$ présentées dans la fig. 3 sont en accord avec les résultats de la solution analytique. La figure montre également une bonne convergence du résultat avec le raffinement du maillage. Même un maillage 2D très grossier donne le bon profil du cisaillement à l'interface, mais pour s'approcher des valeurs de contrainte au bord, un raffinement peut être utile.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Les contraintes généralisées de cisaillement à l’interface de la poutre mixte considérée</figcaption></figure>

## 4.2 Modèle de couche interfaciale élastoplastique parfaite

Le deuxième développement interface réalisé a pour but de modéliser la couche de colle ou résine existante dans les interfaces de composites et dans les joints. L'interface considérée est élasto-plastique parfaite obéissant au critère de Von-Mises. Contrairement à ce qui a été fait au paravant dans [17] où la colle était une simple frontière entre les couches, ici l'interface-colle a une épaisseur physique, afin de proposer son propre comportement mais reste une interface au sens où elle ne travaille qu'hors plan, cisaillements transverses et contrainte normale. L'algorithme de projection est utilisé pour résoudre ce problème. Dans cette section, un modèle 3D élasto-plastique de la colle obéissant à un critère Von-Mises est adapté et implémenté dans le modèle M4-5n. Ce modèle suppose que l'interface-colle ne subit que des contraintes de cisaillement $τ_x$, $τ_y$ et d'arrachement $v$. En introduisant dans le modèle une telle couche-interface, on peut identifier les équations constitutives de cette interface qui relient les déformations généralisées de cisaillement d'interface $\Omega_x^{j,j+1}(x,y)$ aux efforts généralisés d'interface de cisaillement $\tau_x^{j,j+1}(x,y)$ et d'arrachement $v^{j,j+1}(x,y)$ (voir[18]):

$$ 
\begin{aligned}
\Omega_{\alpha}^{j,j+1}(x,y) &= \frac{2e^{j,j+1}(1+v^{j,j+1})}{E^{j,j+1}} \tau_{\alpha}^{j,j+1}(x,y) - g_{\alpha}^{j,j+1}(x,y) + \Omega_{\alpha}^{j,j+1p}(x,y) & \alpha, \beta = 1,2 \\
\Omega_{z}^{j,j+1}(x,y) &= \frac{e^{j,j+1}}{E^{j,j+1}} v^{j,j+1}(x,y) + \Omega_{z}^{j,j+1p}(x,y)
\end{aligned}
\quad (2) $$

où $\Omega_{\alpha}^{j,j+1p}(x,y)$, $\Omega_{z}^{j,j+1p}(x,y)$ sont les glissements dus aux déplacements plastiques de la colle en matériau plastique parfait incompressible obéissant au critère de Von-Mises et,

$$ g_{\alpha}^{j,j+1}(x,y) = e^{j,j+1} \left( \frac{U_{3,\alpha}^{j}(x,y)e^{j+1} + U_{3,\alpha}^{j+1}(x,y)e^{j}}{e^{j} + e^{j+1}} \right) \quad (3) $$

$U_{3,\alpha}^{j}(x,y)$ sont les déplacements généralisés normaux à la plaque de la couche $j$. On trouvera des détails dans [18]. Une adaptation de l'algorithme de retour radial (ou de projection) est implémenté dans le code d'éléments finis MPFEAP pour résoudre le problème.

Un joint à double recouvrement qui subit un cycle de charge-décharge-charge est considéré (Fig.4, cellule hachurée). Afin de tester l'algorithme, le schéma de chargement est sous la forme charge-décharge-charge (fig. 5).

Après une étude de convergence du maillage, les contraintes de cisaillement et d'arrachement et les glissements de l'interface-colle obtenus par MPFEAP sont comparées avec ceux obtenus par le modèle 3D dans ABAQUS Fig.5. Ces contraintes et ces glissements généralisés du modèle
<figure><img src="image_5.png" /><figcaption>Figure 4 – Joint à double recouvrement étudié et maillages 2D utilisés</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5 – Schéma de chargement</figcaption></figure>

correspondent respectivement aux moyennes des contraintes 3D de la colle et aux déplacements relatifs entre les deux surfaces supérieure et inférieure de la colle.

# 5 Conclusion

Ainsi les avantages de MPFEAP sont les suivants :

- Comme on ne maille pas l'épaisseur, le maillage dans MPFEAP est beaucoup plus simple que celui dans ABAQUS et le nombre de degrés de liberté est considérablement diminué. Dans l'application sur joint collé, un rapport de 1/100 est montré entre le nombre de maille MPFEAP et ABAQUS. Bien évidemment, dans le cas d'un multicouche avec plusieurs interfaces, ce rapport diminue encore d'avantage avec le nombre d'interface.

- Le temps de calcul de MPFEAP est donc beaucoup plus faible que celui d'ABAQUS (25 minutes par rapport à 7 heures).

- Les résultats de MPFEAP ne sont pas singuliers aux bords et donc des critères de rupture, délaminage ou seuils plastiques portant sur ces contraintes au bord peuvent être proposés. C'est ce qui a été fait dans [17], ou des confrontations entre une approche analytique du problème de bord libre avec le *M4-5n* et une étude expérimentale de délaminage sur des composites carbone-epoxy a permis de valider des critères en contrainte maximale au bord (cisaillements et arrachement d'interface) et en taux de restitution d'énergie.

L'implémentation de cet élément dans Abaqus est en cours.

# Références

[1] O. Allix and P. Ladevèze, Interlaminar interface modelling for the prediction of laminate delamination. *Compos Struct*, 22, pp. 235-242, 1992.
[2] P. Ladevèze, G. Lubineau, D. Marsal, Towards a bridge between the micro- and mesomechanics of delamination for laminated composites. *Composites Science and Technology*, Volume 66, Issue 6, Pages 698-712, 2006.

[3] D. Leguillon, G. Marion, R. Harry and F. Lécuyer, The onset of delamination at stress-free edges in angle-ply laminates - analysis of two criteria. *Composites Science and Technology*, Volume 61, Issue 3, Pages 377-382, 2001

[4] F. Lebon, R. Rizzoni, S. Ronel-Idrissi, Asymptotic analysis of some non-linear soft thin layers, *Computers and Structures*, 82, 1929-1938, 2004

[5] T. Naciri, A. Chabot, A. Ehrlacher, Interlaminar stress analysis with a new multiparticle modelization of multilayered materials (M4), *Composites Science and Technology*, 58, 337-379, 1998

[6] J.F. Caron., A. Ehrlacher, Modelling the kinetics of transverse cracking in composite laminates. *Composites Science and Technology*, 57, 1261-1330, 1997

[7] N.J. Pagano, Stresses fields in composite laminates. *International Journal of Solids and Structures*, 14, 385-400, 1978.

[8] R.P. Carreira, J.F. Caron, A. Diaz Diaz, Model of multilayered materials for interface stresses estimation and validation by finite element calculations, Mechanics of Materials, 34, 217-246, 2002.

[9] O. Polit, M. Touratier, A multilayered/sandwich triangular finite element applied to linear and non-linear analyses. *Composites tructures*, Vol.58 Issue 1 Pages : 121-128, 2002

[10] E. Carrera, C0 Reissner-Mindlin multilayered plate elements including zig-zag and interlaminar stress continuity. *Int. J. Numer. Meth. Eng.* 39, pp. 1797-1820, 1996

[11] V.T. Nguyen, J.F. Caron, A new finite element for free edge effect analysis in laminated composites. *Computers Structures*, 84, 1538-1546, 2006.

[12] A.K. Noor. Free vibrations of multilayered composite plates. *AIAA Journal*, 11 :1038-1039, 1972.

[13] R.K. Khare, T. Kant, and A.K. Garg. Free vibration of composite and sandwich laminates with a higher-order facet shell element. *Composite Structures*, 65(3-4) :405-418, 2004.

[14] J.N. Reddy and N.D. Phan. Stability and vibration of isotropic, orthotropic and laminated plates according to a higher-order shear deformation theory. *Journal of Sound and Vibration*, 98(2) :157-226, 1985.

[15] K. Karas. Platten unter seitchem stoss. *Ingenieur Archiv*, 10 :237-250, 1939

[16] I.H. Choi and C.H. Lim. Low-velocity impact analysis of composite laminates using linearized contact law. *Composite Structures*, 66(1-4) :125-132, 2004.

[17] J.F. Caron and A. Diaz Diaz and R.P. Carreira and A. Chabot and A. Ehrlacher, Multi-particle modelling for the prediction of delamination in multi-layered materials, *Composites Science and Technology*, 2006, 66(6), p.755-765.

[18] G.S.A.L. Rios and R.C. Balderas and V.A. Duong and J.F. Caron and A. Ehrlacher and G. Foret and A. Diaz Diaz, Laminated plates with plastic interfaces : Modeling and calculation, *Modelling Simul. Mater. Sci. Eng.* 17 (2009)

[19] Z.Q. Cheng and L.H. Hea and S. Kitipornchaib, Influence of imperfect interfaces on bending and vibration of laminated composite shells, *International Journal of Solids and Structures*, 2000, 37(15), p.2127-2150,

[20] Z. Q. Yue and J.H. Yin, Backward Transfer-Matrix Method for Elastic Analysis of Layered Solids with Imperfect Bonding, *emphJournal of Elasticity*, 1998, 50(2), p.109-128.