# Modélisation du comportement inélastique avec endommagement de composites renforcés par des fibres courtes non parfaitement orientées

Delphine Notta-Cuvier, Franck Lauro, Bruno Bennani, Romain Balieu

► To cite this version:

Delphine Notta-Cuvier, Franck Lauro, Bruno Bennani, Romain Balieu. Modélisation du comportement inélastique avec endommagement de composites renforcés par des fibres courtes non parfaitement orientées. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717030

HAL Id: hal-01717030

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

Submitted on 25 Feb 2018

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

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

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Modélisation du comportement inélastique avec endommagement de composites renforcés par des fibres courtes non parfaitement orientées

Delphine NOTTA-CUVIER<fnref n="1" />*, Franck LAURO<fnref n="1" />, Bruno BENNANI<fnref n="1" />, Romain BALIEU<fnref n="1" />

**Résumé** — Le modèle proposé permet une prédiction du comportement de composites renforcés par fibres courtes sans recourir aux techniques complexes d'homogénéisation. Tous les types de comportements élastoplastiques peuvent être considérés pour la matrice, y compris des écoulements plastiques compressibles. Un modèle d'endommagement fortement anisotrope, induit par le renforcement, est également implémenté. Enfin, toutes les distributions d'orientation des fibres dans le matériau, y compris les orientations aléatoires, peuvent être modélisées simplement.

**Mots clés** — Fibres courtes, distributions d'orientation, plasticité non associée, critère de Drucker-Prager, endommagement anisotrope.

## 1. Introduction

Les hautes rigidités des fibres permettent un renforcement efficace de matrices métalliques ou polymères. Le comportement macroscopique d'un matériau composite à fibres courtes dépend fortement de la teneur en fibres et de leur orientation. Dans le cas de procédés par injection, celle-ci est distribuée autour de la direction d'injection [2]. Des orientations aléatoires sont aussi fréquemment rencontrées [5]. Le modèle proposé considère le matériau composite comme la réunion d'un milieu matriciel élastoplastique endommageable et de plusieurs milieux unidimensionnels et élastiques de fibres.

Les fibres courtes sont supposées être uniformément dispersées dans le matériau matriciel. Chaque fibre est orientée selon un vecteur d'orientation exprimé dans le repère global, lié à la matrice. Les fibres ayant le même comportement et le même vecteur d'orientation, $\bar{a}^\alpha$, sont regroupées dans la « famille » $\alpha$. $n_{\text{fam}}$ familles de fibres sont ainsi définies, chacune étant caractérisée par une fraction volumique $v_F^\alpha$, telle que $\sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha = v_F = 1 - v_M$, où $v_F$ et $v_M$ sont les fractions volumiques totales des fibres et de la matrice dans le matériau composite, respectivement. Grâce à cette répartition des fibres en plusieurs familles, des distributions d'orientation complexes sont facilement modélisables en associant à chaque angle d'orientation une proportion donnée de fibres. La matrice a un comportement élastoplastique à endommagement ductile. En l'absence de fibres, celui-ci est isotrope. Dans le matériau composite, l'endommagement peut devenir complètement anisotrope, en fonction des caractéristiques du renforcement et est donc modélisé par un tenseur du 4ème ordre.

Le comportement macroscopique du composite est déterminé grâce à une décomposition additive du potentiel thermodynamique. Cette approche a l'avantage de ne pas nécessiter la mise en œuvre de procédures complexes d'homogénéisation ou de double homogénéisation [4].

<footnote n="1">1 Univ. Valenciennes et du Hainaut-Cambrésis, LAMIH</footnote>

<footnote>* Auteur correspondant : delphine.notta@univ-valenciennes.fr</footnote>
## 2. Modélisation

### 2.1. Endommagement fortement anisotrope de la matrice induit par les fibres

Le modèle d'endommagement du matériau matriciel obéit aux lois de la mécanique de l'endommagement continu et le concept de contrainte effective est utilisé [6]. Un tenseur d'endommagement du 4ème ordre, **D**, est défini pour relier le tenseur des contraintes vraies de la matrice, $\sigma_M$, au tenseur des contraintes effectives, $\tilde{\sigma}_M$ (1). Pour faciliter la notation, le tenseur **M** est également défini (2).

