# Méthode éléments finis robuste pour la simulation industrielle de procédés composites

Arnaud Dereims, Sylvain Drapier, Jean-Michel Bergheau, Patrick De Luca

► To cite this version:

Arnaud Dereims, Sylvain Drapier, Jean-Michel Bergheau, Patrick De Luca. Méthode éléments finis robuste pour la simulation industrielle de procédés composites. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722091

HAL Id: hal-01722091

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

Submitted on 2 Mar 2018

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

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

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Méthode éléments finis robuste pour la simulation industrielle de procédés composites

Arnaud DEREIMS<fnref n="1" />,<fnref n="2" />, Sylvain DRAPIER<fnref n="2" />, Jean-Michel BERGHEAU<fnref n="3" />, Patrick DE LUCA<fnref n="4" />

**Résumé** — Les procédés d'infusion sont une alternative économique aux procédés classiques d'injection pour la réalisation de grandes pièces en composites. Cependant la mauvaise maitrise des propriétés finales de la pièce moulée est un frein à leur démocratisation au niveau industriel. Nous proposons un modèle complet pour la simulation du procédé Liquid Resin Infusion basé sur les travaux de P. Celle. L'approche repose sur un découpage du domaine en trois zones (drainant, préformes imprégnées, préformes sèche), consistant, ainsi, à coupler un écoulement de Stokes dans le drainant à un écoulement de Darcy dans les préformes. De plus, l'influence mutuelle de la résine sur le comportement des préformes et de la déformation des préformes sur la perméabilité sont considérées, à travers la loi de Terzaghi et des lois exprimant la perméabilité en fonction de la fraction de fibres. Enfin, le procédé est découpé en trois phases : compaction des préformes sèches, remplissage et « post-infusion ». Nous proposons une approche numérique applicable à des cas réels d'infusion souvent mis de coté dans les publications récentes, impliquant des perméabilités très faibles, une drainant fin (1mm) et des géométries complexes.

**Mots clés** — Infusion, composite, LRI, post-infusion, Stokes, Darcy, grandes déformations, interaction fluide/structure.

## 1 Introduction

Les procédés de mise en œuvre des matériaux composites par infusion de résine sont développés depuis quelques années pour apporter des solutions économiques à la réalisation de grandes pièces. Ces procédés permettent une réduction des coûts de stockage des matières premières et de fabrication de moule, du temps d'infusion et limitent la formation de vides. Cependant, la mauvaise maitrise des propriétés finales de la pièce, induisant une mise au point du procédé longue et couteuse, réduit considérablement les avantages précités. Il apparaît donc nécessaire de développer des méthodes de simulation numérique permettant d'anticiper ces propriétés.

## 2 Description du procédé LRI

Le procédé LRI compte parmi les procédés d'infusion les plus diffusés au niveau industriel. Dans ce procédé, une couche de résine liquide, recouvrant toute la surface de la pièce à infuser, est formée à l'aide d'un tissu fortement perméable, le drainant, placé au dessus de l'empilement de fibres. C'est ensuite la différence de pression entre l'entrée de résine (pression atmosphérique) et l'évent (pression nulle) qui provoque l'infusion de la résine dans le drainant puis à travers l'épaisseur des préformes (voir figure 1(a)). Après remplissage, l'ensemble subit un cycle de cuisson qui conduit la résine dans un état solide.

La souplesse du sac à vide ne permet pas de maintenir une épaisseur constante pendant toute la durée de l'opération, or de cette épaisseur dépendent la qualité dimensionnelle et les propriétés mécaniques de la pièce fabriquée. Nous présentons dans la suite une approche isotherme pour simuler au mieux ce procédé.

<footnote n="1">¹ ESI Group, Le Récamier, 70 rue Robert, 69458 Lyon Cedex 6, France, arnaud.dereims@esi-group.com</footnote>

<footnote n="2">² École des Mines de Saint-Étienne, centre Science des Matériaux et Structures, Laboratoire Georges Friedel UMR CNRS 5307, 158 cours Fauriel, 42023 Saint-Étienne Cedex 02, France, drapier@emse.fr</footnote>

