# Définition Microstructurale du Mélange : Évaluation de l'Impact des Relations de Fermeture en Écoulement Complexe

Guillaume Maitrejean, Amine Ammar, Francisco Chinesta

► To cite this version:

Guillaume Maitrejean, Amine Ammar, Francisco Chinesta. Définition Microstructurale du Mélange : Évaluation de l'Impact des Relations de Fermeture en Écoulement Complexe. Rhéologie, 2011, 20, pp.44-52. hal-01006718

HAL Id: hal-01006718

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

Submitted on 18 Mar 2017

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
# Définition microstructurale du mélange : évaluation de l'impact des relations de fermeture en écoulement complexe

G. Maîtrejean¹, A. Ammar² et F. Chinesta³

¹Laboratoire de Rhéologie, UJF-INPG-CNRS
1301, Rue de la piscine, BP 53, Domaine universitaire, 38041 Grenoble Cedex 9

²Arts et Métiers ParisTech
2, Boulevard du Ronceray, BP 93525, 49035 Angers Cedex 01

³EADS Corporate International Chair
Ecole Centrale de Nantes, 1, Rue de la Noe, BP 92101 44321 Nantes Cedex 3

**Résumé** : Le mélange de fluides immiscibles est fréquemment utilisé dans de nombreuses applications industrielles. Une approche directe de la représentation de cette interface est impossible du fait de la différence de taille caractéristique entre la microstructure et le procédé lui-même. Pour cette raison, le mélange est décrit à l'aide du tenseur d'aire qui est une variable microstructurale. Cependant, le calcul du tenseur d'aire nécessite l'introduction d'une relation de fermeture dont nous évaluerons l'impact dans un écoulement complexe.

**Mots-clé** : Mélange liquide-liquide, Tenseur d'aire, Proper Generalized Decomposition.

## 1. Introduction

Lorsque deux fluides immiscibles sont mélangés, les propriétés du mélange sont fortement dépendantes de sa microstructure. Cependant, la description du mélange n'est pas possible à l'aide des seules variables hydrodynamiques classiques car le mélange contient une interface qui évolue en temps. Une approche directe de la représentation de cette interface est impossible du fait de la différence de taille caractéristique entre la microstructure et le procédé lui-même. Ainsi, au lieu de formuler le problème en fonction des gouttes, Doi et Ohta [1] ont considéré l'interface entre les deux phases et ont développé une équation cinétique semi-phénoménologique décrivant l'aire d'interface par unité de volume, ainsi que son anisotropie, pour un champ de vitesse associé. Cette première approche dénommée *tenseur d'interface* a ensuite été reprise et améliorée par Wetzel et Tucker [2], définissant alors le *tenseur d'aire*.

Cependant, l'utilisation du tenseur d'aire est sujette à l'introduction d'une relation de fermeture dont l'impact sur la solution ne peut être évalué a priori. Pour éviter l'utilisation d'une relation de fermeture, il est possible de résoudre le tenseur d'aire dans le

cadre de la théorie cinétique [3]. Il a été comparé dans les références [3] puis [4] la précision de deux relations de fermeture en écoulement simple, avec un modèle prenant en compte les différents phénomènes physiques que sont la coalescence, le fractionnement et la relaxation des gouttes, modèle initialement proposé par Lee et Park [5]. Nous proposons ici d'évaluer la précision de ces deux mêmes relations de fermeture pour un écoulement complexe en prenant en compte l'influence de la tension de surface à travers un terme de relaxation. Pour ce faire, nous comparerons les résultats obtenus à partir des relations de fermeture avec ceux issus de la formulation en théorie cinétique. Cette dernière est une approche micro-macro et, bien qu'ayant une formulation pouvant être considérée comme exacte, elle a comme inconvénient majeur d'être un problème à haute dimension. Sa résolution à l'aide de techniques numériques classiques est très coûteuse en temps de calcul et en mémoire. C'est pourquoi nous utiliserons la *Proper Generalized Decomposition* [6-8], qui est une technique de réduction qui permet de s'affranchir du problème de dimensionnalité.
## 2. Définition du tenseur d'aire

### 2.1. Equation avec relation de fermeture

Considérons tous d'abord le problème tel que défini par Doi et Ohta [1] en omettant les effets de la relaxation. Tous les points du domaine d'étude $\Omega$ sont définis par un vecteur position $\mathbf{x}$. Soit un volume représentatif microscopique $V(\mathbf{x})$, centré en un point $\mathbf{x} \in \Omega$, petit par rapport à l'échelle macroscopique mais assez grand par rapport à la micro-structure. Soit $S(\mathbf{x}, t)$ l'interface de $V(\mathbf{x})$ à un instant $t$ et $\mathbf{n}$ le vecteur unitaire normal à l'élément d'interface $dS$. Leur évolution respective au cours du temps est donnée par :

