# Modélisation du procédé de thermoestampage de composites préimprégnés à matrice thermoplastique

Eduardo Guzman Maldonado, Hu Xiong, Nahiene Hamila, Philippe Boisse

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

## To cite this version:

Eduardo Guzman Maldonado, Hu Xiong, Nahiene Hamila, Philippe Boisse. Modélisation du procédé de thermoestampage de composites préimprégnés à matrice thermoplastique. Revue des composites et des matériaux avancés = Journal of Composite and Advanced Materials, 2018, 28, pp.9-33. 10.3166/rcma.28.9-33. hal-04674381

HAL Id: hal-04674381

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

Submitted on 21 Aug 2024

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>
# Modélisation du procédé de thermoestampage de composites préimprégnés à matrice thermoplastique

Eduardo Guzman-Maldonado¹, Hu Xiong², Nahiene Hamila³, Philippe Boisse⁴

1. Université de Lyon, LaMCoS UMR 5259
27 Av. Jean Capelle
69621 Villeurbanne, France
eduardo.guzman-maldonado@insa-lyon.fr

2. Université de Lyon, LaMCoS UMR 5259
27 Av. Jean Capelle
69621 Villeurbanne, France
hu.xiong@insa-lyon.fr

3. Université de Lyon, LaMCoS UMR 5259
27 Av. Jean Capelle
69621 Villeurbanne, France
nahiene.hamila@insa-lyon.fr

4. Université de Lyon, LaMCoS UMR 5259
27 Av. Jean Capelle
69621 Villeurbanne, France
philippe.boisse@insa-lyon.fr

<footnote>RÉSUMÉ. Au cours des dernières années, l'industrie automobile a montré un intérêt croissant pour les procédés de fabrication et de transformation de matériaux composites à matrice thermoplastique. Cela est favorisé par le développement et l'optimisation des procédés de mise en forme tels que le thermostampage. La modélisation et la simulation de ce procédé sont des étapes importantes pour la prédiction des propriétés mécaniques et de la faisabilité technique des pièces à géométrie complexe. Dans ce papier, un modèle viscohyperélastique avec une dépendance thermique est présenté afin de décrire le comportement du composite thermoplastique à l'état fondu. Au cours des simulations, des calculs thermiques et mécaniques sont effectués de manière séquentielle afin d'actualiser les propriétés mécaniques en fonction de la température. L'identification des propriétés thermiques est obtenue par homogénéisation à partir des analyses au niveau mésoscopique du matériau. L'influence de l'évolution du champ de température est analysée.</footnote>
ABSTRACT. Pre-impregnated thermoplastic composites are widely used in the aerospace industry for their excellent mechanical properties, impact resistance and fatigue strength all at lower density than other common materials. In recent years, the automotive industry has shown increasing interest in the manufacturing processes of thermoplastic-matrix composites materials, especially in thermoforming techniques for their rapid cycle times and the possible use of pre-existing equipment.

An important step in the prediction of the mechanical properties and technical feasibility of parts with complex geometry is the use of modelling and numerical simulations of these forming processes which can also be capitalized to optimize manufacturing practices.

This paper aims to present an approach to the simulation of thermoplastic prepreg composites forming. The proposed model is based on convolution integrals defined under the principles of irreversible thermodynamics and within a hyperelastic framework. The hyperelastic potential is built from the contribution of four major deformation modes that assumed independent: the elongation in the warp and weft directions, the in-plane shear, and the bending deformations. Consequently the corresponding strain energy potentials which are uncoupled can be identified by classical tests on textile composite shells, such as tension, in plane-shear and bending tests. The viscoelasticity is introduced exclusively for the in-plane shear mode. The simulation is composed of thermal analyses and forming simulations performed alternatively. This ensures a coupling between the mechanical and thermal simulations. The deformation of a unit cell and the modification of contacts with the tools change the local thermal properties and temperatures. The thermal conductivities are determined by mesoscopic analysis. The woven reinforced prepreg has a periodicity that can be exploited to perform a homogenization simulation and get the macroscopic conductivity from the geometry of the mesostructure and from the thermal properties of the yarns and the polymer.

MOTS-CLÉS : préimprégné, thermoplastique, thermomécanique, viscoélasticité, analyse par éléments finis, mise en forme.

KEYWORDS: prepreg, thermoplastic, thermomechanical, viscoelasticity, finite element analysis, forming.

DOI:10.3166/RCMA.28.9-33 © 2018 Lavoisier

# 1. Introduction

Les matériaux composites à matrice thermoplastique et fibres continues sont relativement jeunes dans le marché industriel comparés aux composites à matrice thermodurcissable. Cependant, les résines thermoplastiques présentent de nombreux avantages : une meilleure stabilité thermique, une meilleure résistance aux chocs et à l'environnement ainsi qu'un cycle de mise en forme plus court. Par conséquent, ce matériau est devenu un candidat privilégié pour la production en série des structures à haute performance.

Depuis quelques années, différents secteurs dans l'industrie, dont le secteur énergétique (Brøndsted et al., 2005), sportif (Easterling, 2012), médical (Ramakrishna et al., 2001) et plus particulièrement le secteur automobile (Verrey et al., 2006; Fuchs et al., 2008), ont montré un intérêt croissant dans leurs processus de fabrication et
de transformation. Parmi ces procédés, le thermoestampage s'avère particulièrement intéressant. Il est basé sur des méthodes similaires à l'emboutissage des tôles métalliques. Cette dernière, étant un procédé facilement automatisable, permet des cycles de production rapides. Afin de prédire les propriétés mécaniques de la pièce, sa morphologie et des éventuels défauts, il est important de prendre en compte l'influence des paramètres de fabrication. Cela est possible à condition de disposer d'un modèle de comportement adapté aux conditions du procédé. Pour les industriels, la méthode d'essai-erreur devient prohibitive pour une large gamme de fabrication de pièces. Afin de contourner ce problème, la simulation se présente comme une très bonne alternative.

Différentes approches ont été développées ces dernières années afin de modéliser la mise en forme des composites. Elles peuvent être classifiées en trois groupes : cinématiques (Van Der Weeën, 1991), discrètes (Cherouat, Billoët, 2001 ; Sze, Liu, 2005 ; Boubaker et al., 2007) et continues.

