# Simulation de grandes déformations de renforts de composites à l'aide de lois de comportement hypoélastiques

Pierre Badel, Emmanuelle Vidal-Salle, Philippe Boisse

► To cite this version:

Pierre Badel, Emmanuelle Vidal-Salle, Philippe Boisse. Simulation de grandes déformations de renforts de composites à l'aide de lois de comportement hypoélastiques. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01504293

HAL Id: hal-01504293

https://hal.science/hal-01504293

Submitted on 9 Apr 2017

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.

Public Domain
Simulation de grandes déformations de
renforts de composites à l'aide de lois de
comportement hypoélastiques

Pierre Badel*, Emmanuelle Vidal-Sallé*, Philippe Boisse*

LaMCoS, CNRS UMR5259
18-20, rue des Sciences
F69621 Villeurbanne Cedex
Philippe.Boisse@insa-lyon.fr

RÉSUMÉ. Les lois hypoélastiques (ou lois en taux) sont des modélisations continues fréquemment utilisées pour les grandes transformations. Dans le cas des milieux fibreux, la plupart des approches utilisent les dérivées objectives classiques de Green Naghdi ou Jaumann et le repère tourné associé. La matrice de comportement dans ce repère est obtenue à partir de celle dans le repère des fibres. Une alternative consiste à définir et utiliser une dérivée objective et un repère tourné à partir de la rotation de la fibre. L'objectif de ce travail est de montrer que seule la seconde approche conduit à des résultats satisfaisants dans certains cas. Des exemples d'application à la déformation de mailles élémentaires en tension biaxiale et en cisaillement seront présentés.

MOTS-CLÉS : Composites, Grandes transformations, Matériau fibreux, Loi de comportement hypoélastique, Dérivée objective, Repère tourné

KEYWORDS: Composites, Finite strain, Fibrous material, Rate constitutive equations, Objective derivative, Rotated frame
# 1. Introduction

