# Modélisation du crash et du thermo-estampage d'une pièce en composite à matrice thermoplastique

Mamadou Adboul Mbacke, Linh Tuan Nguyen, Patrick Rozycki

► To cite this version:

Mamadou Adboul Mbacke, Linh Tuan Nguyen, Patrick Rozycki. Modélisation du crash et du thermo-estampage d'une pièce en composite à matrice thermoplastique. Revue des composites et des matériaux avancés = Journal of Composite and Advanced Materials, 2017, 27 (3-4), pp.357-379. 10.3166/rcma.2017.00023. hal-02443158

HAL Id: hal-02443158

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

Submitted on 13 Jul 2021

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 du crash et du thermo-estampage d'une pièce en composite à matrice thermoplastique

Mamadou Abdoul Mbacke¹, Tuan-Linh NGuyen¹, Patrick Rozycki²

1. IRT Jules Verne, Chemin du Chaffault, 44340 Bouguenais, France {mamadou.mbacke, tuan-linh.nguyen}@irt-jules-verne.fr

2. GeM/École centrale de Nantes, 1, rue de la Noë, 44300 Nantes, France patrick.rozycki@ec-nantes.fr

**RÉSUMÉ.** Ces travaux sont dédiés à la simulation du procédé de thermo-estampage d'une pièce automobile en composite à matrice thermoplastique suivie de la simulation de la tenue mécanique de la pièce fabriquée. La simulation du procédé s'appuie sur un modèle thermo-visco-hyperélastique. La simulation de la tenue mécanique a nécessité d'abord le développement et l'implémentation d'une loi de comportement en crash sous forme de subroutine Abaqus de type VUMAT. Les différents travaux de simulation ont été précédés de campagnes expérimentales pour alimenter les modèles.

**ABSTRACT.** This works consist of the simulation of thermo-stamping process on a thermoplastic composite automotive part followed by the mechanical simulation of the manufactured part. The process simulation is based on a thermo-visco-hyperelastic model. The simulation of mechanical behavior first required the development then the implementation of crash behavior law in Abaqus subroutines VUMAT. The different simulation works have been preceded by experimental tests in order to provide input data for the models.

**MOTS CLÉS :** composite, thermo-estampage, loi de comportement, essais expérimentaux, crash.

**KEYWORDS:** composite, thermo-stamping, constitutive law, experimental tests, crash.

## Extended abstract

IRT Jules Verne, industrial members and academic laboratories have decided to work on R&T projects on TP composite materials, on forming processes and simulation of product / process.
Simulations of the stamping process of thermoplastic matrix composites have been carried out. By using Plasfib with a new thermo-visco-hyperelastic model adapted to thermoplastic composites (Guzman-Maldonado et al., 2015), main results of these simulations are the fiber orientations, shear angles and wrinkles in the manufactured part. Experimental test campaigns (bias-extension (Guzman-Maldonado et al., 2015) and cantilever (Liang et al., 2014)) to determine the material properties have been performed. The influence of various parameters has been demonstrated through sensitivity studies (numerical DOE) to better understand and control thermo-stamping process. Then, coupling between thermostamping and overmolding process could be also studied in order to obtain composite pieces in one-shot.

In conjunction with these process simulations, simulations of the mechanical behavior of a representative automotive part have been carried out. The objective is to work on development of crash behavior law, representative case of automotive design loads. A new damage model for crash behavior has been developed taking into account strain rate developed with GeM (Ecole Centrale de Nantes). This model is based on principles of irreversible thermodynamics and is inspired by P. ROZYCKI and S. MARGUET works (Rozycki, 2000 ; Marguet, 2007).

# 1. Introduction

Dans un objectif d'allègement des véhicules par l'introduction du composite, l'IRT Jules Verne, ses partenaires industriels et académiques cherchent à travers ces travaux à développer une méthodologie riche pour modéliser à la fois le procédé de mise en forme et le comportement mécanique en service d'un démonstrateur automobile. La pièce est conçue pour représenter un certain nombre de contraintes à la fois géométriques et comportementales d'une pièce structurelle. La figure 1 montre la forme de la pièce d'étude. Le matériau utilisé un composite tissé avec une résine thermoplastique.

<figure><img src="image_2.png" /><figcaption>Figure 1. Géométrie de la pièce</figcaption></figure>
Dans un premier temps, les simulations du procédé de mise en forme de cette pièce par thermo-estampage seront tout d'abord présentées. Dans un second temps, le modèle de comportement et la simulation de structure utilisant le modèle développé seront présentés.

## 2. Simulation des procédés de mise en forme

Avant leur utilisation pour assurer la tenue mécanique, les pièces composites sont mises en forme par des divers procédés. La simulation des procédés de mise en forme joue donc un rôle important. Dans le cas de notre pièce, elle a été fabriquée grâce au procédé de thermo-estampage. Dans ce contexte, des simulations numériques du thermo-estampage sont réalisées afin d'aider le développement de ce procédé.

### 2.1. Modèle thermo-visco-hyperélastique

Les simulations du procédé de thermo-estampage ont été réalisées à l'aide du logiciel Plasfib. Plasfib est un logiciel de simulation du drapage de renforts de composites tissés et de préimprégnés développé par le Lamcos-INSA de Lyon, un partenaire de l'IRT JV. Il permet des simulations en fonction de la géométrie, du comportement des renforts textiles secs ou préimprégnés et pour des conditions aux limites données.