Parmi ces familles, les modèles continus sont les plus étudiés. Un des premiers modèles pour la modélisation des préimprégnés a été proposé par Rogers (1989) en par tant des travaux de Spencer (1972 ; 1984). Dans ce modèle, les fibres sont considérées quasi inextensibles et le matériau composite est considéré comme un modèle renforcé idéalisé. Ce modèle continue à être aujourd'hui une référence pour différents travaux (Cherouat, Billoët, 2001 ; Thije et al., 2007 ; Jauffrès et al., 2009 ; Harrison, Clifford, Long, Rudd, 2004) qui proposent des modifications et implémentent des techniques plus complexes, telles que les méthodes d'homogénéisation (Hsiao, Kikuchi, 1999). Le but de cet article est de présenter une approche différente pour la modélisation de la mise en forme de matériaux composites à matrice thermoplastique.

Le modèle proposé est construit dans le cadre de la thermodynamique des processus irréversibles avec variables internes. Il permet de découpler la contribution des principaux modes de déformation macroscopiques du renfort à l'échelle macroscopique. Le comportement viscoélastique est uniquement associé au mode de déformation en cisaillement plan.

## 2. Le procédé de thermostatampage

Le thermoestampage est une solution très intéressante pour la production en masse des structures composites. En effet, étant une adaptation de la mise en forme des tôles métalliques par emboutissage, il existe la possibilité de réutiliser les équipements déjà existants ainsi que son système d'automatisation (Campbell, 2010).

La première étape du thermoformage consiste à chauffer le matériau au-dessus de sa température de fusion. Les principales méthodes utilisées pour le préchauffement sont les lampes infrarouges (IR), les fours à convection et les presses à plateaux chauffants. Après chauffage, le matériau est transporté vers la presse. Le temps de transfert est un paramètre critique. Le matériau doit être transporté le plus rapidement possible afin d'éviter que le matériau ne refroidisse au-dessus de sa température de cristallisation.
sation. Les systèmes d'automatisation actuels tels que les convoyeurs et les robots permettent de réduire ce temps. La perte de température dans cette étape est de l'ordre de 10 %. Le matériau est ensuite mis en forme par application d'une pression entre deux outils complémentaires, empreinte de la structure finale. L'étape suivante est la phase de consolidation, la pièce est maintenue sous pression afin d'enlever les porosités et d'assurer la cohésion entre les différentes couches. La structure est finalement démoulée et poursuit son refroidissement à l'air ambiant l'ensemble des phases sont schématisées sur la figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1. Étapes du procédé de thermoestampage</figcaption></figure>

L'ensemble des étapes sont réalisées en quelques minutes, ce qui fait du thermo-stampage un candidat privilégié pour la production en masse de pièces composites, notamment dans le secteur automobile. Cependant, afin d'optimiser le procédé, les phénomènes physiques en jeu doivent être maîtrisés. La température apparaît comme un facteur essentiel puisqu'elle modifie le comportement de la résine. Les propriétés mécaniques et les mécanismes de déformations particuliers du matériau, de par la nature fibreuse du renfort, évoluent constamment. L'identification, la modélisation et l'évolution des propriétés et des mécanismes de déformation au cours de la mise en forme sont essentielles.

### 3. Modèle visco-hyperélastique pour la mise en forme de préimprégnés

#### 3.1. Modèle hyperélastique pour les renforts tissés

Un renfort tissé peut être considéré comme un matériau orthotrope. Les trois directions principales sont définies respectivement par les directions chaîne ($L_1$), trame
$(L_2)$ et une troisième direction perpendiculaire aux deux précédentes $(L_3)$. Des travaux récents ont permis la mise en place des modèles de comportement du type hyperélastique adaptés aux tissus secs (Aimene et al., 2008; Charmetant et al., 2012). Cette approche, inspirée des travaux de Criscione (Criscione et al., 2001), est basée sur une décomposition multiplicative du gradient de la transformation. Chacun des termes de cette décomposition est caractérisé par un invariant dit « physique », relié aux invariants classiques. L'avantage de ce type de formulation réside dans le découplage des modes de déformation, permettant ainsi de les identifier séparément. Les modes de déformation considérés dans les cas d'un renfort composite 2D initialement orthotrope ainsi que ses respectifs invariants sont :

<table><caption>Tableau 1. Modes de déformation et invariants associés des matériaux composites 2D</caption><thead><tr><th>Mode de déformation</th><th>Invariant</th></tr></thead><tbody><tr><td>Élongation direction chaîne</td><td>$I_{elong1} = \frac{1}{2} \ln (\underline{C} : \underline{L}_{11}) = \ln (\lambda_{chane})$</td></tr><tr><td>Élongation direction trame</td><td>$I_{elong1} = \frac{1}{2} \ln (\underline{C} : \underline{L}_{22}) = \ln (\lambda_{trame})$</td></tr><tr><td>Cisaillement dans le plan</td><td>$I_{sh} = \frac{\underline{C} : \underline{L}_{12}}{\sqrt{(\underline{C} : \underline{L}_{11})(\underline{C} : \underline{L}_{22})}} = \sin (\gamma)$</td></tr></tbody></table>

Avec $\underline{C}$ le tenseur de déformations de Cauchy-Green droit et $\underline{L}_{ij} = L_i \otimes L_j$ le tenseur de structure. L'énergie surfacique de déformation est exprimée en fonction de ces invariants physiques grâce à l'hypothèse de découplage :

$$w = w_{elong1}(I_{elong1}) + w_{elong2}(I_{elong2}) + w_{sh}(I_{sh}) \quad (1)$$

Et la contribution en termes du second tenseur de contraintes de Piola-Kirchhoff, à partir de la relation classique suivante :

$$\underline{S} = 2 \frac{\partial w}{\partial \underline{C}} = 2 \frac{\partial w_k}{\partial I_k} \frac{\partial I_k}{\partial \underline{C}} \quad (2)$$

Avec $w_k$ le potentiel associé à chaque mode de déformation.

## 3.2. Prise en compte du comportement viscoélastique

Afin de prendre en compte l'effet de la résine, il est possible d'introduire sur le modèle précédent une contribution viscoélastique (Guzman-Maldonado et al., 2015). Le premier pas dans la formulation de ce type de modèle consiste à découpler le comportement membranaire et le comportement en flexion du matériau. Le potentiel d'énergie interne (surfacique) est alors défini, comme la somme des potentiels d'énergie associés aux déformations membranaires et à la contribution en flexion.
$$w = w_{mem} + w_{ben} \qquad (3)$$

Le potentiel d'énergie interne membranaire ($w_{mem}$) est exprimé en fonction des invariants associés aux trois modes de déformations cités précédemment. Pour la contribution en flexion ($w_{ben}$), on peut définir également deux invariants $I_{\chi_1}$ et $I_{\chi_2}$ associées à ce mode de déformation dans le sens chaîne et trame, définis d'après (Spencer, Soldatos, 2007) par :

