# Utilisation d'une méthode Level Set pour le couplage Stokes-Darcy : application aux procédés d'élaboration par infusion de résine

Guillaume Pacquaut, Julien Bruchon, Nicolas Moulin, Sylvain Drapier

► To cite this version:

Guillaume Pacquaut, Julien Bruchon, Nicolas Moulin, Sylvain Drapier. Utilisation d'une méthode Level Set pour le couplage Stokes-Darcy : application aux procédés d'élaboration par infusion de résine. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413706

HAL Id: hal-01413706

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

Submitted on 10 Dec 2016

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Utilisation d'une méthode Level Set pour le couplage Stokes-Darcy : application aux procédés d'élaboration par infusion de résine

G. Pacquaut¹, J. Bruchon¹, N. Moulin¹, S. Drapier¹

¹ Ecole Nationale Supérieure des Mines de Saint-Etienne
Centre Sciences des Matériaux et des Structures & LTDS UMR CNRS 5513
158, cours Fauriel
42023 Saint-étienne cedex 2 - France
pacquaut@emse.fr, bruchon@emse.fr, nicolas.moulin@emse.fr, drapier@emse.fr

**Résumé** — Le but de ce travail est de simuler numériquement les procédés d'élaboration de composites structuraux par infusion de résine. Les équations de Stokes-Darcy sont couplées et résolues par une méthode éléments finis mixte, stabilisée, avec un maillage non structuré. Sur l'interface entre les milieux fluides et poreux, représentée par une fonction Level set, une condition de Beaver-Joseph-Saffman est appliquée. L'avancée de la résine est également modélisée par une approche Level Set. Nous présentons dans cet article une validation de la formulation implémentée pour le couplage Stokes-Darcy grâce à une solution analytique ainsi qu'une simulation en 2D du procédé d'infusion.

**Mots clés** — Infusion, Stokes, Darcy, éléments finis, Level Set.

# 1 Position du problème

## 1.1 Contexte de l'étude

Les matériaux composites trouvent de plus en plus d'applications dans de nombreux secteurs industriels tels que le transport, le bâtiment, l'aérospatial ainsi que dans les sports et les loisirs, notamment grâce à des procédés dits par voie sèche apparus ces dernières années. Ces procédés sont utilisés pour élaborer des matériaux composites à matrice organique. Leur particularité est de réaliser le mélange résine/fibres au cours de l'élaboration des structures réduisant ainsi les coûts.

Le but de ce travail est de simuler numériquement les procédés d'élaboration par infusion de résine en utilisant la méthode des éléments finis. Ces procédés peuvent se modéliser schématiquement de la façon suivante : un moule contient initialement de la résine (une partie purement fluide) et des préformes (un milieu poreux). Sous l'effet d'une pression mécanique extérieure, appliquée sur l'ensemble résine+préformes, la résine infuse dans les préformes, elles-mêmes sujettes à de grandes déformations.

Cet article se concentre sur la modélisation et la simulation du couplage de l'écoulement dans une zone purement fluide, décrit par les équations de Stokes (1), et de l'écoulement dans un milieu poreux, décrit par les équations de Darcy (2).
## 1.2 Couplage Stokes-Darcy

Notons $\Omega_s \subset \mathbb{R}^d$ et $\Omega_d \subset \mathbb{R}^d$ ($d=2$ ou $3$) deux domaines de frontières $\partial\Omega_s$ et $\partial\Omega_d$ tels que l'interface $\Gamma$ entre les deux domaines soit définie par : $\Gamma = \partial\Omega_s \cap \partial\Omega_d$. Un fluide Newtonien de viscosité dynamique $\eta$ occupe le domaine purement fluide $\Omega_s$ et s'écoule dans un milieu poreux $\Omega_d$. Dans $\Omega_s$, l'écoulement est gouverné par les équations de Stokes :