$$ \frac{dn}{dt} = -(\nabla \mathbf{v})^T \mathbf{n} + (\nabla \mathbf{v} : (\mathbf{n} \otimes \mathbf{n})) \mathbf{n} \quad (1) $$

$$ \frac{dS}{dt} = -(\nabla \mathbf{v} : (\mathbf{n} \otimes \mathbf{n})) dS \quad (2) $$

où $\frac{d}{dt}$ est la dérivée matérielle, définie par

$\frac{d\cdot}{dt} = \frac{\partial\cdot}{\partial t} + \mathbf{v} \cdot \nabla\cdot$, $\nabla$ est l'opérateur gradient et $\mathbf{v}$ la vitesse du fluide. Le tenseur d'aire est alors défini par [2] :

$$ A(\mathbf{x}, t) = \frac{1}{V(\mathbf{x})} \int_{S(\mathbf{x}, t)} \mathbf{n} \otimes \mathbf{n} \, dS . \quad (3) $$

Dans la suite, par souci de clarté, nous omettrons la dépendance de $A$ en $\mathbf{x}$ et $t$.

L'équation gouvernant l'évolution du tenseur d'aire peut être obtenue en prenant la dérivée temporelle de (3) et en utilisant les Eqs. (1) et (2). Il vient après quelques manipulations :

$$ \frac{dA}{dt} = -(\nabla \mathbf{v})^T \cdot \mathbf{A} - \mathbf{A} \cdot \nabla \mathbf{v} + \nabla \mathbf{v} : \mathbf{A}^4 . \quad (4) $$

où : est le produit doublement contracté et $\mathbf{A}^4$ le tenseur d'ordre quatre, tel que :

$$ \mathbf{A}^4 = \frac{1}{V(\mathbf{x})} \int_{S(\mathbf{x}, t)} \mathbf{n} \otimes \mathbf{n} \otimes \mathbf{n} \otimes \mathbf{n} \, dS . \quad (5) $$

Il n'est pas possible de calculer directement $\mathbf{A}^4$ et l'introduction d'une relation de fermeture permet d'exprimer $\mathbf{A}^4$ en fonction de $\mathbf{A}$. Nous évaluerons la fermeture quadratique, initialement proposée par Doi et Ohta [1] :

$$ \mathbf{A}^4 = \frac{1}{\text{tr}(\mathbf{A})} (\mathbf{A} \otimes \mathbf{A}) \quad (6) $$

ainsi que la fermeture orthotrope (appelée aussi ORE pour Orthotropic Rational Ellispoidal), proposée par Wetzel et Tucker [2, 9].

### 2.2. Propriétés du tenseur d'aire

Le tenseur d'aire précédemment défini est symétrique et comporte des propriétés intéressantes :

- la trace du tenseur d'aire $A$ correspond à l'aire d'interface par volume caractéristique. Cette grandeur, qui est appelée surface spécifique [2], est notée $S_v(\mathbf{x}, t)$ et est une information particulièrement intéressante pour quantifier le mélange.

- si l'on définit la fraction volumique de la phase disperse par $\phi$, que l'on assume uniformément distribuée dans $\Omega$, alors on peut définir la longueur caractéristique $L(\mathbf{x}, t)$ de la microstructure au point $\mathbf{x}$ du volume de la phase discrète : $V_d(\mathbf{x}) = \phi \cdot V(\mathbf{x})$.

- lors d'un mélange passif, le taux de déformation varie linéairement à l'échelle de $V(\mathbf{x})$ de telle manière qu'un domaine discret initialement sphérique sera déformé en un ellipsoïde. On peut alors aisément représenter le tenseur d'aire sous forme d'ellipsoïdes dont les longueurs d'axes sont données par les valeurs propres du tenseur d'aire normalisé et leur orientation par les vecteurs propres associées.

### 2.3. Modèle de théorie cinétique (sans relation de fermeture)

Pour définir un modèle de théorie cinétique, introduisons tout d'abord une fonction de distribution d'aire $\Psi(\mathbf{x}, t, \mathbf{n})$, qui associe à tout point de l'espace physique et quel que soit le temps $t$, l'orientation de l'interface dans la direction $\mathbf{n}$ [3].

Dans le cadre de la théorie cinétique, le tenseur d'aire s'écrit :

$$ A = \int_{B} (\mathbf{n} \otimes \mathbf{n}) \Psi(\mathbf{x}, t, \mathbf{n}) d\mathbf{n} \quad (7) $$