$$I_{\chi_i} = \underline{l}_i \cdot \underline{\chi} \cdot \underline{l}_i = \chi_{ii} \quad i = 1, 2 \qquad (4)$$

où $\underline{\chi}$ représente le tenseur de courbure du renfort dans la configuration actuelle et $\underline{l}_i$ la direction actuelle de la fibre matérielle. Les moments de flexion dans le sens chaîne et trame sont donnés par :

$$M_{ii} = \frac{1}{J_s} \frac{\partial w_{ben}}{\partial I_{\chi}^i} \qquad (5)$$

Avec $J_s$ le rapport de surfaces entre la configuration actuelle et initiale.

La deuxième étape dans cette formulation consiste à considérer que le comportement viscoélastique est uniquement associé au cisaillement plan. Afin de prendre en compte cette contribution, un modèle viscohyperélastique non linéaire est construit à partir de la théorie des processus thermodynamiques irréversibles et le concept de variables internes.

Basé sur les travaux de Simo (1987), ce modèle considère une fonction de l'énergie libre sous la forme suivante :

$$\Psi(\underline{C}, \underline{L}_1, \underline{L}_2) = w_{mem}(\underline{C}, \underline{L}_1, \underline{L}_2) - \frac{1}{2} \sum_{i=1}^{N} (\underline{C} : \underline{Q}_i) + \Xi \left( \sum_{i=1}^{N} \underline{Q}_i \right) \qquad (6)$$

où les $\underline{Q}_i$ représentent un set des variables internes non physiquement mesurables dans la configuration de référence. Ces variables sont associées aux différents mécanismes de relaxation dus aux interactions au sein du matériau (contacts entre les mèches et les fibres).

Le second tenseur de Piola-Kirchoff peut être décomposé suivant les différentes modes de déformation afin exprimer chacune de leurs contributions :

$$\underline{S} = \underline{S}_{elong1} + \underline{S}_{elong2} + \underline{S}_{sh} \qquad (7)$$

où $\underline{S}_{elong1}$ et $\underline{S}_{elong2}$ représentent les contributions purement hyperélastiques liées à l'élongation dans le sens chaîne et trame respectivement.

$$\underline{S}_{elong_i} = 2 \frac{\partial w_{elong_i}}{\partial I_{elong_i}} \frac{\partial I_{elong_i}}{\partial \underline{C}} \quad (i = 1, 2) \qquad (8)$$
Et $\underline{S}_{sh}$ la contribution viscohyperélastique liée au cisaillement plan :

$$ \underline{S}_{sh} = 2 \frac{\partial w_{sh}(\underline{C})}{\partial \underline{C}} - \sum_{i=1}^{N} \underline{Q}_{i} \qquad (9) $$

Les variable internes $\underline{Q}_{i}$ sont déterminées à partir du système d'équations d'évolution suivant :

$$ \left| \begin{array}{l} \dot{\underline{Q}}_{i} + \frac{1}{\tau_{i}} \underline{Q}_{i} = \frac{\gamma_{i}}{\tau_{i}} \left[ 2 \frac{\partial w_{sh}(I_{sh})}{\partial \underline{C}} \right] \\ \lim_{t \to -\infty} \underline{Q}_{i}(t) = 0 \end{array} \right. \qquad (10) $$

où les $\gamma_i \in [0, 1)$ et $\tau_i \ge 0$ représentent les différents paramètres matériaux viscoélastiques et les temps de relaxation associés au matériau.

Finalement, le second tenseur de Piola-Kirchhoff associé à la contribution membranaire viscohyperélastique est donné (sous sa forme de convolution) par :

$$ \underline{S}(t) = \underline{S}_{elong1} + \underline{S}_{elong2} + \int_{-\infty}^{t} G(t-s) \frac{d}{ds} \left[ 2 \frac{\partial w_{sh}}{\partial \underline{C}} \right] ds \qquad (11) $$

Avec $G(t)$ la fonction de relaxation associée :

$$ G(t) = \gamma_{\infty} + \sum_{i=1}^{N} \exp \left( -\frac{t}{\tau_{i}} \right) \gamma_{i} \qquad (12) $$

Le modèle de comportement précédent fait intervenir uniquement des tenseurs lagrangiens qui sont objectifs par définition. L'objectivité de la loi est par conséquent respectée à condition que les potentiels d'énergie de déformation associés à chaque mode de déformation soient admissibles (Itskov, Aksel, 2004).

### 3.3. Identification des paramètres

#### 3.3.1. Matériau

Le matériau considéré dans cette section, est une plaque composite thermoplastique préconsolidée. La plaque est composée de 5 couches Satin-8 en fibre de verre et préimprégnés d'une résine Polyamide 66 (figure 2a). Les principales caractéristiques associées à ce matériau sont détaillées dans le tableau 2. La température de fusion de la résine thermoplastique varie entre 257 °C et 265 °C en fonction de la vitesse de chauffage (figure 2b).
<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Microtomographie du préconsolidé thermoplastique et (b) DSC pour différents vitesses de chauffage/refroidissement</figcaption></figure>

<table><caption>Tableau 2. Principales caractéristiques du préconsolidé</caption><tr><td>Matrice</td><td>PA 66</td></tr><tr><td>Fibres</td><td>E-Glass</td></tr><tr><td>Type d'armure</td><td>Satin-8H</td></tr><tr><td>Volume fibre (%)</td><td>51-53</td></tr><tr><td>Contenu vide (%)</td><td>&lt; 0.2</td></tr><tr><td>Densité (g·cm<sup>-3</sup>)</td><td>1.87</td></tr><tr><td>Épaisseur (mm)</td><td>2</td></tr><tr><td>Morphologie</td><td>5 plis, Satin-8H</td></tr></table>

### 3.3.2. Cisaillement dans le plan

La caractérisation du comportement en cisaillement plan a été effectuée en utilisant le Bias-extension test à hautes températures (Lebrun *et al.*, 2003; Wang *et al.*, 2012; Guzman-Maldonado *et al.*, 2015; Boisse *et al.*, 2016). Cependant, d'autres moyens de caractérisation peuvent être utilisés, comme le Picture Frame (Harrison, Clifford, Long, 2004; Dangora *et al.*, 2015).

Le Bias Extension Test, dans le cas de thermoplastiques mono- et multi-plis, est un essai assez délicat à mettre en œuvre. Afin d'étudier le comportement des préimprégnés pendant la mise en forme, il faut amener l'éprouvette au-dessus du point de fusion. Deux aspects doivent être assurés au cours de l'essai : le respect des conditions limites et le contrôle de la température au sein du matériau.