$$ \left\{ \begin{array}{l} -\overrightarrow{div}(2\eta\dot{\bar{\epsilon}}(\vec{v}_s)) + \overrightarrow{grad} p_s = \vec{0} \\ \operatorname{div} \vec{v}_s = 0 \end{array} \right. \qquad (1) $$

et dans $\Omega_d$, l'écoulement est gouverné par les équations de Darcy :

$$ \left\{ \begin{array}{l} \frac{\eta}{K} \vec{v}_d + \overrightarrow{grad} p_d = \vec{0} \\ \operatorname{div} \vec{v}_d = 0 \end{array} \right. \qquad (2) $$

L'approche adoptée ici est une approche macroscopique. Contrairement aux travaux précédents sur le même sujet [1], l'interface entre le milieu purement fluide et le milieu poreux est décrite de manière Eulérienne, via l'utilisation d'une méthode Level Set.

Des conditions de continuité à l'interface sont prises en compte dans le modèle. Ces conditions sont la continuité de la vitesse due à la conservation de la masse (3a), la continuité de la contrainte normale (3b) et la condition de Beaver-Joseph-Saffman [4] due aux mécanismes de frottement entre les filets fluides (3c). Cette dernière condition impose une vitesse tangentielle non nulle à l'interface avec le milieu poreux :

$$ \vec{v}_s \cdot \vec{n} = \vec{v}_d \cdot \vec{n} \qquad (3a) $$

$$ 2\vec{n} \cdot \dot{\bar{\epsilon}}(\vec{v}_s) \cdot \vec{n} = p_s - p_d \qquad (3b) $$

$$ 2\vec{n} \cdot \dot{\bar{\epsilon}}(\vec{v}_s) \cdot \vec{\tau} = \frac{-\alpha}{\sqrt{K}} \vec{v}_s \cdot \vec{\tau} \qquad (3c) $$

où $v_s$ et $v_d$ sont respectivement la vitesse du fluide dans le domaine de Stokes et dans le domaine de Darcy, $p_s$ et $p_d$ la pression dans le domaine de Stokes et dans le domaine de Darcy, $K$ la perméabilité, $\dot{\bar{\epsilon}}$ le tenseur des taux de déformation, $\alpha$ le coefficient de glissement, $\vec{n}$ la normale extérieure au domaine de Stokes et $\vec{\tau}$ le vecteur tangent à l'interface.

# 2 Résolution numérique

## 2.1 Formulation faible

La formulation faible des équations de Stokes (1) s'écrit : Trouver $\vec{v}_s \in H_{\Gamma_s}^1(\Omega_s)$ et $p_s \in L^2(\Omega_s)$, $\forall \delta\vec{v}_s \in H_{\Gamma_s}^1(\Omega_s)$ et $\forall \delta p_s \in L^2(\Omega_s)$, tel que :

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega_s} 2\eta \dot{\bar{\epsilon}}(\vec{v}_s) : \dot{\bar{\epsilon}}(\delta\vec{v}_s) dv - \int_{\Omega_s} p_s \operatorname{div} \delta\vec{v}_s dv + \underbrace{\displaystyle \int_{\Gamma} \frac{\alpha\eta}{\sqrt{K}} (\vec{v}_s \cdot \vec{\tau})(\delta\vec{v}_s \cdot \vec{\tau}) ds}_{C} + \underbrace{\displaystyle \int_{\Gamma} p_s \delta\vec{v}_s \cdot \vec{n} ds}_{A} = 0 \\ \\ - \int_{\Omega} \operatorname{div} \vec{v}_s \delta p_s dv = 0 \end{array} \right. \qquad (4) $$

avec :

$$ H_{\Gamma_s}^1(\Omega_s) = \{ \mathbf{w} \in H^1(\Omega_s) \mid \mathbf{w} = 0 \text{ sur } \Gamma_s \} \text{ avec } \Gamma_s = \partial\Omega_s \setminus (\Gamma \cup \Gamma_{s,N}) $$

$\Gamma_{s,N}$ : partie de la frontière où une condition de Neumann est appliquée.
La formulation faible des équations de Darcy (2), s'écrit sous la forme suivante : Trouver $\vec{v}_d \in H_{\Gamma_d}(div, \Omega_d)$ et $p_d \in L^2(\Omega_d)$, $\forall \delta\vec{v}_d \in H_{\Gamma_d}(div, \Omega_d)$ et $\forall \delta p_d \in L^2(\Omega_d)$, tel que :

$$ egin{cases} 
\displaystyle \int_{\Omega_d} \frac{\eta}{K} \vec{v}_d \cdot \delta \vec{v}_d \, dv - \int_{\Omega_d} p_d \, div \, \delta \vec{v}_d \, dv - \underbrace{\int_{\Gamma} p_d \, \delta \vec{v}_s \cdot \vec{n} \, ds}_{B} = 0 \\ \\ - \int_{\Omega} div \, \vec{v}_d \, \delta p_d \, dv = 0 
\end{cases} 
\quad (5) $$

avec :

$$ H_{\Gamma_d}(div, \Omega_d) = \{ w \in H(div, \Omega_d) \mid w \cdot n = 0 \text{ sur } \Gamma_d \} \text{ avec } \Gamma_d = \partial \Omega_d \setminus \Gamma $$

## 2.2 Discrétisation

Nous nous donnons une discrétisation du domaine de calcul, à savoir un maillage $\mathcal{T}(\Omega)$, constitué de simplexes K (triangles en 2D et tétraèdres en 3D). L'ensemble du domaine de calcul est donc discrétisé à l'aide d'un maillage non structuré.

Lors de la résolution du système couplé nous considérons un champ de vitesse $\vec{v}$ unique et continu sur $\Omega$ et un champ de pression $p$ unique et continu sur $\Omega$ tel que $\vec{v}|_{\Omega_i} = \vec{v}_i$ avec $i = s$ ou $d$. Ces deux champs sont solutions d'une formulation mixte en vitesse-pression définie sur tout le domaine obtenue en ajoutant les équations (4) et (5). Dans cette formulation les termes A et B s'annulent. Afin que la résolution de ce système devienne effective, les intégrales sur $\Omega_s$ et $\Omega_d$ doivent être re-définie sur tout $\Omega$. Pour cela on introduit une fonction Heaviside H (définie au paragraphe 2.3) tel que : $\int_{\Omega_i} (...) d\Omega_i = \int_{\Omega} (...) * H_i d\Omega$ ($i = s$ ou $d$). Dans Stokes, nous utilisons l'élément P1+/P1. Dans Darcy, cet élément n'étant pas stable pour la formulation utilisée [2], nous adoptons une formulation P1/P1 stabilisée à l'aide d'une formulation multi-échelles [3].

## 2.3 Représentation de l'interface entre le milieu purement fluide et le milieu poreux

Nous avons utilisé une fonction Level Set pour décrire l'interface entre le milieu purement fluide et le milieu poreux. L'interface est donc représentée par la ligne de niveau zéro d'une fonction régulière $\phi$ (fonction distance signée). Cette fonction détermine pour chaque nœud du maillage la distance de ce nœud à l'interface. Le signe placé devant la distance permet de savoir si ce nœud se situe sur l'interface, ou d'un côté ou de l'autre de l'interface. La fonction distance nous permet de calculer la fonction Heaviside H (définie au paragraphe 2.2) continue sur $\Omega$, valant 1 sur un domaine, 0 sur l'autre et variant linéairement sur une épaisseur $\varepsilon$ à l'interface.

L'interface n'est pas définie de manière explicite et traverse donc les éléments du maillage. Ceci implique que la condition de Beaver-Joseph-Saffman, appliquée sur l'interface (terme C dans l'équation variationnelle (4)) doit être exprimée à l'aide d'une intégrale de volume grâce à la fonction Level Set $\phi$. Pour effectuer cela, on applique une approximation volumique d'une mesure surfacique qui consiste à utiliser une fonction une fois dérivable $\zeta(r)$, à support dans [-1,1], tel que $r \to \frac{1}{\varepsilon}\zeta(\frac{r}{\varepsilon})$ converge vers $\delta_0$ (fonction Dirac) au sens des distributions et $\int_{\mathbb{R}} \zeta = 1$. On a lorsque $\varepsilon \to 0$ :