où $B$ représente la surface de la sphère unitaire. Par souci de clarté, on omettra par la suite la dépendance de $\Psi(\mathbf{x}, t, \mathbf{n})$ en $\mathbf{x}$, $t$ et $\mathbf{n}$. L'évolution temporelle de $\Psi$ est alors déduite de :

$$ \frac{dA}{dt} = \int_{B} \left( \frac{dn}{dt} \otimes \mathbf{n} \right) \Psi dn + \int_{B} \left( \mathbf{n} \otimes \frac{dn}{dt} \right) \Psi dn + \int_{B} (\mathbf{n} \otimes \mathbf{n}) \frac{d\Psi}{dt} dn \quad (8) $$

En utilisant l'équation (1), il vient :
$$ \int_B \left( \frac{dn}{dt} \otimes n \right) \Psi dn = -(\nabla v)^T \cdot A + \nabla v : A^4 \quad (9) $$

$$ \int_B \left( n \otimes \frac{dn}{dt} \right) \Psi dn = \int_B \left( n \otimes n \left( \nabla v : (n \otimes n) \right) \right) \Psi dn - A \cdot \nabla v \quad (10) $$

D'où :

$$ \frac{d\Psi}{dt} = -(\nabla v : (n \otimes n))\Psi \quad (11) $$

De plus :

$$ \frac{d\Psi}{dt} = \frac{\partial \Psi}{\partial t} + \frac{\partial}{\partial n} \left( \frac{dn}{dt} \Psi \right) \quad (12) $$

Ainsi, l'évolution temporelle de $\Psi(x,t,n)$ est donnée par :

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial n}\left(\frac{dn}{dt}\Psi\right) - (\nabla v : (n \otimes n))\Psi. \quad (13) $$

La trace du tenseur d'aire, $S_v$, s'écrit :

$$ \mathrm{Tr}(A) = S_v = \int_B \mathrm{Tr}(n \otimes n) \Psi dn = \int_B \Psi dn \quad (14) $$

Si la formulation dans le cadre de la théorie cinétique a pour avantage de ne pas introduire d'erreur à travers l'utilisation d'une relation de fermeture, elle a cependant l'inconvénient majeur d'être définie sur chaque point de l'espace physique. Le calcul du tenseur d'aire dans ce cadre est alors très délicat sans l'aide de techniques numériques spécifiques. Nous présenterons dans la section 3 une technique permettant de s'affranchir de cette difficulté.

## 2.4. Tension de surface

Dans la définition du tenseur d'aire présentée précédemment, nous n'avons pas pris en compte l'influence de la tension de surface sur l'évolution du mélange. Or, de nombreux phénomènes apparaissent lors de la prise en compte de ce terme, notamment des phénomènes de coalescence et de fractionnement des gouttes de la phase discrète. Ceux-ci influencent fortement l'évolution du mélange, ainsi que l'écoulement fluide. Ainsi, de nombreux auteurs [1, 10, 11] ont proposés des modèles de relaxation prenant en compte les phénomènes liés à la tension de surface. Nous choisissons d'introduire un phénomène de relaxation de forme, identique à celui exposé dans les références [1] et [10], dans les équation d'évolution du tenseur d'aire (Eqs. (4) et (13)). Ce terme correspond à la diminution de l'anisotropie du tenseur d'aire au

cours du temps, c'est-à-dire que les gouttes de la phase dispersée tendent vers un forme sphérique. Le terme de relaxation s'écrit :

$$ \left. \frac{dA}{dt} \right|_{\text{relaxation}} = -c S_v \left( A - \frac{S_v}{3} I \right) \quad (15) $$

$$ \text{où } c = a \frac{\Gamma}{\eta} \quad (16) $$

avec $a$ un coefficient positif gérant l'intensité de la relaxation, $\Gamma$ la tension de surface et $\eta$ la viscosité (les deux phases sont considérées equi-visqueuses). L'équation d'évolution du tenseur d'aire prenant en compte le phénomène de relaxation s'écrit alors :

$$ \frac{dA}{dt} = \left. \frac{dA}{dt} \right|_{\text{écoulement}} + \left. \frac{dA}{dt} \right|_{\text{relaxation}} \quad (17) $$

dont le terme $\left. \frac{dA}{dt} \right|_{\text{écoulement}}$ est défini Eq. (4).

Dans le cadre de la théorie cinétique, un terme de relaxation strictement équivalent est introduit à travers le terme de diffusion [12], tel que :

$$ \left. \frac{\partial \Psi}{\partial t} \right|_{\text{relaxation}} = \frac{\partial}{\partial n} \left( D \frac{\partial \Psi}{\partial n} \right) \quad (18) $$

$$ \text{où } D = \frac{c}{6} S_v \text{ et } c \text{ est défini Eq. (16).} $$

Finalement, l'équation d'évolution temporelle de $\Psi$ est donnée par :

