# Intégration Numérique de lois constitutives pour le pli élémentaire

Philippe Jetteur

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

## To cite this version:

Philippe Jetteur. Intégration Numérique de lois constitutives pour le pli élémentaire. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01416074

HAL Id: hal-01416074

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

Submitted on 14 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Intégration Numérique de lois constitutives pour le pli élémentaire

Ph. Jetteur

SAMTECH

Rue des Chasseurs Ardennais 8, B-4031Liège, Belgique

Philippe.Jetteur@samcef.com

http://www.samcef.com

**Résumé** – Ce papier décrit l’intégration des lois constitutives composites. Le modèle retenu est un des modèles proposés par Cachan, il fait intervenir à la fois du dommage et de la plasticité. Le problème global est analysé de manière implicite. Au niveau d’un pas de calcul, l’accroissement de déformation peut donc devenir grand. La méthode choisie intègre en même temps l’endommagement et la plasticité et permet d’obtenir, au niveau du matériau, une loi tangente consistante avec le schéma d’intégration.

**Mots clés** – Matériau composite, intégration numérique, schéma implicite

## 1 Introduction

On ne rappelle ici ni l'intérêt des matériaux composites ni l'historique des différents modèles disponibles pour représenter leur comportement. Le but du papier n'est pas de présenter un nouveau modèle mais bien de montrer comment un de ces modèles peut être introduit dans le code de calcul implicite SAMCEF [5]. Les détails portent donc sur l'intégration numérique de la loi constitutive.

On commence par présenter le modèle retenu qui est dérivé du modèle initial de Ladevèze- Le Dantec [1]. Ensuite, on décrit l'intégration numérique ainsi que l'obtention d'une matrice tangente consistante. Quelques exemples sont présentés. On termine par les perspectives futures.

## 2 Modèle matériau

Comme dit plus haut, le modèle a comme point de départ le modèle initial de Ladevèze- Le Dantec [1]. Des enrichissements à ce modèle ont été présentés dans [2], ils concernent :

- la prise en compte d'un dommage associé à la rupture des fibres

- la possibilité de prendre en compte les contraintes transverses

- l'introduction d'un effet retard

Il s'agit d'un modèle choisi à un moment donné, il ne s'agit pas de la dernière version du modèle de Cachan.

### 2.1 Dommage

On définit un potentiel fonction des contraintes :
$$
\begin{align*}
e_d &= \frac{\sigma_{11}^2}{2(1-d_{11})E_1^0} - \frac{\nu_{12}^0}{E_1^0} \sigma_{11} \sigma_{22} - \frac{\nu_{13}^0}{E_1^0} \sigma_{11} \sigma_{33} + \\
    &\quad \frac{\langle \sigma_{22} \rangle_+^2}{2(1-d_{22})E_2^0} + \frac{\langle \sigma_{22} \rangle_-^2}{2E_2^0} + \frac{\langle \sigma_{33} \rangle_+^2}{2(1-d_{22})E_3^0} + \frac{\langle \sigma_{33} \rangle_-^2}{2E_3^0} - \frac{\nu_{23}^0}{E_2^0} \sigma_{22} \sigma_{33} + \\
    &\quad \frac{\sigma_{12}^2}{2(1-d_{12})G_{12}^0} + \frac{\sigma_{13}^2}{2(1-d_{12})G_{13}^0} + \frac{\sigma_{23}^2}{2(1-d_{23})G_{23}^0}
\end{align*}
$$

Les déformations élastiques sont égales à la dérivée de ce potentiel par rapport aux contraintes. En pratique, on prend $d_{23} = d_{22}$. On définit ensuite les forces thermodynamiques qui vont gouverner l'évolution du dommage. Dans la formulation actuelle, on prend des valeurs locales, on ne fait pas de moyenne sur l'épaisseur du pli. Si nécessaire, il est possible de ne prendre qu'un point d'intégration sur l'épaisseur du pli. On a par exemple pour la force thermodynamique associée au dommage $d_{22}$ :

$$Y_{22} = \frac{\partial e_d}{\partial d_{22}} = \frac{\langle \sigma_{22} \rangle_+^2}{2(1-d_{22})^2 E_2^0} + \frac{\langle \sigma_{33} \rangle_+^2}{2(1-d_{22})^2 E_3^0} + \frac{\sigma_{23}^2}{2(1-d_{22})^2 G_{23}^0}$$

On introduit une force équivalente : $Y = \max_{\tau \le t} [Y_{12} + b_2 Y_{22}]$. On définit un dommage « statique » en fonction des forces thermodynamiques :