$$ \frac{1}{\varepsilon}\zeta(\frac{\phi}{\varepsilon})||\vec{\nabla}\phi|| \to \delta_{\{\phi=0\}} $$

Et donc pour une fonction régulière $f$, on fait l'approximation suivante lorsque $\varepsilon \to 0$ :

$$ \int_{\Gamma} f(x) \, ds \approx \int_{\Omega} \frac{1}{\varepsilon} \zeta(\frac{\phi}{\varepsilon}) ||\vec{\nabla} \phi|| f(x) \, dv $$

avec : $f(x) = \frac{\alpha\eta}{\sqrt{K}}(\vec{v}_s \cdot \vec{\tau})(\delta\vec{v}_s \cdot \vec{\tau})$ et $\zeta(\frac{\phi}{\varepsilon}) = \begin{cases} \frac{1}{2}(1 + \cos(\pi\frac{\phi}{\varepsilon})) & \text{si } -\varepsilon < \phi < \varepsilon (\varepsilon = 1,5 \times \text{taille de maille}) \\ 0 & \text{sinon} \end{cases}$
## 2.4 Représentation du front fluide

Le front fluide est également réprésenté par une fonction Level Set $\phi_f$. Cette fonction est transportée par le champ de vitesse de la résine, grâce à la résolution d'une équation de transport (6).