$$ \sigma_M = \mathbf{D}\tilde{\sigma}_M, \text{ i.e. } \sigma_{Mij} = \sum_{k,l=1}^{3} \mathbf{D}_{ijkl}\tilde{\sigma}_{Mkl} \quad (1) $$

$$ \tilde{\sigma}_M = \mathbf{M}\sigma_M, \text{ i.e. } \tilde{\sigma}_{Mij} = \sum_{k,l=1}^{3} \mathbf{M}_{ijkl}\sigma_{Mkl} \quad (2) $$

Les tenseurs **D** et **M** sont construits à partir des caractéristiques du renforcement. Plus précisément, ils résultent de l'assemblage de l'endommagement induit par chaque famille de fibres $\alpha$ sur une portion volumique $v_F^\alpha$ de la matrice. Les contraintes vraies, $\sigma_M^{\alpha 0}$, et effectives, $\tilde{\sigma}_M^{\alpha 0}$, sont définies pour un matériau fictif constitué de la matrice renforcée uniquement par la famille de fibres $\alpha$. L'exposant 0 indique que ces tenseurs sont définies dans le repère lié aux fibres (i.e. de premier vecteur $\bar{\alpha}^\alpha$). $\sigma_M^{\alpha 0}$ et $\tilde{\sigma}_M^{\alpha 0}$ sont reliées par le tenseur d'endommagement $\mathbf{D}^\alpha$ et par son « inverse » $\mathbf{M}^\alpha$, de façon similaire aux relations (1) et (2). Leurs composantes sont définies en supposant que les fibres empêchent l'endommagement de la matrice dans leur direction d'orientation et en introduisant deux variables scalaires d'endommagement : $D_n$ qui agit sur la contrainte hydrostatique et $D_s$ qui agit sur la contrainte déviatorique (3). Enfin, une portion $v_M$ de la matrice présente un endommagement isotrope, non influencé par la présence des fibres. Finalement, l'assemblage de ces différentes contributions permet d'exprimer les tenseurs **D** et **M** dans le repère global (4-5).

