Etude numérique de mélanges visco-élastiques

Laurent CHUPIN

Mathématiques Appliquées de Bordeaux, UMR 5466, Université Bordeaux 1, 351 cours de la libération, 33405 Talence cedex

Franck BOYER

Laboratoire d’Analyse, Topologie et Probabilités, UMR 6632 CNRS/Université de Provence, 39 rue Joliot-Curie, 13453 Marseille cedex 13

Mots-clefs : Modèles de Cahn-Hilliard, Navier-Stokes, Oldroyd; simulations numériques

Introduction

Suite à de nombreux travaux dans ce domaine [1, 6, 9, 18] nous proposons dans [2] un modèle physique pour l’étude numérique des écoulements de mélanges de deux phases non-newtoniennes. Les fluides dits visco-élastiques possèdent une mémoire des déformations antérieures, leur comportement ne peut donc pas être prédit par les lois usuelles des fluides newtoniens. Par conséquent, une loi constitutive de type Oldroyd sera ajoutée aux classiques équations de Navier-Stokes pour modéliser l’écoulement. Dans le modèle considéré, afin de prendre en compte le caractère diphasique de l’écoulement, un champ scalaire décrivant la concentration d’un fluide dans le mélange est introduit. L’interface entre les deux fluides est supposée diffuse de sorte que les deux phases peuvent se mélanger. Il en résulte un système couplant les équations de Navier-Stokes, la loi constituitive de type Oldroyd et une équation de Cahn-Hilliard. Dans un tel modèle, les deux phases sont incompressibles mais peuvent avoir des densités différentes.

Dans la première partie, nous donnons une dérivation du système couplé Cahn-Hilliard/Navier-Stokes/Oldroyd. Une mise sous forme adimensionelle permet d’introduire quatre grandeurs indépendantes, caractéristiques du modèle.

Une approche par “splitting” pour la résolution numérique est présentée dans la partie 2. Ce schéma est une extension de celui proposé dans [1] dans le cas de fluides newtoniens. Nous nous intéressons donc ici plus précisemment à la discrétisation de la loi constitutive et nous proposons un schéma de type exponentiel pour la discétisation en temps. En conclusion de cette partie, nous décrivons la discrétisation spatiale utilisée, et en particulier concernant les termes de convection pour lesquels nous avons implémenté la méthode de Runge-Kutta.

Une dernière partie est consacrée aux simulations numériques. Les premiers résultats présentés concernent les cisaillements à la fois monophasiques et diphasiques. Dans le premier cas, la différence de comportement entre les fluides newtoniens et visco-élastiques est mise en évidence. Dans le second, nous observons la décomposition spinodale et la structuration des phases sous cisaillement dans les cas newtoniens et non-newtoniens.

Une des raisons de l’intérêt pour les écoulements visco-élastiques est lié aux nombreuses avancées technologiques concernant les mélanges de gouttelettes, les sprays, les élongations de fibres... Nous présentons deux applications numériques allant dans ce sens. Tout d’abord un test d’allongement de fibres réalisé en collaboration avec A. Colin et O. Greffier au Centre de Recherche Paul Pascal (Bordeaux) puis des expériences de remplissage de cuves mettant en évidence la formation plus importante de bulles dans le cas visco-élastique que dans le cas newtonien.

1 Dérivation du modèle

1.1 Loi constitutive pour le tenseur des contraintes : cas monophasique

Les lois de la mécanique pour un fluide newtonien stipulent que la contrainte appliquée à un élément de fluide est proportionnelle à la déformation que subit cet élément. Pour des fluides plus complexes tels les fluides visco-élastiques, nous devons définir une loi plus générale liant le tenseur des déformations et le tenseur des contraintes. Cette relation peut être vue soit sous forme intégrale, mettant en avant les effets “mémoires” d’un fluide, soit sous forme différentielle. Puisque nous voulons obtenir en dernier lieu une équation différentielle, c’est cette dernière option que nous allons choisir. Il faut signaler que de très nombreuses lois existent et chacune d’elle est une loi empirique convenant à diverses situations physiques (voir par exemple [15]).

La loi constitutive que nous allons utiliser par la suite est due à J.G. Oldroyd [17]. Elle peut être vue comme une généralisation d’un modèle mécanique simple monodimensionnel (modèle de Jeffrey, [14]). Si on note τ le tenseur des contraintes et D(v) celui des déformations, le modèle d’Oldroyd reliant ces deux quantités s’exprime ainsi