$$ \frac{\partial \phi_f}{\partial t} + \vec{v} \cdot \overrightarrow{\text{grad} \phi_f} = 0 \qquad (6) $$

Cette équation de transport permet de faire évoluer la fonction $\phi_f$ et donc de connaître l'évolution du front de résine en fonction du temps. Cette équation a été résolue par une méthode éléments finis continue stabilisée (SUPG). Une étape de réinitialisation est nécessaire pour éviter que la solution se dégrade. Elle permet de récupérer la propriété de distance signée ($\|\overrightarrow{\text{grad} \phi_f}\| = 1$) qui garantit la régularité de la fonction $\phi_f$.

# 3 Validation de la formulation implémentée

## 3.1 Simulation d'un écoulement horizontal

Afin de valider la formulation implémentée sous ZéBuLoN, nous avons effectué des simulations d'un écoulement horizontal. La figure 1 montre les résultats obtenus. Pour ce cas test, un gradient de pression est appliqué suivant l'axe $\vec{x}$ sur toute la hauteur de la pièce. Les autres conditions aux limites sont une vitesse nulle sur le bord supérieur de la pièce et une vitesse normale nulle sur le bord inférieur. Ce cas test permet de vérifier que la vitesse tangentielle à l'interface entre le

<figure><img src="image_2.png" /><figcaption>Figure 1 – Ecoulement horizontal : (a) résultats obtenus pour la pression et (b) résultats obtenus pour la vitesse. La position de l'interface se situe en y=0.</figcaption></figure>

milieu purement fluide et le milieu poreux n'est pas nulle (d'après la condition de Beaver-Joseph-Saffman). Nous avons également étudié l'influence de certains paramètres comme la perméabilité K et le coefficient de glissement $\alpha$. La figure 2 (a) montre que plus la perméabilité augmente et plus la vitesse de la résine augmente. Ce qui est normal puisque la perméabilité correspond schématiquement au diamètre moyen des canaux dans lesquels la résine avance lors de l'infusion et donc plus la perméabilité est élevée et plus ces canaux seront importants et faciliteront l'écoulement de la résine. La figure 2 (b) montre l'influence du coefficient de glissement $\alpha$ à l'interface. Plus la valeur de ce coefficient augmente et plus la vitesse de la résine diminue.

Le cas test de l'écoulement horizontal a une solution analytique qui permet de valider l'implémentation numérique puisque : $v_x = -\frac{K}{2\eta}(\frac{\sigma^2+2\alpha\sigma}{1+\alpha\sigma})\frac{dP}{dx}(1+\frac{\alpha}{\sqrt{K}}y)+\frac{1}{2\eta}(y^2+2\alpha y\sqrt{K})\frac{dP}{dx}$ avec $\sigma = \frac{h}{\sqrt{K}}$ où h est la hauteur du milieu purement fluide.