Pour les matériaux composites thermoplastiques tissés, Plasfib dispose d'une loi de comportement thermo-visco-hyperélastique non-linéaire développé par GUZMAN et al. (Guzman-Maldonado et al., 2015). Les détails du modèle sont amplement présentés dans l'article de GUZMAN et al. (Guzman-Maldonado et al., 2015). Pour décrire le comportement de renfort textile 2D, quatre modes de déformation indépendants sont pris en compte :

- élongation dans le sens chaîne ;

- élongation dans le sens trame ;

- cisaillement dans le plan ;

- flexion.

Basée sur une approche continue, le modèle hyperélastique s'appuie donc sur un découplage entre le comportement de membrane et le comportement en flexion. La fonction de l'énergie potentielle w peut s'écrire donc sous la forme :

$$w = w_{\text{membrane}} + w_{\text{flexion}} = w_1(I_1) + w_2(I_2) + w_{12}(I_{12}) + w_f(I_{f1}, I_{f2}) \quad (1)$$

où $w_1(w_2)$ est l'énergie potentielle de déformation associée à l'élongation dans le sens chaîne (trame), $w_{12}$ celle associée au cisaillement dans le plan et $w_f$ celle associée à la flexion. Ces énergies potentielles sont exprimées à travers les invariants physiques $I_1, I_2, I_{12}, I_{f1}, I_{f2}$.
Le tenseur de contrainte $\underline{S}$ (second tenseur de Piola-Kirchoff) est ensuite calculé :

$$ \underline{S} = 2 \frac{\partial w}{\partial \underline{C}} \qquad (2) $$

où $\underline{C}$ est le tenseur de déformation de Cauchy-Green droit.

Le comportement visqueux est lié uniquement au cisaillement dans le plan. L'énergie potentielle de déformation associée au cisaillement $w_{12}$ dans le plan étant supposée sous la forme d'un polynôme.

$$ w_{12}(I_{12}) = \sum_{k=1}^{n} c_k (I_{12})^k \qquad (3) $$

En se basant sur le modèle de Maxwell généralisé où la contrainte totale est la somme des contraintes de chaque élément du système amortisseur-ressort, le tenseur de contrainte lié à la déformation en cisaillement dans le plan $\underline{S}_{12}$ peut s'écrire sous la forme :

$$ \underline{s}_{12} = 2 \frac{\partial w_{12}}{\partial \underline{C}} - \sum_{i=1}^{N} \underline{Q}_{i} \qquad (4) $$

où $\underline{Q}_i$ représente un jeu de variables internes.

Ce tenseur peut être représenté sous la forme d'intégrale de convolution suivante :

$$ \underline{s}_{12} = \int_{-\infty}^{t} g(t-s) \frac{d}{ds} \left[ 2 \frac{\partial w_{12}}{\partial \underline{C}} \right] ds \qquad (5) $$

où $g(t)$ est la fonction de relaxation.

Elle s'écrit sous la forme :

$$ g(t) = \gamma_{\infty} + \sum_{i=1}^{N} \gamma_i \exp\left(-\frac{t}{\tau_i}\right) \qquad (6) $$