$$d_{11}^s = \begin{cases} 0 & \text{si } Y_{11} \le Y_{11}^+ \\ 1 & \text{si } Y_{11} > Y_{11}^+ \text{, } \sigma_{11} > 0 \text{ ou } Y_{11} > Y_{11}^- \end{cases} \quad \text{et} \quad Y_{11} \le Y_{11}^-$$

$$d_{12}^s = \begin{cases} \frac{\sqrt{Y} - \sqrt{Y_{12}^0}}{\sqrt{Y_{12}^c} - \sqrt{Y_{12}^0}} & \text{si } d_{12} < 1, d_{22} < 1, Y_{12} < Y_{12}^F \\ 1 & \text{sin on} \end{cases}$$

$$d_{22}^* = \begin{cases} b_3 d_{12} & \text{si } d_{12} < 1, d_{22} < 1, Y_{22} < Y_{22}^F \\ 1 & \text{sin on} \end{cases}$$

On remarque une évolution d'abord progressive du dommage transverse (endommagement diffus) suivi d'une augmentation brutale jusqu'à la valeur unitaire (endommagement fragile). Les « vrais » dommages vérifient les équations suivantes :

$$\dot{d}_{ij} = \frac{1}{\tau} \{1 - \exp[-a_c \langle d_{ij}^s - d_{ij} \rangle_+]\}$$

$\tau$ représente un temps caractéristique. Si la vitesse de chargement est lente, c'est à dire si le temps de chargement est grand par rapport à $\tau$, la différence entre le « vrai » dommage et le dommage « statique » est négligeable. Pour finir, dans la direction des fibres, un comportement élastique non-linéaire est introduit.

## 2.2 Plasticité

On introduit la notion de contrainte effective. Par exemple, pour la direction 22, on a:

$$\tilde{\sigma}_{22} = \frac{\langle \sigma_{22} \rangle_+}{1 - d_{22}} + \langle \sigma_{22} \rangle_-$$

Le critère de plasticité porte sur les contraintes effectives :

$$f(\tilde{\sigma}, p) = \sqrt{\tilde{\sigma}_{12}^2 + (\tilde{\sigma}_{13}^2 + \tilde{\sigma}_{23}^2) + a^2 (\tilde{\sigma}_{22}^2 + \tilde{\sigma}_{33}^2)} - R_0 - R(p) \le 0$$

$$R(p) = K p^\gamma$$

Les vitesses de déformations plastiques conjuguées aux contraintes effectives se calculent par :
$$
\dot{\tilde{\epsilon}}^p = \dot{p} \frac{\partial f}{\partial \tilde{\sigma}}
$$

Pour calculer les vitesses de déformations plastiques conjuguées aux contraintes, on
utilise : $Tr(\tilde{\sigma} \dot{\tilde{\epsilon}}^p) = Tr(\sigma \dot{\epsilon}^p)$

3 Intégration numérique

Le modèle est basé sur un potentiel en contraintes. Ceci permet de distinguer facilement les états de traction et de compression. Par contre, d'un point de vue numérique, ce qui est connu lors de l'intégration des lois constitutives, ce sont les déformations. L'intégration numérique n'est dès lors pas immédiate.

Une première solution consiste à calculer le dommage au pas de calcul *n*+1 en fonction des contraintes au pas *n*. Il s'agit d'une intégration explicite. L'avantage d'une telle solution est que le dommage est constant lors d'un pas de calcul et perturbe peu la convergence. L'inconvénient est qu'il faut limiter la taille des pas et avoir des critères de contrôle/erreur de l'intégration de la loi constitutive. Pour un code implicite, la limitation sur la taille des pas peut être gênante. Nous ne retenons pas cette solution.

Les équations à intégrer sont des équations différentielles, nous utilisons un schéma d'Euler complètement implicite (α = 1). Comme ces équations sont non-linéaires, un schéma local de Newton est utilisé. Il est possible de faire une double boucle sur les itérations. Dans la boucle intérieure, on résout les équations à plasticité (endommagement) constante. Une fois la convergence obtenue, la boucle extérieure met à jour la plasticité (le dommage). Nous préférons une solution où on met à jour à la fois la plasticité et le dommage. De cette manière, il est simple de trouver une loi constitutive tangente.

Il faut définir les inconnues à mettre dans le système de Newton-Raphson, le choix n'est pas unique. On peut éliminer analytiquement une variable en fonction des autres pour avoir un système plus petit à résoudre ou, au contraire, on peut ajouter des variables pour essayer de rendre le système moins non linéaire. La convergence du Newton-Raphson dépendra du choix de ces inconnues et de la manière dont les équations sont écrites.

