# Modélisation tridimensionnelle pour le calcul de singularités à la jonction de matériaux anisotropes

Vincent Magnier, Géry de Saxcé

► To cite this version:

Vincent Magnier, Géry de Saxcé. Modélisation tridimensionnelle pour le calcul de singularités à la jonction de matériaux anisotropes. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01412122

HAL Id: hal-01412122

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

Submitted on 8 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
# Modélisation tridimensionnelle pour le calcul de singularités à la jonction de matériaux anisotropes

V. Magnier¹, G. de Saxcé¹

¹ Laboratoire de Mécanique de Lille (USTL, Centrale Lille, ENSAM, CNRS)
Boulevard Paul Langevin, 59655 VILLENEUVE D'ASCQ
magnier.vincent@gmail.com ; gery.desaxce@univ-lille1.fr

## Résumé —

Dans l'étude de l'endommagement des multicouches, des phénomènes de "surcontraintes" apparaissent. Néanmoins, leur modélisation échappe largement au traitement numérique par éléments finis réguliers à moins de mailler extrêmement finement. Pour palier cette difficulté, nous proposons une méthode numérique qui s'appuie sur une extension de la mécanique de la rupture aux cas des matériaux anisotropes. Nous avons développé une méthode d'éléments finis mixtes enrichis. Après avoir expliqué la méthode, nous proposons un cas original qui est celui d'un insert plongé dans un quadri-couches.

**Mots clés** — Éléments finis mixtes enrichis, singularité, matériaux anisotropes, rupture fragile.

## 1 Introduction

De nos jours, les matériaux composites prennent une place de plus en plus conséquente en raison de leur légèreté et de leur rigidité à l'effort. Cependant, le comportement des matériaux composites est difficile à prédire à cause de discontinuités géométriques et matérielles. En effet, ces discontinuités engendrent un état tridimensionnel local de contrainte pouvant être singulier dans le sens où ces contraintes tendent localement vers l'infini comme le prévoit la mécanique élastique linéaire de la rupture.

Ces contraintes, aussi appelées "surcontraintes", sont souvent néfastes pour la rigidité de la structure provoquant généralement des sources d'amorçage locale de fissure. De par sa complexité, ce type de problème reste aujourd'hui l'un des défis à la conception des structures composites stratifiées. Sa compréhension et sa maîtrise suscitent toujours un vif intérêt autant chez les industriels et que pour les universitaires.

Le problème de surcontraintes ne peut être traité en pratique par des méthodes numériques telles que celles des éléments finis car il échappe largement au traitement numérique par éléments finis réguliers tridimensionnels.

Pour ce faire, nous proposons de déterminer précisément les champs mécaniques en nous basant sur une extension de la mécanique élastique linéaire de la rupture (ou mécanique de la rupture fragile) aux cas des multi-matériaux anisotropes. Cette théorie prévoit que quel que soit le corps Ω considéré contenant une ligne singulière, appelée Γf, alors le champ des contraintes solution du problème d'élasticité linéaire posé sur ce corps s'écrit en fond d'entaille :

$$ \sigma \sim r^{\alpha} K_{\alpha} \text{ lorsque } r \text{ est voisin de } 0 $$

r est la distance d'un point appartenant au matériau au front de fissure. (1)

$K_\alpha$ sont les facteurs d'intensité des contraintes.

$\alpha$, appelé "exposant de singularité", est de partie réelle négative.
La figure (1) présente quelques exemples de structures propices à l'amorçage local de fissure.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Configurations de composite présentant potentiellement une surcontrainte</figcaption></figure>

Il est à noter que l'exposant de singularité $\alpha$ ne dépend que de la géométrie de l'entaille et des propriétés des matériaux avoisinants ; il est en particulier indépendant de la géométrie du corps et du chargement. En revanche les fonctions $K_\alpha$ dépendent de la géométrie du corps et du chargement appliqué.

La première section de ce papier a pour objectif de présenter une méthode permettant d'extraire l'exposant de singularité $\alpha$. Nous développerons une méthode semi-analytique qui s'appuie sur une étude asymptotique d'un voisinage de ligne singulière couvrant une large gamme de situations rencontrées en pratique. La valeur de $\alpha$ sera obtenue par la résolution d'un problème aux valeurs propres dont le système traduit l'expression des conditions aux limites.

Les facteurs d'intensité des contraintes $K_\alpha$ résulteront d'un problème par éléments finis enrichis. La méthode sera détaillée dans la deuxième section. Pour cela, nous utiliserons une formulation d'élément fini hybride, appelé "métis" initialement proposée par Nguyen [13].

Enfin dans la dernière section, nous présenterons un exemple complet d'un quadri-couches avec un insert en traction.

# 2 Calcul de la singularité $\alpha$

## 2.1 Introduction