Les milieux composés de fibres (de verre, de carbone, d'aramide en particulier) sont utilisés dans des applications à fortes sollicitations mécaniques où la masse doit être minimale. Ils constituent notamment une grande partie des renforts des matériaux composites. L'analyse de la mise en forme des renforts de matériaux composites nécessite la modélisation du comportement mécanique de ces milieux fibreux en grandes transformations. Celui-ci est spécifique car des glissements entre les fibres sont possibles. Afin de simplifier la présentation, le milieu fibreux étudié ici ne présentera qu'une seule direction de fibre. C'est le cas d'un pli unidirectionnel ou encore d'une analyse mésoscopique à l'échelle des mèches constituées de milliers de fibres. Pour la simulation de la mise en forme (drapage) les approches continues sont les plus utilisées actuellement en raison du grand nombre de fibres.

Une grande partie des logiciels s'appuie sur des lois de comportement de type hypoélastique (ou loi en taux) :

$$ \underline{\underline{\sigma}}^{\nabla} = \underline{\underline{C}} : \underline{\underline{D}} \quad [1] $$

g et D sont respectivement le tenseur des contraintes de Cauchy et le tenseur des taux de déformations. Ces deux tenseurs duaux sont classiquement admis comme variables eulériennes représentant la vitesse de déformation et la contrainte dans un milieu continu. g$^{\nabla}$ est une dérivée objective (ici en rotation) de g.

$$ \underline{\underline{\sigma}}^{\nabla} = \underline{\underline{Q}} \cdot \left( \frac{d}{dt} \left( \underline{\underline{Q}}^{\mathrm{T}} \cdot \underline{\underline{\sigma}} \cdot \underline{\underline{Q}} \right) \right) \cdot \underline{\underline{Q}}^{\mathrm{T}} = \dot{\underline{\underline{\sigma}}} + \underline{\underline{\sigma}} \cdot \underline{\underline{\Omega}} + \underline{\underline{\Omega}} \cdot \underline{\underline{\sigma}} \quad [2] $$

où Q est une rotation qui suit au mieux la matière et Ω le spin correspondant défini par $\underline{\underline{\Omega}} = \dot{\underline{\underline{Q}}} \cdot \underline{\underline{Q}}^{\mathrm{T}}$. Les logiciels de simulation utilisent classiquement la rotation de la décomposition polaire ou celle du repère corotationnel (où le taux de rotation de la matière est nul à chaque instant). Dans le premier cas la dérivée objective est celle de Green Naghdi, dans le second celle de Jaumann. Cette dérivée objective peut être vue comme la vitesse pour un observateur lié au repère tourné par Q. On appelle ce repère « repère tourné ». Cette approche est très utilisée dans les codes de calculs en grandes transformations.

Pour un milieu fibreux le tenseur de comportement est orienté par **f** la direction de la fibre qui est une direction matérielle. **f** n'est en général pas fixe dans le repère {eᵢ} formé des vecteurs de base initiaux tournés par **Q** car il s'agit d'une direction matérielle. La direction initiale de la fibre **f**⁰ est transformée par **E**, le tenseur gradient, et non pas par **Q** (figure 1). Pour résoudre ce problème, deux approches sont développées pour les milieux fibreux, dans le cadre hypoélastique présenté ci-
dessus. La première approche qui est la plus courante (Dong et. al., 2001, Peng et. al., 2005, Yu et. al., 2002) utilise la dérivée classique de Green Naghdi (ou de Jaumann) et effectue l'ensemble des calculs dans le repère tourné correspondant. La deuxième approche (Boisse et. al., 2005, Hagege, 2004) consiste à considérer une dérivée objective définie à partir de la rotation des fibres. L'objectif de la communication proposée est de montrer que seule la seconde approche est correcte pour la simulation des grandes déformations des milieux fibreux.

<figure><img src="image_1.png" /><figcaption>Figure 1. Repères de base {$e_i^0$}, tourné {$e_i$} et lié aux fibres {$f_i$}</figcaption></figure>

## 2. Ecriture incrémentale et calcul des contraintes

Dans les codes de simulation, le comportement mécanique joue un rôle majeur lors du calcul des contraintes. L'objectif est de calculer $\underline{\underline{\underline{g}}}^{n+1}$ connaissant $\underline{\underline{\underline{g}}}^n$ et les champs de déplacements à $t^n$ et $t^{n+1}$.

Notons d'un indice A les quantités non tournées telles que $\underline{\underline{\underline{S}}}_A=\underline{\underline{\underline{Q}}}^T \underline{\underline{\underline{S}}} \underline{\underline{\underline{Q}}}$ et $\underline{\underline{\underline{G}}}_A : \underline{\underline{\underline{S}}}_A = \underline{\underline{\underline{Q}}}^T (\underline{\underline{\underline{G}}}_A : \underline{\underline{\underline{S}}}_A) \underline{\underline{\underline{Q}}}$ pour $\underline{\underline{\underline{S}}}_A$ et $\underline{\underline{\underline{G}}}_A$ tenseurs d'ordre 2 et 4. $\underline{\underline{\underline{S}}}_A$ peut aussi être considéré comme $\underline{\underline{\underline{S}}}_A$ pour un observateur fixe dans le repère tourné. [1] devient :

$$ (\underline{\underline{\sigma}})^V_A = \underline{\underline{\underline{C}}}_A : \underline{\underline{\underline{D}}}_A \quad \text{ou} \quad \underline{\underline{\underline{\sigma}}}_A = \underline{\underline{\underline{C}}}_A : \underline{\underline{\underline{D}}}_A \qquad [3] $$

Par intégration de [3], et en notant que les composantes de $\underline{\underline{\underline{S}}}_A$ dans {$\underline{\underline{\underline{e}}}_i^0$} sont les mêmes que $\underline{\underline{\underline{S}}}_A$ dans le repère tourné {$\underline{\underline{\underline{e}}}_i$}, on obtient (Hughes et. al., 1980) :

$$ \begin{aligned} [\underline{\underline{\sigma}}^{n+1}]_{e_i^{n+1}} &= [\underline{\underline{\sigma}}^n]_{e_i^n} + [\underline{\underline{C}}^{n+1/2}]_{e_i^{n+1/2}} [\underline{\underline{D}}^{n+1/2}]_{e_i^{n+1/2}} \Delta t \\ [\underline{\underline{\sigma}}^{n+1}]_{e_i^{n+1}} &= [\underline{\underline{\sigma}}^n]_{e_i^n} + [\underline{\underline{C}}^{n+1/2}]_{e_i^{n+1/2}} [\underline{\underline{\Delta\epsilon}}]_{e_i^{n+1/2}} \end{aligned} \qquad [4] $$
où $[\mathbf{S}]_{e_i}^n$ représente la matrice des composantes de $\mathbf{S}$ dans la base $\mathbf{e}_i \otimes \mathbf{e}_j \otimes \mathbf{e}_k$ au temps $t^n$. La notation de Voigt est utilisée : les composantes de tenseurs de second ordre symétriques sont organisées en colonne.

## 2.1 Approches par repères tournés classiques

Dans ces approches, la rotation $\mathbf{Q}$ du repère tourné et de la dérivée objective [2] est la rotation de la décomposition polaire du tenseur gradient $\mathbf{F}=\mathbf{R}.\mathbf{U}$ : $\mathbf{Q}=\mathbf{R}$ comme il est souvent le cas dans les codes éléments finis (une alternative fréquente est le repère corotationnel). Ces codes fournissent à chaque incrément les contraintes en début d'incrément $[\boldsymbol{\sigma}^n]_{e_i}^n$ et l'incrément de déformation $[\Delta\boldsymbol{\varepsilon}]_{e_i}^{n+1/2}$. Le seul calcul à effectuer pour utiliser [4] dans le cas d'un milieu fibreux est celui de la matrice de comportement. Cette dernière est connue dans $\{\mathbf{f}_i\}$ où $\mathbf{f}=\mathbf{f}_1$ est la direction des fibres. Ainsi faut il réaliser un changement de base entre $\{\mathbf{f}_i^{n+1/2}\}$ et $\{\mathbf{e}_i^{n+1/2}\}$.

Les deux repères sont supposés confondus à l'instant initial. $\{\mathbf{e}_i\}$ est tourné par $\mathbf{R}$ alors que $\mathbf{f}_1$ est obtenue à partir du tenseur gradient.

$$ \mathbf{e}_i^n = \mathbf{R}^n \cdot \mathbf{e}_i^0 \quad \text{et} \quad \mathbf{f}_1^n = \frac{\underline{\underline{\mathbf{F}}}^n \cdot \mathbf{e}_1^0}{\|\underline{\underline{\mathbf{F}}}^n \cdot \mathbf{e}_1^0\|} \qquad [5] $$

$\mathbf{R}$ et $\mathbf{F}$ sont généralement calculés dans les codes de calcul. La transformation par $\underline{\underline{\mathbf{F}}}$ du vecteur $\mathbf{e}_2^0$ est utilisée pour déterminer $\mathbf{f}_2$ orthogonal à $\mathbf{f}_1$ :

$$ \tilde{\mathbf{f}}_2 = \underline{\underline{\mathbf{F}}} \cdot \mathbf{e}_2^0 \quad \text{puis} \quad \mathbf{f}_2 = \frac{\tilde{\mathbf{f}}_2 - (\tilde{\mathbf{f}}_2 \cdot \mathbf{f}_1)}{\|\tilde{\mathbf{f}}_2 - (\tilde{\mathbf{f}}_2 \cdot \mathbf{f}_1)\|} \quad \text{et} \quad \mathbf{f}_3 = \mathbf{f}_1 \wedge \mathbf{f}_2 \qquad [6] $$

On note $\mathbf{Q}$ la rotation entre $\{\mathbf{e}_i\}$ et $\{\mathbf{f}_i\}$ (orienté par $\mathbf{f}_1$). $\mathbf{Q} = \mathbf{f}_1 \otimes \mathbf{e}_i = (\mathbf{f}_j \cdot \mathbf{e}_i) \mathbf{e}_i \otimes \mathbf{e}_j$.

La matrice de comportement peut alors être exprimée dans le repère tourné $\{\mathbf{e}_i\}$ et [4] permet d'actualiser les contraintes.

## 2.2 Approche par le repère orienté par la direction de fibres

Dans cette approche, le repère tourné est le repère $\{\mathbf{f}_i\}$ défini par la direction de la fibre $\mathbf{f}_1$. La rotation $\mathbf{Q}$ de l'équation [2] est celle de la fibre, notée $\mathbf{\Phi}$. Une nouvelle dérivée objective en rotation est ainsi définie. Cette dérivée est la dérivée temporelle pour un observateur lié au repère de la fibre. On peut montrer le
caractère objectif de cette dérivée. L'équation [4] de calcul incrémental des contraintes devient alors :

$$[\boldsymbol{\sigma}^{n+1}]_{f_i^{n+1}} = [\boldsymbol{\sigma}^n]_{f_i^n} + [\mathbf{C}^{n+1/2}]_{f_i^{n+1/2}} [\Delta\boldsymbol{\varepsilon}]_{f_i^{n+1/2}} \quad [7]$$

Le repère courant {$\mathbf{f}_i$} est calculé à partir de [5]-[6] et permet alors de déterminer
$$\mathbf{\Phi} = \mathbf{f}_i \otimes \mathbf{e}_i^0 = (\mathbf{f}_j)_i^0 \mathbf{e}_i^0 \otimes \mathbf{e}_j^0 = (\mathbf{f}_j \cdot \mathbf{e}_i^0) \mathbf{e}_i^0 \otimes \mathbf{e}_j^0.$$

Le principal intérêt de cette approche réside dans le fait que la matrice de comportement apparaît dans le repère des fibres {$\mathbf{f}_i^{n+1/2}$} où elle a sa forme spécifique aux matériaux fibreux dont les fibres sont dirigées par $\mathbf{f}_1$.

Compte tenu que le terme de contraintes initiales $[\boldsymbol{\sigma}^n]_{f_i^n}$ a été calculé au pas précédent, seul le terme $[\Delta\boldsymbol{\varepsilon}]_{f_i^{n+1/2}}$ est nécessaire pour utiliser [7]. L'incrément $\Delta\boldsymbol{\varepsilon}$ est fourni dans le repère de travail du code. Abaqus, par exemple, fournit $\Delta\boldsymbol{\varepsilon}$ dans le repère de Green Naghdi {$\mathbf{e}_i$}. Connaissant la rotation $\mathbf{\Phi}$ de la fibre, et connaissant le repère où est donné $\Delta\boldsymbol{\varepsilon}$, la matrice de changement de base entre {$\mathbf{e}_i$} et {$\mathbf{f}_i$} peut être calculée. Quand le repère de travail du code est celui de Green Naghdi, il s'agit du changement de base associé à la rotation $\mathbf{\Phi} \cdot \mathbf{R}^T$.

L'équation [7] peut alors être utilisée et l'on prendra soin de retourner $\mathbf{g}^{n+1}$ dans le repère du code de calcul par un changement de base inverse.

### 3. Tests et comparaison des approches

On présente ici les résultats de deux tests réalisés sur un cube de côté unitaire dont le matériau fibreux présente une seule direction de fibres orientée par $\mathbf{e}_1^0$ à l'instant initial. Le seul paramètre matériau non nul est le module de tension dans la direction des fibres : $E_1= 35400$ MPa. Pour la suite, RC désigne la première approche présentée, RF la seconde.

- Test de cisaillement simple à 45° :

RC : $\forall t, [\boldsymbol{\sigma}]_{f_i} = [0]$ ; RF : $\forall t, [\boldsymbol{\sigma}]_{f_i} = [0]$

- Test de cisaillement simple à 45° sous tension de 24.5E+03 N :

RC : $[\boldsymbol{\sigma}]_{f_i} = \begin{bmatrix} \boldsymbol{\sigma}_{11} \\ \boldsymbol{\sigma}_{22} \\ \boldsymbol{\sigma}_{12} \end{bmatrix} = \begin{bmatrix} 22.9\,\mathrm{E} + 03 \\ 1.6\,\mathrm{E} + 03 \\ -6\,\mathrm{E} + 03 \end{bmatrix}$ ; RF : $[\boldsymbol{\sigma}]_{f_i} = \begin{bmatrix} 24.5\,\mathrm{E} + 03 \\ 0 \\ 0 \end{bmatrix}$
## 4. Analyse et conclusion

Les tests précédents montrent clairement que l'approche RC basée sur l'utilisation du repère tourné de Green Naghdi pourtant largement utilisée conduit à des résultats erronés dans certains cas. Les résultats seraient équivalents avec la dérivée de Jaumann. On ne peut donc pas utiliser les repères tournés classiques pour la simulation des grandes déformations de matériaux fibreux et il convient de faire jouer un rôle majeur à la direction de la fibre dans l'approche hypoélastique.

Les erreurs mises en évidence proviennent du terme $[\sigma^n]_{e_i^n}$ dans l'équation [4].

Sous cisaillement et incrément de contrainte nul, ce terme reste constant alors que la fibre tourne dans le repère {$e_i$}. En fait, la rotation polaire qui définit le repère de calcul ne représente qu'une moyenne de la rotation matérielle. Ainsi le repère de calcul ne suit il pas les directions matérielles. L'approche RF ne présente pas cet inconvénient et permet d'éviter les problèmes évoqués.

Une application concerne les simulations de renforts de composites à échelle mésoscopique où un modèle du comportement des fils constitués de milliers de fibres est nécessaire.

<figure><img src="image_2.png" /><figcaption>**Figure 2.** *Simulation mésoscopique du cisaillement d'un taffetas*</figcaption></figure>

Boisse P., Gasser A., Hagege B., Billoet J.L., « Analysis of the mechanical behaviour of woven fibrous material using virtual tests at the unit cell level », *Int. Journal of Material Science*, vol. 40, 2005, p. 5955-5962

Dong L., Lekakou C. and Bader M.G., « Processing of Composites: Simulations of the Draping of Fabrics with Updated Material Behaviour Law », *Journal of Composite Materials*, vol. 35, n°2, 2001, p. 138-163

Hagège B., « Simulation du comportement mécanique des milieux fibreux en grandes transformations : application aux renforts tricotés », Thèse de doctorat, ENSAM Paris, 2004

Hughes T.J.R, Winget J., « Finite rotation effects in numerical integration of rate constitutive equations arising in large deformation analysis », *Int. J. Numer. Methods Eng.*, vol. 15, 1980, p.1862-1867

Peng X., Cao J., « A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics», *Composites: Part A* vol. 36, 2005, p. 859–874

Yu W.R., Pourboghrat F., Chung K., Zamploni M. and Kang T.J., « Non-orthogonal Constitutive Equation for Woven Fabric Reinforced Thermoplastic Composites », *Composites Part A*, vol. 33, 2002, p. 1095-1105