où $\tau_i > 0$ représentent les temps de relaxation (liés à la viscosité $\eta_i$ et modules d'Young $E_i : \tau_i \frac{\eta_i}{E_i}$) et $\gamma_i$ sont les paramètres matériaux adimensionnés liés au module d'Young $\gamma_i = E_i/E_0 (E_0 = E_{\infty} + \sum_{i=1}^{N} \underline{E}_i^j, \text{ donc } \gamma_i \in [0, 1])$.

Ces différents coefficients sont obtenus grâce à l'essai de bias-extension présenté ci-après (figure 4) pour une température donnée.

La thermo-dépendance est prise en compte par des coefficients d'interpolation qui modifient les propriétés matériaux en fonction de la température. Il est important de préciser que le modèle est valable uniquement dans la plage de température de mise en forme étudiée (supérieure à la température de fusion de la résine).
<figure><img src="image_3.png" /><figcaption>Figure 2. Élément triangulaire 123 et ses voisins (136, 235, 124) (Hamila et al., 2009)</figcaption></figure>

Compte-tenu de la faible épaisseur des renforts textiles et pré-imprégnés TP, les éléments de type coque ont été utilisés. Par souci d'efficacité numérique, des éléments finis spécifiques de coques sans degrés de liberté (ddl) en rotation ont été développés par Hamila *et al.* (2009).

Effectivement, les degrés de liberté en rotations servent à prendre en compte la rigidité en flexion. Elle est calculée dans ce cas à l'aide des courbures des éléments. La courbure des éléments est calculée à partir des positions et des déplacements des éléments voisins (V. figure 2). Comme il n'y a que des ddl en déplacements aux nœuds, le temps de calcul est réduit.

## 2.2. Campagnes de caractérisation

Afin d'alimenter cette loi de comportement, des essais de caractérisation des propriétés en cisaillement (bias-extension test) et en flexion (cantilever test) ont été réalisés à différentes températures afin de prendre en compte la thermo-dépendance.

Il est à remarquer que l'histoire thermique (cycle de chauffage-refroidissement du procédé du thermo estampage) doit être respectée pour ces essais.

### 2.2.1. Essai de bias-extension

L'essai de bias-extension (Guzman-Maldonado et al., 2015) permet de caractériser le comportement en cisaillement plan qui est le principal mécanisme de déformation du renfort dans le procédé thermo-estampage. Il s'agit d'un essai de traction sur une plaque à l'orientation de fibres +/-45° par rapport au chargement dans une enceinte thermique afin de pouvoir la tester à plusieurs températures. On impose un chargement sous forme de palier et on récupère l'effort généré.
<figure><img src="image_4.png" /><figcaption>Figure 3. Résultat d'un bias-extension test-courbe effort/temps</figcaption></figure>

La figure 3 montre un exemple de l'effort expérimental obtenu avec cet essai à une vitesse de déformation (Gassoumi et al., 2017). La relaxation de contrainte liée à la nature d'un matériau thermoplastique a été mise en évidence grâce à cet essai. L'essai permettra par la suite d'alimenter le comportement viscoélastique lié au cisaillement.

Une campagne de post-traitement de cet essai est ensuite réalisée afin de déterminer la carte matériau de Plasfib (les paramètres $τ_i, E_i, C_k$ évoqués précédemment) à travers un recalage entre l'effort théorique et l'effort expérimental (comme montre la figure 3).

Afin de valider la carte matériau obtenue, les simulations de cet essai de bias-extension test peuvent être réalisées. La figure 4 (à droite) montre un exemple du résultat de simulation où l'on voit les différentes zones de cisaillement (rouge, blanc, bleu). Ces résultats sont en bonne corrélation avec l'essai de bias-extension test réalisé par le Lamcos (à gauche).

### 2.2.2. Essai de cantilever

Normalement, les propriétés en flexion sont obtenues par les essais classiques comme la flexion 3 points. Néanmoins, dans le cas de la mise en forme des composites TP, ces essais ne sont pas adaptés car la rigidité en flexion est très faible (surtout en T°). C'est la raison pour laquelle le cantilever test est proposé dans les travaux de Liang et al. (2014). Il a pour but d'obtenir l'évolution du moment de flexion en fonction de la courbure à chaque température.

Le principe du cantilever test est la flexion de l'éprouvette due à son propre poids. La figure 5 montre un exemple de cet essai à deux températures différentes. On voit qu'à une température plus élevée, la flexion est plus importante.
<figure><img src="image_5.png" /><figcaption>Figure 4. Essai (à gauche) et Simulation (à droite) du bias-extension test</figcaption></figure>

Comme pour le bias-extension test, une campagne de post-traitement sera ensuite réalisée afin de déterminer la carte matériau de Plasfib.

Un exemple de simulation de validation de cet essai est présenté dans la figure 6. On voit également que la courbure est plus marquée à une température plus élevée (blanc).

Dans les deux cas, une bonne corrélation a été obtenue, permettant ainsi de valider la carte matériau pour simuler par la suite la pièce démonstrateur (figure 1).

<figure><img src="image_6.png" /><figcaption>Figure 5. Essai de cantilever à T° faible (à gauche) et T° élevée (à droite)</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Simulation de cantilever à T° faible (rouge) et T° élevée (blanc)</figcaption></figure>

## 2.3. Résultat de simulation

Dans ce paragraphe, des simulations du thermo-estampage sur pièce présentée précédemment (figure 1) sont présentées.

La figure 7 en haut montre l'angle de cisaillement et la forme finale de la pièce après le thermo-estampage. On identifie quatre zones de fort cisaillement (en rouge).

<figure><img src="image_8.png" /><figcaption>Figure 7. Angle de cisaillement-Pièce finale (haut) et Orientation de fibres : sens chaîne-violet et sens trame - blanc (bas)</figcaption></figure>
Ces simulations (résultats) serviront à placer de manière optimale les capteurs sur les zones d'intérêt pour le suivi du process. Quelques plissements ont été également constatés dans ces zones. Les bords de la pièce ont été modifiés suite à l'avalement de matière.

Il est à noter qu'en fonction des matériaux utilisés, les zones de cisaillements et de plissements sont différentes.

La figure 7 en bas montre l'orientation de fibres sens chaîne et sens trame dans une
zone de la pièce. Ces résultats peuvent être utilisés pour alimenter les logiciels de
simulation de structure par un transfert de données.

## 2.4. Etude de sensibilité et l'interaction entre paramètres procédé

L'interaction entre les paramètres procédés peut être mise en évidence grâce à un Plan d'expériences (PEX). Dans cette partie, un exemple d'étude de sensibilité entre la vitesse de mise en forme (v) et la température de la plaque d'organosheet (T) sera présenté à l'aide d'un PEX très simple: un PEX factoriel de 2 facteurs (v et T) à deux niveaux.

Une analyse de ce PEX sur l'angle de cisaillement maximal a été réalisée en utilisant
le logiciel Minitab©. La figure 8 montre l'évolution de la valeur moyenne de l'angle de
cisaillement max dans le domaine des facteurs (v et T) étudiés.

<figure><img src="image_9.png" /><figcaption>Figure 8. Graphique des effets principaux pour l'Angle de cisaillement maximal</figcaption></figure>
En tenant compte de ces analyses du PEX, on peut tirer les conclusions suivantes pour réduire l'angle de cisaillement max lors de la mise au point du procédé dans la fenêtre procédé étudiée :

- la température de mise en forme doit être la plus faible possible ;

- la vitesse de mise en forme doit être la plus importante possible.

## 2.5. Détermination de la forme développée des plaques

Lors de la conception, la forme finale de la pièce est définie afin de concevoir les moules. Par contre, la forme développée, c'est-à-dire la forme initiale de l'organosheet en vue de l'obtention de la forme finale après la mise en forme « net-shape », est souvent inconnue. De plus, la connaissance de la forme développée de la pièce est importante pour la conception et la mise au point du procédé de thermo-estampage tels que la disposition d'outils de maintien et la réduction de l'encombrement des moules.

La figure 9 montre le résultat de la forme développée obtenue pour notre pièce basée sur une méthode itérative : essais-erreurs numériques. L'idée est de remonter aux positions initiales des nœuds du contour de la pièce finale.

En utilisant cette forme développée comme la géométrie initiale de la plaque, la forme finale obtenue comme montre que la figure 10 coïncide avec la forme finale souhaitée (représentée ici par les points blancs).

Il est à remarquer que la robustesse de la simulation joue un rôle très important car elle influence la forme finale de la pièce. Il est évident également que les résultats obtenus dépendent des paramètres utilisés (conditions de mise en forme, matériaux...). C'est la raison pour laquelle la forme développée obtenue correspond uniquement à un jeu de conditions de mise en forme.

<figure><img src="image_10.png" /><figcaption>Figure 9. Forme développée de la pièce initiale</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. Forme finale obtenue utilisant la forme développée précédente</figcaption></figure>

## 2.6. Influence du champ thermique

Dans le procédé de thermo-estampage, la plaque composite sera en contact avec les moules lors de l'étape de mise en forme. Au cours du procédé, avec l'avancement du poinçon et la déformation de la plaque, ces zones de contact évoluent. Le transfert de chaleur entre la plaque et le moule sera donc influencé.

Toutefois, le couplage thermique (transitoire)-mécanique n'est pas encore disponible dans la version actuelle de Plasfib, le champ thermique de la plaque est donc considéré comme constant lors de la fermeture du moule (à noter que ce temps reste très court, dans l'ordre de quelques secondes). Afin d'analyser l'influence du champ de température de la plaque, un cas test utilisant un champ de température non homogène mais constant au cours du procédé seront présentés (figure 11). Toutefois, il est à noter que le problème de transfert thermique entre les moules et la plaque étant un