Le banc expérimental est présenté dans la figure 3. L'éprouvette est placée dans un enceinte thermique et installée dans une machine à traction universelle. Le déplacement est mesuré par le capteur interne de la machine et l'effort à l'aide d'un capteur 500 N (± 0.3 %). Le maintien en position est assuré par des vis à bout téton qui blo-
queront l'éprouvette à travers l'épaisseur, ainsi qu'à l'aide d'un système de ressorts dans le mors permettant d'assurer un effort de serrage constant.

<figure><img src="image_5.png" /><figcaption>Figure 3. Montage expérimental et détail de conception des mors</figcaption></figure>

Afin de mettre en évidence et de caractériser le comportement viscoélastique de ce mode de formation, le test est effectué avec un déplacement par paliers. La géométrie déformée de l'éprouvette et les résultats obtenus pour une température de 27 °C sont présentés dans la figure 4.

L'algorithme d'optimisation de Levenberg-Marquardt (Schnur, Zabaras, 1992) est utilisé pour identifier, à partir des résultats expérimentaux, les différents paramètres matériaux $\gamma_i$, les temps de relaxation $\tau_i$ ainsi que les coefficients $p_{sh}$ du potentiel $w_{sh}$, définie par :

$$w_{sh} = \sum_{i=1}^{3} p_{sh} (I_{sh})^{2i} \qquad (13)$$

L'effort externe théorique $F_{ext}$ est déduit à partir de la relation entre la puissance interne développée par les couples/efforts en cisaillement et la puissance externe :

$$P_{int} = P_{ext} \Rightarrow \int_{\Omega_0} \frac{1}{2} S : \dot{\underline{C}} dV = F_{ext} \dot{\underline{d}} \qquad (14)$$

L'identification est effectuée à partir de quatre températures entre 260 °C et 300 °C. Les paramètres intermédiaires peuvent être interpolés linéairement en fonction de la température (Eq. 15).
<figure><img src="image_6.png" /><figcaption>Figure 4. (a) Bias-extension test : éprouvette en fin d'essai. (b) Résultats expérimentaux pour une température de 270 °C</figcaption></figure>

$$w_{sh}(I_{sh}, T) = \begin{cases} (\alpha_1^w T + \beta_1^w) \sum_{i=1}^{3} p_{sh_i}(I_{sh})^{2i} & 260 \le T \le 270 \\ (\alpha_2^w T + \beta_2^w) \sum_{i=1}^{3} p_{sh_i}(I_{sh})^{2i} & 270 < T \le 300 \end{cases} \quad (15)$$

$$\gamma_1(T) = \alpha_1^\gamma T + \beta_1^\gamma$$

$$\gamma_2(T) = \alpha_2^\gamma T + \beta_2^\gamma$$

$$\tau_1(T) = \alpha_1^\tau T + \beta_1^\tau$$

$$\tau_2(T) = \alpha_2^\tau T + \beta_2^\tau$$

### 3.3.3. Flexion hors plan

La caractérisation en flexion a été effectuée en utilisant un test du type cantilever à hautes températures (Bilbao et al., 2010; Liang et al., 2014). Le principe de l'essai consiste à faire fléchir un échantillon, en porte-à-faux, soumis à son propre poids et pour une température donnée. Le profil obtenu est enregistré par une capture optique et pour différentes températures. Après dépouillement, le moment de flexion est tracé en fonction de la courbure puis une rigidité en fonction de la température peut être associée.

Le dispositif expérimental est placé dans une enceinte thermique comme illustré figure 5, le support est équipé d'un système de chauffage indépendant. En effet, le temps de chauffage du massif métallique est plus important que celui de l'éprouvette, ce qui génère un gradient de température.
Afin d'assurer une répartition homogène de la température, une série de thermo-couples sont placés aux voisinages de l'éprouvette et dans le dispositif.

<figure><img src="image_7.png" /><figcaption>Figure 5. Essai cantilever à hautes températures (Liang et al., 2014)</figcaption></figure>

Une série de tests a été effectuée dans le sens chaîne et trame sur des éprouvettes de 80 mm de longueur et sur un intervalle de température [260 – 320 °C]. Les profils déformés mesurés ainsi que la courbe moment de flexion en fonction de la courbure, sont montrés dans la figure 6. La déformée maximale est obtenue pour une température de 300 °C, au-dessus de cette température la rigidité en flexion reste constante et est voisine à celle d'un tissu sec.

<figure><img src="image_8.png" /><figcaption>Figure 6. Mesure de la rigidité en flexion du préimprégné à différentes températures
(Liang et al., 2014)</figcaption></figure>
La rigidité en flexion dans le cas des thermoplastiques dépend fortement de la température. Le potentiel associé à ce mode de déformation peut être représenté par une fonction quadratique :

$$w_{\text{ben}}(I_{\chi_i}) = p_{\chi_i}(T) (I_{\chi_i})^2 \qquad (16)$$

où $p_{\chi_i}$ est une fonction linéaire thermo-dépendante.

$$p_{\chi_i}(T) = \alpha_i^{\text{ben}} \cdot T + \beta_i^{\text{ben}} \qquad (17)$$

### 3.3.4. Élongation dans le sens des fibres

Dans une première approche, il est considéré que la rigidité en élongation du matériau dépend uniquement de la rigidité des fibres $E_f$. La réponse en traction d'une mèche relève d'un comportement quasi linéaire (Dumont, 2003). Afin de reproduire ce comportement, le potentiel de déformation sera défini par :

$$w_{\text{elongi}} = p_{\text{elongi}}(I_{\text{elongi}})^2 \quad i = 1, 2 \qquad (18)$$

La valeur de $p_{\text{elongi}}$ peut être approximée à partir du taux de fibres dans les directions chaîne et trame ($v_{f_i}$). Pour une sollicitation uni-axiale on peut montrer approximer cette valeur par :

$$p_{\text{elongi}} \approx \frac{E_f \cdot v_{f_i}}{2} \qquad (19)$$

Le taux de fibres est estimé, par exemple, à partir des images issues de la microtomographie X du matériau (Naouar et al., 2014).

### 3.4. Détermination des propriétés thermiques par homogénéisation

Avant et durant l'étape de mise en forme, des pertes de chaleur par conduction ont lieu suite au contact du matériau avec les outils. Or, comme mentionné précédemment, le comportement mécanique du matériau varie de manière importante dans l'intervalle de température de la mise en forme. Par conséquent, il est important de prendre en compte ce champ de température hétérogène induit par le contact ainsi que son évolution au cours de la mise en forme (Guzman-Maldonado et al., 2016).