$$ \left\{ \begin{aligned} \mathbf{D}_{ijkl}^{\alpha} &= \delta_{ik} \delta_{jl} [1 - D_n \delta_{ij} (1 - \delta_{il}) - D_s (1 - \delta_{ij})] \\ \mathbf{M}_{ijkl}^{\alpha} &= \delta_{ik} \delta_{jl} / [1 - D_n \delta_{ij} (1 - \delta_{il}) - D_s (1 - \delta_{ij})] \end{aligned} \right. \quad \forall \alpha \quad (3) $$

$$ \left\{ \begin{aligned} \mathbf{D}_{ijkl}^{\alpha} &= \delta_{ik} \delta_{jl} [1 - D_n \delta_{ij} (1 - \delta_{il}) - D_s (1 - \delta_{ij})] \\ \mathbf{M}_{ijkl}^{\alpha} &= \delta_{ik} \delta_{jl} / [1 - D_n \delta_{ij} (1 - \delta_{il}) - D_s (1 - \delta_{ij})] \end{aligned} \right. \quad \forall \alpha \quad (3) $$

$$ \mathbf{D}_{ijkl} = v_M (1 - D_n) \delta_{ik} \delta_{jl} + \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha \sum_{a,b=1}^{3} T_{ia}^\alpha T_{bj}^{\alpha-1} T_{ak}^{\alpha-1} T_{lb}^\alpha \mathbf{D}_{abab}^\alpha \quad \forall \alpha \quad (4) $$

$$ \mathbf{M}_{ijkl} = [v_M / (1 - D_n)] \delta_{ik} \delta_{jl} + \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha \sum_{a,b=1}^{3} T_{ia}^\alpha T_{bj}^{\alpha-1} T_{ak}^{\alpha-1} T_{lb}^\alpha \mathbf{M}_{abab}^\alpha \quad \forall \alpha \quad (5) $$

Les matrices $T^\alpha$ sont les matrices de passage du repère lié aux fibres $\alpha$ vers le repère global.

L'endommagement se développe dans la matrice à partir d'un seuil de déformation plastique cumulée ($p>p_D$). L'évolution des variables scalaires d'endommagement, $D_n$ et $D_s$, est déterminée dans le cadre de la plasticité non-associée couplée à l'endommagement. Pour cela, le potentiel de dissipation du matériau matriciel, $\psi$, est défini, tel que $\psi = \psi^P + \psi^D(D_n, D_s)$. $\psi^P$ est la partie élastoplastique du potentiel qui dépend des variables internes $\epsilon$ (tenseur des déformations toales), $\epsilon^P$ (tenseur des déformations plastiques, en supposant la partition des déformations) et $\epsilon^R$ (variable d'écrouissage isotrope), pour une transformation isotherme et en l'absence d'écrouissage cinématique. Les lois de normalité (6) permettent alors de relier les variables internes à leur variable associée, avec $\dot{\Lambda}$ le multiplicateur plastique. $Y_n$ et $Y_s$ sont les taux de densité de relaxation d'énergie associés aux variables d'endommagement $D_n$ et $D_s$, respectivement.
$$
\left\{
\begin{array}{l}
\dot{D}_n = \dot{\Lambda} \frac{\partial \psi}{\partial Y_n} = \dot{\Lambda} \frac{\partial \psi^D}{\partial Y_n} \qquad \dot{D}_s = \dot{\Lambda} \frac{\partial \psi}{\partial Y_s} = \dot{\Lambda} \frac{\partial \psi^D}{\partial Y_s} \\
\dot{\epsilon}^P = \dot{\Lambda} \frac{\partial \psi}{\partial \sigma_M} = \dot{\Lambda} \frac{\partial \psi^P}{\partial \sigma_M} \qquad \dot{r} = -\dot{\Lambda} \frac{\partial \psi}{\partial R} = -\dot{\Lambda} \frac{\partial \psi^P}{\partial R}
\end{array}
\right.
\qquad (6)
$$

La partie du potentiel liée à l'endommagement, $\psi^D$, est une fonction de $Y_n^2/(1-D_n)$ et de $Y_s^2/(1-D_s)$ [6]. $Y_n$ et $Y_s$ peuvent être exprimés à partir de la densité d'énergie complémentaire, définie par $w_e^* = 1/2 \epsilon^e : \sigma_M$, selon $Y_n = \partial w_e^* / \partial D_n$ et $Y_s = \partial w_e^* / \partial D_s$ [6]. En supposant que le comportement élastique de la matrice est isotrope linéaire et suit la loi de Hooke (paramètres $E_M$ et $v_M$), les lois d'évolution des variables d'endommagement $D_n$ et $D_s$ sont donc connues (7). L'expression de $\dot{D}_s$ est identique à celle de $\dot{D}_n$ en remplaçant $D_n$ par $D_s$ et $S_n$ par $S_s$. $S_n$ et $S_s$ sont deux paramètres matériaux.

$$
\dot{D}_n = \frac{\dot{\Lambda}}{2E_M S_n} \sum_{i,j=1}^{3} \left[ \begin{array}{l} (1 + v_M) \sigma_{Mij} \sum_{k,l=1}^{3} \frac{\partial \mathbf{M}_{ijkl}}{\partial D_n} \sigma_{Mkl} \\ - \left( v_M \text{tr}(\sigma_M) \sum_{k,l=1}^{3} \frac{\partial \mathbf{M}_{ijkl}}{\partial D_n} \sigma_{Mkl} \right) \delta_{ij} \end{array} \right]
$$

$$
\text{où } \frac{\partial \mathbf{M}_{ijkl}}{\partial D_n} = \frac{\left( v_M \delta_{ik} \delta_{jl} + \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha \sum_{a=2}^{3} T_{ia}^\alpha T_{aj}^{\alpha-1} T_{ak}^{\alpha-1} T_{la}^\alpha \right)}{(1 - D_n)^2} \quad \text{et} \quad \frac{\partial \mathbf{M}_{ijkl}}{\partial D_s} = \frac{\sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha \sum_{\substack{a,b=1 \\ a \neq b}}^{3} T_{ia}^\alpha T_{bj}^{\alpha-1} T_{ak}^{\alpha-1} T_{lb}^\alpha}{(1 - D_s)^2}.
$$

La modélisation du comportement élastoplastique de la matrice repose sur des principes bien connus et ne sera pas exposée en détail ci. Toutefois, un point important est que le potentiel élastoplastique $\psi^P$ est sensible à la pression, avec $\psi^P = \tilde{J}_2 + \gamma\tilde{\sigma}_H$. En l'absence d'écrouissage cinématique, $\tilde{J}_2$ est la contrainte équivalente de von Mises effective, définie par $\tilde{J}_2 = \sqrt{\frac{3}{2}\tilde{S}_M : \tilde{S}_M}$.

$\tilde{S}_M$ est la partie déviatorique de $\tilde{\sigma}_M$ et $\tilde{\sigma}_H$ est la contrainte hydrostatique effective associée. $\gamma$ est un paramètre matériau qui traduit la variation de volume au cours de l'écoulement plastique. $\gamma$ peut donc aussi dépendre de la nature du chargement appliqué (traction/compression notamment). Le critère de plasticité de Drucker-Prager, sensible en pression, est utilisé ce qui permet de considérer des limites d'élasticité différentes en compression et en traction, en particulier. Enfin, toutes les lois d'écrouissage, linéaires ou non, peuvent être considérées sans que cela ne nécessite de modification du cœur de l'implémentation, grâce au traitement en préambule de la loi d'écrouissage (approximation linéaire par morceaux).

Un schéma classique prédiction élastique + correction plastique [7] est utilisé pour la résolution du comportement élastoplastique matriciel. Les relations de mises à jour principales sont données par l'équation (8). L'exposant (T) désigne les grandeurs définies à l'état test, c'est-à-dire en supposant une évolution purement élastique de l'incrément n (déterminé) à l'incrément n+1. Dans les relations (8), **M** est défini à partir des variables $D_n^{(T)} = D_n^{(n)}$ et $D_s^{(T)} = D_s^{(n)}$. Les contraintes vraies sont ensuite calculées (1) en utilisant $D_n^{(T)} = D_n^{(n)}$ et $D_s^{(T)} = D_s^{(n)}$ dans l'expression de **D** (4). Les variables d'endommagement sont à leur tour mises à jour selon le schéma semi-couplé $\Delta D_n = g_n(\sigma_M^{(n+1)}(D_n^{(T)}, D_s^{(T)}), D_n^{(T)})$ et $\Delta D_s = g_s(\sigma_M^{(n+1)}(D_n^{(T)}, D_s^{(T)}), D_s^{(T)})$. Les expressions des fonctions $g_n$ et $g_s$ proviennent des équations (7).
$$
\left\{
\begin{aligned}
\epsilon^{p(n+1)} &= \epsilon^{p(T)} + \Delta\Lambda \left[ \frac{3 \mathbf{M} \tilde{\mathbf{S}}_{\mathrm{M}}^{(T)}}{2 \tilde{\mathbf{J}}_{2}^{(T)}} + \frac{\gamma}{3} \mathbf{M} \mathbf{I} \right] \\
\epsilon^{e(n+1)} &= \epsilon^{e(T)} - \Delta\Lambda \left[ \frac{3 \mathbf{M} \tilde{\mathbf{S}}_{\mathrm{M}}^{(T)}}{2 \tilde{\mathbf{J}}_{2}^{(T)}} + \frac{\gamma}{3} \mathbf{M} \mathbf{I} \right] \\
\tilde{\sigma}_{\mathrm{M}}^{(n+1)} &= \tilde{\sigma}_{\mathrm{M}}^{(T)} - \Delta\Lambda \left[ 3\mu_{\mathrm{M}} \frac{\mathbf{M} \tilde{\mathbf{S}}_{\mathrm{M}}^{(T)}}{\tilde{\mathbf{J}}_{2}^{(T)}} + \frac{2}{3}\mu_{\mathrm{M}} \gamma \mathbf{M} \mathbf{I} + \lambda_{\mathrm{M}} \left( \frac{3}{2} \frac{\mathrm{tr}(\mathbf{M} \tilde{\mathbf{S}}_{\mathrm{M}}^{(T)})}{\tilde{\mathbf{J}}_{2}^{(T)}} + \frac{\gamma}{3} \mathrm{tr}(\mathbf{M} \mathbf{I}) \right) \right] \mathbf{I}
\end{aligned}
\right.
\quad (8)
$$

2.2. Comportement élastique des fibres

Une hypothèse fondamentale est que les fibres ne se déforment que suivant leur axe. Elles ont un comportement linéaire élastique, de module d'Young E<sub>F</sub><sup>α</sup> dans la famille α. Le tenseur des gradients de transformation appliquée à chaque famille de fibres, F<sub>F</sub><sup>α</sup>, est la projection du tenseur des gradients de transformation appliquée au matériau composite, F, selon la direction d'orientation des fibres (9). Le tenseur F est intégralement transmis à la matrice. La relation (9) permet de relier les tenseurs droits de Cauchy-Green, C et C<sub>F</sub><sup>α</sup>, définis pour la matrice et chaque famille de fibres, respectivement (10).

$$
F_F^\alpha = F A^\alpha \quad \forall \alpha \quad \text{avec} \quad A^\alpha = \bar{a} \otimes \bar{a} \tag{9}
$$

$$
C_F^\alpha = A^\alpha C A^\alpha \quad \forall \alpha \tag{10}
$$

Par construction, chaque tenseur C<sub>F</sub><sup>α</sup> possède une unique valeur propre, ζ<sup>α</sup>, associée au vecteur
propre a<sup>α</sup>. D'après les propriétés des C<sub>F</sub><sup>α</sup>, ζ<sup>α</sup> est donc le carré de l'élongation axiale des fibres
regroupées dans la famille α. La déformation 1D de Hencky (petites déformations) des fibres est alors
donnée par ε<sub>F</sub><sup>α0</sup> = 1/2 ln(ζ<sup>α</sup>). Elle est localement commune à la fibre et à la matrice et est considérée
pour définir l'effort qui est réellement transmis de la matrice vers chaque fibre de la famille α. Pour
cela, la théorie de Bowyer et Bader [1] est adaptée aux cas 3D et à des orientations complexes de
fibres. Une correction supplémentaire est apportée en considérant que les parties endommagées de la
matrice ne peuvent transmettre aucun effort. Au final, la contrainte 1D des fibres dans chaque famille
α est donnée par la relation (11).

$$
\left\{
\begin{array}{ll}
\sigma_F^{\alpha 0} = \frac{1}{2}(1 - D_n)\ln(\zeta^\alpha)\left(1 - \frac{E_F^\alpha \ln(\zeta^\alpha) r^\alpha}{4L^\alpha \tau^\alpha}\right) E_F^\alpha & \text{si } L^\alpha > L_\epsilon^\alpha \\
\\
\sigma_F^{\alpha 0} = \frac{1}{2}(1 - D_n)\frac{L^\alpha \tau^\alpha}{2r^\alpha} E_F^\alpha & \text{sinon}
\end{array}
\right.
\quad \forall \alpha \qquad (11)
$$

avec L<sup>α</sup>, r<sup>α</sup> et τ<sup>α</sup> la longueur, le rayon et la contrainte interfaciale de cisaillement (IFSS) pour les
fibres dans la famille α, respectivement. L<sub>ε</sub><sup>α</sup> est une longueur critique de fibre définie pour chaque
valeur de déformation dans la fibre par L<sub>ε</sub><sup>α</sup> = (E<sub>F</sub><sup>α</sup>r<sup>α</sup>) / (τ<sup>α</sup>ε<sub>F</sub><sup>α0</sup>) [1]. La longueur des fibres, de même que
l'IFSS ont donc une influence importante sur l'état de contrainte des fibres et donc sur le
comportement macroscopique du composite. On peut souligner que l'IFSS dépend bien sûr de la
nature des fibres et de la matrice (e.g., 45 MPa environ pour des cas PA6,6 fibres de verre [1] ou 18.2
MPA pour des cas PP/fibres de carbone [3]), mais aussi du diamètre des fibres et de leur densité dans
le composite [8][9].

Les tenseurs des contraintes des fibres dans le repère global sont finalement donnés par la relation
(12).
$$ \sigma_F^\alpha = T^\alpha \begin{bmatrix} \sigma_F^{\alpha 0} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} T^{\alpha-1} \quad \forall \alpha \qquad (12) $$

## 2.3. Comportement macroscopique du composite

Le potentiel thermodynamique du matériau composite, $\Phi$, est défini comme la somme des potentiels des différents milieux (13), de densités constantes $\rho, \rho_M$ et $\rho_F^\alpha$, pondérée de leur fraction volumique. La dérivée temporelle du potentiel thermodynamique de la matrice peut être exprimée par la forme aux dérivées partielles (14). La seule variable d'état de chaque famille de fibres est $\varepsilon_F^{\alpha 0}$. Toutefois, on peut choisir de considérer la déformation tensorielle $\varepsilon_F^\alpha$, exprimée dans le repère global, pour plus de commodité. La dérivée temporelle du potentiel thermodynamique des fibres est alors donnée par $\frac{d\Phi_F^\alpha}{dt} = \frac{\partial\Phi_F^\alpha}{\partial\varepsilon_F^\alpha} : \dot{\varepsilon}_F^\alpha, \forall \alpha$. En petites déformations, les déformations de Hencky, $\varepsilon_F^\alpha$, peuvent être assimilées aux déformations de Green, $E_F^\alpha = \frac{1}{2}(C_F^\alpha - I)$. On démontre alors rapidement que $\dot{\varepsilon}_F^\alpha \approx A^\alpha \dot{\varepsilon} A^\alpha, \forall \alpha$. Enfin, en remarquant que $X : A^\alpha Y A^\alpha = A^\alpha X A^\alpha : Y$, pour tous tenseurs X et Y, on aboutit à la relation (15).

$$ \rho\Phi = v_M\rho_M\Phi_M + \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha\rho_F^\alpha\Phi_F^\alpha \qquad (13) $$

$$ \frac{d\Phi_M}{dt} = \frac{\partial\Phi_M}{\partial\varepsilon} : \dot{\varepsilon} + \frac{\partial\Phi_M}{\partial\varepsilon^p} : \dot{\varepsilon}^p + \frac{\partial\Phi_M}{\partial r} \cdot \dot{r} + \frac{\partial\Phi_M}{\partial D_n} \cdot \dot{D}_n + \frac{\partial\Phi_M}{\partial D_s} \cdot \dot{D}_s \qquad (14) $$

$$ \frac{d\Phi_F^\alpha}{dt} = A^\alpha \frac{\partial\Phi_F^\alpha}{\partial\varepsilon_F^\alpha} A^\alpha : \dot{\varepsilon}, \forall \alpha \qquad (15) $$

Les expressions (13) à (15) sont injectées dans l'inégalité de Clausius-Duhem, formulée pour une transformation isotherme, pour aboutir à l'inégalité (16), où $\sigma$ est le tenseur des contraintes du matériau composite. L'inégalité (16) doit être vérifiée pour toute valeur de $\dot{\varepsilon}$. Le système (S) (17) est donc une solution admissible. La considération des lois d'état $\rho_M \frac{\partial\Phi_M}{\partial\varepsilon} = \sigma_M$ et $\rho_F^\alpha \frac{\partial\Phi_F^\alpha}{\partial\varepsilon_F^\alpha} = \sigma_F^\alpha, \forall \alpha$ permet ensuite d'exprimer l'état de contrainte du matériau composite par la relation (18).

$$ \left[ \sigma - v_M \rho_M \frac{\partial \Phi_M}{\partial \varepsilon} - \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha \rho_F^\alpha A^\alpha \frac{\partial \Phi_F^\alpha}{\partial \varepsilon_F^\alpha} A^\alpha \right] : \dot{\varepsilon} - v_M \rho_M \left[ \frac{\partial \Phi_M}{\partial \varepsilon^p} : \dot{\varepsilon}^p + \frac{\partial \Phi_M}{\partial r} \cdot \dot{r} + \frac{\partial \Phi_M}{\partial D_n} \cdot \dot{D}_n + \frac{\partial \Phi_M}{\partial D_s} \cdot \dot{D}_s \right] \ge 0 \qquad (16) $$

$$ (S) \left\{ \begin{array}{l} \sigma = v_M \rho_M \frac{\partial \Phi_M}{\partial \varepsilon} + \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha \rho_F^\alpha A^\alpha \frac{\partial \Phi_F^\alpha}{\partial \varepsilon_F^\alpha} A^\alpha \\ v_M \rho_M \left[ \frac{\partial \Phi_M}{\partial \varepsilon^p} : \dot{\varepsilon}^p + \frac{\partial \Phi_M}{\partial r} \cdot \dot{r} + \frac{\partial \Phi_M}{\partial D_n} \cdot \dot{D}_n + \frac{\partial \Phi_M}{\partial D_s} \cdot \dot{D}_s \right] \le 0 \end{array} \right. \qquad (17) $$

$$ \sigma = v_M \sigma_M + \sum_{\alpha=1}^{n_{\text{fam}}} v_F^\alpha A^\alpha \sigma_F^\alpha A^\alpha \qquad (18) $$
### 3. Application : composite PA6,6 et courtes fibres de verre à orientation distribuée

Le modèle présenté dans la section précédente est appliqué au cas d'un composite à matrice polyamide 6,6 renforcée par des fibres de verre, formé par injection. Les résultats expérimentaux sont issus de la bibliographie [4]. Le matériau est soumis à des essais de traction uniaxiale, selon plusieurs directions de chargement caractérisées par un angle θ par rapport à la direction d'injection. Les fibres ont une longueur, supposée constante, de 0.24 mm, un rayon moyen de 5 µm et un module d'Young de 76 MPa. Elles présentent une orientation distribuée autour de la direction d'injection qui est modélisée grâce à une répartition des fibres en 10 familles, avec un angle d'orientation par rapport à la direction de chargement variant de 0 à 90°, par pas de 10°. Les fractions volumiques associées sont calculées à partir de la distribution d'orientation présentée par Kammoun et al. [4]. Le comportement élastoplastique de la matrice est modélisé en fonction des données fournis par ces auteurs [4], avec une limite initiale d'élasticité, σ₀, de 25 MPa (égale en traction et en compression). La loi d'écrouissage isotrope est linéaire-exponentielle et l'écoulement plastique est incompressible (γ=0). Enfin, le seuil d'endommagement est arbitrairement fixé à 0 et Sₙ=Sₛ=0.1 MPa, sauf mention contraire. L'IFSS dans le composite décroit à mesure que la teneur en fibres augmente. Dans l'exemple présenté ici, vF=10%, ce qui correspond à une IFSS de 45 MPa.

Le modèle est implémenté sous Abaqus Explicit (VUMAT). La simulation EF consiste en un unique élément cubique de 1mm subissant un déplacement imposé de 0.03 mm.

Les figures 1 à 3 montrent des exemples de comportements en traction du composite calculés pour des angles θ entre la direction de chargement et la direction d'injection (i.e. direction d'orientation préférentielle des fibres) de 0, 15 et 30°. Les résultats sont en très bon accord avec les données expérimentales, ainsi qu'avec les résultats de la procédure à double homogénéisation, proposée par Kammoun et al. [4], à l'exception d'une légère sous-estimation des niveaux de contrainte axiale pour le cas θ=30°. Cette première évaluation du modèle est donc tout à fait satisfaisante.

<figure><img src="image_2.png" /><figcaption>Fig. 1. θ=0°.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. θ=15°.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. θ=30°.</figcaption></figure>

L'influence d'une variation des paramètres d'endommagement de la matrice sur le comportement macroscopique du composite est illustrée par la Figure 4, pour le cas θ=0° (S=S_n=S_s). Dans cette configuration, les contraintes calculées pour les fibres sont très largement prépondérantes devant celles de la matrice en raison de la grande différence de rigidité des deux matériaux (76GPa pour les fibres de verre contre 3.1GPa pour le PA6,6). Pour des valeurs de S relativement élevées, l'endommagement de la matrice a alors une influence modérée sur le comportement du composite. En revanche, une diminution drastique de la valeur de S accélère fortement le développement de l'endommagement (Figure 5). L'adoucissement du composite est alors beaucoup plus marqué, d'autant plus que la transmission de l'effort de la matrice vers les fibres est très rapidement affectée (11). De plus, les niveaux de contrainte et de déformation à rupture (ici caractérisée par D_n=1) peuvent être fortement impactés. Sur la Figure 5(b), on peut constater un développement de l'endommagement dit « de cisaillement », c'est-à-dire affectant la partie déviatorique des contraintes σ_F^{α0} (cf section 2.1). Celui-ci est dû à l'orientation parfaite des fibres. En effet, dans le cas de fibres parfaitement alignées selon la direction d'injection (ici confondue avec la direction de chargement, puisque θ=0°), la variable D_s reste nulle.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Influence de l'endommagement de la matrice sur le comportement du composite (θ=0°).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Influence des paramètres d'endommagement sur les variables d'endommagement ($\theta=0°$).</figcaption></figure>

## 4. Conclusion

Ce papier présente une approche originale pour la modélisation du comportement de matériaux élastoplastiques endommageables renforcés par des fibres courtes présentant des orientations complexes. L'anisotropie de l'endommagement de la matrice est directement gouvernée par les caractéristiques du renforcement. L'endommagement affecte également l'effort transmis de la matrice vers les fibres.

La comparaison des résultats de simulation à des données expérimentales issues de la bibliographie démontre la pertinence de cette modélisation.

## Remerciements

Ces travaux ont bénéficié du support du CISIT (www.cisit.org), de la Région Nord-Pas-de-Calais, de la Communauté Européenne, de la Délégation Régionale pour la Recherche et la Technologie, du Ministère de l'Enseignement Supérieur et de la Recherche et du CNRS.

## Références

[1] W.H. Bowyer, M.G. Bader. On the re-inforcement of thermoplastics by imperfectly aligned discontinuous fibres, J. Materials Science, 7, 1315-1321, 1972.

[2] I. Doghri, L. Tinel. Micromechanics of inelastic composites with misaligned inclusions : numerical treatment of orientation, Comp. Meth. Appl. Mech. Engnrg, 195, 1387-1406, 2006.

[3] S.Y. Fu, B. Lauke, E. Mäder, C.Y. Hue, X. Hu. Tensile properties of short-glass-fiber- and short-carbon-fiber-reinforced polypropylene composites, Composites: Part A, 31, 1117-1125, 2000.

[4] S. Kammoun, I. Doghri, L. Adam, G. Robert, L. Delannay. First pseudo-grain failure model at finite elastic-plastic strains, Composites Part A, 42, 1892-1902, 2011.

[5] G.Z. Kang, C. Yang, J.X. Zhang. Tensile properties of randomly oriented short $\delta$-Al$_2$O$_3$ fiber reinforced aluminium alloy composites. I. Microstructure characteristics, fracture mechanisms and strength prediction, Composites: Part A, 33, 647-656, 2002.

[6] J. Lemaitre. A course on damage mechanics, 2nd Edition, Springer Berlin, 1996.

[7] J.C. Simo, T.J.R. Hughes. Computational Inelasticity, vol 7 Interdisciplinary Applied Mechanics, 2nd Edition, Springer Berlin, 1998.

[8] J.L. Thomason. Structure-property relationships in glass-reinforced polyamide, part 1: the effects of fibre content, Polymer Composites, 27, 552-562, 2006.

[9] J.L. Thomason. Structure-property relationships in glass-reinforced polyamide, part 2: The effects of average fiber diameter and diameter distribution, Polymer Composites, 28, 331-343, 2007.