En pratique, on prend comme inconnues les contraintes effectives et la variation de déformation plastique équivalente sur le pas. On exprime pendant les itérations le dommage en fonction des contraintes et le tenseur de déformations plastiques en fonction des inconnues. Les 7 équations qui sont résolues sont :

$$
\begin{align*}
\underline{R}_{\sigma} &= \underline{\epsilon}_{el} + \underline{\epsilon}_{pl} + \underline{\epsilon}_{therm} - \underline{\epsilon} = \frac{\partial e_d}{\partial \underline{\sigma}} + \Delta \underline{\epsilon}_{pl} + \underline{\epsilon}_{pl}^{n} + \underline{\epsilon}_{therm} - \underline{\epsilon} = 0 \\
R_p &= f(\tilde{\sigma}, p) = 0
\end{align*}
$$

Pendant les itérations, on calcule une matrice tangente :

$$
\left|
\begin{array}{cc}
\displaystyle \frac{\partial R_{\sigma}}{\partial \tilde{\sigma}} & \displaystyle \frac{\partial R_{\sigma}}{\partial p} \\
\displaystyle \frac{\partial R_{p}}{\partial \tilde{\sigma}} & \displaystyle \frac{\partial R_{p}}{\partial p}
\end{array}
\right|
$$

A la convergence du processus de Newton, on peut exprimer les variations des inconnues en fonction de la variation des déformations. On utilise :

$$
\left| \frac{d\tilde{\sigma}}{dp} \right| = - \frac{\left| \frac{\partial R_{\sigma}}{\partial \tilde{\sigma}} \quad \frac{\partial R_{\sigma}}{\partial p} \right|^{-1} \left| \frac{\partial R_{\sigma}}{\partial \varepsilon} \quad \frac{\partial R_{\sigma}}{\partial \varepsilon} \right|}{\left| \frac{\partial R_{p}}{\partial \tilde{\sigma}} \quad \frac{\partial R_{p}}{\partial p} \right|}
$$

La contrainte s'exprime en fonction des contraintes effectives. On peut en déduire la matrice tangente, consistante avec le schéma d'intégration :
$$ 
\begin{gathered}
\underline{\sigma} = S(\tilde{\sigma}, p, \underline{\varepsilon}) \\
d\underline{\sigma} = \frac{\partial S}{\partial \tilde{\sigma}} d\tilde{\sigma} + \frac{\partial S}{\partial p} d p + \frac{\partial S}{\partial \underline{\varepsilon}} d\underline{\varepsilon}
\end{gathered}
$$

Avant d'activer ce mécanisme de Newton-Raphson, il y a une phase de prédiction élastique des contraintes, calculées avec le dommage en début de pas. Si toutes les équations sont vérifiées, le pas est élastique.

La procédure d'intégration correspond à une procédure de type prédicteur élastique, correcteur plastique et dommage.

D'un point de vue pratique, pour éviter les problèmes numériques, on introduit une valeur limite du dommage légèrement inférieure à 1. Pour pouvoir intégrer les équations au cours du temps, il faut garder des variables d'état. On prend comme variables d'état :

- tenseur déformation plastique

- déformation plastique équivalente $p$

- dommage

- sup $Y$

Le principe de résolution se résume à :

Boucle sur les itérations ( locales au point de gauss )

Entrée : déformation et estimation de $\tilde{\sigma} p$ à l'itération en cours

Calcul du dommage, avec ou sans effet retard, selon les données

Calcul des 7 équations et de leurs dérivées par rapport aux inconnues

Résolution

Mise à jour de $\tilde{\sigma} p$

Fin de la boucle

Calcul des contraintes et de la matrice tangente consistante

Le modèle possède une forte discontinuité quand $Y$ atteint $Y^F$. Le dommage subit un saut. Avec effet retard et un petit pas de temps, ce saut peut être limité. Si pendant les itérations locales, on s'aperçoit que l'on cycle entre le dommage limite et une valeur finie, on arrête les itérations, on fixe la composante du dommage qui oscille à sa valeur limite, on reprend les itérations avec la valeur limite du dommage jusqu'à atteindre la convergence.

# 4 Exemples

## 4.1 Plaque trouée