τ = 2η(1 − r)D(v) + σ, \frac{D σ}{Dt} + \frac{σ}{λ} = \frac{2ηr}{λ} D(v).  (1)

λ étant le temps de relaxation, r le paramètre de retard, η étant la viscosité du fluide. Bien entendu, la relation obtenu dans le cadre monodimensionnel est formulé dans les variables lagrangiennes, c’est pour cette raison que nous avons introduit l’opérateur \frac{D}{Dt} qui correspond à une dérivée convective en temps. Plusieurs modèles sont possibles pour définir cette dérivée, la seule contrainte étant l’invariance par toute transformation galiléenne (voir [14, 17]). Dans la littérature, les modèles classiques ont la forme suivante :

\frac{D M}{Dt} = \frac{∂ M}{∂t} + v.∇M − W (v).M + M.W (v) − a(D(v).M + M.D(v)),  (2)

le coefficient a ∈ [−1, 1] est un paramètre rhéologique, D(v) et W (v) désignant respectivement les parties symétrique et anti-symétrique du gradient de vitesse.

1.2 Conservation de la masse

Nous allons maintenant prendre en compte la présence des deux phases dans l’écoulement. Nous avons choisi de considérer les modèles à interface diffuse (en opposition aux modèles à interface singulière). Ce type de modèle a une origine thermodynamique et permet de décrire le système final avec seulement des quantités continues (variant très rapidement à l’interface). Les deux fluides incompressibles considérés seront par la suite distingués par les indices 1 et 2. Pour décrire le mélange de façon continue, nous introduisons la fraction volumique d’un des deux fluides dans le mélange (voir [1]) :

Φ = Φ(t, x) = \frac{dV_{1}}{dV} ,

où dV_{1} est le volume occupé par la phase 1 dans le volume élémentaire dV. Pour des raisons de compatibilité avec les modèles fréquemment étudiés [1, 5, 10], le paramètre d’ordre utilisé est une renormalisation de Φ donnée par φ = 2Φ − 1.

Les densités de chacune des phases au sein du mélange sont notées ρ_{1} et ρ_{2} alors que les densités des deux fluides seront définies par ρ_{1}^{0} et ρ_{2}^{0}. Les vitesses propres à chaque phase étant notées v_{1} et v_{2}, on peut alors définir des quantités moyennées :

ρ = ρ_{1} + ρ_{2}, ρ_{1} = \frac{1 + φ}{2}ρ_{1}^{0}, ρ_{2} = \frac{1 − φ}{2}ρ_{2}^{0},

v = \frac{1 + φ}{2}v_{1} + \frac{1 − φ}{2}v_{2}.

Ainsi, en accord avec les calculs effectués dans [1, 6] les équations de conservation de la masse pour chacune des phases deviennent

div (v) = 0,

\frac{∂ φ}{∂t} + v.∇φ + div (\frac{1 − φ^{2}}{2} w) = 0,

où w = v_{1} − v_{2}. Cette dernière équation d’évolution sur le paramètre d’ordre contient deux termes de natures très différentes. Le premier est un terme de transport de l’interface à la vitesse moyenne du mélange et le second est un terme de diffusion à travers l’interface, le flux de cette diffusion étant proportionnel à la vitesse relative des deux phases.

1.3 Conservation de la quantité de mouvement

Afin de modéliser les propriétés hydrodynamiques du mélange, nou avons besoin d’une équation d’évolution sur le champ de vitesse v. Dans ce but, nous écrivons les équations de conservation de la quantité de mouvement pour chaque phase :

\frac{∂ (ρ_{i}v_{i})}{∂t} − div (Σ_{i}) + div (ρ_{i}v_{i} ⊗ v_{i}) = F_{i}, i = 1, 2

où F_{i} désigne les forces volumiques extérieures appliquées à chaque fluides et Σ_{i} le tenseur des contraintes pour chaque phase.

Forces intérieures : Le tenseur des contraintes totales s’écrit Σ_{i} = −p_{i}Id + τ_{i} où la pression correspond à la partie diagonale. Nous supposons qu’il existe des quantités moyennes p et τ définie par (voir [18]) :