<figure><img src="image_12.png" /><figcaption>Figure 11. Champ thermique hétérogène testé</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12. Champ thermique homogène (gauche) et champ thermique hétérogène 2D (droite)</figcaption></figure>

problème très complexe (choc thermique, convection, cristallisation, gradient thermique...), un calcul thermique beaucoup plus fin devrait être réalisé si l'on veut obtenir un champ thermique exact de la plaque.

La figure 12 montre le résultat avec un champ thermique homogène (à gauche) et le champ thermique hétérogène testé. On voit que l'apparition de plissements varie en fonction du champ thermique utilisé. Néanmoins, cette différence n'est pas très marquée car l'hétérogénéité n'est pas très importante.

## 2.7. Vers un couplage thermo-estampage/surmoulage

Le procédé de thermo-estampage peut être couplé avec le procédé de surmoulage (injection sur l'insert estampé) afin d'obtenir les pièces one-shot. Il est donc intéressant de réaliser la simulation du procédé de surmoulage. L'objectif est d'étudier l'influence de la présence de l'insert (sa température, sa forme...) sur l'écoulement de surmoulage.
Inversement, cet écoulement peut également influencer l'insert (déformation, fusion de la résine sur les surfaces). Il s'agit donc d'une étude de l'interaction fluide-structure. L'échange de chaleur entre la résine et l'insert est également important à étudier.

Normalement, la simulation du procédé de surmoulage est réalisée en utilisant un insert idéal sans défauts issu de la conception. Toutefois, l'insert qu'on obtient à la suite d'un procédé de thermo-estampage peut parfois présenter des défauts comme une géométrie déformée, des plissements, une variation d'épaisseur... Il est donc intéressant d'étudier l'impact de ces défauts sur l'écoulement du surmoulage. On s'intéressera donc au couplage thermo-estampage/surmoulage. Cette thématique sera abordée dans les futurs travaux.

# 3. Modélisation du comportement mécanique de la structure

L'un des objectifs de ces travaux est le développement et l'implémentation d'une loi de comportement suffisamment représentative pour prédire la tenue mécanique de structures en composite tissé, particulièrement sollicitées en grande vitesse (crash).

Par la nature du matériau et par la nature des sollicitations, la loi de comportement devra permettre de modéliser l'endommagement progressif du matériau jusqu'à rupture et le comportement élastoplastique principalement gouverné par la nature de la résine et de prendre en compte la sensibilité éventuelle du matériau aux effets de vitesse.

## 3.1. Analyse expérimentale

Des essais mécaniques ont été réalisés à différentes conditions d'hygrométrie, de température et de vitesse pour comprendre le comportement macroscopique du matériau étudié. Les éprouvettes ont été conditionnées dans une enceinte climatique avec trois niveaux d'humidité relative (HR) : 0%, 50% et 85%. De même, trois niveaux de température ont été testés : -40 °C, 23 °C et 80 °C. Seule une partie des résultats expérimentaux est présentée dans cet article.

Il a été ainsi réalisé des essais de traction monotone sur des éprouvettes à 0° et 90°, découpées suivant les directions sens chaîne et sens trame pour identifier les modules d'élasticité et les grandeurs à rupture. Les résultats d'essais confirment l'équilibre du matériau dans le sens chaîne et le sens trame. On observe un comportement à rupture presque fragile. Les niveaux d'endommagement avant rupture sont très faibles. Les propriétés mécaniques restent très peu sensibles à l'effet d'hygrométrie et de température.

Parallèlement, des essais cycliques hors axe ont été réalisés sur des éprouvettes à 45° pour analyser le comportement en cisaillement. On constate une forte non-linéarité du matériau due à la fois à la plasticité et à l'endommagement. En effet, le comportement de la matrice thermoplastique contribue grandement au comportement global du matériau ; la matrice gouverne notamment la plasticité observée en cisaillement dans le plan 1-2. Au-delà du module de cisaillement et des grandeurs à rupture, ces essais cycliques
<figure><img src="image_14.png" /><figcaption>Figure 13. Influence de l'hygrométrie relative sur l'évolution des contraintes en fonction de la déformation</figcaption></figure>

permettent de relever l'évolution de l'endommagement et de définir la loi de plasticité correspondante.

On observe ainsi une forte influence de l'hygrométrie sur le comportement mécanique du matériau matérialisé par la baisse des rigidités et des résistances. Comme l'ont montré Launay *et al*. (2013) pour une matrice composite à matrice polyamide 66 (PA66) renforcée de fibres de verre courtes, la température et l'hygrométrie affecte le comportement mécanique du matériau à travers une modification de sa température de transition vitreuse (Tg) du matériau (figure 13).

Parallèlement, des essais dynamiques ont été réalisés pour analyser la sensibilité du matériau aux effets de vitesse. On observe une évolution des propriétés mécaniques en fonction de la vitesse de sollicitation. Les courbes ci-dessous présentent les résultats expérimentaux pour deux niveaux de vitesse.

Par ailleurs, on constate que l'influence de la vitesse de déformation est conditionnée par le couplage de la température et de l'hygrométrie, matérialisé soit par un
adoucissement soit par un durcissement du matériau. La compréhension de ces phénomènes observés nécessitera forcément l'analyse de l'évolution de la Tg en fonction des conditions d'essais.

## 3.2. Loi de comportement élastoplastique endommageable

### 3.2.1. Présentation du modèle implémenté

À la suite d'une étude bibliographique sur les modèles présents dans la littérature et sur les modèles disponibles dans les codes de calcul commerciaux et d'une analyse des premiers essais expérimentaux, il a été décidé de s'inspirer des travaux de Johnson et al. (2001) et des travaux de thèse de Rozycki (2000). Ce modèle a également été repris dans le cadre des travaux de thèse de Marguet (2007).

Il s'agit d'un modèle initialement proposé pour les tissus verre/époxy, dont l'approche adoptée pour décrire le comportement repose sur les nombreux travaux menés par Ladeveze et al. (Ladevèze, 1992 ; Ladevèze et al., 2000) sur les composites unidirectionnels. Le modèle de comportement proposé permet de décrire les comportements suivant les différentes directions du composite et est formulé en contraintes planes à l'échelle du pli, dans son repère d'orthotropie. La formulation théorique du modèle en quasi statique se divise en :

- un comportement élastique endommageable fragile dans le sens chaîne et le sens trame du renfort tissé du composite ;

- un comportement élastoplastique endommageable dans le plan 12 (cisaillement).

Le matériau composite étudié a un comportement équilibré dans les directions sens chaîne et sens trame. Plutôt que d'utiliser une formulation en contrainte comme les auteurs précédents, nous avons choisi d'utiliser une formulation en déformation, plus adaptée à la finalité de la loi dans un code de calcul explicite, comme le proposent les auteurs dans [3,4]. Cela se traduit par l'introduction, dans la matrice de rigidité $C̄$, de variables dites d'endommagement ($d_1$, $d_2$ et $d_{12}$) qui évoluent de 0 (matériau non endommagé) à 1 (matériau rompu).

$$C = \begin{pmatrix} C_{11}^0 (1 - d_1) & v_{12} C_{22}^0 & 0 \\ v_{21} C_{11}^0 & C_{22}^0 (1 - d_2) & 0 \\ 0 & 0 & G_{12}^0 (1 - d_{12}) \end{pmatrix} \quad (7)$$

$C_{11}^0$ et $C_{22}^0$ sont les coefficients de la matrice de rigidité du composite non endommagé. $G_{12}^0$ est le module de cisaillement du composite non endommagé.

Pour ce faire, les équations décrivant les évolutions de l'endommagement sont introduites et relient les paramètres de dégradation qui sont des variables d'état à l'énergie libérée au sein du pli tissé.
Ainsi, l'énergie de déformation élastique s'exprime comme suit :

$$2W_e^a = C_{11}^0(1-d_1)(\epsilon_{11}^e)^z + C_{22}^0(1-d_2)(\epsilon_{22}^e)^z + 2\nu_{21}^0 C_{11}^0 \epsilon_{11}^e \epsilon_{22}^e + G_{12}^0(1-d_{12})^2 \quad (8)$$

Les variables thermodynamiques sont déduites de cette énergie par dérivation par rapport aux variables d'endommagement :

$$\left\{ \begin{array}{l} \gamma_{11} = \frac{1}{2} c_{11}^0 (\epsilon_{11}^e)^2 \\ \gamma_{22} = \frac{1}{2} c_{22}^0 (\epsilon_{22}^e)^2 \\ \gamma_{12} = \frac{1}{2} G_{12}^0 (2\epsilon_{12}^e)^2 \end{array} \right. \qquad (9)$$

Les lois d'évolution des variables d'endommagement sont définies en fonction des variables thermodynamiques :

$$\left\{ \begin{array}{l} d_1 = f_{11}(\gamma_{11}, \gamma_{110}, \gamma_{11R}) \\ d_2 = f_{22}(\gamma_{22}, \gamma_{220}, \gamma_{22R}) \\ d_{12} = f_{12}(\gamma_{12}, \gamma_{120}, \gamma_{12R}) \end{array} \right. \qquad (10)$$

Où :

– les grandeurs $\gamma_{ij}(avec \ ij = 11, 22, 12)$ correspondent aux racines carrées des variables thermodynamiques tels que : $\gamma_{ij} = \max(\sqrt{\gamma_{ij}(t)}); \ ij = (11, 22, 12)$;

– les grandeurs $\gamma_{ij0}(avec \ ij = 11, 22, 12)$ et $\gamma_{ijR}(avec \ ij = 11, 22, 12)$ correspondent respectivement aux seuils inférieurs et supérieurs ;

– les fonctions $f_{ij}(avec \ ij = 11, 22, 12)$ sont des fonctions d'approximation adaptées aux observations expérimentales.

Quant au comportement en cisaillement, la réponse du matériau est, en grande partie, gouvernée par la nature de la matrice thermoplastique utilisée. Elle se traduit, dans le plan 12, par la présence prédominante d'une déformation non réversible en plus de la déformation élastique $\epsilon_{12}^f = \epsilon_{12}^e + \epsilon_{12}^p$. La modélisation du comportement anélastique du matériau est modélisée en introduisant une fonction de charge $f_p(\tilde{\sigma}_{12}, R)$, définie à l'aide de la variable d'endommagement $d_{12}$ et d'une loi d'écrouissage $R(p)$ de type isotrope tel que :

$$f_p(\tilde{\sigma}_{12}, R) = |\tilde{\sigma}_{12}| - R_{(p)} - R_0 \qquad (11)$$

$$R_{(p)} = Qp^8 \qquad (12)$$

$$\tilde{\sigma}_{12} = \frac{\sigma_{12}}{1 - d_{12}} \qquad (13)$$

$\sigma_{12}$ et $\tilde{\sigma}_{12}$ sont respectivement la contrainte de cisaillement dans le plan et la contrainte effective associée.
Les variables $p$ et $R_0$ sont respectivement la déformation plastique cumulée et la contrainte seuil de plasticité. $Q$ et $\beta$ sont les paramètres matériaux de la loi d'écrouissage. Enfin, une loi complémentaire d'écoulement y est associée en introduisant le multiplicateur de plasticité $\dot{\lambda}$ :

$$f_p(\bar{\sigma}_{12}, R) = \dot{f}_p(\bar{\sigma}_{12}, R) = 0 \qquad (14)$$

$$ \left\{ \begin{array}{l} \dot{e}_{12}^p = \dot{\lambda} \frac{\partial f_p}{\partial \sigma_{12}} \\ \\ \dot{p} = \dot{\lambda} \frac{\partial f_p}{\partial R} \end{array} \right. \quad \text{avec } \dot{\lambda} \ge 0 \qquad (15) $$

On suppose que l'incrément de déformation plastique est normal à la fonction de domaine élastique.

<figure><img src="image_15.png" /><figcaption>Figure 14. Influence de la vitesse de sollicitation sur l'évolution des contraintes en fonction de la déformation</figcaption></figure>
Cette loi de comportement est, par la suite, implémentée dans ABAQUS/Explicit
sous la forme d'une subroutine de type VUMAT. Ce choix est justifié par la perspective
de l'intégration des effets de vitesse dans le modèle. Le langage utilisé est fortran 77.

Les données d'entrée de la VUMAT sont :

- les propriétés élastiques du composite $E_1^0, E_2^0, v_{12}, v_{21}$ et $G_{12}^0$ ;

- les variables thermodynamiques seuils $\Upsilon_{ij0}(\text{avec } ij = 11, 22, 12)$ et $\Upsilon_{ijR}(\text{avec } ij = 11, 22, 12)$ ;

- la limite d'élasticité et les paramètres matériaux de la loi d'écrouissage $R_0$, $Q$ et $\beta$ ;

- la densité du matériau.

3.2.2. Validation sur un élément

Le programme implémenté, un essai de traction uniaxial suivant la direction sens chaîne (figure 15) et un essai de cisaillement (figure 16) ont été simulés pour valider l'implémentation.

On obtient une excellente corrélation entre modèle et essai pour le comportement
dans le sens trame et sens chaîne. Le modèle prédit un très faible taux
d'endommagement dans ces directions.

Quant au comportement en cisaillement, la corrélation des contraintes entre modèle
et essai est moins satisfaisante que celle obtenue pour le comportement dans le sens
trame et chaîne. Une des problématiques qui n'est pas prise en compte est la

<figure><img src="image_16.png" /><figcaption>Figure 15. Essai de traction sens chaîne (direction 1) : Courbe $\sigma_1$ - $\varepsilon_1$</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16. Essai de cisaillement : Courbe $\sigma_{12}$ - $\gamma_{12}$</figcaption></figure>

réorientation des fibres au sein de la microstructure. Cette réorientation des fibres sens trame et sens chaîne impose une évolution « non conventionnelle » et pourrait remettre en cause l'hypothèse d'orthotropie. Toutefois, on obtient des énergies de déformations presque équivalentes. En perspective, le modèle fait actuellement l'objet d'une étude d'optimisation des fonctions d'approximation. En effet, les solutions proposées et déjà implémentées pour améliorer la prédiction des niveaux de contraintes reposent sur une meilleure approximation de l'endommagement en utilisant deux fonctions couplées. Les premiers résultats sont très satisfaisants.

### 3.3. Intégration de l'effet de vitesse

#### 3.3.1. Présentation du modèle implémenté

Comme le montre les résultats d'essais (figure 14), le matériau composite à atrice thermoplastique étudié a un comportement sensible à la vitesse de déformation. En effet, certaines grandeurs physiques comme la contrainte et les modules d'élasticité croissent quand on augmente la vitesse de sollicitation $\dot{c}$ (figure 14).

Afin de décrire les effets de la vitesse de déformation, une extension du modèle au cas dynamique est réalisée. Pour ce faire, le calcul de la vitesse de déformation est défini et basé sur le second invariant du tenseur des déformations élastiques :

$$ \dot{c} = \frac{1}{\Delta t} \sqrt{\frac{4}{3} \left( (\varepsilon_{11}^e)^2 + (\varepsilon_{22}^e)^2 + (\varepsilon_{11}^e \varepsilon_{22}^e) + (\varepsilon_{12}^e)^2 \right)} \quad (16) $$

Des observations expérimentales ont permis d'identifier les paramètres du modèle qui sont sensibles aux effets de la vitesse. Ainsi, un paramétrage des modules d'élasticité, de la contrainte d'écoulement initiale et des taux de restitution d'énergie d'initiation et de rupture en fonction de la vitesse de déformation a été réalisé à la
manière d'un modèle de Johnson Cook (Johnson et Cook, 1985) :

$$ \gamma(\dot{e}) = \gamma(\dot{e}_{ref}) \left[ 1 + D_{\gamma} \ln \left( \frac{\dot{e}}{\dot{e}_{ref}} \right) + \ln(m\gamma) \right]^{n_{\gamma}} \quad (17) $$

- η représente n'importe lequel des grandeurs citées ci-dessus ;

- D, m, η, n, η, sont des coefficients propres à chaque η ;

- ε<sub>ref</sub> correspond à la vitesse de déformation de référence, généralement celle utilisée en quasi-statique.

### 3.3.2. Validation sur un élément

Des essais de cisaillement ont été simulés à différentes vitesses de sollicitation. Le modèle permet ainsi de reproduire l'évolution des propriétés mécaniques en fonction de la vitesse de sollicitation du matériau comme le montre la courbe ci-dessous (figure 17).

Une optimisation du modèle est menée afin de décroître les temps de calcul CPU et d'y intégrer la non-dépendance au maillage (Germain, 2006). Enfin, le modèle sera validé avec les essais dynamiques réalisés dans le projet.

## 3.4. Simulation du comportement de la pièce automobile

Après mise en forme de la pièce (figure 1), elle est percée au milieu pour introduire de nouvelles singularités géométrie afin de pouvoir ultérieurement évaluer la stabilité du modèle développé.

<figure><img src="image_18.png" /><figcaption>Figure 17. Résultats numériques - Évolution des contraintes en fonction de la vitesse de déformation</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18. Essais d'impact : flexion trois points (a) et compression (b)</figcaption></figure>

Le modèle développé et implémenté est utilisé pour simuler des prototypes automobiles. Ces simulations numériques seront, par la suite, corrélées avec des essais réalisés sur le banc multivérin de l'IRT Jules Verne et avec d'autres essais réalisés chez les partenaires du projet (RENAULT, PSA et FAURECIA).

Deux types de modèles ont été créés pour valider le modèle sur le démonstrateur automobile retenu. Il s'agit d'essais d'impact pour une sollicitation en flexion trois points et pour une sollicitation en compression.

Pour avoir un calcul de référence tant au niveau de la prédictibilité des résultats que du temps de temps CPU, ces deux essais ont été d'abord simulés en supposant un comportement élastique linéaire. Les résultats sont présentés dans les figures ci-dessous (figure 18).

Un critère de rupture de type Tsai Wu est utilisé en post-traitement des contraintes pour observer les zones où la rupture est plus probable. Cette approche permet, malheureusement, d'avoir juste une analyse macroscopique de la tenue mécanique de la structure.

Parallèlement, les premières simulations réalisées avec le modèle en cours d'optimisation permettent une analyse plus fine et donc plus prédictive. Elles fournissent des informations sur la nature des déformations observées, notamment les déformations élastiques et des déformations irréversibles. En même temps, le modèle permet de visualiser l'évolution d'endommagement et sa localisation au sein de la pièce.

## 4. Conclusion et perspectives

Dans un premier temps, le procédé de thermo-estampage de la pièce a été simulé en utilisant une nouvelle loi de comportement thermo-visco-hyperélastique développée par le Lamcos, une loi dont les paramètres ont été déterminés grâce à deux essais de caractérisation (essai de bias-extension et essai de cantilever). Ces calculs ont permis d'optimiser le procédé thermo-estampage : identification des zones à risques (fort cisaillement, plissements...), influence des paramètres réglables du process, détermination de la forme développée, impact de transfert de chaleur dans les
outillages,... Dans l'avenir, ces travaux pourront être couplés avec la simulation du surmoulage afin de développer le procédé thermo-estampage/surmoulage en vue de la fabrication des pièces composites en one-shot.

Dans un deuxième temps, la tenue mécanique de la même pièce a été analysée. Ainsi, il a été proposé et implémenté un modèle de comportement prenant en compte l'endommagement progressif du matériau, la déformation plastique et la sensibilité aux effets de vitesse du matériau composite étudié. Le choix du modèle est le fruit de l'analyse des essais expérimentaux réalisés sous différentes conditions de température et d'hygrométrie et sous différentes vitesses de sollicitation. La première phase de validation du modèle est satisfaisante. Cette version du modèle a été testée sur une structure réelle pour analyser sa tenue mécanique en simulant deux essais d'impact. La loi est en cours d'optimisation pour améliorer son niveau de prédicibilité et le temps de calcul CPU.

Enfin, une réflexion sur un chaînage numérique permettant le couplage de simulation produit-process est entamée. Avec ces futurs travaux, l'IRT Jules Verne proposera un outil de simulation permettant un dimensionnement optimal en intégrant les effets du process dans les calculs de structures en composites.

## Remerciements

Ce travail est réalisé dans le cadre des projets IRT Jules Verne :

- COMMANDO_STAMP (IRT JV, CEMCAT, PSA, SAFRAN, SOLVAY, 3SR, GeM, LaMCoS, LTN) ;

- COMPOSTAMP (IRT JV, RENAULT, PSA, FAURECIA, CETIM, CEMCAT, AIRBUS, DAHER, COMPOSE, DEDIENNE) ;

- COPERSIM (IRT JV, PSA, RENAULT, FAURECIA, SOLVAY, FLEX |N| GATE, Ecole Centrale de Nantes).

## Bibliographie

Gassoumi M., Rolland Du Roscat S., Casari P., Dumont P.J.J., Orgeas L., Jacquemin F. (2017). *Étude du comportement et des mésostructures de préimprégnés thermoplastiques à renforts tissés au cours du thermoestampage*. JNC 20, Champs-sur-Marne.

Germain N. (2006). *Modélisation non locale de l'endommagement dans les structures composites*. Thèse de doctorat, ENSMP – Centre des matériaux.

Guzman-Maldonado E., Hamila N., Boisse P., Bikard J. (2015). Thermomechanical analysis, modelling and simulation of the forming of pre-impregnated thermoplastics composites. *Composites: Part A*, vol. 78, p. 211-222.

Hamila N., Boisse P., Sabourin F., Brunet M. (2009). A semi-discrete shell finite element for textile composite reinforcement forming simulation. *International Journal of Numerical Methods in Engineering*, vol. 79, p. 1443-1466.
Johnson G.R., Cook W.H. (1985). Fracture characteristics of three metal subjected to various strains, strain rates, temperature and pressures. *Engineering Fracture Mechanics*, vol. 21, p. 31-48.

Johnson A.F., Pickett A.K., Rozycki P. (2001). Computational methods for predicting impact damage in composite structures. *Composites Sciences and Technology*, vol. 61, p. 2183-2192.

Ladevèze P. (1992). A damage computational method for composite structures. *Computers & Structures*, vol. 44, p. 79-87.

Ladevèze P., Allix O., Deu J.F., Lévêque D. (2000). A mesomodel for localisation and damage computation in laminates. *Computer Methods in Applied Mechanics and Engineering*, vol. 183, p. 105-122.

Launay A., Marco Y., Maitournam M.H., Raoult I. (2013). Modelling of the influence of temperature and relative humidity on the time-dependent mechanical behaviour of a short glass fibre reinforced polyamide. *Mechanics of Materials*, p. 1-10.

Liang B., Hamila N., Peillon M., Boisse P. (2014). Analysis of thermoplastic prepreg bending stiffness during manufacturing and of its influence on wrinkling simulations. *Composites: Part A*, vol. 67, p. 111-122.

Marguet S. (2007). *Contribution à la modélisation du comportement mécanique des structures sandwichs soumises à l'impact*. Thèse de doctorat, Université de Nantes.

Rozycki P. (2000). *Contribution au développement de lois de comportement pour matériaux composites soumis à l'impact*. Thèse de doctorat, Université de Valencienne et du Hainaut-Cambrésis.