$$ \frac{\partial \Psi}{\partial t} = \left. \frac{\partial \Psi}{\partial t} \right|_{\text{écoulement}} + \left. \frac{\partial \Psi}{\partial t} \right|_{\text{relaxation}} \quad (19) $$

avec $\left. \frac{\partial \Psi}{\partial t} \right|_{\text{écoulement}}$ défini par l'Eq. (13).

# 3. Proper Generalized Decomposition

## 3.1. Présentation générale

Le tenseur d'aire dans le cadre de la théorie cinétique est défini à la fois sur l'espace physique et sur l'espace des configurations, qui peuvent être respectivement de dimensions trois et deux. Si nous prenons en compte le caractère instationnaire, le problème peut être défini sur un espace de dimension six. Dans de tels cas, les Eléments Finis (EF) classiques (ou tout autre méthode numérique classique, de type Volumes Finis, Différences Finies...) ne peuvent être utilisés qu'à condition d'avoir une grande puissance de calcul et une importante capacité de stockage. En effet, ces
techniques sont sujettes à la malédiction de la dimensionnalité : la complexité du problème évolue exponentiellement en fonction du nombre de dimensions. Une technique de réduction de dimensions précédemment développée, la Proper Generalized Decomposition (PGD) [6-8], nous permet, si ce n'est de s'affranchir de ce problème, du moins d'en diminuer l'importance : dans le cadre de cette méthode, la complexité du problème évolue linéairement en fonction du nombre de dimensions.

L'idée principale de la PGD est de trouver une solution pouvant s'écrire sous forme séparée telle que :

$$ \Psi(\mathbf{x}, \mathbf{n}) = \sum_{i} \alpha_{i} F_{i}(\mathbf{x}) G_{i}(\mathbf{n}). \quad (20) $$

où F et G sont des fonctions de bases appartenant respectivement à l'espace physique et à l'espace des configurations. Il est intéressant de noter que le temps peut aussi être considéré comme une dimension supplémentaire dont la discrétisation est identique aux autres espaces. Cependant, la dimension temporelle sera traitée ici de manière classique.

Le fait d'écrire $\Psi$ sous forme séparée permet une intégration séparée sur chaque espace. D'un point de vue numérique, cela signifie qu'en discrétisant chaque espace avec respectivement $N_x$ et $N_n$ nœuds, la discrétisation totale de $\Psi$ sera de $N_x + N_n$ degrés de libertés, au lieu de $N_x \times N_n$ lors d'une discrétisation éléments finis classique. L'utilisation de cette méthode permet alors de réduire fortement le temps de calcul et la quantité de mémoire nécessaires. Le schéma numérique consiste en une procédure itérative qui résout à chaque itération les trois étapes suivantes :

- projection de la solution dans une base discrète en calculant le coefficient $\alpha$ de l'équation (20);

- vérification du critère de convergence;

- enrichissement de la base en ajoutant de nouvelles fonctions F et G.

Le processus d'enrichissement est stoppé lorsque le critère de convergence précédemment défini est atteint ou lorsque le nombre de fonctions de base est arbitrairement atteint.

## 3.2. Forme tensorielle

L'équation d'évolution de $\Psi$ (13) peut se réécrire sous la forme simplifiée :

$$ \frac{\partial \Psi}{\partial t} = \frac{\partial \Psi}{\partial t} + \mathbf{v} \cdot \nabla \Psi = -E_0 \Psi - \frac{\partial}{\partial \mathbf{n}} (E_1 \Psi) + \frac{\partial}{\partial \mathbf{n}} \left( E_2 \frac{\partial \Psi}{\partial \mathbf{n}} \right) \quad (21) $$

$$ \text{où : } E_0(\mathbf{x}, \mathbf{n}) = \nabla \mathbf{v} : (\mathbf{n} \otimes \mathbf{n}) \quad (22) $$

$$ E_1(\mathbf{x}, \mathbf{n}) = \frac{dn}{dt} \quad (23) $$

$$ E_2 = D \quad (24) $$

Les coefficients précédents sont définis à la fois sur l'espace physique et l'espace des configurations, mais leur dépendance sur l'espace physique réside seulement à travers le gradient de vitesse. Pour réécrire les coefficients précédents sous forme séparée, nous choisissons d'exprimer leur contribution dans la base canonique du gradient de vitesse. Le coefficient $E_0(\mathbf{x}, \mathbf{n})$ se réécrit alors :

$$ E_0(\mathbf{x}, \mathbf{n}) = \sum_{i=1}^{3} \sum_{j=1}^{3} E_0^{ij}(\mathbf{n}) \frac{dv_i}{dx_j}(\mathbf{x}) \quad (25) $$