p_{1} = \frac{1 + φ}{2}p, p_{2} = \frac{1 − φ}{2}p, τ_{1} = \frac{1 + φ}{2}τ , τ_{2} = \frac{1 − φ}{2}τ .

Là aussi, comme dans le cas monophasique (voir partie 1.1), pour dériver une loi constitutive pour le mélange de deux phases nous supposons que le tenseur des contraintes τ se décompose en une partie purement newtonienne et une partie purement élastique. Plus précisement, nous faisons l’hypothèse que τ a un comportement rhéologique identique au cas monophasique (voir équation (1)), les coefficients η, λ et r apparaisant sont choisis en fonction des coefficients η_{i}, λ_{i} et r_{i} de chacune des phases de la façon suivante (voir [7, 18]) :

η_{1}(φ) = \frac{1 + φ}{2}η(φ), η_{2}(φ) = \frac{1 − φ}{2}η(φ).

λ = \frac{1 + φ}{2}λ_{1} + \frac{1 − φ}{2}λ_{2}, r = \frac{1 + φ}{2}r_{1} + \frac{1 − φ}{2}r_{2}.

Il est important de noter que le modèle obtenu permet de traiter aussi bien les cas monophasiques (en choisissant par exemple φ = 1) que les cas diphasiques.

Forces extérieures : Elles sont de trois types différents. Les forces de gravité, les forces de frottement entre les deux phases et les forces dues aux interactions chimiques :

F_{1} = ρ_{1}g − ξ(φ)w − ρ_{1}∇μ_{1} , F_{2} = ρ_{2}g + ξ(φ)w − ρ_{2}∇μ_{2}.

A. Onuki [18] suggère d’écrire les potentiels chimiques sous la forme suivante

μ_{1} = μ_{0}(φ) + \frac{1 − φ}{2ρ(φ)} \frac{∂ F}{∂φ} , μ_{2} = μ_{0}(φ) − \frac{1 + φ}{2ρ(φ)} \frac{∂ F}{∂φ} ,

où μ_{0}(φ) représente le potentiel d’un mélange uniforme de composition φ. Les autres termes décrivant la non-uniformité du mélange, prenant en compte l’energie libre de Gibbs données par une expression du type

\frac{∂ F}{∂φ} = −AΔφ + EF′(φ).

Les constantes physiques A et E décrivent les propriétés de l’interface (voir par exemple [1, 8]). La fonction F est le potentiel de Cahn-Hilliard généralement donné par une expression de type polynomiale [4] :

F (φ) = C_{1}φ^{4} − C_{2}φ^{2}, C_{1} > 0, C_{2} > 0.

1.4 Modèle final

Suivant la philosophie des modèles à interface diffuse, les deux équations de conservation de la quantité de mouvement obtenues permettent d’obtenir deux équations d’évolution, l’une sur la vitesse moyenne v, l’autre sur la vitesse relative w. Nous ne donnons pas les détails de ces calculs ni de l’analyse dimensionelle du système qui permet de négliger certains termes devant d’autres dans des considérations physiques raisonnables (ces calculs sont entièrement décrits dans [1, 6]).

Pour l’adimensionnement du système, il est usuel d’introduire les grandeurs caractéristiques suivantes (toutes les quantités surmontées d’une barre correspondant à des grandeurs de référence : L̄ pour la longueur du domaine, T̄ pour le temps...)

Re = ρ̄v̄L̄/η̄, le nombre de Reynolds,
K = E/ρ̄v̄^{2}, le nombre capillaire,
Pe = ρ_{2}^{0}/βρ_{1}^{0}, le nombre de Peclet,
We = We(φ) = λ(φ)/T̄, le nombre de Weissenberg.

Finalement, le système adimensionné couplant à la fois le paramètre d’ordre φ, la vitesse v et le tenseur des contraintes σ a la forme suivante :