<footnote n="3">³ Université de Lyon, École Nationale d'Ingénieurs de Saint-Étienne, LTDS UMR CNRS 5513, 58 rue Jean Parot, 42023 Saint-Étienne Cedex 02, France, bergheau@enise.fr</footnote>

<footnote n="4">⁴ ESI Group, Aero Business Center - Aeroparc, 25 Rue Marcel Issartier, BP 20005, 33700 Mérignac, patrick.de.luca@esi-group.com</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 – Découpage spatial du domaine d'étude</figcaption></figure>

## 3 Modélisation physique et aspects numériques

En se référant à la description du procédé, on peut remarquer, la présence de deux milieux différents, le drainant et les renforts, partiellement imprégnés par un fluide, la résine. Ainsi, Pierre Celle propose de découper le domaine d'étude en trois zones distinctes (figure 1(b)), dont les frontières vont évoluer au cours du temps [1, 2] :

- une zone d'écoulement rapide constituée de la résine et du drainant (Remplissage), que nous appelons la « zone de Stokes » ;

- une zone d'écoulement plus lent constituée du renfort imprégné de résine (Remplissage et Post-Filling), que nous appelons « zone de Darcy » ;

- une zone constituée des préformes non encore imprégnées de résines (Compression initiale, Remplissage), que nous appelons « zone sèche ».

En se basant sur ce découpage, nous allons maintenant nous attacher à présenter les modèles physiques, ainsi que les méthodes numériques mises en place pour les résoudre.

### 3.1 Compression initiale

Initialement, il n'y a pas de résine dans le système, seule la préforme sèche est considérée. Lorsque le vide est tiré, la préforme est soumise à la pression atmosphérique (1 bar) et se déforme. Les contraintes et déformations ainsi induites sont à l'origine des variations de perméabilités à prendre en compte lors de la phase de remplissage. En négligeant les effets d'inertie, nous pouvons écrire l'équation de conservation de la quantité de mouvement comme suit :

$$ \mathrm{div}\underline{\underline{\sigma}}(\underline{u}) = \underline{f}_v \quad (1) $$

où $\underline{\underline{\sigma}}$ est le tenseur des contraintes de Cauchy, $\underline{u}$ le champ de déplacement et $\underline{f}_v$ les forces de volume. Numériquement, nous utilisons une formulation lagrangienne réactualisée à chaque itération (2). Cette approche aussi appelée formulation quasi eulérienne fait apparaître les grandeurs eulériennes et nous permet, lors du couplage fluide / structure, d'imposer directement la pression hydrostatique de résine comme chargement interne sans la transporter sur la configuration de référence.