La même approche est effectuée pour le coefficient $E_1(\mathbf{x}, \mathbf{n})$ afin d'obtenir une forme séparée des coefficients de l'équation d'évolution de $\Psi$. Le coefficient $E_2(\mathbf{x}, \mathbf{n})$ fait quant à lui intervenir $\Psi$ qui est déjà écrit sous forme séparée.

La forme faible associée à l'équation (21) s'écrit :

$$ \begin{aligned} & \iint_{\Omega B} \Psi^* \left( \frac{\partial \Psi}{\partial t} + \mathbf{v} \cdot \nabla \Psi \right) dx dn \\ &= \iint_{\Omega B} \Psi^* \left( -E_0 \Psi - \frac{\partial}{\partial \mathbf{n}} (E_1 \Psi) + \frac{\partial}{\partial \mathbf{n}} \left( E_2 \frac{\partial \Psi}{\partial \mathbf{n}} \right) \right) dx dn \end{aligned} \quad (26) $$

où $\Psi^*$ est la fonction test.

L'équation précédente est linéarisée en exprimant $E_1$ et $E_2$ en fonction de $\Psi$ obtenu à l'itération précédente. Après discrétisation du terme temporel selon le schéma d'Euler implicite, il vient :

$$ \begin{aligned} & \iint_{\Omega B} \Psi^* (1 + \Delta t (\mathbf{v} \cdot \nabla) + \Delta t \cdot E_0) \Psi^{t+1} dx dn \\ & + \iint_{\Omega B} \Psi^* \Delta t \left( E_1 \frac{\partial}{\partial \mathbf{n}} - E_2 \frac{\partial^2}{\partial \mathbf{n}^2} \right) \Psi^{t+1} dx dn = \iint_{\Omega B} \Psi^* \Psi^t dx dn \end{aligned} \quad (27) $$

L'équation précédente est alors discrétisée en utilisant une approximation éléments finis classique sur chaque espace. Pour ce faire, on introduit les vecteurs $\mathbf{N}$ et $\mathbf{M}$ qui contiennent les fonctions de
forme associées respectivement à l'espace de coordonnées **x** et **n**.

Ainsi, la forme discrète nécessite l'introduction des matrices suivantes :

$$A_1^1 = \int_{\Omega} \mathbf{N} \mathbf{N}^{\mathrm{T}} d\mathbf{x}$$

$$A_2^1 = \int_B \mathbf{M} \mathbf{M}^{\mathrm{T}} d\mathbf{n}$$

$$A_1^2 = \Delta t \int_{\Omega} \mathbf{N} \mathbf{v} \frac{d\mathbf{N}^{\mathrm{T}}}{d\mathbf{x}} d\mathbf{x}$$

$$A_2^2 = \int_B \mathbf{M} \mathbf{M}^{\mathrm{T}} d\mathbf{n}$$

$$A_1^3 = \Delta t \int_{\Omega} \mathbf{N} \Psi_x' \mathbf{N}^{\mathrm{T}} d\mathbf{x}$$

