Validation d'un modèle cohésif code_aster / MFront pour le délaminage des composites stratifiés (DCB)

1. Objectif

Le délaminage est le mode de ruine dominant des composites stratifiés, et sa simulation prédictive est un enjeu industriel majeur : elle permet de dimensionner par exemple des assemblages collés, des raidisseurs d'aile d'avion ou des structures sous impact sans multiplier les essais. Cet article montre que de telles simulations sont aujourd'hui accessibles avec des outils open-source : une loi cohésive écrite en MFront [1] est intégrée dans code_aster via des éléments finis d'interface.

Le cas d'étude retenu est le benchmark d'Alfano & Crisfield [2] : un essai Double Cantilever Beam (DCB) [3] sur une éprouvette composite stratifiée à fibres de carbone. Son intérêt est de posséder une solution analytique de la réponse force-déplacement pendant la propagation de la fissure, prédite par la mécanique linéaire de la rupture (LEFM) [4]. L'objectif est donc de comparer cette référence à la simulation code_aster + MFront, et de montrer qu'elle la reproduit fidèlement.

2. Le problème Double Cantilever Beam (DCB) [3]

2.1 Géométrie et données

L'éprouvette est un stratifié unidirectionnel (selon $x$) pré-fissuré sur une longueur initiale $a_0$. Chaque bras est chargé par un déplacement d'ouverture imposé $\pm u$, auquel correspond une réaction $P$.

Grandeur	Symbole	Valeur	Unité
Géométrie
Longueur totale	$L$	100	mm
Largeur	$b$	20	mm
Épaisseur d'un pli	$t$	1.5	mm
Épaisseur totale	$h=2t$	3	mm
Longueur de pré-fissure	$a_0$	30	mm
Matériau (pli)
Module longitudinal	$E_{11}$	135 300	MPa
Module de cisaillement	$G_{12}$	5 200	MPa
Coefficient de Poisson	$\nu_{12}$	0.24	-
Matériau (interface)
Taux de restitution d'énergie critique (mode I)	$G_c$	0.28	N/mm
Contrainte à rupture	$\sigma_c$	1.7	MPa

3. La solution analytique de référence (LEFM)

La rupture des matériaux fragiles et quasi-fragiles se décrit depuis Griffith [4] par un bilan d'énergie : une fissure ne se propage que si l'énergie élastique libérée par son avancée suffit à créer les nouvelles surfaces. Irwin [5] a ensuite reformulé ce critère sous une forme directement exploitable en ingénierie, en introduisant le taux de restitution d'énergie $G$, l'énergie disponible par unité de surface fissurée, et son seuil critique $G_c$ au-delà duquel la fissure se propage.

L'essai DCB possède une solution analytique fermée fondée sur cette approche énergétique, ainsi $G$ prend la forme :

$$G = \frac{P^2}{2b}\,\frac{\mathrm{d}C}{\mathrm{d}a}\tag{1}$$

Où $C = 2u/P=\delta/P$ est la souplesse (inverse de la rigidité) directement mesurable expérimentalement.

Une fissure assouplit toujours la structure ($\mathrm{d}C/\mathrm{d}a > 0$), donc $G$ est toujours positif : de l'énergie est toujours disponible pour fissurer, la question étant seulement de savoir s'il y en a assez pour atteindre $G_c$.

Propagation de la fissure. En combinant le critère de Griffith $G = G_c$ (énergie empiriquement déterminée) avec la théorie des poutres appliquée aux deux bras, la réaction $P$ prend la forme:

$$\boxed{\; P(u) = \frac{\lambda}{\sqrt{u}}, \qquad \lambda(G_c) = 62.1\ \mathrm{N\,mm^{1/2}} \;}\tag{2}$$

Physiquement, à mesure que la fissure avance, la structure s'assouplit : la charge nécessaire décroît alors que l'ouverture imposée augmente, d'où une réponse qui chute continûment. La théorie situe l'amorçage à $R_i = 68.8$ N pour $u_i = 0.81$ mm, et la sortie de fissure ($a=L$) pour un déplacement des bras théorique $u_{max} = 9.04$ mm.