\{ \begin{matrix}ρ(φ) = 1 − (\frac{ρ_{2}^{0} − ρ_{1}^{0}}{ρ̄} )\frac{1 − φ}{2}, \\ \frac{∂ φ}{∂t} + v.∇φ − \frac{1}{Pe} div (\frac{B(φ)}{ρ}∇ (\frac{μ}{ρ}) ) = 0,  μ = −α^{2}Δφ + F′(φ), \\ ρ(\frac{∂ v}{∂t} + v.∇v) − \frac{1}{Re} div (2η(φ)(1 − r(φ))D(v)) + ∇p − div σ \\ = ρg + K(\frac{ρ_{2}^{0} − ρ_{1}^{0}}{ρ̄} )\frac{1 − φ^{2}}{2} ∇ (\frac{μ}{ρ}) + Kμ.∇φ, \\ div (v) = 0, \\ \frac{D σ}{Dt} + \frac{σ}{We(φ)} = \frac{2η(φ)r(φ)}{Re We(φ)}D(v), \end{matrix}  (3)

où la dérivée \frac{D}{Dt} est donnée par (2), B(φ) = (1−φ^{2})^{2}/8ξ(φ) (coefficient de mobilité) et α mesure l’épaiseur de l’interface.

La structure de ce système est maintenant claire : il est composé d’une équation de Cahn-Hilliard comprenant un terme de convection, d’une équation de Navier-Stokes avec deux termes sources (capillaire et visco-élastique) et d’une équation d’Oldroyd donnant l’évolution de la contrainte visco-élastique.

1.5 Conditions aux limites

Pour compléter cette description physique, nous spécifions les conditions au bord du domaine. Concernant la vitesse, nous choisirons des conditions de type Dirichlet ou Neumann selon les tests réalisés (couplant parfois ces deux types de conditions, voir partie 3) :

v = v_{b}  ou  \frac{∂ v}{∂n} = 0.

Si on impose la non-diffusion à travers les parois ainsi que l’angle de contact de l’interface avec cette paroi (ici en supposant que cette angle est un angle droit), nous avons

\frac{∂ μ}{∂n} = 0  et  \frac{∂ φ}{∂n} = 0.

Concernant les conditions aux limites pour la contrainte σ, il suffit d’imposer la contrainte aux endroits où les lignes de courants entrent dans le domaine [2].

2 Schéma numérique

Au cours de cette partie, nous proposons une approximation numérique du système précédent dans le cas bi-dimensionnel. Dans toutes les applications, le domaine d’étude sera rectangulaire et la seule différence sera la donnée de la condition initiale et des conditions aux limites.

2.1 Discrétisation en temps

Nous utilisons une méthode de pas fractionnaires. Les trois étapes correspondent à la résolution des trois équations principales du système (3).

Equation de Cahn-Hilliard : Connaissant φ^{n}, μ^{n} et v^{n} (les valeurs de φ, μ et v à l’instant t^{n}), la première étape consiste à résoudre l’équation de Cahn-Hilliard sans le terme de transport. Le principe est d’effectuer une θ-méthode avec θ > 0.5 pour assurer à la fois une bonne précision et la stabilité asymptotique. Une méthode de point fixe est utilisée pour résoudre le système non linéaire obtenu [1]. Pour les termes de convection, nous utilisons un schéma de type Runge-Kutta d’ordre 3 : par exemple sur une équation du type ∂_{t}φ + v.∇φ = 0, nous proposons

φ^{n+1} = φ^{n} − dtK(φ^{n}) + \frac{1}{2}dt^{2}K^{2}(φ^{n}) − \frac{1}{6}dt^{3}K^{3}(φ^{n}),

où K(f) = u.∇f. Une discrétisation spatiale de cet opérateur est proposée dans le paragraphe 2.2.2.

Equations de Navier-Stokes : Puisque nous considèrons des tests à vitesse relativement petite, le problème n’est pas raide et les termes non linéaires ne génèrent pas de problème. La méthode utilisée est implicite pour la pression et les termes visqueux, explicite pour les termes d’inertie et les termes couplés :

ρ^{n+1}(\frac{v^{n+1} − v^{n}}{dt} ) − \frac{1}{Re} div (2η^{n+1}(1 − r^{n+1})D(v^{n+1})) + ∇p^{n+1} = div σ^{n}

− ρ^{n+1}(v.∇v)^{n} + ρ^{n+1}g + K\frac{ρ_{2}^{0} − ρ_{1}^{0}}{ρ̄} \frac{1 − (φ^{n+1})^{2}}{2} ∇ (\frac{μ^{n+1}}{ρ^{n+1}} ) + Kμ^{n+1}.∇φ^{n+1},

div (v^{n+1}) = 0,

où nous avons introduit ρ^{n+1}, η^{n+1} et r^{n+1} pour noter ρ(φ^{n+1}), η(φ^{n+1}) et r(φ^{n+1}). Finalement, une méthode de lagrangien augmenté permet de résoudre ce problème de Stokes.

Equation d’Oldroyd : La dernière étape est celle permettant de déterminer σ^{n+1} en connaissant v, φ et σ à l’instant t^{n}. Pour cela, il est utile de récrire l’équation d’évolution en la contrainte sous la forme suivante (identifier les tenseurs d’ordre 2 et les vecteurs de ℝ^{4}) :

\frac{∂ σ}{∂t} + (\frac{1}{We(φ)}Id − M_{a}).σ = \frac{2η(φ)r(φ)}{We(φ)}D(v).

La matrice M_{a} étant clairement déterminée par la donnée de la vitesse v. Nous trouvons de manière exacte la solution de cette équation linéaire que l’on approche par (voir [2, 11]) :

σ^{n+1} = e^{−dt(\frac{1}{We(φ^{n+1})} Id−M_{a}^{n+1})}.σ^{n} + \frac{2η^{n+1}r^{n+1}}{We(φ^{n+1})} dt e^{−dt(\frac{1}{We(φ^{n+1})} Id−M_{a}^{n+1})}.D(v^{n+1}).

La dernière étape consiste à déterminer l’exponentielle de la matrice M_{a}. Cette matrice étant de taille 4 ce calcul est réalisé directement en déterminant les valeurs propres et vecteurs propres de M_{a}. Nous ne donnons pas ici tous les détails (voir [2]) mais une remarque essentielle est que la propriété de symétrie du tenseur est préservée. Plus précisemment, si la donnée initiale σ(0) est symétrique alors σ^{n} sera symétrique pour tout n ∈ ℕ.

2.2 Discrétisation en espace

2.2.1 Traitement des conditions au bord

Le domaine est supposé rectangulaire et nous utilisons la méthode des différences finies. Le maillage (de type MAC) est formé d’une grille de N × M cellules uniformes où les inconnues sont prises soit au centre de la cellule (pour la pression, les termes de contrainte et le paramètre d’ordre) soit sur les milieux des cotés (pour la vitesse). De plus, pour adapter très facilement la méthode aux nombreuses conditions au bord envisagées (voir la partie 1.5), nous avons introduit des inconnues artificielles autour du domaine physique.

2.2.2 Schéma pour les termes de convection

En ce qui concerne les termes de convection (aussi bien sur φ que sur v ou σ) nous avons proposé une discrétisation en temps de type Runge-Kutta (voir partie 2.1). Une discrétisation spatiale de l’opérateur K(f) = u.∇f est nécessaire. De plus il est important que cette discrétisation soit L^{∞}-stable et positive (physiquement, le paramètre d’ordre φ doit toujours rester entre 1 et −1). Dans [2], nous proposons un schéma centré en espace muni de limiteur de flux de façon à vérifier de telles propriétés. Plus précisement, nous prouvons la

Proposition 2.1 Sous la C.F.L. suivante

\frac{dt}{δx} max_{i,j} (|u_{i+1/2,j}| + |u_{i−1/2,j}|) + \frac{dt}{δy} max_{i,j} (|v_{i,j+1/2}| + |v_{i,j−1/2}|) ≤ 1.  (4)

le schéma construit est L^{p}-stable pour tout p ≥ 1.

2.3 Propriété de stabilité du schéma exponentiel

Nous montrons aussi dans [2] que le schéma introduit pour discrétiser la loi constitutive (paragraphe 2.1) est stable sous certaines conditions. Le résultat est prouvé dans le cas du choix d’une dérivée corotationnelle pour la contrainte (a = 0) et dans le cas monophasique (φ = 1) :

Proposition 2.2 Le schéma exponentiel est L^{2}-stable sous les conditions (4) et

dt ≤ \frac{(1 − r)We^{2}}{4ηr^{2}} .  (5)

La condition de stabilité (5) n’est en général pas très contraignante pour des valeurs physiques de We et r. De plus, pour un fluide newtonien (r = 0) cette condition devient dt ≤ +∞.

3 Simulations numériques

3.1 Canal sous cisaillement

Nous présentons ici un exemple de test dans lequel le comportement d’un fluide visco-élastique est très différent de celui d’un fluide newtonien. Nous considérons un écoulement dans une cellule de Couette (ramenée ici à un canal bidimensionnel avec des conditions de périodicité sur deux des cotés : voir la figure 1).

Figure 1: Conditions au bord pour la cellule de Couette

Dans un premier temps, nous nous sommes intéressés au contexte monophasique pour lequel il est connu qu’un régime stationnaire est rapidement atteint (en initialisant avec une vitesse nulle). Pour ce régime stationnaire, la vitesse ne dépend que de la variable transverse y et cette vitesse est linéaire en y (voir [1]). Dans le cas visco-élastique, si on cherche des solutions stationnaires du même type (c’est-à-dire ne dépendant que de la variable transverse) nous montrons qu’il peut exister, selon les valeurs des paramètres r, a et We plusieurs solutions (toutes ont des profils affines par morceaux). Ces résultats mathématiques [2, 16] ont été vérifiés numériquement. Le test proposé ici correspond aux valeurs r = 0.9, We = 1 et a = 0. Nous montrons (figure 2) les profils de vitesse, la pression ainsi que le premier coefficient de la contrainte au temps t = 200. Un autre résultat mathématique que nous avons pu vérifier concerne l’instabilité du profil linéaire en vitesse. En effet, C. Guillopé et J.C. Saut [13] montrent que pour certaines valeurs de r, a et We le profil linéaire est instable. En prenant de telles valeurs, nous avons obtenu une déstabilisation du profil linéaire [2].

Une étude de la décomposition spinodale sous cisaillement a ensuite été réalisée. Cette expérience consiste à placer deux fluides (de densités identiques) dans un système de Couette-Taylor en imposant des conditions de cisaillement sur le champ de vitesse. La température est initialement portée sous la température critique ce qui impose la séparation des deux phases. Dans le cas newtonien, une fois les deux phases séparées, l’influence du cisaillement se fait sentir et les deux

Profil vitesse  Pression  Contrainte

Figure 2: Fluide visco-élastique sous cisaillement

phases s’organisent sous forme de bandes parallèles [1, 3].

Nous avons réalisé des tests identiques avec un mélange visco-élastique/newtonien. Initialement, les deux phases sont au repos et mélangées, les paramètres de visco-élasticité étant r = 0.9, We = 10 et a = 0. Il apparait que le temps de formation de bandes parallèles est nettement plus long dans ce cas que dans le cas purement newtonien. La figure 3 montre les résultats aux instants t = 10 et t = 100. Sur la gauche, nous montrons le paramètre d’ordre, le fluide “blanc” étant le fluide newtonien alors que le “noir” est le fluide visco-élastique. On voit alors que les bandes sont plus larges dans le cas visco-élastiques et que le profil de vitesse est non linéaire.

Mélange newtonien  Mélange visco-élastique

t=10

t=100

Figure 3: Décomposition spinodale sous cisaillement

3.2 Etirements de fibres

Cette partie a été réalisée en collaboration avec A. Colin et O. Greffier au C.R.P.P. (Bordeaux). L’expérience physique consiste à injecter un fluide dans un autre à l’aide d’une buse tout en tirant cette buse vers l’arrière (voir figure 4). Nous avons négligé l’effet de la gravité, les deux fluides ont des densités supposées égales alors que la viscosité du fluide injecté est 10 fois plus importante que celle de l’autre fluide [7]. Deux cas de figure peuvent être envisagés : soit le jet se brise en goutelettes, soit le jet produit un long filament. Les deux modes obtenus dépendent uniquement des deux vitesses imposées, la vitesse d’injection et la vitesse de déplacement de la buse. L’apparition de ces deux régimes peut être résumé dans un diagramme de phase. En ce qui concerne le cas newtonien (lorsque les deux fluides sont newtoniens) le diagramme de phase est le même dans le cas des expérimentations physiques et numériques. Nous avons obtenu le diagramme de phase pour des fluides non-newtoniens et ainsi comparé les résultats (figure 5).

Figure 4: Processus expérimental

Figure 5: Diagrammes de phase newtonien (gauche) et visco-élastique (droite)

D’autres aspects dus aux effets visco-élastiques ont été observés, par exemple la différence de forme du jet à la sortie de la buse dans le cas newtonien et dans le cas visco-élastique [14], l’apparition d’instabilités de Rayleigh à des vitesses relativements petites [2]...

3.3 Remplissage de cuve

Dans cette dernière expérience, un fluide (noir) est injecté dans une cuve contenant initialement l’autre fluide (blanc). La gravité est ici prise en compte et le fluide injecté est 1000 fois plus dense que l’autre fluide (rapport eau/air). Nous remplissons la cuve jusqu’au trois cinquièmes avant d’arrêter l’injection et d’observer la stabilisation.

Plusieurs tests ont été réalisés (voir [2]) : newtonien ou visco-élastique. Les résultats que nous présentons ici (figures 7 et 8) correspondent au cas newtonien et totalement visco-élastique (les deux fluides ont des propriétés d’élasticités), la vitesse d’injection est égale à 3 (nous stoppons donc le jet au temps t = 6). Dans le cas non-newtonien, de plus nombreuses bulles sont formées ce qui perturbe la stabilisation. Ce phénomène est bien connu des physisiens et est la cause de très nombreux problèmes dans l’industrie.

t=1  t=1.4  t=1.8  t=2.2  t=2.5

Figure 6: Comportement du jet visco-élastique dans le cas v_{inj} = v_{dis} = 1

t=3  t=6  t=7  t=9

Figure 7: Remplissage de cuve : cas newtonien

t=3  t=6  t=7  t=9

Figure 8: Remplissage de cuve : cas visco-élastique

Références

[1] F. Boyer, A theoretical and numerical model for the study of incompressible mixture flows, Computers and Fluids Vol. 31 no 1 (2002)

[2] F. Boyer, L. Chupin, P. Fabrie Numerical study of viscoelastic mixtures through a Cahn-Hilliard flow model, en préparation

[3] J.W. Cahn On spinodal decomposition, Acta. Metall., 9, pp 795-801 (1961)

[4] J.W. Cahn, J.E. Hilliard, Free energy of non-uniform system. I. Interfacial free energy, J. Chem. Phys. 28, pp 258-267 (1958)

[5] R. Chella, J. Vinals, Mixing of two-phase fluids by a cavity flow, Physical Review E 53, pp 3832-3840 (1996)

[6] L. Chupin, Existence result for a mixture of non Newtonian flows with stress diffusion using the Cahn-Hilliard formulation, Discrete and Continuous Dynamical Systems B, vol. 3, 1 , pp 45-68 (2003)

[7] A. Colin, O. Greffier et S. Manneville, communications personelles

[8] M. Doi, Dynamics of domains and textures, Theoretical Challenges in the dynamics of complex fluids, pp 293-314 (1997)

[9] J. Eggers, Nonlinear dynamics and breakup of free-surface flows, Rev. of Mod. Physics, Vol.60, No.3, pp 865-929 (1997)

[10] C. M. Elliott, H. Garcke, On the Cahn-Hilliard Equation with degenerate mobility, Siam J. Math. Anal. 27, n 1 , pp 404-423 (1996)

[11] P. Fabrie, Contribution à l’étude de la convection naturelle en milieu poreux, Thèse d’état, Université Bordeaux 1 (1987)

[12] M. Fortin, R. Glowinski, Augmented lagrangian methods: applications to the numerical solution of boundary-value problems, Vol 15 of Studies in Mathematics and its applications. North-Holland Publishing Comp. (1983)

[13] C. Guillopé, J.C. Saut, Global existence and one-dimensional nonlinear stability of shearing motions of viscoelastic fluids of Oldroyd type, Rairo model Math. Anal. Num., Vol. 24 n.3, pp 369-401 (1990)

[14] D. D. Joseph, Fluid Dynamics of Viscoelastic Liquids, Springer (1990)

[15] R. Keunings, Computer Modeling for Polymer Processing, C.L.Tucker III, Hanser Publishers (1989)

[16] H. Le Meur, Non-uniqueness and linear stability of the one-dimensional flow of multiple viscoelastic fluids, Math. Model. and Num. Anal., Vol. 31, 2, pp 185-212 (1997)

[17] J.G. Oldroyd, On the formation of rheological equations of state, Proc. R. Soc. Lond. A200, pp 523-541 (1950)

[18] A. Onuki, Phase transitions of fluids in shear flow, J. Phys. Condens. Matter 9, pp 6119-6157 (1997)