Le premier exemple concerne une plaque trouée. Un seul pli est utilisé, les fibres sont parallèles au grand coté. Différents maillages sont utilisés pour tester la sensibilité au maillage, le raffinement se fait en surface et dans l'épaisseur. La figure 1 montre le logarithme du dommage dans la direction perpendiculaire aux fibres en fin de chargement. Une valeur nulle correspond à un dommage unitaire. Sans effet retard, la convergence n'est pas atteinte, le dommage fragile pose des problèmes. Avec effet retard, la convergence ne pose pas trop de problème, le relâchement de la contrainte lors de l'apparition de l'endommagement fragile peut se faire progressivement, sur quelques pas de temps. Les zones endommagées ne sont pas trop sensibles au maillage.
<figure><img src="image_3.png" /><figcaption>Figure 1 – distribution du logarithme du dommage dans la direction transverse</figcaption></figure>

## 4.2 Voilement d'un panneau cylindrique

Le deuxième exemple correspond à un panneau cylindrique raidi soumis à une compression, c'est un problème traité lors du projet européen COCOMAT [4]. Des coques multicouches sont utilisées. En plus du dommage dans les plis, on prend en compte le décollement entre la peau et les raidisseurs en introduisant des éléments cohésifs, la loi utilisée dans ces éléments est celle décrite dans [3]. La figure 2 montre la distribution du dommage maximum dans les coques en fin d'analyse, le maximum correspond au maximum sur les plis et composantes. On remarque que les raidisseurs sont complètement décollés au centre du panneau en fin d'analyse. La partie droite de la figure 2 montre l'évolution de la charge axiale en fonction du déplacement axial. Les deux courbes correspondent à deux valeurs de la constante de temps de l'effet retard pour le pli (0.002 et 0.02 pour le pli, 0.001 pour l'interface). La sensibilité est faible, seule la décharge est affectée. Dans ce type de calcul, l'avantage de l'effet retard est de stabiliser le calcul et de faciliter la convergence. La valeur utilisée pour la constante de temps (0.02 ou 0.002) est petite vis-à-vis du temps de chargement (4) mais est grande par rapport aux valeurs utilisées en dynamique explicite.

# 5 Conclusion et perspective

L'implantation faite correspond à un modèle local. Les premiers résultats sont jugés encourageants, notamment en ce qui concerne la convergence en analyse implicite. Le travail n'est cependant pas terminé.

Dans [1], le dommage est supposé constant sur l'épaisseur du pli. La première amélioration consiste à prendre en compte le caractère moyen du dommage. Si on suppose un élément sur l'épaisseur du pli, cela implique de traiter en même temps tous les points d'intégration se situant sur l'épaisseur. En particulier, cela signifie que l'état de contrainte au point I dépend des déformations en I mais aussi en J (et de variables d'état). Dans la matrice tangente du matériau, la variation de contrainte au point I dépendra de la variation de déformation au point J. Ce sont des termes qui ne sont généralement pas pris en compte dans le développement des éléments. Cela signifie, que si on veut avoir une matrice tangente au niveau de l'élément, il faut revoir
aussi l'élément. Il ne suffit pas de modifier l'intégration de la loi constitutive.

Le modèle intégré correspond à un modèle à un moment donné, une version actualisée du modèle sera introduite très prochainement.

Comme signalé dans le deuxième exemple, un modèle pour le pli et un modèle pour l'interface sont disponibles. Pour le moment, ces modèles sont découplés. Il faudrait pouvoir les coupler. Le challenge pour obtenir une matrice tangente dans ce cas paraît grand.

Les tests présentés ont été pilotés en déplacement. Ceci permet de voir ce qui se passe après le maximum de la charge. Selon le type de chargement, ce type de pilotage n'est pas toujours possible. Des travaux ont été publiés pour adapter les méthodes de continuation de type Riks à ce type de problématique, en remplaçant l'équation supplémentaire globale par une équation plus locale. Il faudra voir si ces méthodes sont générales et non spécifiques à un type de comportement. En particulier quand on voit l'apparition « brutale » du dommage fragile dans le modèle traité.

La sensibilité au maillage devra être approfondie.

<figure><img src="image_4.png" /><figcaption>Figure 2 – panneau cylindrique</figcaption></figure>

Références

[1] Ladeveze P., Le Dantec E. : Damage modelling of the elementary ply for laminated composites, *Composite Science and Technology*, 43, pp.257-267, 1992.

[2] Gallucio A.C. : Modèle de comportement du pli – implémentation dans SAMCEF, Technical report 2007-27158/1-CTO/IW/SP/ME, EADS IW, 2007

[3] Allix O., Ladevèze P.: Interlaminar interface modelling for the prediction of laminate delamination, *Comp. Structures*, 22, pp.235-242, 1992.

[4] http://www.cocomat.de/

[5] http://www.samcef.com