$$ \left\{
\begin{array}{l}
\text{On cherche } \underline{u} \in E_{\underline{u}} / \forall \underline{v}^* \in E_{\underline{v}^*} \\
\displaystyle\int_{\Omega_t} \underline{D}^* : \underline{\underline{\sigma}}(\underline{u}) \, dv = \displaystyle\int_{\partial\Omega_t^N} \underline{v}^* \underline{\underline{\sigma}}(\underline{u}) \, n \, ds + \displaystyle\int_{\Omega_t} \underline{v}^* \cdot \underline{f}_v \, dv \quad \forall \underline{x} \in \Omega_t \, dv \\
\\
\underline{\underline{\sigma}}\underline{n} = \underline{t} \qquad \forall \underline{x} \in \partial\Omega_t^N \\
\\
\underline{E}_{\underline{u}} = \{\underline{u} \in H^1(\Omega_t) \mid \underline{u} = \underline{u}^d \quad \forall \underline{x} \in \partial\Omega_t^D\} \\
\\
\underline{E}_{\underline{v}^*} = \{\underline{v}^* \in H^1(\Omega_t) \mid \underline{v}^* = \underline{0} \quad \forall \underline{x} \in \partial\Omega_t^D\}
\end{array}
\right.
\quad (2) $$

Dans l'équation (2), $\underline{D}^*$ est le tenseur des taux de déformations virtuelles eulérien, $\underline{v}^*$ un champ de vitesse virtuelle cinématiquement admissible, $\underline{\Omega}_t$ la dernière configuration connue, $\partial\underline{\Omega}_t^N$ la frontière de normale sortante $\underline{n}$ sur laquelle est imposé un chargement et $\partial\underline{\Omega}_t^D$ la frontière sur laquelle est imposé un déplacement.

De plus, on considère la préforme comme un milieu homogène constitué de fibres indéformables. Cette hypothèse permet d'assimiler les variations de volume observées à l'échelle macroscopique à des
réarrangement de fibres à l'échelle microscopique et, ainsi, de relier directement les variations du jacobien de la transformation (J) aux variations de porosités (φ) selon l'équation suivante :

$$J(\underline{x}, t + \Delta t)(1 - \phi(\underline{x}, t + \Delta t)) = J(\underline{x}, t)(1 - \phi(\underline{x}, t)). \quad (3)$$

Enfin la préforme pouvant être soumise à de grandes déformations et son comportement étant généralement isotrope transverse élastique non-linéaire, un algorithme de Newton-Raphson est utilisé pour traiter les non-linéarité matérielles et géométriques.

## 3.2 Remplissage : une histoire de couplage

### 3.2.1 Écoulements de la résine : couplage Stokes/Darcy

Comme décrit plus haut, la résine s'écoule dans un tissu drainant puis dans les préformes sous l'action du vide. Dans la zone dite de Stokes, la très forte perméabilité du drainant, relativement à celle des fibres, nous permet de considérer, en première approximation, un écoulement incompressible régi par les équations de Stokes (4),

$$\left| \begin{array}{l} \mathrm{div}(2\mu \underline{D}(\underline{v})) - \underline{\nabla} p = 0 \\ \mathrm{div} \underline{v} = 0 \end{array} \right. \qquad (4)$$

où $\mu$ est la viscosité de la résine, $\underline{v}$ la vitesse de la résine, $p$ la pression de la résine et $\underline{D}(\underline{v}) = \frac{1}{2}(\underline{\nabla}\underline{v} + \underline{\nabla}^T \underline{v})$ le tenseur des taux de déformation eulérien associé au champ de vitesse $\underline{v}$. Tandis que dans la zone de Darcy, on représente un fluide Newtonien incompressible s'écoulant dans un milieu poreux à faible perméabilité (de l'ordre de $10^{-12} m^2$ pouvant descendre jusqu'à $10^{-14} m^2$) régi par la loi de Darcy (5),

$$\left| \begin{array}{l} \underline{v} = -\frac{1}{\mu} \underline{K} \underline{\nabla} p \\ \mathrm{div} \underline{v} = 0 \end{array} \right. \qquad (5)$$

avec $\underline{v}$ la vitesse de Darcy ou la vitesse relative moyenne de la résine par rapport aux préformes, $\underline{K}$ le tenseur de perméabilité, $\mu$ la viscosité de la résine et $p$ la pression interstitielle. Les deux domaines sont séparés par une interface $\Gamma = \Omega_S \cap \Omega_D$ à travers laquelle la conservation de la masse et la continuité de la contrainte normale doivent être respectées. De plus nous choisissons de contrôler la vitesse tangente sur l'interface à l'aide d'une condition de type Beavers-Joseph-Saffman en introduisant un coefficient de glissement ($\alpha$). La loi de Darcy impliquant que $\underline{\sigma}_{\underline{D}} = -p_D \underline{I}$, la continuité de la contrainte normale est équivalente à une condition de Dirichlet sur la pression hydrostatique, conduisant au système de conditions suivant sur $\Gamma$ :

$$\left| \begin{array}{l} \underline{v}_S \cdot \underline{n} = \underline{v}_D \cdot \underline{n} \\ \underline{v}_S \cdot \underline{\tau} = -\frac{1}{2\mu} \underline{K} \left( \frac{\lambda^2 + 2\alpha\lambda}{1 + \alpha\lambda} \right) \underline{\nabla} p \cdot \underline{\tau} \\ p_D = p_S \end{array} \right. \qquad (6)$$

où $\underline{n} = \underline{n}_S = -\underline{n}_D$ est la normale sortante de l'interface Stokes/Darcy, $\underline{\tau}$ le vecteur unitaire tangent à l'interface et $\lambda = \frac{H}{\sqrt{\underline{K}}}$ avec $H$ l'épaisseur de la zone de Stokes ($\Omega_S$). Dans la pratique, la vitesse normale prise dans la zone de Darcy, $\underline{v}_D \cdot \underline{n}$, est imposée par une méthode de pénalité dans les équations de Stokes, et la pression hydrostatique prise dans la zone de Stokes, $p_S$, est imposée dans les équations de Darcy.

Du point de vue numérique, même si la forme de l'équation de Darcy nous permet de résoudre séparément les champs de vitesse et de pression, afin, d'être conforme à la formulation utilisée pour le problème de Stokes en raison du couplage, on choisit de résoudre une formulation mixte en vitesse-pression (7) et (8) par une méthode éléments finis stabilisée à l'aide des mini-éléments P1+/P1 [2]. Les mini-éléments P1+/P1 s'appuient sur l'enrichissement de l'approximation du champs de vitesse par l'introduction d'un nœud supplémentaire situé au centre de l'élément (voir figure 2). Le degré d'interpolation du champs de vitesse est, par conséquent, supérieur à celui du champs de pression, ce qui permet de respecter la condition de Brezzi-Babuska et d'éviter l'apparition de modes de pression parasites [3].

La formulation mixte des équations de Stokes et Darcy peut prendre deux formes différentes qui font apparaître les conditions limites sur les pressions en tant que conditions de Dirichlet ou faiblement. Ces formulations dites primale ou duale prennent les formes suivantes.
<figure><img src="image_3.png" /><figcaption>Fig. 2 – mini-élément P1-bulle/P1 : (a) élément triangle, (b) représentation graphique de la fonction de forme associée à la bulle [1, 2]</figcaption></figure>

La formulation duale mixte en vitesse pression du système de Stokes s'écrit :

$$
\left|
\begin{array}{l}
\displaystyle \text{soit } (\underline{v}, p) \in E_{\underline{v}} \times L^2(\Omega) / \forall (\underline{v}^*, p^*) \in E_{\underline{v}^*} \times L^2(\Omega), \\
\displaystyle \int_{\Omega} 2\mu \underline{D}(\underline{v}^*) : \underline{D}(\underline{v}) \, dv - \int_{\Omega} p \, \text{div} \, \underline{v}^* \, dv = \underbrace{\int_{\partial \Omega_t} \mu \, \underline{v}^* \, \underline{\nabla} \underline{v} \, n \, ds - \int_{\partial \Omega_t} p \, \underline{v}^* \, \cdot n \, ds}_{\displaystyle \int_{\partial \Omega_t} \underline{v}^* \underline{\sigma} n \, ds} \\
\\
\displaystyle \int_{\Omega} p^* \, \text{div} \, \underline{v} \, dv = 0 \\
\\
\displaystyle E_{\underline{v}^*} = \{ \underline{v}^* \in H^1(\Omega) / \underline{v}^* = 0 \, \text{sur } \partial \Omega_{\underline{v}} \} \\
\displaystyle E_{\underline{v}} = \{ \underline{v} \in H^1(\Omega) / \underline{v} = \underline{v}^d \, \text{sur } \partial \Omega_{\underline{v}} \}
\end{array}
\right.
\qquad (7)
$$

La formulation primale mixte en vitesse pression du système de Darcy s'écrit :

$$
\left|
\begin{array}{l}
\text{soit } (\underline{v}, p) \in L^2(\Omega) \times E_p / \forall (\underline{v}^*, p^*) \in L^2(\Omega) \times E_{p^*}, \\
\displaystyle \int_{\Omega} \underline{v}^* \mu \underline{K}^{-1} \underline{v} \, dv + \int_{\Omega} \underline{v}^* \cdot \underline{\nabla} p \, dv = 0 \\
\\
\displaystyle \int_{\Omega} \underline{v} \cdot \underline{\nabla} p^* \, dv = \int_{\partial \Omega_v} p^* \, \underline{v} \, n \cdot ds \\
\\
\displaystyle E_{p^*} = \{ p^* \in H^1(\Omega) / p^* = 0 \, \text{sur } \partial \Omega_p \} \\
\displaystyle E_p = \{ p \in H^1(\Omega) / p = p^d \, \text{sur } \partial \Omega_p \}
\end{array}
\right.
\qquad (8)
$$

Le choix des formulations (7) et (8) s'explique pour deux raisons. D'une part, la formulation duale en vitesse-pression de Darcy n'est pas stable avec les éléments P1+/P1 et nécessite l'implémentation d'une méthode de stabilisation telle que la Hughes Varational Methods (HVM) [4], [5] et [6]. D'autre part, la formulation primale en vitesse-pression du problème de Stokes n'est pas adaptée à l'imposition d'une vitesse, or lors du couplage nous faisons le choix d'imposer la vitesse normale sur l'interface.

Numériquement, l'imposition des conditions cinématiques nécessite quelques précautions. D'une part, les très faibles perméabilités rencontrées dans la simulation de procédés composites, impliquent de très faibles vitesses de résine qui posent des problèmes numériques en ce qui concerne la formulation du problème purement fluide. En effet, la méthode de stabilisation employée (mini-éléments P1+/P1) repose sur la relaxation de la contrainte d'incompressibilité et est équivalente à l'équation suivante :

$$
\mathrm{div}\vec{v} + \beta p = 0 \tag{9}
$$

où $\beta$ est un terme de stabilisation proportionnel à $\frac{1}{\mu}$ pour Stokes et $\frac{1}{\mu}\underline{K}$ pour Darcy. Ainsi, la résine est transformer en un fluide quasi-incompressible et doit le rester afin de respecter la physique du problème traité. Cela signifie que $\beta p$ doit rester négligeable devant la vitesse et la pression, ce qui peut ne pas être le cas lors d'un couplage avec un écoulement très lent. Dans la pratique, on constate qu'en introduisant une viscosité « artificielle » dix fois supérieure à la viscosité réelle de la résine dans le terme de stabilisation ($\beta$), on s'assure de la quasi-incompressibilité jusqu'à des perméabilités de l'ordre de $10^{-15}m^2$ sans répercussion sur la qualité de la stabilisation ou les performances du solveur. D'autre part, la continuité de la vitesse étant imposée par une méthode de pénalité, il est judicieux de dimensionner le coefficient de pénalité en fonction de la perméabilité afin d'assurer la bonne pénalisation du système.
### 3.2.2 Validation du couplage Stokes/Darcy par la méthodes solutions manufacturées (MMS)

La bonne implémentation de l'algorithme de couplage Stokes/Darcy a été vérifiée à l'aide de la méthode des solutions manufacturées [8, 9]. La méthode consiste à construire une solution analytique du problème qui est ensuite introduite dans le système d'équations à résoudre pour en déduire un terme source. À partir de ce terme source le problème éléments finis est résolu, enfin on calcule l'erreur de la solution discrète par rapport à la solution analytique initiale.

Soit un domaine $\Omega = (0; 1) \times (0; 2)$ divisé en deux sous domaines $\Omega_S = (0; 1) \times (1; 2)$ régi par les équations de Stokes et $\Omega_D = (0; 1) \times (0; 1)$ régi par les équations de Darcy en contact sur une interface $\Gamma = \Omega_S \cap \Omega_D$. On choisi la solution analytique proposée par Discacciati [10] :

$$ 
\begin{cases}
\vec{v}_S \cdot \vec{x} = -\cos(\pi/2y) \sin(\pi/2x) \\
\vec{v}_S \cdot \vec{y} = \cos(\pi/2x) \sin(\pi/2y) - 1 + x \\
p_S = 1 - x
\end{cases}
\quad (10) $$

$$ 
\begin{cases}
\vec{v}_D \cdot \vec{x} = \sin(\pi/2x) \cos(\pi/2y) + y \\
\vec{v}_D \cdot \vec{y} = \cos(\pi/2x) \sin(\pi/2y) + x - 1 \\
p_D = 2/\pi \cos(\pi/2x) \cos(\pi/2y) - y(x-1)
\end{cases}
\quad (11) $$

Nous considérons des maillages réguliers construit à partir de carrés de dimensions $h \times h$ divisés en deux triangles. Chaque calcul est réalisé sur un maillage de plus en plus fin (h divisé par deux) en commençant avec $h = 0.1$ jusqu'à $h = 0.0125$. Pour simplifier on choisi $\mu = 1 Pa.s$, $K = 1 m^2$ comme proposé dans [10].

La convergence de la solution a été étudiée, la figure 3 montre que la pression et la vitesse convergent respectivement avec des taux de 2 et 1.4 dans le domaine Darcy et des taux de 1.8 et 2.2 dans le domaine Stokes. Ces résultats sont en accord avec les taux de convergence théoriques déterminés dans la littérature [10, 9].

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Taux de convergence dans les domaines de Stokes et Darcy : (a) vitesse, (b) pression</figcaption></figure>

### 3.2.3 Couplage fluide / Solide

Le couplage fluide/structure est traité, dans une approche quasi-statique, par un couplage fort itératif. En effet, le problème fluide couplé Stokes/Darcy est résolu sur un maillage fixe (i.e. pour des préformes indéformables). L'influence de la résine sur les fibres est prise en compte à travers la pression hydrostaticque selon la loi de Terzaghi (12),

$$ \underline{\underline{\sigma_{tot}}} = \underline{\underline{\sigma_{eff}}} + p\underline{\underline{I}} \quad (12) $$

avec $\underline{\underline{\sigma_{tot}}}$ la contrainte totale appliquée au système, $\underline{\underline{\sigma_{eff}}}$ la contrainte effective dans les préformes et $\underline{\underline{I}}$ le tenseur unité. Le problème non linéaire en grandes déformations est résolu pour un champs de pression fixée (et donc un comportement équivalent au sens de Terzaghi), puis, l'influence des déformations se traduit par l'évolution de la porosité qui influe sur la perméabilité du milieu, déterminée selon une loi
explicitant la perméabilité en fonction de la porosité, on peut par exemple citer la loi de Carman-Kozeny (13) :

$$ \underline{K} = \frac{d_f^2}{16 k_h} \frac{\phi^3}{(1-\phi)^2} \qquad (13) $$

avec $d_f$ le diamètre moyen de fibre, $k_h$ constantes de Kozeny et $\phi$ la porosité. Finalement, la convergence est vérifiée par un calcul d'erreur relative sur les champs de pression, vitesse et déformations entre deux itérations (Figure 4).

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Algorithme de couplage pour la simulation de l'infusion</figcaption></figure>

### 3.2.4 Avancée du front de fluide

L'évolution du front de fluide est gérée par l'algorithme de remplissage de PAM-RTM™. Cet algorithme repose sur un découpage du phénomène transitoire en une succession d'états quasi-statiques. Les nouveaux éléments remplis entre deux états ou temps de remplissage sont déterminés à partir du champ de vitesse calculé au pas de temps précédent. Un facteur de remplissage compris entre 0 et 1 est associé à chaque élément.

## 3.3 Écoulements post-infusion

À la fin du remplissage, le système n'est pas équilibré, la différence de pression entre l'évent et l'injection induit des écoulements que nous qualifierons d'écoulements post-infusion. Ces écoulements, peu étudiés dans la littérature, influent pourtant beaucoup sur la qualité finale de la pièce infusée. En effet, le gradient de pression final, associé à la souplesse du sac à vide, implique une répartition de la résine et une épaisseur non-homogènes. La résine va donc migrer petit-à-petit des zones à faible fraction de fibres vers les zones à forte fraction de fibres afin d'atteindre un équilibre mécanique.

La simulation de cette phase met en jeux les mêmes problèmes physiques que la phase de remplissage et repose sur le même genre de couplage. Cependant, les états quasi-statiques successifs sont plus délicats à déterminer, en effet, la pièce étant entièrement remplie, il n'y a plus de front de fluide à faire avancer. Ainsi, nous proposons de déterminer le nouvelle fraction volumique de résine ($\phi$) à partir d'un bilan de masse sur chaque élément fini, comme suit :

$$ \phi_{t+\Delta t} (V_t^e + \Delta t \sum S_i^e \vec{v}_{ti} \cdot \vec{n}_i^e) = V_t^e \phi_t + \Delta t \sum S_i^e \vec{v}_{ti} \cdot \vec{n}_i^e \qquad (14) $$

où $V_t^e$ est le volume de l'élément au temps $t$, $S_i^e \vec{v}_{ti} \cdot \vec{n}_i^e$ correspond au débit de résine passant à travers la face $i$ de l'élément au temps $t$. Connaissant la nouvelle répartition de résine dans le domaine, on peut déduire le jacobien de la transformation à $t + \Delta t$ à partir de l'équation (3). La variation de volume de l'élément n'étant due qu'à la variation du volume de résine contenu dans l'élément, nous faisons l'hypothèse que la déformation associée n'est due qu'à la variation de la pression interstitielle de résine et qu'elle est donc hydrostatique. Ceci nous permet, avec l'équation (15), de déterminer un incrément de déformation équivalent $\Delta\epsilon_{ii}$ permettant de déduire l'incrément de pression associé à la migration de la résine dans la pièce.

$$ \prod (1 + \epsilon_{ii} + \Delta\epsilon_{ii}) = J_{t+\Delta t} \qquad (15) $$
avec $\epsilon_{ii}$ l'état déformation courant (connu). Finalement, en appliquant la loi de Darcy (5), la nouvelle vitesse de la résine est calculée. Le processus est réitéré jusqu'à ce qu'il n'y ai plus d'écoulement résiduel ($v = 0$).

## 4 Application

Afin d'illustrer le potentiel de l'outil développé dans le cadre de nos travaux, nous présentons, un cas d'infusion d'une pièce complexe en trois dimensions qui correspond à une géométrie typique de pièces aéronautiques telles que des pales de rotor.

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Exemple 3D : maillage et délimitation des zones de Stokes et de Darcy</figcaption></figure>

La préforme est supposée être constituée d'un empilement d'unidirectionels de type UD G1157 produit par Hexcel reinforcement, dont le comportement a été étudié et caractérisé dans les travaux de P. Wang [7]. La géométrie correspond à un secteur angulaire d'épaisseur initiale 30mm.

La perméabilité isotrope du matériau peut être modélisée par la loi de Carman-Kozeny :

$$K_{xx} = K_{yy} = K_{zz} = 1.5625 \cdot 10^{-13} (1 - V_f)^3 / V_f^2 \quad (16)$$

avec $V_f$ la fraction de fibres, la viscosité de la résine RTM6 est de 0.058 Pa.s à 120°C. L'épaisseur du drainant est de 1 mm (en vert sur la figure 5). Les conditions limites sont décrites sur la figure 5.

La simulation a été réalisée avec un maillage de 4 497 éléments, 76 pas de temps ont été nécessaires pour remplir la préforme. Le calcul, réalisé sur un processeur, a duré 40 minutes en tout et a consommé 60Mo de RAM. La figure 6(a) montre le taux de remplissage finale de la pièce, le temps de remplissage obtenu est de 1 092 seconds. On remarquera la présence d'une zone sèche en fin de remplissage due à la compaction initiale importante dans les zones de faibles courbures. La figure 6(b), montre les variations d'épaisseur en fonction du temps, on peut identifier trois points clefs :

- Au point ①, la résine atteint les points de mesures, l'épaisseur commence à augmenter,

- Au point ②, après une phase d'augmentation linéaire de l'épaisseur, la résine atteint le bas de la préforme et la pression augmente brutalement impliquant une forte augmentation de l'épaisseur,

- Au point ③, la résine termine d'imprégner la préforme petit-à-petit et la pression interstitielle augmente doucement jusqu'à la fin du remplissage.

## 5 Conclusion et travail futur

Nous avons présenté un modèle innovant pour la simulation des procédés d'élaboration de pièces composites par infusion de résine. L'originalité de cette approche tient en sa capacité à traiter des cas réels d'infusion souvent mis de coté dans les travaux récents. Pour cela nous avons couplé les équations de Stokes et de Darcy à la mécanique du solide en grandes déformations. Les conditions de couplages spécifiques à la prise en compte de paramètres physiques réels ont été présentés et validés. Enfin, Une validation expérimentale complète est actuellement menée avec des partenaires industriels.
<figure><img src="image_7.png" /><figcaption>Fig. 6 – Exemple de simulation 3D : (a) Fraction volumique de résine en fin de remplissage, (b) Épaisseur en fonction du temps</figcaption></figure>

Cependant, pour répondre pleinement aux besoins des industriels des améliorations doivent être apportées. Ainsi, un couplage avec la thermochimie doit être introduit afin de simuler les procédés non-isothermes. De plus d'après les expériences réalisées sur les préformes, un comportement élasto-plastique du milieu homogène équivalent semble plus réaliste.

# 6 remerciement

Les auteurs remercient l'union européenne pour le soutien du projet INFUCOMP (European Community's Seventh Framework Program FP7/2009-2013 under grant agreement n°233926) qui finance une partie des recherches présentées ici.