La mécanique élastique linéaire de la rupture prévoit que dans un matériau homogène, l'exposant de singularité $\alpha$ est égal à $-\frac{1}{2}$. Cependant, ce n'est pas forcément le cas dans les matériaux hétérogènes où celui-ci peut être complexe tel que sa partie réelle est comprise strictement entre 0 et -1 (pour des raisons d'énergie finie et de déplacement fini).

Ainsi de nombreux auteurs se sont intéressés à cette problématique. L'une des premières tentatives pour résoudre le problème de bord libre en élasticité anisotrope fut proposée par Pipes et Pagano [18] par une méthode de différences finies. Par la suite, d'autres méthodes purement numériques (éléments finis, éléments finis modifiés, éléments de frontière...) ont été utilisées pour résoudre ce genre de problème.

Parallèlement, des méthodes analytiques ont permis de résoudre des cas relativement simples comme les bi-matéraux isotropes, des multi-matéraux isotropes... Cependant, l'extension est difficile pour les matériaux anisotropes à cause du couplage entre les modes élémentaires de rupture (Modes I, II et III).

Pour palier ces difficultés, des méthodes semi-analytiques, donc des méthodes offrant un cadre
théorique avec une résolution numérique, sont apparues (Ting [20]...). C'est dans ce dernier cadre que nous nous placerons.

## 2.2 Principe mixte d'Hellinger-Reissner

Une idée pour caractériser les surcontraintes est définir ces celles-ci comme étant la solution du principe mixte d'Hellinger-Reissner :

$$ \begin{cases} \text{Trouver } (u, \underline{\sigma}) \text{ tel que :} \\ \forall \hat{u} : \langle \underline{\sigma}, \partial_M \hat{u} \rangle = \{f, \hat{u}\} \\ \forall \hat{\underline{\sigma}} : \langle \underline{\hat{\sigma}}, \partial_M \underline{u} \rangle = \langle \underline{\hat{\sigma}}, \underline{S} \cdot \underline{\sigma} \rangle \end{cases} \quad (2) $$

f et $\underline{\underline{S}}$ représentent respectivement les forces volumiques et la matrice de souplesse. Dans un second temps, nous décomposerons le problème en trois parties, respectivement plane, transverse et normale dans le repère local (P;e₁,e₂,e₃) associé au front de fissure (voir Fig. (2)) :

$$ \underline{\underline{\sigma}} = \begin{pmatrix} \underline{\underline{\sigma}}_T & \underline{\underline{\sigma}}_S \\ \underline{\underline{\sigma}}_S^\top & \underline{\underline{\sigma}}_N \end{pmatrix} \qquad \underline{\underline{\varepsilon}} = \begin{pmatrix} \underline{\underline{\varepsilon}}_T & \underline{\underline{\varepsilon}}_S \\ \underline{\underline{\varepsilon}}_S^\top & \underline{\underline{\varepsilon}}_N \end{pmatrix} \quad (3) $$

Puis nous effectuons une approche asymptotique de type "couche-limite" sur le principe d'Hellinger-Reissner en faisant tendre le paramètre r vers zéro pour obtenir le comportement dominant comme l'illustre la Figure (2).

<figure><img src="image_3.png" /><figcaption>Figure 2 – Repère local et repère global</figcaption></figure>

## 2.3 Équations de compatibilité

Le problème nous conduit aux équations de Beltrami-Mitchell dans lequel nous introduisons les potentiels d'Airy $\phi$ et de Prandtl $\psi$. Nous obtenons finalement un système dans lequel il faut déterminer $\phi$ et $\psi$ :

$$ \left[ \begin{matrix} \mathcal{D}_4 & \mathcal{D}_3 \\ \mathcal{D}_3 & \mathcal{D}_2 \end{matrix} \right] \left\{ \begin{matrix} \phi \\ \psi \end{matrix} \right\} = 0 \qquad (4) \\ \| \mathcal{D}_i \text{ est un opérateur d'ordre } i $$
## 2.4 Solutions singulières de S.S. Wang

L'expression de $\phi$ fut donnée par S.S. Wang [21] sous la forme d'un quotient où l'on voit apparaître explicitement $\alpha$ qui est l'exposant de singularité à déterminer :

$$ \phi(Z_k) = \frac{C_k Z_k^{\alpha+2}}{(\alpha+2)(\alpha+1)} \qquad (5) $$

$\alpha$ est l'exposant de singularité à déterminer
$C_k$ est une constante à déterminer
$Z_k = \xi^1 + \bar{\omega}_k \xi^2$ avec $\bar{\omega}_k$ sont les racines de l'équation caractéristique $(4) \in \mathbb{C}$

Finalement, la forme explicite du couple $(\underline{\sigma}^0, u^0)$ à l'ordre dominant est :

$$ \underline{\sigma}^0 = \begin{pmatrix} \underline{\sigma}_T^0 & \underline{\sigma}_S^0 \\ \underline{\sigma}_S^{0\top} & \underline{\sigma}_N^0 \end{pmatrix} \Rightarrow \left\{ \begin{array}{l} \underline{\sigma}_T^0 = \sum_{k=1}^3 \mathcal{R}e(C_k \underline{\tau}_T^k Z_k^\alpha + C_{k+3} \overline{\underline{\tau}_T^k} \overline{Z_k}^\alpha) \\ \underline{\sigma}_S^0 = \sum_{k=1}^3 \mathcal{R}e(C_k \underline{\tau}_S^k Z_k^\alpha + C_{k+3} \overline{\underline{\tau}_S^k} \overline{Z_k}^\alpha) \\ \underline{\sigma}_n^0 = \frac{-S_{nt'} \underline{\sigma}_{t'} - S_{ns'} \underline{\sigma}_{s'}}{S_{nn}} \end{array} \right. \qquad (6) $$
$$ \begin{array}{l} \underline{\tau}_T^k = \bar{\omega}_k^2 e_1 \cdot e_1^\top + e_2 \cdot e_2^\top - \bar{\omega}_k (e_1 \cdot e_2^\top + e_2 \cdot e_1^\top) \\ \underline{\tau}_S^k = \eta_k \bar{\omega}_k e_1 - \eta_k e_2 \end{array} $$

$$ \mathbf{u}^0 = \begin{pmatrix} \mathbf{u}_t^0 \\ \mathbf{u}_3^0 \end{pmatrix} \Rightarrow \left\{ \begin{array}{l} u_t^0 = \sum_{k=1}^3 \mathcal{R}e[C_k V_t^k \frac{Z_k^{\alpha+1}}{\alpha+1} + C_{k+3} \overline{V_t^k} \frac{\bar{Z}_k^{\alpha+1}}{\alpha+1}]] \\ u_3^0 = \sum_{k=1}^3 \mathcal{R}e[C_k V_3^k \frac{Z_k^{\alpha+1}}{\alpha+1} + C_{k+3} \overline{V_3^k} \frac{\bar{Z}_k^{\alpha+1}}{\alpha+1}]] \end{array} \right. $$
$$ \begin{array}{l} V_t^k = (\tilde{S}_{11}\bar{\omega}_k^2 + \tilde{S}_{12} - \tilde{S}_{16}\bar{\omega}_k + \tilde{S}_{15}\bar{\omega}_k\eta_k - \tilde{S}_{14}\eta_k)e_1 \\ \phantom{V_t^k = (} + (\tilde{S}_{21}\bar{\omega}_k + \tilde{S}_{22}/\bar{\omega}_k - \tilde{S}_{26} - \tilde{S}_{24}\eta_k/\bar{\omega}_k + \tilde{S}_{25}\eta_k)e_2 \\ V_3^k = \tilde{S}_{14}\bar{\omega}_k + \tilde{S}_{42}/\bar{\omega}_k - \tilde{S}_{44}\eta_k/\bar{\omega}_k + \tilde{S}_{45}\eta_k - \tilde{S}_{46} \end{array} $$

Où $\tilde{S}_{ij}$ sont les composantes de la matrice de souplesse réduite.

## 2.5 Condition aux limites

Afin de déterminer les inconnues $C_k$ et l'exposant de singularité $\alpha$ introduits dans la fonction (5), il faut préciser les conditions aux limites.

Considérons le cas général où la ligne singulière se situe à la jonction de $n$ secteurs angulaires $\omega_k$. Chaque secteur est constitué d'un matériau comme il est indiqué sur la figure (3). Il est plus commode alors d'introduire un repère pseudo-cylindrique $(r, \rho_k, \theta)$ caractérisé par :

$$ Z_k = \xi^1 + \bar{\omega}_k(\xi^3)\xi^2 = r\rho_k \quad ; \quad \rho_k = \cos\theta + \bar{\omega}_k(\xi^3)\sin\theta \qquad (7) $$

Dorénavant, la notation contractée sera utilisée et les composantes du tenseur des contraintes $\underline{\sigma}^0$ seront stockées dans un vecteur $\underline{\sigma}^0$ :

$$ \left\{ \begin{array}{l} \underline{\sigma}^0 = \mathcal{R}e \sum_{\alpha} r^{\alpha} (\underline{\tau} \cdot \underline{\rho}^{\alpha} \cdot \mathbf{A} + \bar{\underline{\tau}} \cdot \bar{\rho}^{\alpha} \cdot \mathbf{B}) \\ \mathbf{u}^0 = \mathcal{R}e \sum_{\alpha} \frac{r^{\alpha+1}}{\alpha+1} (\underline{V} \cdot \underline{\rho}^{\alpha+1} \cdot \mathbf{A} + \bar{\underline{V}} \cdot \bar{\rho}^{\alpha+1} \cdot \mathbf{B}) \end{array} \right. \qquad (8) $$

avec les matrices suivantes :

$$ \begin{array}{l} \underline{\tau} = [\tau^1, \tau^2, \tau^3] \quad \text{avec} \quad \tau^k = [\tau_1^k, \tau_2^k, \mathbf{0}, \tau_4^k, \tau_5^k, \tau_6^k]^\top \\ \underline{V} = [V^1, V^2, V^3] \quad \text{avec} \quad V^k = \mathbf{V}_t^k + V_3^k \cdot e_3 \\ \underline{\rho} = \text{Diag}(\rho_1, \rho_2, \rho_3) \\ \mathbf{A}^\top = (C_1, C_2, C_3) \quad ; \quad \mathbf{B}^\top = (C_4, C_5, C_6) \end{array} $$
<figure><img src="image_4.png" /><figcaption>Figure 3 – Voisinage de la ligne singulière $\Gamma_f$ à la jonction de $n$ secteurs angulaires $\omega_k$</figcaption></figure>

Dans les sections qui suivent, nous allons présenter différentes conditions aux limites (bord libre et continuité en contraintes et déplacements) que nous pouvons rencontrer sur une interface donnée. Chacune d'entre elles impose aux champs des déplacements et/ou des contraintes de respecter certaines conditions.

### 2.5.1 Condition de bord libre

Ce cas signifie que la surface est libre de contrainte. Ceci se traduit par $\sigma^0.n = 0$, où $n$ est le vecteur unitaire normal à la surface. Étant donnée que cette condition est valable quel que soit $r$, pour une surface donnée, nous imposons aux équations (6) qu'elles respectent :

$$ \left[ \begin{matrix} \underline{t} \cdot \underline{\rho}^{\alpha} & \underline{t} \cdot \bar{\underline{\rho}}^{\alpha} \\ \underline{0} & \underline{0} \end{matrix} \right] \left\{ \begin{matrix} A \\ B \end{matrix} \right\} = \underline{0} \qquad (9) \\ \| \underline{t} = [t^1, t^2, t^3] \text{ tel que } t^k = \underline{\tau}^k \cdot n $$

### 2.5.2 Condition de continuité

Cette condition, établie sur une interface entre deux secteurs adjacents $\omega_k$ et $\omega_{k+1}$, conduit à la continuité des champs des déplacements $[[u^0]] = u_{k+1}^0 - u_k^0 = 0$ et à celle des champs des contraintes $[[\sigma^0.n]] = \sigma_{k+1}^0 - \sigma_k^0 = 0$, $n$ étant le vecteur unitaire normal à l'interface comme il est illustré sur la figure (3). Cette condition s'écrit :

$$ \left( \left[ \begin{matrix} \underline{t} \cdot \underline{\rho}^{\alpha} & \underline{t} \cdot \bar{\underline{\rho}}^{\alpha} \\ \underline{V} \cdot \underline{\rho}^{\alpha+1} & \underline{V} \cdot \bar{\underline{\rho}}^{\alpha+1} \end{matrix} \right] \left\{ \begin{matrix} A \\ B \end{matrix} \right\} \right) = \underline{0} \qquad (10) \\ \text{où } \underline{t} = [t^1, t^2, t^3] \text{ tel que } t^k = \underline{\tau}^k \cdot n $$

## 2.6 Résolution du problème aux valeurs propres

La méthode se résume finalement à un problème aux valeurs propres où l'exposant de singularité $\alpha$ apparaît explicitement :

$$ \underline{\Lambda}(\alpha) \cdot \underline{X} = \underline{0} $$

La construction de $\underline{\Lambda}(\alpha)$ et de $\underline{X}$ se fait par un assemblage de toutes les conditions entre les différents matériaux autour d'une ligne singulière. Les solutions sont localisées par la méthode de
la fausse position (Regula-Falsi). Une des difficultés de la compilation est la présence de valeurs propres multiples. Ainsi, pour les détecter, nous utiliserons le critère de Rouché.

$$ \text{Multiplicité de } \alpha = \frac{1}{2\pi i} \oint_C \frac{\partial_\alpha[\det(\Delta)]}{\det(\Delta)} dC = \frac{1}{2\pi} \oint_C \text{Argument}[\det(\Delta)] dC \quad (11) \\ \Vert C \text{ est un lacet entourant la solution } \alpha $$

Finalement, l'expression des déplacements et des contraintes près de la ligne singulière $\Gamma_f$ s'écrit :

$$ \begin{aligned} \boldsymbol{\sigma} &= \sum_{\alpha} \mathcal{R}e[K_{\alpha} r^{\alpha} \tilde{\sigma}] \\ \boldsymbol{u} &= \sum_{\alpha} \mathcal{R}e[K_{\alpha} r^{\alpha+1} \tilde{u}] \end{aligned} \quad (12) $$

L'intégralité de la méthode est expliquée dans [10].

Cependant, l'importance d'une singularité est liée non seulement à son exposant mais aussi aux facteurs d'intensité des contraintes associés. Une méthode d'évaluation de ces paramètres font l'objet de la prochaine section.

# 3 Calcul des facteurs d'intensité de contraintes $K_\alpha$

L'objet de cette section est de déterminer les facteurs d'intensité des contraintes $K_\alpha$ introduits dans l'équation (12).

## 3.1 Position du problème

Généralement, les facteurs d'intensité des contraintes sont obtenus numériquement par des méthodes par post-traitement (citons à titre d'exemples l'intégrale-J [19], la $\theta$-méthode [7]...). Cependant, avec ses méthodes pré-citées, la difficulté essentielle réside dans l'estimation précise des champs des contraintes et des déformations qui doivent être singuliers comme le prédit la théorie de la mécanique linéaire de la rupture. Elle s'explique par l'utilisation de fonctions de type polynômial peu propices à représenter des singularités dans la formulation des éléments finis classiques. Pour palier ces difficultés, une des méthodes proposées est d'utiliser un maillage fin local près des lignes singulières. Cependant, on comprend aisément que dans le but d'obtenir une économie de coût numérique, une limitation du nombre d'éléments finis est souhaitable tout en préservant une précision satisfaisante des résultats. Ceci a conduit à l'élaboration de nouvelles méthodes d'éléments finis traitant de problèmes fissurés. Une avancée significative (et l'une des plus connues) dans l'utilisation des éléments finis pour la mécanique élastique linéaire de la rupture fut l'apparition de la technique du "nœud de quart" introduite par Henshell [9] et Barsoum [3]. Ces chercheurs ont montré que les déplacements en fond de fissure, ainsi que les champs des contraintes et des déformations, peuvent être modélisés en manipulant certains nœuds dans un élément quadratique. La technique consiste à déplacer les nœuds du milieu, B et C, au quart du côté de l'élément, B' et C', en se rapprochant du nœud singulier A (voir Fig. (4)). Ce processus engendre alors une singularité dans la transformation entre les coordonnées d'espace paramétrique de l'élément et d'espace cartésien permettant aux champs des contraintes et des déformations d'être inversement proportionnels à la racine carrée de la distance d'un point au front de fissure. Cependant la technique de Barsoum nécessite tout de même un maillage fin local. L'autre convénient de la méthode est la génération d'une singularité $\alpha$ égale à -0.5 peu convenable aux matériaux hétérogènes.

Une autre possibilité de déterminer les facteurs d'intensité des contraintes est de passer par la méthode des éléments finis enrichis. Cette méthode sera décrite dans les sections suivantes avec une formulation hybride nommée "Métis".
<figure><img src="image_5.png" /><figcaption>Figure 4 – Élément quadratique et élément quadratique de Barsoum centré sur le point A</figcaption></figure>

## 3.2 L'élément fini "Métis"

### 3.2.1 Définition

Le modèle Métis est une classe particulière de modèle hydride qui constitue un pont entre les modèles classiques à un champ d'inconnues et les modèles à deux champs d'inconnues. Pour définir ce modèle Métis, nous supposons que le champ des déplacements **u** est cinématiquement admissible sur le bord de l'élément, tandis que le champ des contraintes **σ** est un champ arbitraire en équilibre dans Ω. En outre, nous imposons la condition d'unisolvance sur le champ des déplacements **u**. Dans ce cas, **u** est défini et continu à l'intérieur de chaque élément et doit être conforme entre éléments (voir Fig. (5)).

La fonctionnelle de ce modèle s'écrit alors :

<figure><img src="image_6.png" /><figcaption>Figure 5 – Conditions intra et inter-élément pour la formulation métis</figcaption></figure>

$$ \Pi(u, \underline{\sigma}) = \int_{\Omega} \frac{1}{2} \underline{\sigma}^{\top} \underline{S} \cdot \underline{\sigma} d\Omega - \int_{\Omega_u \cup \Omega_{\sigma}} u^{\top} \underline{\sigma} \cdot n dS + \int_{\Omega_{\sigma}} u^{\top} \cdot \mathbf{T}^{\text{imposé}} d\Omega_{\sigma} \quad (13) $$

### 3.2.2 Avantages de la formulation

Le premier intérêt de la formulation métis est que la convergence monotone décroissante alors qu'elle est oscillante pour un élément hybride quelconque. En outre, pour un maillage fixé, la convergence est plus rapide que celle de l'élément déplacement pur (voir Nguyen Dang [13] qui le montre grâce à l'analyse convexe) comme illustré sur la figure (6). Le deuxième point d'intérêt est que le Lagrangien peut s'écrire en fonction d'intégrales de volume. En effet, la fonctionnelle d'énergie gouvernant l'élément fini Métis pour les matériaux anisotropes, s'écrit avec les notations indicielles :

$$ \Pi(\underline{\sigma}, u) = \int_{\Omega} \frac{1}{2} \underline{\sigma}^{\top} \underline{S} \cdot \underline{\sigma} d\Omega - \int_{\Omega_u \cup \Omega_{\sigma}} u^{\top} \underline{\sigma} \cdot n dS + \int_{\Omega_{\sigma}} u^{\top} \cdot \mathbf{T}^{\text{imposé}} dS \\ \left\| \mathbf{T}^{\text{imposé}} \text{ sont les tractions imposées sur } \Omega_{\sigma} \partial \Omega_{e} \right. \\ \left. \mathbf{n} \text{ est le vecteur unitaire tangent à la surface } \Omega_u \cup \Omega_{\sigma} \right. \quad (14) $$
<figure><img src="image_7.png" /><figcaption>Figure 6 – Convergence vers la solution de différents formulations d'éléments finis en fonction du degré de liberté</figcaption></figure>

Cette fonctionnelle (14) peut être décomposée en une somme de fonctionnelles sur les éléments (ou fonctionnelles élémentaires) :

$$
\begin{gathered}
\Pi(\underline{\sigma}, \mathbf{u}) = \sum_{e} \Pi_{e}(\underline{\sigma}, \mathbf{u}) \\
\left\|
\begin{aligned}
\Pi_{e}(\underline{\sigma}, \mathbf{u}) &= \int_{\Omega_{e}} \frac{1}{2} \underline{\sigma}^{\top} \underline{S} \underline{\sigma} d\Omega_{e} - \int_{\partial\Omega_{e}} \mathbf{u}^{\top} \underline{\sigma} \cdot \mathbf{ndS} + \int_{\partial\Omega_{e}} \mathbf{u}^{\top} \mathbf{T} dS \\
\Omega_{e} &\text{ représente le domaine d'un élément fini} \\
\partial\Omega_{e} &\text{ est la frontière du domaine de l'élément fini } \Omega_{e}
\end{aligned}
\right\| \tag{15}
\end{gathered}
$$

Comme nous l'avons écrit précédemment, le champ des contraintes **σ** doit être statiquement ad-
missible à l'intérieur de l'élément Ω<sub>e</sub> et le champ des déplacements **u** doit l'être cinématiquement
sur les frontières de celui-ci. Cependant, en plus de ces conditions, nous imposons au champ des
déplacements la condition d'unisolvance, ce qui signifie que le champ des déplacements est pro-
longé à l'intérieur de chaque élément par un polynôme complet.

Ce principe permet de transformer de manière élégante la fonctionnelle (15) en une intégrale de
volume sans termes de bord :

$$
\int_{\partial \Omega_e} \mathbf{u}^{\top} \underline{\boldsymbol{\sigma}} \cdot \mathbf{ndS} = \int_{\Omega_e} \left[ \mathrm{Tr}(\underline{\boldsymbol{\sigma}}^{\top} \cdot \mathrm{Sym}(\partial_M \mathbf{u})) + \mathbf{u}^{\top} \cdot \mathrm{div}_M(\underline{\boldsymbol{\sigma}}) \right] d\Omega_e \quad (16)
$$

Étant donné que le champ des contraintes est statiquement admissible à l'intérieur de chaque élément fini, le Lagrangien élémentaire $\Pi_e(\underline{\sigma}, u)$ prend alors la forme :

$$
\Pi_e(\underline{\sigma}, \mathbf{u}) = \int_{\Omega_e} \left[ \frac{1}{2} \underline{\sigma}^{\top} \underline{S} \underline{\sigma} - \mathrm{Tr}(\underline{\sigma}^{\top} \cdot \mathrm{Sym}(\partial_M \mathbf{u})) \right] d\Omega_e + \int_{\partial\Omega_e} \mathbf{u}^{\top} \cdot \mathbf{T}^{\text{imposé}} dS \quad (17)
$$

Comme nous le voyons, la fonctionnelle ne comporte que des intégrales de volume ; cela a pour
conséquence d'obtenir une évaluation numérique plus précise, surtout quand on traite de matériaux
anisotropes. En effet, il est bien connu que les calculs d'intégrales de volume sont plus robustes,
plus stables et plus faciles à appliquer dans un code que les calculs d'intégrales sur un contour,
qui sont fortement "numerically path-dependent" (le lecteur intéressé pourra aller voir page 300
de Nguyen QS [15] et le papier de Banks-Sills [2]).

Le troisième point d’intérêt concerne la bonne tenue des éléments Métis vis-à-vis des distorsions géométriques. En traitant une poutre de Saint-Venant soumise à la traction, Crépin [6] montre que ces éléments se comportent très bien quel que soit l’angle de distorsion des éléments. Ceci est directement lié au fait que le théorème de Brezzi [4] a servi à la construction de la formulation discrète, de sorte que cette dernière possède les bonnes propriétés de la formulation continue ; celles-ci sont trop souvent absentes dans les formulations mixtes d’éléments tridimensionnels comme dans l’élément hybride de contraintes de [17].
## 3.3 L'élément fini "Métis" enrichi

### 3.3.1 Discrétisation en éléments finis "réguliers" et "singuliers"

Nous allons enrichir les éléments finis près de la ligne singulière (éléments singuliers) en incorporant la solution asymptotique (12) trouvée dans la première section. En supposant que ces surcontraintes sont confinées au voisinage immédiat de la ligne singulière, il est raisonnable de les tronquer en les incorporant uniquement au maillage voisin des bords de la ligne singulière $\Gamma_f$. Les autres éléments de la structure sont considérés comme réguliers (voir Fig. (7)).

<figure><img src="image_8.png" /><figcaption>Figure 7 – Maillage d'une structure $\Omega$ contenant une ligne singulière $\Gamma_f$ par éléments finis singuliers et réguliers</figcaption></figure>

Dans notre étude, nous supposerons que la fissure ne coïncide avec les nœuds de certains éléments finis.

Ce type de méthode simple et générale autorise ainsi des maillages relativement lâches sans perte de précision significative sur les singularités de contrainte.

### 3.3.2 Éléments finis "réguliers"

Le champ des déplacements est discrétisé par une approximation nodale $u = N(x)q$ telle que $N(x)$ représente la matrice contenant les fonctions d'interpolation $N_i(x)$ et $q$ est le vecteur des déplacements nodaux.

Les nœuds géométriques (correspondant aux coins de l'élément) sont confondus avec les nœuds d'interpolation. Cet élément comporte donc 24 variables.

Les fonctions de forme $N_i(x)$ sont celles de Lagrange pour un élément hexaédrique linéaire :

$$ N_i(x) = \frac{1}{8}(1+x_ix)(1+y_iy)(1+z_iz) \quad \text{avec } i \in \{1;..;8\} \quad (18) \\ \| \quad x_i, y_i \text{ et } z_i \text{ représentent les coordonnées des nœuds d'interpolation} $$

Dans la suite, la notation contractée sera utilisée. Donc, le champ des déformations est interpolé comme suit :

$$ \varepsilon = \underline{B} \cdot q \\ \| \underline{B} = \partial_{D_f} \cdot N(x) \quad \text{tel que} \quad \partial_{D_f} = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \\ 0 & \frac{\partial}{2\partial z} & \frac{\partial}{2\partial x} \\ \frac{\partial}{2\partial z} & 0 & \frac{\partial}{2\partial x} \\ \frac{\partial}{2\partial y} & \frac{\partial}{2\partial x} & 0 \end{bmatrix} \quad (19) $$
En ce qui concerne les contraintes régulières, elles sont exprimées de la façon suivante :

$$
\boldsymbol{\sigma}^R = \{\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5, \sigma_6\} = \underline{\mathbf{R}}_R \cdot \mathbf{h}_R \qquad (20)
$$

Le vecteur $h_R$ contient les paramètres des contraintes régulières. La matrice d'interpolation $\underline{\mathbf{R}}_R$ est construite de sorte que les deux conditions du théorème de Brezzi [4] soient satisfaites. De façon plus précise, la matrice $\underline{\mathbf{R}}_R$ est construite comme ceci :

$$
\underline{\mathbf{R}}_R =
\begin{bmatrix}
y & z & xz & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & x & z & xz & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & y & x & xy & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & x \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & y & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & z & 0
\end{bmatrix}
\quad (21)
$$

Discrétisons le Lagrangien élémentaire (17) dans le cas des éléments réguliers :

$$
\begin{align}
\Pi_e(\mathbf{h}_R, q) &= \frac{1}{2} \mathbf{h}_R^\top \underline{\mathbf{F}}_R \cdot \mathbf{h}_R - \mathbf{h}_R^\top \underline{\mathbf{G}}_R^\top \cdot q + \mathbf{L}_e^\top \cdot q \nonumber \\
&\qquad \left\| \underline{\mathbf{F}}_R = \int_{\Omega_e} \underline{\mathbf{R}}_R^\top \cdot \underline{\mathbf{S}} \cdot \underline{\mathbf{R}}_R d\Omega_e \right. \nonumber \\
&\qquad \left. \underline{\mathbf{G}}_R = \int_{\Omega_e} \underline{\mathbf{B}}^\top \cdot \underline{\mathbf{R}}_R d\Omega_e \right. \nonumber \\
&\qquad \left. L_e = \int_{\partial \Omega_e} \underline{\mathbf{N}}^\top \cdot T dS \right. \tag{22}
\end{align}
$$

La dérivée du Lagrangien élémentaire (22) par rapport aux paramètres des contraintes $\mathbf{h}_R$ donne :

$$
\underline{\mathbf{F}}_R \cdot \mathbf{h}_R - \underline{\mathbf{G}}_R^\top \cdot \mathbf{q} = \mathbf{0} \qquad (23)
$$

ce qui est une forme discrétisée des équations de compatibilité pour un élément.
D'autre part, en considérant la variation du Lagrangien élémentaire $\Pi_e(\mathbf{q}, \mathbf{h}_R)$ par rapport à $\mathbf{q}$, nous
obtenons :

$$
\underline{\mathbf{G}}_R \cdot \mathbf{h}_R - \mathbf{L}_e = \mathbf{0} \qquad (24)
$$

ce qui est une forme discrétisée des équations d'équilibre.
Au vu des équations (24) et (23), nous proposons d'utiliser la méthode de condensation statique
ou de sous-structuration statique introduite par Imbert [8]. Cette méthode est une technique nu-
mérique qui permet de diminuer le nombre d'inconnues des modèles éléments finis ; les degrés
de liberté du modèle initial sont séparés en degrés de liberté de jonction (ou maîtres) et degrés de
libertés internes (ou esclaves) qui vont être éliminés. Donc, en proposant de considérer $\mathbf{h}_R$ comme
degrés de libertés internes, $q$ devient par conséquent les degrés de liberté de jonction et le problème
à résoudre pour les éléments réguliers se résume à :

$$
\underline{\mathbf{G}}_R \cdot \underline{\mathbf{F}}_R^{-1} \cdot \underline{\mathbf{G}}_R^\top \cdot \mathbf{q} = \mathbf{L}_e \qquad (25)
$$

**3.3.3 Éléments finis "singuliers"**

Dans cette approche, on s'intéresse aux éléments touchant le front de fissure.

Étant donné que l'élément Métis est un élément hybride, donc à deux champs de variables (u, σ),
l'enrichissement est double, dans le sens où chaque variable sera enrichie. Nguyen Dang [14] ap-
pelle ceci la "double singularité"; une singularité de contrainte et une singularité de déplacement.
Antérieurement, différents auteurs ont travaillé avec l'élément métis enrichi. Tout d'abord, Nguyen
Dang qui, dans sa thèse, propose de générer les singularités, dans les éléments singuliers, en utilisant la technique de Barsoum et en décomposant les contraintes en une partie régulière et singulière. La partie singulière est directement déduite de l'approche asymptotique de la première partie. Cependant, l'ensemble de la méthode nécessite d'utiliser des éléments finis quadratiques. De plus, comme nous l'avions vu, la technique de Barsoum génère automatiquement un exposant de singularité égal à -0.5, ce qui n'est pas adapté aux problèmes d'entaille à la jonction de multi-matériaux anisotropes.

Ainsi, nous avons développé une nouvelle méthode en nous appuyant sur la méthode dite de la partition de l'unité (Melenk [11]).

L'approximation du champ des déplacements classique $u^R$, issue du calcul par éléments finis, est enrichie par une correction $u^S$ sous la forme :

$$ \mathbf{u} = u^R + u^S \qquad (26) $$

Cette approximation, inspirée des travaux de Moës [12], ajoute des degrés de liberté dans l'élément, dans le but de mieux représenter le phénomène physique à évaluer pour tout type de problème.

$$ \mathbf{u}^S = N(x).[\gamma(\xi^i, \mu_k, \alpha + 1).h_S] \\ \text{où} \quad \begin{array}{l} \mathbf{\gamma} = [\mathbf{u}_1^0, \dots, \mathbf{u}_k^0, \dots, \mathbf{u}_m^0], \quad k \text{ est le numéro du nœud.} \\ \bar{\mathbf{u}}_k^0 \text{ est le déplacement singulier au nœud } k \text{ décrit dans l'équation (12).} \end{array} \qquad (27) $$

La construction du champs des contraintes se fait sur le même principe que pour les déplacements ; ces champs sont décomposés en une partie régulière $\sigma^R$ et une partie singulière $\sigma^S$ :

$$ \sigma = \sigma^R + \sigma^S \qquad (28) $$

La partie régulière est évaluée exactement de la même manière que pour les éléments réguliers. Les fonctions additionnelles pour les contraintes seront celles décrites dans le premier chapitre (12).

$$ \begin{gathered} \sigma^S = \underline{R}_S \cdot h_S \\ \underline{R}_S = [\sigma_1^0, \dots, \sigma_m^0], m \text{ étant le numéro du nœud} \\ \sigma_k^0 \text{ est le vecteur des contraintes singulières au nœud } k \text{ décrit dans (12)} \end{gathered} \qquad (29) $$

La forme discrétisée du Lagrangien est :

$$ \begin{gathered} \Pi_e(\mathbf{h}_R, h_S, q) = \frac{1}{2}h_R^\top \cdot \underline{\mathbf{F}}_R \cdot h_R + h_S^\top \cdot \underline{\mathbf{F}}_{SR} \cdot h_R + \frac{1}{2}\mathbf{h}_S^\top \cdot \underline{\mathbf{F}}_S \cdot h_S + L_e^\top \cdot q \\ - \mathbf{h}_R^\top \cdot \underline{\mathbf{G}}_R^\top \cdot q - \mathbf{h}_S^\top \cdot \underline{\mathbf{G}}_S^\top \cdot q - \mathbf{h}_R^\top \cdot \underline{\mathbf{H}}_{RS} \cdot h_S - \mathbf{h}_S^\top \cdot \underline{\mathbf{H}}_S \cdot \mathbf{h}_R \\ \underline{\mathbf{F}}_{SR} = \int_{\Omega_e} \underline{\mathbf{R}}_S^\top \cdot \underline{\mathbf{S}} \cdot \underline{\mathbf{R}}_S d\Omega_e \qquad \underline{\mathbf{G}}_S = \int_{\Omega_e} \underline{\mathbf{B}}^\top \cdot \underline{\mathbf{R}}_S d\Omega_e \\ \underline{\mathbf{F}}_S = \int_{\Omega_e} \underline{\mathbf{R}}_S^\top \cdot \underline{\mathbf{S}} \cdot \underline{\mathbf{R}}_S d\Omega_e \qquad \underline{\mathbf{H}}_S = \int_{\Omega_e} \underline{\mathbf{R}}_S^\top \cdot \underline{\mathbf{A}}_S d\Omega_e \\ \underline{\mathbf{H}}_{RS} = \int_{\Omega_e} \underline{\mathbf{R}}_R^\top \cdot \underline{\mathbf{A}}_S d\Omega_e \end{gathered} \qquad (30) $$

L'annulation de la variation première de la fonctionnelle (30) par rapport aux paramètres du champ des déplacements $q$ conduit aux équations d'équilibre :

$$ \frac{\partial \Pi_e(\mathbf{h}_R, h_S, q)}{\partial q} = \mathbf{L}_e - \underline{\mathbf{G}}_R \cdot \mathbf{h}_R - \underline{\mathbf{G}}_S \cdot \mathbf{h}_S = \mathbf{0} \qquad (31) $$

L'annulation de la dérivée première de l'équation (30) par rapport aux paramètres de contrainte $h_R$ et $h_S$ conduit aux conditions de stationnarité :

$$ \begin{gathered} \frac{\partial \Pi_e(\mathbf{h}_R, h_S, q)}{\partial h_R} = \underline{\mathbf{F}}_R \cdot h_R + \underline{\mathbf{F}}_{SR}^\top \cdot h_S - \underline{\mathbf{G}}_R^\top \cdot q - \underline{\mathbf{H}}_{RS} \cdot h_S = \mathbf{0} \\ \frac{\partial \Pi_e(h_R, h_S, q)}{\partial h_S} = \underline{\mathbf{F}}_{SR} \cdot h_R + \underline{\mathbf{F}}_S \cdot h_S - \underline{\mathbf{H}}_{RS}^\top \cdot h_R - \underline{\mathbf{G}}_S^\top \cdot q - \underline{\mathbf{H}}_S \cdot h_S - \underline{\mathbf{H}}_S^\top \cdot h_S = \mathbf{0} \end{gathered} \qquad (32) $$
<table><caption>Tableau 2 – Valeurs des intégrales en fonction du nombre de points de Gauss et erreurs relatives avec le schéma de Gauss-Radau.</caption><thead><tr><th colspan="2">Gauss-Radau</th></tr><tr><th>Fonction</th><th>18 Points de Gauss</th></tr></thead><tbody><tr><td>F<sub>1</sub><br>(erreur relative)</td><td>1%</td></tr><tr><td>F<sub>2</sub><br>(erreur relative)</td><td>8%</td></tr><tr><td>F<sub>3</sub><br>(erreur relative)</td><td>8%</td></tr></tbody></table>

<table><caption>Tableau 3 – Valeurs des intégrales en fonction du nombre de points de Gauss et erreurs relatives avec le schéma de Gauss-Jacobi.</caption><thead><tr><th colspan="2">Gauss-Jacobi</th></tr></thead><tbody><tr><td>F<sub>1</sub><br>(erreur relative)</td><td>1.2 * 10<sup>-4</sup>%</td></tr><tr><td>F<sub>2</sub><br>(erreur relative)</td><td>1.5 * 10<sup>-4</sup>%</td></tr><tr><td>F<sub>3</sub><br>(erreur relative)</td><td>2.1 * 10<sup>-4</sup>%</td></tr></tbody></table>

Au vu des précisions sur les résultats, le schéma de Gauss-Radau et celui de Gauss-Jacobi se distinguent très clairement par leur précision. En effet, l'erreur maximale pour ce type de fonctions singulières testes est de 8%. Ces deux schémas ont été privilégiés pour le calcul d'intégrales dans les éléments singuliers.

# 4 Exemple du rivet dans un quadri-couches orthotropes

## 4.1 Présentation du problème

Les inserts dans les structures stratifiées sont souvent sources d'amorçage de fissure. Cependant, le problème des singularités introduites est difficile à résoudre en raison de son caractère tridimensionnel. Néanmoins, grâce à la formulation présentée précédemment, la résolution de ce type de problème, présentant une courbure, nous est facilement accessible. Nous proposons donc l'exemple d'un insert en aluminium plongé dans un quadri-couches de graphite/epoxy. Les propriétés élastiques de chaque constituant sont résumées dans le tableau (4). Il est à noter que les propriétés mécaniques de chaque couche du carbone/epoxyde sont obtenues en utilisant la "loi des mélanges". Ici, $E_{ii}$, $G_{ij}$ et $v_{ij}$ représentent respectivement les modules de Young, les modules

<table><caption>Tableau 4 – Propriétés élastiques d'une couche de Carbone/Epoxyde et de l'Aluminium.</caption><thead><tr><th></th><th>E<sub>11</sub><br>(GPa)</th><th>E<sub>22</sub><br>(GPa)</th><th>E<sub>33</sub><br>(GPa)</th><th>G<sub>12</sub><br>(GPa)</th><th>G<sub>13</sub><br>(GPa)</th><th>G<sub>23</sub><br>(GPa)</th><th>v<sub>12</sub></th><th>v<sub>13</sub></th><th>v<sub>23</sub></th></tr></thead><tbody><tr><td>Graphite/Epoxyde</td><td>137.9</td><td>14.48</td><td>14.48</td><td>5.86</td><td>5.86</td><td>5.86</td><td>0.21</td><td>0.21</td><td>0.21</td></tr><tr><td>Aluminium</td><td></td><td>72.</td><td></td><td></td><td>27.07</td><td></td><td></td><td>0.33</td><td></td></tr></tbody></table>

de cisaillement et les coefficients de Poisson dans le repère matériau. Les indices 1, 2 et 3 du repère se réfèrent respectivement à la direction longitudinale des fibres, à la direction transverse et à l'épaisseur d'une couche.

La stratification du quadri-couches sera fixée par la suite à $[45°/ - 45°]_S$. Sur cette éprouvette, nous effectuerons un essai de traction (suivant $X_1$) d'une intensité de 10MPa.
<figure><img src="image_9.png" /><figcaption>Figure 8 – Rivet en aluminium plongé dans un quadri-couches de graphite/epoxyde ayant une stratification $[45^{\circ}/-45^{\circ}]_s$</figcaption></figure>

## 4.2 Calcul de la singularité α

<figure><img src="image_10.png" /><figcaption>Figure 9 – Numérotation des lignes singulières</figcaption></figure>

Six lignes devraient être, *a priori*, sujettes à un comportement singulier, ce qui sera confirmé ou non par notre code de calcul. La première ligne $\Gamma_1$ se situe au bord de la tête du rivet reposant sur la surface extérieure du composite. La seconde ligne $\Gamma_2$ se trouve à l'intersection du pli supérieur avec le rivet au niveau du trou. Les troisième $\Gamma_3$ et quatrième ligne $\Gamma_4$ concerne une discontinuité de propriétés mécaniques. Et finalement, la dernière ligne $\Gamma_5$ est un bord droit entre deux matériaux différents. Étant donné que le rivet au niveau des lignes singulières a une forme cylindrique de révolution, ces cinq lignes singulières sont des cercles. Un paramètre $\theta$ est introduit pour désigner la position sur chaque cercle (voir Fig. (9)).

**Notre but est de déterminer l'évolution de $\alpha$ en fonction de l'angle $\theta$ sur chaque ligne singulière.**
Il est évident que les configurations des lignes singulières n°3 et n°4 sont les mêmes, par conséquent nous ne tracerons qu'une courbe sur les deux. La ligne $\Gamma_1$ présente trois exposants propres dont chacun ont une périodicité de phase $\pi$. Ceci est la conséquence d'une symétrie "miroir" du problème physique et des propriétés mécaniques. La ligne $\Gamma_2$, quand à elle, n'est pas représentée graphiquement car elle n'est pas sujette à un comportement singulier. Ceci est directement aux valeurs des propriétés mécaniques car pour la même géométrie mais avec un différentiel plus grand des modules élastiques, nous avons constaté la présence de valeur propre. Le problème lié aux lignes $\Gamma_3$ et $\Gamma_4$ présente qu'une valeur propre possédant une périodicité de $\frac{\pi}{2}$. Cette périodicité s'explique par le fait que si nous nous plaçons sur un point du cercle $\Gamma_3$ ou $\Gamma_4$ à un angle $\theta$ à
<figure><img src="image_11.png" /><figcaption>Figure 10 – Exposant de singularité pour un rivet universel plongé dans un quadri-couches de graphite/epoxyde ayant une stratification $[45^{\circ}/-45^{\circ}]_{s}$</figcaption></figure>

qui correspond une stratification $[\phi^1/\phi^2]$, alors à un angle $\theta + \pi/2$ correspondra la stratification $[\phi^2/\phi^1]$, expliquant la périodicité d'angle $\pi/2$. Enfin la ligne $\Gamma_5$ est une alternance de comportement singulier et régulier. Il n'y a qu'un exposant propre de périodicité $\pi$. Il est à noter que la résolution de problème global des singularités prend moins d'une seconde sur un PC cadencé à 1Ghz.

## 4.3 Calcul des Facteurs d'Intensité des Contraintes $K_{\alpha}$

Pour cela, nous avons utilisé un maillage constitué de 8394 éléments comme l'illustre la figure (11).

<figure><img src="image_12.png" /><figcaption>Figure 11 – Maillage utilisé pour la détermination des Facteurs d'Intensité des Contraintes K</figcaption></figure>

Notre but est de déterminer les facteurs d'intensité des contraintes $K_{\alpha}$. Sur les figures (12) et (13) sont représentées les facteurs d'intensité des contraintes associés aux lignes $\Gamma_1, \Gamma_3, \Gamma_4$ et $\Gamma_5$.
<figure><img src="image_13.png" /><figcaption>(a) FIC pour la ligne Γ₁</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) FIC pour les lignes Γ₃ et Γ₄</figcaption></figure>

Figure 12 – Facteurs d’intensité des contraintes d’un un quadri-couches de graphite/epoxyde pos-sédant un rivet soumis à de la traction

Il est à noter que les facteurs d'intensité des contraintes associés aux lignes Γ₃ et Γ₄ sont les mêmes. Comme nous pouvons le constater, les facteurs d'intensité des contraintes associé à la ligne singulière Γ₁ sont les plus importants. Remarquons aussi qu'une périodicité de phase Π existe étant donné les propriétés du problème traité.
<figure><img src="image_15.png" /><figcaption>(a) FIC pour la ligne Γ<sub>5</sub></figcaption></figure>

Figure 13 – Facteurs d’intensité des contraintes d’un un quadri-couches de graphite/epoxyde pos-
sédant un rivet soumis à de la traction

**5 Conclusion**

Dans ce papier, nous nous sommes intéressés au développement d'éléments finis enrichis pour
caractériser des paramètres de la mécanique de la rupture dans un contexte tridimensionnel au
sein de multi-matériaux anisotropes. Pour cela, nous avons proposé une méthode se décompos-
sant en deux étapes. La première consiste à déterminer l'exposant de singularité en s'inspirant de
modèle type "couche-limite". Dans une seconde étape, nous évaluons les facteurs d'intensité des
contraintes (F.I.C) en passant par un calcul par éléments finis.

L'intégralité de la méthode a été montrée à travers l'exemple de l'insert plongé dans un quadri-couches. Il est intéressant de remarquer que la méthode présentée ne requiert qu'un faible temps CPU.

Références

[1] JE. Akin. The generation of element with singularities. *International Journal for Numerical Methods in Engineering*, 10, 1249-1259, 1976.

[2] L. Banks-Sills. A conservative integral for determining stress intensity factors of a bimaterial strip. *International Journal of Fracture*, 86, 385-398, 1997.

[3] RS. Barsoum. On the use of isoparametric finite elements in linear fracture mechanics. *International Journal for Numerical Methods in Engineering*, 10, 25-37, 1976.

[4] F. Brezzi, M. Fortin. *Mixed and hybrid finite element methods*, Springer Verlag, Berlin, 1991.

[5] SK. Chan, IS. Tuba, WK. Wilson. On the finite element method in linear fracture mechanics. *Engineering Fracture Mechanics*, 2, 1-17, 1970.

[6] D. Crepin. Thèse de doctorat. Éléments finis mixtes singuliers tridimensionnels pour l'étude des surcontraintes dans les structures en matériaux composites, Université des Sciences et Technologies de Lille 1, 2001.
[7] P. Destuynder, M. Djaoua. Sur une interprétation mathématique de l’intégrale de Rice en théorie de la rupture fragile. *Mathematical Methods in the Applied Sciences*, Éditeur, 3, 70-87, 1981.

[8] JF. Imbert. *Analyse des structures par éléments finis*, Cepadues, 1991.

[9] RD. Henshell, KG. Shaw. Crack tip finite elements are unnecessary. *International Journal for Numerical Methods in Engineering*, 9, 495-507, 1975.

[10] V. Magnier. Thèse de doctorat. *Modélisation tridimensionnelle par éléments finis enrichis pour le calcul de singularités de délaminage et à la jonction de matériaux anisotropes*, Université des Sciences et Technologies de Lille 1, 2008.

[11] JM. Melenk, I. Babuska. The partition of unity finite element method : basic theory and applications. *Computer Methods in Applied Mechanics and Engineering*, 139, 289-314, 1996.

[12] N. Moes, J. Dolbow, T. Belytschko. A finite element method for crack growth without reme-shing. *International Journal for Numerical Methods in Engineering*, 139, 289-314, 1999.

[13] H. Nguyen Dang. On the monotony and the convergence of a special class of hybrid finite element : the mongrel element. *Variational Methods in Mechanics of Solids*, 1, 208-213, 1978.

[14] H. Nguyen Dang, G. de Saxce, K. Che-Hang. The computation of 2-D stress intensity factors using hybrid mongrel displacement finite elements. *Engineering Fracture Mechanics*, 38, 197-205, 1991.

[15] QS. Nguyen. *Stability and Non-Linear Solid Mechanics*, Wiley, New-York, 2000.

[16] D. Nguyen Tien, H. Nguyen Dang. Regular and singular metis finite element models for delamination in composite laminates. *Finite Elements in Analysis and Design*, 42, 650-659, 2006.

[17] THH. Pian. Derivation of element stiffness matrices by assumed stress distributions. *AIAA Journal*, 2, 1333-1336, 1964.

[18] RB. Pipes, NJ. Pagano. Interlaminar stresses in composite laminates under uniform axial extension. *Journal of Composite Materials*, 39, 538-548, 1970.

[19] JR. Rice. *Mathematical analysis in the mechanics of fracture*, Fracture : a advanced treatise, II, Liebowitz H., Academic Press, 1968.

[20] TCT. Ting, SC. Cou. Edge singularities in anisotropic composites. *International Journal of Solids and Structures*, 17, 1057-1068, 1981.

[21] SS. Wang, I. Choi. Boundary layer effect in composite laminates, part I : Free-edge stress singularities ; part II : free-edge stress solutions and characteristics. *Journal of Applied Mechanics*, 49, 541-550, 1982.