Cette solution analytique est comparée à la solution numérique. La figure 3 montre une bonne corrélation entre les résultats analytiques et numériques dans le domaine purement fluide.
<figure><img src="image_3.png" /><figcaption>Figure 2 – (a) Amplitude de la vitesse dans le milieu purement fluide et dans le milieu poreux pour différentes valeurs de la perméabilité $K$ avec $\eta = 1$ Pa.s et $\alpha = 1$. (b) Amplitude de la vitesse dans le milieu purement fluide et dans le milieu poreux pour différentes valeurs de $\alpha$ avec $\eta = 1$ Pa.s et $K = 10^{-2} m^2$. La position de l’interface est située en $y=0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3 – Comparaison des résultats analytiques et numériques pour la vitesse suivant la direction x. La solution analytique est représentée en trait fin et la solution numérique en pointillé.</figcaption></figure>

## 3.2 Simulation d'un écoulement vertical

Nous avons également effectué des simulations d'un écoulement vertical. La figure 4 montre les résultats obtenus. Pour ce cas test, le gradient de pression est appliqué suivant l'axe $\bar{y}$. Les autres conditions aux limites sont une vitesse normale nulle sur les bords gauche et droit de la pièce. On a vérifié sur ce cas test la conservation du débit. En effet, la vitesse dans le domaine de Darcy est plus élevée que dans le domaine de Stokes ($v_d = v_s S_s / S_d$ avec : $S_d$ la surface dans Darcy, $S_s$ la surface dans Stokes et $S_d < S_s$).

<figure><img src="image_5.png" /><figcaption>Figure 4 – Écoulement vertical avec $K = 10^{-2} m^2$, $\eta = 1$ Pa.s et $\alpha = 0,1$. (a) résultats obtenus pour la pression et (b) résultats obtenus pour la vitesse. La position de l'interface se situe en $y=0$.</figcaption></figure>
# 4 Simulation numérique de l'infusion en 2D

Nous présentons ici les premiers résultats obtenus sur un cas réel d'infusion de résine (figure 5). Pour cette simulation, on impose une pression en entrée (au niveau du canal d'injection) et une vitesse normale nulle sur les autres bords du domaine. La simulation montre que l'infusion dans le milieu poreux se produit bien avant le remplissage complet du domaine fluide (figure 5(b)) et montre également des phénomènes d'emprisonnement de l'air initialement présent dans le domaine.

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

Figure 5 – Simulation sous ZéBuLoN du procédé d'infusion de résine (la ligne rouge correspond à l'interface entre le milieu purement fluide (haut) et le milieu poreux (bas) et la ligne noire correspond au front fluide. Les isovaleurs représentées sont celles de la vitesse.

# 5 Conclusions et perspectives

L'objectif de cet article était de présenter une utilisation d'une méthode Level Set pour le couplage Stokes-Darcy afin de modéliser l'interface entre un domaine purement fluide et un milieu poreux. Elle a également été utilisée pour suivre l'évolution du front fluide au cours du temps grâce à la résolution d'une équation de transport.

Concernant les perspectives à donner à ce travail, nous souhaitons prendre en compte la déformation élastique des préformes dont le comportement dépend de la résine (modèle de Terzaghi).

# Références

[1] P. Celle, S. Drapier, J.M. Bergheau Numerical modelling of liquid infusion into fibrous media undergoing compaction. *European Journal of Mechanics A/Solids*, 647–661, 2008.

[2] A. Masud, T.J.R. Hugues A stabilised mixed finite element method for Darcy flow. *Computer Methods in Applied Mechanics and Engineering*, 4341–4370, 2002.

[3] K.B. Nakshatrala, D.Z. Turner, K.D. Hjelmstad, A. Masud. A stabilized mixed finite element method for Darcy flow based on a multiscale decomposition of the solution. *Computer Methods in Applied Mechanics and Engineering*, 195:4036–4049, 2006.

[4] J.M. Urquiza, D. N'Dri, A. Garon, M.C. Delfour. Coupling Stokes and Darcy equations. *Applied Numerical Mathematics*, 58 (1):525–538, 2008.