Le comportement thermique des matériaux composites préimprégnés est complexe à décrire de manière analytique. En effet, la structure du matériau est fortement hé térogène et anisotrope. Les propriétés thermiques du matériau dépendent à la fois de la taille et de la répartition spatiale des constituants (fibre + matrice), ainsi que des directions privilégiées du matériau. Si en outre, les propriétés du matériau sont
thermo-dépendantes le problème devient non linéaire. Tel est le cas de la conductivité thermique.

Dans ce papier, les différentes propriétés thermiques du matériau ont été déterminées par une méthode d'homogénéisation (voir par exemple (Ozdemir et al., 2008)). En effet, les renforts tissés ont naturellement une propriété de périodicité qui peut être exploitée pour réaliser des calculs d'homogénéisation à l'échelle mésoscopique. Le tenseur de conductivité thermique macroscopique peut être estimé à partir de la géométrie, et les propriétés thermiques des constituants à l'échelle du VER (Volume Élémentaire Représentatif).

### 3.4.1. Schéma d'homogénéisation

Le début de ce schéma consiste à définir le problème à l'échelle mésoscopique du matériau ($\Omega_m$), à savoir, à l'échelle du VER. Le problème est traité en régime stationnaire. En absence des sources internes l'équilibre thermique est donné par :

$$ \underline{\nabla}_m \cdot \underline{q}_m(\underline{x}) = 0 \qquad (20) $$

où $q_m$ représente le flux de chaleur volumique dans tout point $\underline{x}$ à l'échelle mésoscopique.

La résolution du système (20) nécessite une définition correcte des conditions aux limites dans le VER. Ces conditions aux limites sont imposées à partir du gradient de température macroscopique $\underline{\nabla}_M T_M$.

Les relations de transition entre l'échelle macroscopique et mésoscopique définies d'après (Ozdemir et al., 2008) par :

$$ \frac{1}{V} \int_{\Omega_m} \underline{\nabla}_m T_m \, dV = \underline{\nabla}_M T_M + \frac{1}{V} \int_{\partial \Omega_m} T_f \underline{n} \, dS \qquad (21) $$

où V représente le volume du domaine $\Omega_m$ de frontière $\partial\Omega_m$, $\underline{\nabla}_m T_m$ représente le gradient de température à l'échelle mésoscopique et $T_f(\underline{x})$ représente une fluctuation thermique (variation de deuxième ordre) (Geers et al., 2003; Ozdemir et al., 2008).

Afin d'imposer l'égalité entre le gradient macroscopique et le gradient mésoscopique moyenné, la condition à imposer sur les fluctuations aux limites du domaine doit être de la forme :

$$ \frac{1}{V} \int_{\partial \Omega_m} T_f \underline{n} \, dS = 0 \qquad (22) $$

Sans perte de généralités, considérons un VER hexaédrique cisaillé de volume V avec un repérage défini sur la figure 7.
<figure><img src="image_9.png" /><figcaption>Figure 7. Domaine d'étude et paramétrisation géométrique du VER</figcaption></figure>

Les conditions suffisantes à imposer peuvent être exprimées en termes des températures à la surface du domaine et du gradient thermique macroscopique sous la forme :

$$T_{mB} - T_{mA} = \underline{\nabla}_M T_M \cdot (\underline{x}_B - \underline{x}_A) \quad (23a)$$

$$T_{mD} - T_{mC} = \underline{\nabla}_M T_M \cdot (\underline{x}_D - \underline{x}_C) \quad (23b)$$

$$T_{mF} - T_{mE} = \underline{\nabla}_M T_M \cdot (\underline{x}_F - \underline{x}_E) \quad (23c)$$

Avec $\underline{x}_I \in \partial\Omega_I$, ($I = A, B, C, D, E, F$) l'ensemble de points matérielles sur la surface $I$ du domaine. Les points $\underline{x}_B, \underline{x}_D$ et $\underline{x}_F$ sont l'image des points par projection périodique sur les surfaces opposées $\underline{x}_A, \underline{x}_C$ et $\underline{x}_E$. Dans le cas d'un VER non cisaillé, ils correspondent aux points vis-à-vis sur chaque facette. Dans le cas général la cinématique de déformation doit être prise en compte afin de respecter cette correspondance spatiale. Ces conditions sont connues sous le nom de conditions aux limites périodiques et conduisent naturellement à des conditions antipériodiques en termes des flux de chaleur aux limites du domaine :

$$\underline{q}_{mB} \cdot \underline{n}_1 = -\underline{q}_{mA} \cdot \underline{n}_1 \quad (24a)$$

$$\underline{q}_{mD} \cdot \underline{n}_2 = -\underline{q}_{mC} \cdot \underline{n}_2 \quad (24b)$$

$$\underline{q}_{mF} \cdot \underline{n}_3 = -\underline{q}_{mE} \cdot \underline{n}_3 \quad (24c)$$

avec $\underline{n}_1, \underline{n}_2$ et $\underline{n}_3$ les vecteurs unitaire normaux aux facettes B, D et C aux points $\underline{x}_B, \underline{x}_D$ et $\underline{x}_F$.

Le flux macroscopique est également déterminé à partir de la moyenne volumique du flux microscopique.

$$\frac{1}{V} \int_{\Omega_m} \underline{q}_m dV = \underline{q}_M \quad (25)$$
Il a été montré dans (Ozdemir et al., 2008) que cette expression permet d'assurer la consistance de l'entropie entre l'échelle macroscopique et l'échelle mésoscopique :

$$ \frac{1}{V} \int_{\Omega_m} \nabla_m T_m \cdot q_m \, dV = \nabla_M T_M \cdot q_M \qquad (26) $$

### 3.4.2. Calcul du tenseur de conductivité équivalente

Afin de calculer le tenseur de conductivité équivalente, on considère le problème de la conduction macroscopique en régime stationnaire.

$$ q_M = \underline{K}_M \nabla_M T_M \qquad (27) $$

ou bien sous forme matricielle :

$$ \left\{ \begin{matrix} q_{M1} \\ q_{M2} \\ q_{M3} \end{matrix} \right\} = \begin{bmatrix} k_{11} & k_{12} & k_{13} \\ k_{21} & k_{22} & k_{23} \\ k_{31} & k_{32} & k_{33} \end{bmatrix} \left\{ \begin{matrix} (\nabla_M T_M)_1 \\ (\nabla_M T_M)_2 \\ (\nabla_M T_M)_3 \end{matrix} \right\} \qquad (28) $$

Le tenseur $\underline{K}_M$ peut être entièrement déterminé en imposant trois gradients macroscopiques unitaires dans les directions principales de la base considérée :