$$A_2^3 = \int_B \frac{d\mathbf{N}}{d\mathbf{n}} \left( \frac{c}{6} \int_B \Psi_n' d\mathbf{n} \right) \frac{d\mathbf{N}^{\mathrm{T}}}{d\mathbf{n}} d\mathbf{n}$$

$$A_1^{3(i-1)+j+3} = \Delta t \int_{\Omega} \mathbf{N} \frac{d\mathbf{v}_i}{d\mathbf{x}_j} \mathbf{N}^{\mathrm{T}} d\mathbf{x}$$

$$A_2^{3(i-1)+j+3} = \int_B \left( \mathbf{M} \mathbf{E}_0^{ij} \mathbf{M}^{\mathrm{T}} + \mathbf{M} \mathbf{E}_1^{ij} \frac{d\mathbf{M}^{\mathrm{T}}}{d\mathbf{n}} \right) d\mathbf{n}$$

$$B_1^1 = \int_{\Omega} \mathbf{N} \mathbf{N}^{\mathrm{T}} d\mathbf{x} \cdot \Psi_x'$$

$$B_2^1 = \int_{\Omega} \mathbf{M} \mathbf{M}^{\mathrm{T}} d\mathbf{n} \cdot \Psi_n'$$

avec $i, j = 1, 2, 3$ et $\Psi_x'$, $\Psi_n'$ les décompositions de la fonction $\Psi$, respectivement sur l'espace physique et l'espace des configurations.

Ainsi, la forme discrète s'écrit :

$$A\Psi = B \qquad (28)$$

avec :

$$A = \sum_{i=1}^{11} A_1^i \otimes A_2^i \qquad (29)$$

$$B = B_1^1 \otimes B_2^1 \qquad (30)$$

La forme tensorielle définie dans l'équation (28) est alors résolue grâce à l'algorithme défini dans [6-8] et nous invitons le lecteur intéressé à s'y référer pour plus d'informations.

# 4. Résultats numériques

## 4.1. Validation en écoulement simple

Nous nous intéressons tout d'abord aux résultats issus de la PGD pour des écoulements simples. Dans ce cas l'opérateur défini dans l'équation (29) devient :

$$A = 1 \otimes \int_B \mathbf{M} \mathbf{M}^{\mathrm{T}} d\mathbf{n} + 1 \otimes \Delta t \int_B \frac{d\mathbf{N}}{d\mathbf{n}} \left( \frac{c}{6} \int_B \Psi_n' d\mathbf{n} \right) \frac{d\mathbf{N}^{\mathrm{T}}}{d\mathbf{n}} d\mathbf{n} \\ + \Delta t \sum_{i=1}^{9} \sum_{j=1}^{9} \frac{d\mathbf{v}_i}{d\mathbf{x}_j} \otimes \int_B \left( \mathbf{M} \mathbf{E}_0^{ij} \mathbf{M}^{\mathrm{T}} + \mathbf{M} \mathbf{E}_1^{ij} \frac{d\mathbf{M}^{\mathrm{T}}}{d\mathbf{n}} \right) d\mathbf{n} \qquad (31)$$

et le second membre de l'équation (28) s'écrit :

$$B = 1 \otimes \int_{\Omega} \mathbf{M} \mathbf{M}^{\mathrm{T}} d\mathbf{n} \cdot \Psi_{n}' \qquad (32)$$

Dans ce cas, l'espace physique est réduit à un seul point et il est possible de traiter le problème indifféremment avec la méthode EF classique ou la méthode PGD. Nous comparons pour chaque méthode l'évolution au cours du temps de la surface spécifique $S_V$ (trace du tenseur **A**) pour un cisaillement simple et une elongation bi-axiale. Initialement, la phase discrète du mélange est sous forme de gouttes sphériques, ce qui signifie que le tenseur d'aire s'écrit :

$$A = \begin{pmatrix} \frac{1}{3} & 0 & 0 \\ 0 & \frac{1}{3} & 0 \\ 0 & 0 & \frac{1}{3} \end{pmatrix} \qquad (33)$$

Le gradient de vitesse est défini pour l'écoulement de cisaillement par :

$$\nabla \mathbf{v} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \qquad (34)$$

et pour l'élongation bi-axiale par :

$$\nabla \mathbf{v} = \begin{pmatrix} -1 & 0 & 0 \\ 0 & \frac{1}{2} & 0 \\ 0 & 0 & \frac{1}{2} \end{pmatrix} \qquad (35)$$

L'espace des configurations, représenté par la sphère unité, est discrétisé à l'aide plus de 5000 éléments
triangles à trois nœuds pour un total de près de 2500 nœuds. Pour les deux méthodes, le schéma temporel est identique (schéma d'Euler implicite).

On trace sur les Figures 1 et 2 l'évolution de la surface spécifique au cours du temps pour trois paramètres de relaxation différents. Les résultats sont identiques entre les deux méthodes de calcul. La méthode PGD est une méthode itérative dont la précision dépend du nombre de fonctions de bases, définies par *F* et *G* dans la section précédente. Pour les écoulements simples, la convergence de la méthode est très forte et est atteinte avec l'ajout d'un seul produit de fonctions à chaque itération.

<figure><img src="image_2.png" /><figcaption>Figure 1. Comparaison de l'évolution de la surface spécifique $S_V$ en fonction du temps et de trois valeurs du paramètre $c$ (en m.s⁻¹) pour un écoulement de cisaillement simple.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Comparaison de l'évolution de la surface spécifique $S_V$ en fonction du temps et de trois valeurs du paramètre $c$ (en m.s⁻¹) pour un écoulement d'élongation bi-axiale.</figcaption></figure>

D'autre part, le terme de relaxation introduit dans la section 2.3 influence l'évolution de $S_V$ comme on pouvait s'y attendre : plus le coefficient $c$ est important

et plus l'évolution de la surface spécifique est affectée et sa croissance faible.

Les figures 3 et 4 représentent $\Psi$ dans l'espace des configurations, c'est-à-dire la sphère unité, respectivement pour l'écoulement de cisaillement simple et d'élongation bi-axiale, pour $c = 0 \text{ m.s}^{-1}$. L'orientation de l'aire d'interface est représentée par la variation de $\Psi$: plus $\Psi$ est grand et plus l'orientation dans la direction considérée est prononcée.

Note : $\Psi$ est identique pour les deux méthodes, à l'erreur machine près.

<figure><img src="image_4.png" /><figcaption>Figure 3. Représentation de l'espace des configurations pour un cisaillement simple à $t = 0,1s$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Représentation de l'espace des configurations pour une elongation bi-axiale à $t = 0,1s$.</figcaption></figure>

## 4.2. En écoulement complexe

Dans la section précédente, nous avons montré que les résultats issus de la PGD pour des écoulements simples étaient en tout point comparables avec ceux issus de la méthode EF. Or, dans ce cadre, nous n'introduisons pas de relation de fermeture et pouvons ainsi considérer les résultats comme exacts.
C'est pourquoi nous prenons les résultats PGD en écoulement complexe comme référence par rapport aux relations de fermetures quadratique et orthotrope.

L'écoulement est une contraction 2:1 telle que présentée à la Figure 5, de largeur maximum 2 m et de hauteur maximum 0,6 m, et dont le bord inférieur représente un axe de symétrie. On impose en entrée un écoulement de Poiseuille, défini par :

$$v = \begin{pmatrix} 0.5 \left( 1 - \frac{y^2}{y_{\max}^2} \right) \\ 0 \end{pmatrix} \qquad (36)$$

avec $v$ la vitesse en m.s⁻¹.

<figure><img src="image_6.png" /><figcaption>Figure 5. Champ de vitesse fluide pour l'écoulement complexe. Noter que seule une moitié de la contraction est représentée : le bord inférieur est un axe de symétrie.</figcaption></figure>

Les fluides sont incompressibles, newtoniens, de même densité et immiscibles. Les contraintes induites par le mélange ne sont pas prises en compte, ce qui signifie que le champ de vitesse reste identique tout au long de la simulation et son état stationnaire est obtenu par résolution de l'équation de Navier-Stokes. De plus, on résout l'équation (4) avec les relations de fermeture ainsi que le système tel que définis dans l'Eq. (28) pour la résolution PGD. Celles-ci prennent en compte l'effet de la tension de surface à travers le terme de relaxation de forme où $c = 3 \text{ m.s}^{-1}$.

Deux conditions aux limites sont imposées : la microstructure est imposée sphérique en entrée et sur le bord supérieur, c'est-à-dire que le tenseur d'aire s'écrit sous la forme donnée équation (33). Ces conditions correspondent à un fluide qui n'a subi aucun mélange avant d'entrer dans la contraction et qui est parfaitement sphérique au contact de la paroi. On simule l'évolution du mélange pour les trois méthodes jusqu'à atteindre l'état stationnaire.

Nous choisissons de représenter dans les Figures 6 à 8 l'orientation de la microstructure sous forme d'ellipsoïdes et de quantifier le mélange en fonction de la surface spécifique et plus exactement du rapport $1/S_v$.

<figure><img src="image_7.png" /><figcaption>Figure 6. Résultats obtenus avec la relation de fermeture quadratique. L'orientation de la microstructure est donnée par les ellipsoïdes et les isocouleurs représentent $1/S_v$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 3. Résultats obtenus avec la relation de fermeture orthotrope. L'orientation de la microstructure est donnée par les ellipsoïdes et les isocouleurs représentent $1/S_v$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4. Résultats obtenus sans relation de fermeture (théorie cinétique). L'orientation de la microstructure est donnée par les ellipsoïdes et les isocouleurs représentent $1/S_v$.</figcaption></figure>

Afin de pouvoir comparer les résultats obtenus entre les différentes simulations, nous résumons Tableau 1 la surface spécifique totale dans la géométrie, ainsi que les maxima atteints.

<table><caption>Tableau 1. Comparaison des résultats issus en écoulement complexe en fonction de la méthode de résolution.</caption><thead><tr><th></th><th>Quadratique</th><th>Orthotrope</th><th>Théorie Cinétique</th></tr></thead><tbody><tr><td>$S_v$ totale</td><td>1,73</td><td>1,55</td><td>1,58</td></tr><tr><td>$S_v$ min</td><td>1</td><td>1</td><td>1</td></tr><tr><td>$S_v$ max</td><td>20,06</td><td>15,75</td><td>18,30</td></tr></tbody></table>

Les Figures 6 à 8 présentent les résultats obtenus en état stationnaire pour les trois simulations différentes. Ces résultats sont relativement proches. En effet, la surface spécifique totale, qui correspond à l'intégration de la surface spécifique sur toute la géométrie (voir Tableau 1), est du même ordre. Cependant, l'erreur relative de la relation orthotrope est bien inférieure à celle de la relation quadratique, puisqu'elle est de l'ordre de 2% contre 9% respectivement. Cela est tout à fait en accord avec les résultats de la littérature [2-4], où est mise en
avant la plus grande précision de la fermeture
orthotrope, proposée par Wetzel et Tucker [2].

D'un point de vue numérique, les simulations avec
relation de fermeture correspondent à un problème
défini par près de 1000 degrés de liberté, alors que la
définition dans le cadre de la théorie cinétique
entraîne 2500 fois plus d'inconnues : à chaque point
de l'espace physique est associé un espace de
configuration. Or, l'utilisation de la PGD permet de
s'affranchir des difficultés qu'aurait engendrée la
résolution EF d'un tel problème. En effet, la
résolution d'un problème de plus de 2 500 000
degrés de liberté dans le cadre de méthode
numérique classique nécessite le stockage et
l'inversion d'une matrice composées de 2 500 000 x
2 500 000 termes, ce qui pose d'importants problè-
mes de mémoire et de temps de calcul (et cela même
lors de l'utilisation de matrice creuse et de résolution
itératives). Il a ainsi été possible de résoudre le
problème couplé micro-macro sur un PC équipé
d'un processeur Intel Core i7 en de deux heures dans
le cadre de la méthode PGD.

Les simulations avec relations de fermeture restent néanmoins bien plus efficaces en termes de temps de calcul (mais au détriment d'une perte de précision) puisque celui-ci ne dépasse pas quelques minutes. Dans les cas où une grande précision n'est pas requise, l'erreur introduite par les relations de fermetures peut être acceptable. Ainsi, la fermeture quadratique, qui est cinq fois plus rapide que la fermeture orthotrope, a une utilisation tout à fait justifiée dans le cadre d'écoulements complexes [4].

**5. Conclusion**

Le tenseur d'aire est un outil efficace pour évaluer
quantitativement et qualitativement le mélange de
fluides immiscibles. Cependant, les équations
gouvernant l'évolution du tenseur d'aire nécessitent
l'introduction d'une relation de fermeture dont
l'erreur ne peut être évaluée a priori. L'évaluation de
cette erreur dans un écoulement complexe a pu se
faire grâce à l'utilisation de la PGD, technique de
réduction de dimension, permettant de résoudre une
approche micro-macro équivalente ne nécessitant
pas de relation de fermeture. Cependant, bien que la
PGD permette de résoudre sans approximation
l'évolution du tenseur d'aire en écoulement
complexe, son coût en terme de calcul reste bien

supérieur à celui d'une résolution avec relation de
fermeture. La comparaison des résultats permet de
justifier l'utilisation de relations de fermeture en
fonction du degré de précision désiré : la fermeture
quadratique, plus simple à mettre en œuvre par
rapport à la fermeture orthotrope, plus précise.

**6. Références**

[1] M. Doi, T. Ohta. Dynamics and rheology of complex interfaces. I. *J. Chem. Phys.*, 95, 15 (1991).

[2] E. D. Wetzel, C. L. Tucker III. Area tensors for modelling microstructure during laminar liquid-liquid mixing. *Int. J. Multiph. Flow*, 25, 35-61 (1999).

[3] F. Chinesta, M. R. Mackley. Microstructure evolution during liquid-liquid laminar mixing: A kinetic theory approach. *Int. J. Mater. Form.*, 1, 47-55 (2008).

[4] G. Maîtrejean, A. Ammar, F. Chinesta, Simulating microstructural evolution during passive mixing, *Int. J. Mater. Form.*, In press (2011).

[5] H. M. Lee, O. Park. Rheology and dynamics of immiscible polymer blends. *Science*, 38, (1994).

[6] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A new family of solvers for some classes of multi-dimensional partial differential equations encountered in kinetic theory modeling of complex fluids *J. Non-Newt. Fluid Mech.*, 139, 153-176 (2006).

[7] A. Ammar, E. Pruliere, F. Chinesta, M. Laso, Reduced numerical modeling of flows involving liquid-crystalline polymers. *J. Non-Newt. Fluid Mech.*, 160, 140-156 (2009).

[8] A. Ammar, D. Ryckelynck, F. Chinesta, R. Keunings, On the reduction of kinetic theory models related to finitely extensible dumbbells. *J. Non-Newt. Fluid Mech.*, 134, 136-147 (2006).

[9] J.S. Cintra Jr, C.L. Tucker III. Orthotropic closure approximations for flow-induced fiber orientation. *J. Rheol.*, 39, 1085-1122 (1995).

[10] H. Lee, O. Park. Rheology and dynamics of immiscible polymer blends. *J. Rheol.*, 38, 1405-1425 (1994).

[11] M. Grmela, A. Ait-Kadi, L.A. Utracki, Blend of two immiscible and rheologically different fluids *J. Non-Newt. Fluid Mech.*, 77, 253-259 (1998).

[12] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, *Dynamics of Polymeric Liquids*, John Wiley, New York (1987).