La sollicitation est appliquée en déplacement imposé : $G \propto \delta^2/a^4$ décroît alors fortement avec $a$, si bien que la propagation s'arrête d'elle-même tant que l'on n'écarte pas davantage les bras. Ce pilotage rend donc la propagation intrinsèquement stable et dispense de méthodes numériques plus avancées (longueur d'arc ou intégration temporelle dynamique).

4. Le modèle de zone cohésive

Là où la théorie LEFM suppose une fissure préexistante et une contrainte infinie en pointe, le modèle de zone cohésive, introduit par Barenblatt [6] et Dugdale [7], place sur le trajet de la fissure une loi reliant la traction à l'ouverture des lèvres de fissure. Cette dernière s'exprime avec deux paramètres, chacun rattaché à un critère différent: une résistance $\sigma_c$ qui gouverne l'amorçage, et une énergie $G_c$ qui gouverne la propagation et redonne le critère de Griffith. Dans la loi bilinéaire retenue, la traction croît linéairement avec l'ouverture jusqu'au pic $\sigma_c$, puis décroît jusqu'à annulation ; l'aire sous cette courbe traction-séparation vaut exactement $G_c$ et le modèle retrouve la théorie LEFM lorsque $\sigma_c \to \infty$.

$$G_c = \int_{0}^{\delta_c} \sigma(\delta)\text{d}\delta = \tfrac{1}{2}\,\sigma_c\,\delta_c \tag{3}$$

L'interface est discrétisée par des éléments d'épaisseur nulle reliant deux nappes de nœuds coïncidents mais distincts. La cinématique n'est pas une déformation mais un saut de déplacement de part et d'autre de l'interface:

$$\delta = u^+-u^-\tag{4}$$

et la variable duale n'est pas un tenseur des contraintes mais un vecteur traction :

$$\mathbf{\sigma} = t_n \tag{5}$$

Avec $t_n$ la composante d'effort de traction normale (mode I), on note $\sigma_c$ la contrainte critique, $\delta_0$ son ouverture associée et $\delta_c$ l'ouverture critique. La loi normale comporte quatre branches :

$$\sigma(\delta) = \begin{cases} \alpha K\,\delta & \delta < 0 & \text{(compression, pénalisation)}\\[4pt] K\,\delta & 0 \le \delta \le \delta_0 & \text{(élastique)}\\[4pt] (1-d)\,K\,\delta & \delta_0 < \delta^{\max} < \delta_c & \text{(endommagement)}\\[4pt] 0 & \delta^{\max} \ge \delta_c & \text{(rupture)} \end{cases} \tag{6}$$

On introduit $d$ la variable d'endommagement, qui représente ici la dégradation progessive de l'interface avant sa rupture complète. Cette variable varie continûment entre $[0,1]$: $0$ interface saine, $1$ interface rompue et vient dégrader la cohésion de l'interface suivant:

$$(1-d)\,K\,\delta \tag{7}$$

avec

$${\ d(\delta) = \frac{1}{\eta}\left(1 - \frac{\delta_0}{\delta}\right)\ } \quad \text{avec} \quad \eta \equiv 1 - \frac{\delta_0}{\delta_c} \tag{8}$$

De plus, l'irréversibilité de l'endommagement impose que $d$ soit une fonction monotone croissante de l'ouverture maximale atteinte $\delta^{\max}$ ainsi numériquement nous imposons:

$$d \leftarrow \min\!\Big(\max\big(d^{\,\text{prev}},\, d(\delta)\big),\ 1\Big) \tag{9}$$

Un paramètre gouverne l'équilibre numérique, la longueur de zone cohésive [6,7] $\ell_{cz} \sim E_{11} \,G_c/\sigma_c^2$, que le maillage doit résoudre. Comme $\ell_{cz} \propto 1/\sigma_c^2$, le choix de $\sigma_c$ arbitre entre fidélité physique (résistance élevée, zone de fissuration courte) et robustesse numérique (résistance faible, zone plus étendue). La valeur retenue ici, $\sigma_c = 1.7$ MPa, privilégie la robustesse.

Les aspects numériques, d'intégration de la loi matériau MFront, de l'écriture et du calcul de la matrice tangente ne seront pas abordés ici, faisant d'ores et déjà l'objet d'articles approfondis ailleurs.

5. Mise en œuvre dans code_aster avec MFront

Le workflow complet repose sur des outils libres : CAO et maillage sous Salome, loi de matériau écrite en MFront, calcul sous code_aster, post-traitement sous ParaVis (module adaptant Paraview dans Salome).

Loi cohésive: La loi bilinéaire du modèle de zone cohésive est écrite en MFront, compilée à la volée et intégrée dans code_aster comme comportement du joint. Écrire la loi séparément permet de la modifier sans recompiler code_aster.

Maillage: L'interface à la zone de fissuration est discrétisée par des éléments d'épaisseur nulle reliant deux nappes de nœuds dédoublés et confondus (modélisation 3D_JOINT). Les bras sont quant à eux modélisés avec un modèle 3D.

Conditions aux limites : C'est le point qui conditionne la réussite du calcul. Un DCB est une structure souple, tenue uniquement par le joint intact en avant de la fissure. On impose donc uniquement un déplacement d'ouverture selon $z$ sur les fibres de valeur respective $+u$ et $-u$, tout en laissant les sections libres de tourner. On ne bloque que quelques nœuds ponctuels pour lever les modes de corps rigide restants.

Résolution. Le calcul est mené en analyse non-linéaire avec incréments et itérations adaptatifs, en pilotage par déplacement imposé.

6. Résultats et validation

Avec $\sigma_c = 1.7$ MPa, on compare la courbe issue de la simulation avec la référence analytique $P=62.1/\sqrt{u}$.

L'accord sur la branche de propagation est excellent (0.29 %) : c'est la signature de la cohérence énergétique du modèle, qui dissipe exactement $G_c$. La phase initiale plus souple et le pic abaissé ($\sigma_c = 1.7$ MPa) ne sont pas des défauts d'implémentation, mais la conséquence directe de $\ell_{cz} \gg L$ : toute l'interface s'endommage dès le début du chargement, sans pointe de fissure localisée. C'est le comportement documenté dans [2] pour les faibles résistances. Un $\sigma_c$ plus élevé rétablirait le pic, au prix d'une convergence plus délicate.

7. Conclusion

Un modèle éléments finis 3D de délaminage en mode I, entièrement construit avec des outils libres (Salome-Meca, MFront, code_aster, ParaVis), reproduit quantitativement le benchmark d'Alfano & Crisfield [1], avec un accord de l'ordre du pour-cent sur la propagation de la fissure par rapport à la solution analytique. Cet accord repose sur trois ingrédients bien identifiés : des conditions aux limites préservant la souplesse de la structure, une loi cohésive relâchant complètement la traction à rupture, et une stratégie de résolution adaptée à la branche d'endommagement adoucissante.

Ce résultat montre que code_aster couvre un besoin habituellement réservé aux codes commerciaux : la simulation prédictive du délaminage. La même démarche s'étend naturellement au mode mixte, au multi-délaminage et aux structures collées réelles.

Références

[1] T. Helfer et al., Introducing the open-source MFront code generator, Computers & Mathematics with Applications, 70 (2015) 994–1023.

[2] G. Alfano, M. A. Crisfield, Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues, International Journal for Numerical Methods in Engineering, 50 (2001) 1701–1736.

[3] ASTM D5528, Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites.

[4] A. A. Griffith, The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society A, 221 (1921) 163–198.

[5] G. R. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate, Journal of Applied Mechanics, 24 (1957) 361–364.

[6] G. I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, Advances in Applied Mechanics, 7 (1962) 55–129.

[7] D. S. Dugdale, Yielding of steel sheets containing slits, Journal of the Mechanics and Physics of Solids, 8 (1960) 100–104.