$$ \underline{\nabla}_M T_M = [1 \ 0 \ 0]^T \quad \underline{\nabla}_M T_M = [0 \ 1 \ 0]^T \quad \underline{\nabla}_M T_M = [0 \ 0 \ 1]^T \qquad (29) $$

### 3.4.3. Application

Le modèle numérique du VER comportant les 5 plis du matériau a été créé à l'aide du logiciel Texgen, développé à l'université de Nottingham (Sherburn, 2007). Les différents paramètres géométriques ont été mesurés directement à partir des images de la microtomographie du matériau préconsolidé (figure 2a). Les directions locales des fibres sont également récupérées et exportées afin de créer des repères locaux associés aux mèches.

La méthode d'homogénéisation choisie implique que les propriétés thermiques et le comportement des constituants soient connus. On suppose que le comportement de la résine polymère (PA 66) et des fibres, est isotrope. La conductivité thermique du PA 66 ($k_p$), dont l'évolution est donnée figure 8 sera considérée thermo-dépendante. Concernant la conductivité des fibres ($k_f$), elle sera supposée constante et égale à $1 \, W \cdot m^{-1} \cdot K^{-1}$ (Kutz, 2002; Thomas et al., 2012).

La conductivité thermique de la mèche est du type isotrope traverse. La conductivité dans le sens longitudinal (i.e. parallèle à la direction des fibres) et la conductivité transversale (i.e. perpendiculaire à la direction des fibres) dans une section de la mèche sont différentes.
<figure><img src="image_10.png" /><figcaption>Figure 8. Conductivité thermique du PA 66 en fonction de la température
(Mark, 2007)</figcaption></figure>

La mèche peut être modélisée localement, comme une série d'inclusions cylin-
driques isotropes (les fibres) noyées dans un milieu isotrope (la matrice). Dans cette
configuration, la conductivité thermique dans le sens longitudinal $k_{lon}^m$ peut être ap-
proximée par une loi de mélange classique de la forme :

$$
k_{lon}^{m} = k_{f} \cdot v_{f}^{m} + k_{p} \cdot (1 - v_{f}^{m}) \quad (30)
$$

où $v_f^m$ représente le taux volumique des fibres dans les mèches.

Dans le cas de la conductivité transversale $k_{tra}^m$, Hasselman (Hasselman *et al.*, 1993) propose une expression analytique dérivée de la théorie de Maxwell et Rayleigh, prenant en compte la conductance entre la fibre et la matrice :

$$
k_{tra}^{m} = k_{p} \cdot \frac{\left( \frac{k_{f}}{k_{m}} - 1 - \frac{k_{f}}{a h_{c}} \right) v_{f}^{m} + \left( 1 + \frac{k_{f}}{a h_{c}} + \frac{k_{f}}{k_{m}} \right)}{\left( 1 + \frac{k_{f}}{a h_{c}} - \frac{k_{f}}{k_{m}} \right) v_{f}^{m} + \left( 1 + \frac{k_{f}}{a h_{c}} + \frac{k_{f}}{k_{m}} \right)} \quad (31)
$$

où $h_c$ représente la conductance à l'interface du contact fibre-matrice et $a$ le rayon de la fibre. Pour cette étude, on considère le cas d'un contact parfait ($h \to \infty$).
$$k_{tra}^{m} = k_{p} \cdot \frac{\left(\frac{k_{f}}{k_{m}} - 1\right) v_{f}^{m} + \left(1 + \frac{k_{f}}{k_{m}}\right)}{\left(1 - \frac{k_{f}}{k_{m}}\right) v_{f}^{m} + \left(1 + \frac{k_{f}}{k_{m}}\right)} \quad (32)$$

Les relations (30) et (32) font intervenir le taux de fibres dans les mèches $v_f^m$. Cette information peut être calculée à partir de la connaissance du taux de fibres dans le matériau $v_f$ (généralement connu) et de l'estimation du taux de mèches $v_m$ dans le VER, soit :

$$v_f^m = \frac{v_f}{v_m} \qquad (33)$$

Le taux de mèches du matériau $v_m$ peut être obtenu, par exemple, à partir de l'analyse des images de la tomographie.

Un modèle numérique voxelisé est considéré afin de simplifier l'application des conditions aux limites périodiques. Le maillage sur les différentes faces du domaine, est par défaut coïncident, indépendamment du type de renfort considéré.

Une série de simulations a été effectuée pour différents angles de cisaillement afin de construire une cartographie de l'évolution de la conductivité thermique en fonction de la température mais également en fonction de l'angle de cisaillement du renfort. La figure 9 illustre le gradient de température macroscopique imposé pour différentes configurations matérielles.

<figure><img src="image_11.png" /><figcaption>Figure 9. Gradient thermique macroscopique imposé dans la direction $e_1$ pour différents angles de cisaillement (à gauche) 0°, (au centre) 20° et (à droite) 40°</figcaption></figure>

Les différentes composantes du tenseur de conductivité sont présentées dans la figure 10 en fonction de la température et de l'angle de cisaillement ($\gamma$) dans la base orthonormée ($e_1, e_2, e_3$), où $e_1$ est colinéaire avec la direction des chaînes.
<figure><img src="image_12.png" /><figcaption>Figure 10. Composantes du tenseur de conductivité en fonction de la température et de l'angle de cisaillement</figcaption></figure>

4. Simulation numérique

Le modèle de comportement décrit précédemment a été implémenté dans un code par éléments finis en dynamique explicite (Plasfib). La plaque préconsolidée est modélisée par un élément de coque triangulaire à trois degrés de liberté où les contributions membranaires et en flexion sont découplées. La méthode employée pour le calcul de la courbure est reprise des travaux de Hamila et Sabourin sur les éléments du type rotational free (Sabourin, Brunet, 2006; Hamila et al., 2009). La courbure dans la direction des fibres est calculée uniquement à partir de la position des nœuds voisins à l'élément fini. Ceci permet de diminuer le nombre global de degrés de liberté et par conséquent le temps de calcul.

Le modèle numérique comporte un quart de la géométrie réelle et est composé par un poinçon, une matrice et une plaque rectangulaire de dimensions 270 x 190 mm maintenue en contact avec la matrice par 6 serre flans (figure 11). La simulation a été effectuée pour deux orientations matérielles : 0°/90° et -45°/45°. Concernant les conditions initiales et les conditions aux limites, il est considéré que la température avant la mise en forme est connue et égale à 300°. La température des outils est fixé à 250° et la vitesse du déplacement du poinçon est fixé à 240 mm/min pour un déplacement total de 60 mm.
<figure><img src="image_13.png" /><figcaption>Figure 11. Modèle numérique - Double Dome</figcaption></figure>