# Références

[1] P. Celle. *Couplage fluide / milieu poreux en grandes déformations pour la modélisation des procédés d'élabo-ration par infusion*, Thèse de doctorat, École Nationale des Mines de Saint-Étienne, 2006.

[2] P. Celle, S. Drapier, J-M. Bergheau. *Numerical modelling of liquid infusion into fibrous media undergoing compaction*, European Journal of Mechanics - A/Solids, 27(4): 647-661, 2008.

[3] D.N. Arnold, F. Brezzi et M. Fortin. *A stable finite element for the Stokes equations*, Calcolo, 21(4): 337-344, 1984.

[4] A. Massud. *A stabilized mixed finite element method for Darcy-Stokes flow*, International journal for numerical methods in fluids, 54: 665-681, 2007.

[5] G. Pacquaut. *Couplage Stokes/Darcy dans un cadre Level-set en grandes déformations pour la simulation des procédés d'élaboration par infusion de résine*, Thèse de doctorat, École Nationale des Mines de Saint-Étienne, 2010.

[6] G. Pacquaut, J. Bruchon, N. Moulin et S. Drapier. *Combining a level set method and a mixed stabilized P1/P1 formulation for coupling Stokes-Darcy flows*, International Journal for Numerical Methods in Fluids, accepté pour publication.

[7] P. Wang, S. Drapier, J. Molimard, A. Vautrin et J-C. Minni. *Characterization of Liquid Resin Infusion (LRI) filling by fringe pattern projection and in situ thermocouples*, Composites Part A : Applied Science and Manufacturing, 41(1): 36-44, 2010.

[8] P. Knupp, K. Salari. *Verification of Computer Codes in Computational Science and Engineering*, Chapman & Hall/CRC, 2003.

[9] R. Troian, S. Drapier. *Efficient decoupled approach of stokes-Darcy layers for LCM process modelling*, Soumis dans International Journal of Material Forming.

[10] M. Discacciati. *Domain decomposition method for the coupling of surface and groundwater flows*, Thèse de doctorat, École Polytechnique Fédéral de Lausane (EPFL), 2004.