Dans un premier temps, la température est maintenue constante pendant la mise en forme pour les deux configurations matérielles. Ces résultats sont ensuite comparés avec une deuxième série de simulations où la réactualisation de la carte de température est effectuée. Un couplage faible est considéré entre le problème mécanique et thermique. Le champ de température à l'instant *t* est obtenu en résolvant l'équation de la chaleur dans le cas d'un transfert transitoire par conduction entre l'instant *t* et (*t* − Δ*t*<sub>th</sub>), où Δ*t*<sub>th</sub> représente le pas de réactualisation thermique. Ce calcul tient compte de l'orientation locale de fibres et de l'angle de cisaillement donné par le calcul mécanique à l'instant *t*. Cela permet de modifier les valeurs locales du tenseur de conductivité thermique du matériau. Cette nouvelle carte de température est ensuite considérée pour la détermination des nouveaux paramètres matériaux mécaniques et le processus de réactualisation thermique continue jusqu'à la fin de la mise en forme.

L'ensemble des calculs thermiques ont été effectués sur le logiciel commercial ABAQUS. La mise en données du calcul thermique : définition des contacts, cartes matériaux, orientations locales ; ainsi que le lancement des calculs et la récupération des données sont pilotés par le logiciel Plasfib au cours de la simulation. La figure 12 montre l'évolution du champ de température réel (moyennée dans l'épaisseur) au cours de la simulation.

La figure 13 montre la distribution de l'angle de cisaillement obtenu dans la pièce finale pour deux orientations initiales : 0°/90° (figures 13a,c) et -45°/45° (figures 13b,d), et deux champs de température : un champ homogène constante (figures 13a,b), et un champ hétérogène réactualisé 13c,d)). En comparant le champ de cisaillement obtenu dans le cas de la réactualisation de la température avec le cas homogène (voir par exemple figure 13b vs figure 13d), on constate un écart de la valeur de l'angle de cisaillement de l'ordre de 3° à 5°, ainsi qu'une légère modification de la forme du champ de cisaillement. L'angle de cisaillement dans le cas d'un champ de température homogène et constante est sous-évalué.
<figure><img src="image_14.png" /><figcaption>Figure 12. Évolution du champ de température au cours de la mise en forme</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13. Comparaison du champ de cisaillement en fonction du type de champ de température considéré (a),(b) : Homogène et constante et (b),(c) : Hétérogène et réactualisé</figcaption></figure>

Finalement, la figure 14 montre le gradient de température dans l'épaisseur obtenu à la fin du procédé de mise en forme. Cette information est accessible au cours de la simulation, puisque les calculs sont effectués sur un élément de coque à 5 points d'intégration dans l'épaisseur sur Abaqus.

<figure><img src="image_16.png" /><figcaption>Figure 14. Gradient de température dans l'épaisseur pour $t = 16$ s</figcaption></figure>
**5. Conclusions**

Une nouvelle approche pour la simulation de la mise en forme des composites préimprégnés a été présentée. Le modèle de comportement est construit dans le cadre de la thermodynamique des processus irréversibles. Le modèle permet de découpler la contribution des principaux modes de déformation macroscopiques du renfort à l'échelle macroscopique. Cela grâce à la définition des invariants dits *physiques*, appliqués dans la formulation d'un modèle hyperélastique. Le comportement viscoélastique est uniquement associé au mode de déformation en cisaillement plan. Les propriétés mécaniques inhérentes au modèle sont identifiées à partir de tests simples à hautes températures.

Au cours des simulations, des calculs thermiques et mécaniques sont effectués de manière séquentielle afin d'actualiser les propriétés mécaniques de la loi à partir de l'évolution du champ de température. Le calcul thermique prend en compte l'évolution locale des directions de fibres, ce qui permet d'établir un couplage faible entre le problème mécanique et thermique.

Les simulations ont permis de mettre en évidence l'influence de la température sur l'orientation des fibres dans la pièce. L'évolution de la température doit être prise en compte afin d'avoir une prédiction plus réaliste des propriétés mécaniques et géométries du produit final.

*Remerciements*

Ce travail a été effectué dans le cadre du projet Composite Cab. Ce projet collabo-
ratif a réuni 8 partenaires industriels de la région Rhône-Alpes : Solvay, Saertex,
Reanult Trucks, Plastic Omnium, Mecacorp, Segula, Addiplast, Altair ; et deux uni-
versités : l'Université de Bourgogne et l'Institut National de Sciences Appliquées
de Lyon.

**Bibliographie**

Aimene Y., Hagege B., Sidoroff F., Vidal-Sallé E., Boisse P., Dridi S. (2008). Hyperelastic Approach for Composite Reinforcement Forming Simulations. *International Journal of Material Forming*, vol. 1, n° S1, p. 811-814.

Bilbao E. d., Soulat D., Hivet G., Gasser A. (2010). Experimental study of bending behaviour of reinforcements. *Experimental Mechanics*, vol. 50, p. 333–351.

Boisse P., Hamila N., Guzman-Maldonado E., Madeo A., Hivet G., dell'Isola F. (2016). The bias-extension test for the analysis of in-plane shear properties of textile composite reinforcements and prepregs: a review. *International Journal of Material Forming*, p. 1–20.

Boubaker B. B., Haussy B., Ganghoffer J. F. (2007). Discrete models of woven structures. Macroscopic approach. *Composites Part B: Engineering*, vol. 38, n° 4, p. 498–505.

Brøndsted P., Lilholt H., Lystrup A. (2005). Composite materials for wind power turbine blades. *Annu. Rev. Mater. Res.*, vol. 35, p. 505–538.

Campbell F. C. (2010). *Structural Composite Materials*. ASM International.
Charmetant a., Orliac J. G., Vidal-Sallé E., Boisse P. (2012). Hyperelastic model for large deformation analyses of 3d interlock composite preforms. *Composites Science and Technology*, vol. 72, n° 12, p. 1352-1360.

Cherouat a., Billoët J. L. (2001). Mechanical and numerical modelling of composite manufacturing processes deep-drawing and laying-up of thin pre-impregnated woven fabrics. *Journal of Materials Processing Technology*, vol. 118, n° 1-3, p. 460-471.

Criscione J. C., Douglas A. S., Hunter W. C. (2001). Physically based strain invariant set for materials exhibiting transversely isotropic behavior. *Journal of the Mechanics and Physics of Solids*, vol. 49, p. 871-897.

Dangora L. M., Hansen C. J., Mitchell C. J., Sherwood J. A., Parker J. C. (2015). Challenges associated with shear characterization of a cross-ply thermoplastic lamina using picture frame tests. *Composites Part A: Applied Science and Manufacturing*, vol. 78, p. 181-190.

Dumont F. (2003). *Contribution à l'expérimentation et à la modélisation du comportement mécanique de renforts de composites tissés*. Thèse de doctorat non publiée, Université Paris 6.

Easterling E. (2012). *Advanced materials for sports equipment: how advanced materials help optimize sporting performance and make sport safer*. Springer Science & Business Media.

Fuchs E. R. H., Field F. R., Roth R., Kirchain R. E. (2008). Strategic materials selection in the automobile body: Economic opportunities for polymer composite design. *Composites Science and Technology*, vol. 68, n° 9, p. 1989-2002.

Geers M. G. D., Kouznetsova V. G., Brekelmans W. A. M. (2003). Multi-scale modelling : Computational homogenization in solid mechanics. *Composites Science and Technology*, vol. 1, n° 4, p. 1235-1260.

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: Applied Science and Manufacturing*, vol. 78, p. 211-222.

Guzman-Maldonado E., Hamila N., Naouar N., Moulin G., Boisse P. (2016). Simulation of thermoplastic prepreg thermoforming based on a visco-hyperelastic model and a thermal homogenization. *Materials & Design*, vol. 93, p. 431-442.

Hamila N., Boisse P., Sabourin F., Brunet M. (2009). A semi-discrete shell finite element for textile composite reinforcement forming simulation. *Int. J. Numer. Methods Eng.*, vol. 79, p. 1443-1466.

Harrison P., Clifford M. J., Long A. C. (2004). Shear characterisation of viscous woven textile composites: a comparison between picture frame and bias extension experiments. *Composites Science and Technology*, vol. 64, p. 1453-1465.

Harrison P., Clifford M. J., Long A. C., Rudd C. D. (2004). A constituent-based predictive approach to modelling the rheology of viscous textile composites. *Composites Part A: Applied Science and Manufacturing*, vol. 35, n° 7-8, p. 915-931.

Hasselman D., Donaldson K., Thomas J. (1993). Effective Thermal Conductivity of Uniaxial Composite with Cylindrically Orthotropic Carbon Fibers and Interfacial Thermal Barrier. *Journal of Composite Materials*, vol. 27, n° 6, p. 637-644.
Hsiao S.-W., Kikuchi N. (1999). Numerical analysis and optimal design of composite thermoforming process. *Computer Methods in Applied Mechanics and Engineering*, vol. 177, p. 1–34.

Itskov M., Aksel N. (2004). A class of orthotropic and transversely isotropic hyperelastic constitutive models based on a polyconvex strain energy function. *International Journal of Solids and Structures*, vol. 41, p. 3833–3848.

Jauffrès D., Morris C. D., Sherwood J. a., Chen J. (2009). Simulation of the therostamping of woven composites: Mesoscopic modelling using explicit FEA codes. *International Journal of Material Forming*, vol. 2, n° SUPPL. 1, p. 173–176.

Kutz M. (2002). *Handbook of materials selection*. John Wiley & Sons.

Lebrun G., Bureau M. N., Denault J. (2003). Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics. *Composite Structures*, vol. 61, n° 4, p. 341–352.

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.

Mark J. E. (2007). *Physical Properties of Polymers Handbook*. Springer.

Naouar N., Vidal-Salle E., Schneider J., Maire E., Boisse P. (2014). Meso-scale FE analyses of textile composite reinforcement deformation based on X-Ray computed tomography. *Composite Structures*, vol. 116, n° 0, p. 165–176.

Ozdemir I., Brekelmans W. A. M., Geers M. G. D. (2008). Computational homogenization for heat conduction in heterogeneous solids. *International Journal for Numerical Methods in Engineering*, vol. 73, n° 2, p. 185–204.

Ramakrishna S., Mayer J., Wintermantel E., Leong K. W. (2001). Biomedical applications of polymer-composite materials: a review. *Composites science and technology*, vol. 61, n° 9, p. 1189–1224.

Rogers T. G. (1989). Rheological characterization of anisotropic materials. *Composites*, vol. 20, n° 1, p. 21–27.

Sabourin F., Brunet M. (2006). Detailed formulation of the rotation-free triangular element S3 for general purpose shell analysis. *Engineering Computations*, vol. 23, p. 469–502.

Schnur D. S., Zabaras N. (1992). An inverse method for determining elastic material properties and a material interface. *International Journal for Numerical Methods in Engineering*, vol. 33, n° 10, p. 2039–2057.

Sherburn M. (2007). *Geometric and mechanical modelling of textiles*. Thèse de doctorat non publiée, Université de Nottingham.

Simo J. C. (1987). On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects. *Computer Methods in Applied Mechanics and Engineering*, vol. 60, n° 2, p. 153–173.

Spencer A. J. M. (1972). *Deformations of Fibre Reinforced Materials*. Oxford University Press, London.

Spencer A. J. M. (1984). *Continuum Theory of the Mechanics of Fibre-Reinforced Composites*. Springer Vienna.
Spencer A. J. M., Soldatos K. P. (2007). Finite deformations of fibre-reinforced elastic solids with fibre bending stiffness. *International Journal of Non-Linear Mechanics*, vol. 42, p. 355-368.

Sze K. Y., Liu X. H. (2005). A new skeletal model for fabric drapes. *International Journal of Mechanics and Materials in Design*, vol. 2, no 3-4, p. 225-243.

Thije R. H. W. ten, Akkerman R., Huétink J. (2007). Large deformation simulation of anisotropic material using an updated Lagrangian finite element method. *Computer Methods in Applied Mechanics and Engineering*, vol. 196, no 33-34, p. 3141-3150.

Thomas S., Joseph K., Malhotra S., Goda K., Sreekala M. (2012). *Polymer Composites, Macro-and Microcomposites*. Wiley-VCH: Weinheim, Germany.

Van Der Weeën F. (1991). Algorithms for draping fabrics on doubly-curved surfaces. *International journal for numerical methods in engineering*, vol. 31, no 7, p. 1415-1426.

Verrey J., Wakeman M., Michaud V., M\a anson J.-A. (2006). Manufacturing cost comparison of thermoplastic and thermoset RTM for an automotive floor pan. *Composites Part A: Applied Science and Manufacturing*, vol. 37, no 1, p. 9-22.

Wang P., Hamila N., Pineau P., Boisse P. (2012). Thermomechanical analysis of thermoplastic composite prepregs using bias-extension test. *Journal of Thermoplastic Composite Materials*, vol. 27, no 5, p. 679-698.
<figure><img src="image_17.png" /></figure>