Simulation directe de l’écoulement dans un Homogénéisateur Haute Pression

 

 

BARRET Luc

Encadré par : ABBAS Micheline
 


Simulation de l'écoulement au sein d'un homogénéisateur haute pression. L'étude à pour but une meilleure compréhension et une quantification des mécanismes à l'origine de phénomènes de ruptures de gouttes dans un écoulement diphasique. L'enjeu est une optimisation des moyens de production d'émulsions, ici dans le domaine de la santé.


    

Présentation du projet

Le système industriel

L'installation que nous étudierons au long de cette étude est un homogénéisateur haute pression utilisé pour produire une émulsion.

Une émulsion est un mélange de deux substances liquides non-miscibles, macroscopiquement homogène mais microscopiquement inhomogène. L'un des deux fluides est ainsi présent en phase dispersée sous la forme de micro-gouttelettes au sein de l'autre fluide constituant une phase continue. Les émulsions sont utilisées dans de nombreux domaines tels que la cosmétique, l'agroalimentaire et, comme c'est le cas ici, le domaine de la santé. Elles permettent l'inclusion d'une grande variété de composants au sein d'une même solution.

En entrant dans le système, la solution comprend déjà une phase dispersée, et l'objectif sera ici de réduire la taille des gouttes présentes par rupture en générant un écoulement turbulent. Pour cela, on amène le fluide à s'écouler entre deux plaques très rapprochées avant d'entrer brusquement dans une cavité beaucoup plus large. Cet élargissement entraîne alors une forte augmentation du nombre de Reynolds lié à l'écoulement. La géométrie du système est la suivante (géométrie axisymétrique) :

 

 

Enjeux pour l'industriel

Si le diamètre initial et le diamètre final des gouttes sont connus, les phénomènes qui ont lieu au sein du système restent difficiles à quantifier, notamment les mécanismes de rupture et de coalescence des gouttes et leur influence respective sur les dimensions de ces dernières en sortie du système. 

Obtenir des gouttes suffisamment fines au sein de l'émulsion représente pour l'industriel un enjeux d'une grande importance. En effet, la taille des gouttes est déterminante pour des raisons d'efficacité de l'émulsion, mais aussi pour des raisons de stabilité dans le temps.

La simulation numérique est un outil privilégié pour étudier l'influence de différents paramètres sur un écoulement. Il est donc particulièrement adapter aux études d'optimisation de systèmes industriels complexes. 

 

Nos objectifs

Le but de notre étude est de caractériser l'écoulement au niveau de l'élargissement brusque et en aval de celui-ci. Nous nous intéresserons particulièrement aux grandeurs caractéristiques de la turbulence, celles-ci étant à l'origine des phénomènes de rupture des gouttes. Ce projet s'inscrit dans une démarche de plus grande envergure dont l'objectif est de déterminer l'influence de paramètres physiques et géométriques sur le processus d'émulsification en vu de son optimisation. Nous poserons ici les bases de ce travail en caractérisant les phénomènes observés dans l'écoulement pour un jeu de paramètres donné sans observer l'effet d'une variation de ces paramètres. 

Etude théorique préliminaire

Phénomènes à l'origine de la rupture

Les différents phénomènes intervenant dans la rupture

De manière générale, on peut distinguer plusieurs phénomènes intervenant dans la rupture de gouttes au sein d'un écoulement. La turbulence, ainsi que des effets visqueux et inertiels peuvent intervenir avec une plus ou moins grande importance selon la nature de l'écoulement. Les forces alors engendrées sont a comparer avec la tension de surface à l'interface entre la goutte et la phase porteuse qui tend à rétablir et conserver la goutte dans une forme sphérique et qui s'oppose ainsi à la déformation et à la rupture. Une telle comparaison permet de déterminer une taille de goutte critique.

 

Influence de la turbulence

Une approche basée exclusivement sur l'intensité de la turbulence peut donner de très bon résultats dans la prédiction des phénomènes de rupture dans la cas où celle-ci est suffisamment importante pour causer un large taux de rupture.

La grandeur permettant de quantifier l'efficacité de la turbulence dans le mécanisme de rupture de gouttes est le nombre de Weber, comparant les forces de pression statistiquement moyennées à l'échelle de la goutte et les forces de tension de surface. On peut exprimer le Weber turbulent de la manière suivante :

\begin{equation} We=\frac{\rho_{c} \overline{\delta u^{2} (x,d) }d}{\sigma} \end{equation}

où $\overline{\delta u^{2} (x,d)}$ est la moyenne du carré de la différence des fluctuations de vitesse sur l'échelle d'une goutte. On le calcul de la manière suivante :

\begin{equation} \overline{\delta u^{2} (x,d_{\alpha})}=\overline{max_\alpha [u_{x}(\textbf{x}+\textbf{d}_{\alpha}/2)-u_{x}(\textbf{x}-\textbf{d}_{\alpha}/2)]^{2} +[u_{y}(\textbf{x}+\textbf{d}_{\alpha}/2)-u_{y}(\textbf{x}-\textbf{d}_{\alpha}/2)]^{2}}\end{equation}

Avec $\textbf{d}_{\alpha}$ de norme d, orienté selon x ou y.

 

Autres mécanismes physiques intervenant dans la rupture

Une autre contribution au phénomène de rupture vient de la décélération du fluide dans la direction de l'écoulement. Le Weber associé à cette contrainte inertielle peut s'écrire :

\begin{equation} We_{I}=\rho (\frac{\partial U_{y}}{\partial y})^{2} \frac{d^{3}}{4 \sigma} \end{equation}

Le jet entraîne de forts gradients de vitesse moyenne qui peuvent participer à la rupture des gouttes. Cette contribution des effets visqueux s'exprime à l'aide du nombre capillaire :

\begin{equation} Ca=\mu_{c} \frac{\partial U_{y}}{\partial x}\frac{d}{2 \sigma} \end{equation}

avec $\mu_{c}$ viscosité de la phase porteuse.

Dans le cadre de notre étude, nous nous intéresserons particulièrement à l'influence de la turbulence et de la viscosité.

 

Étude dynamique de la goutte.

Le modèle d'étude des phénomènes de rupture de gouttes présenté précédemment ne tient pas  compte de la dynamique de la goutte. En effet, lorsque celle-ci est soumise à une perturbation extérieure, elle peut entrer en résonance et se déformer jusqu'à une amplitude qui l'entraîne à se rompre. Il s'agit donc pour mener une étude dynamique d'étudier les temps caractéristiques de la goutte (fréquence propre d'oscillation et taux de d'amortissement) ainsi que les champs de contraintes instantanés de l'écoulement.

Cette étude ne sera pas menée dans le cadre de notre travail.

Etude de la turbulence

La turbulence est un phénomène fondamentale au sein d'un processus d'émulsification. La connaissance des échelles de la turbulence renseigne sur l'efficacité de cette dernière à venir briser les gouttes présentes dans la phase dispersée. Il s'agit donc ici de proposer une étude de la turbulence engendrée par le jet en sortie de restriction dans l'optique d'en comparer les échelles caractéristiques avec les dimensions des gouttes.

 

Estimation des échelles spatiales

L'estimation des échelles caractéristiques de la turbulence se fait de proche en proche à partir de grandeurs liées à l'énergie plus facilement approchables. Ainsi, le flux d'énergie nous est donné par la relation :

\begin{equation} E=\frac{\rho U_{0}^{2}}{2}.Q=\frac{\rho U_{0}^{2}}{2}.U_{0}.h.B \end{equation}

avec :

  • U0 : vitesse moyenne d'entrée
  • Q : débit d'entrée
  • h : hauteur de la cavité d'entrée
  • B : largeur de la cavité d'entrée (ici arbitraire)

On a de plus pour le taux de dissipation par unité de masse la relation suivante :

\begin{equation} \epsilon=\frac{E}{V_{diss}.\rho}=\frac{U_{0}^{2}.h.B}{2.A_{diss}.B} \end{equation}

avec :

  • Vdiss : volume de dissipation d'énergie
  • Adiss : aire de dissipation d'énergie

Par ailleurs, on peut estimer la grandeur Adiss en estimant les longueurs selon x et y sur lesquelles le fluide perd la majeur partie de sa vitesse, et donc de son énergie cinétique.

On peut alors estimer l'échelle de Kolmogorov à l'aide de l'expression :

\begin{equation}\eta=(\frac{\nu^{3}}{\epsilon})^{\frac{1}{4}}\end{equation}

et les grandes échelles de la turbulence à l'aide de l'équation :

\begin{equation}l_{0}=\eta.Re^{\frac{3}{4}}\end{equation}

avec $Re=\frac{\rho.h.U_{0}}{\mu}$

Simulations

Présentation du logiciel Jadim

Jadim est un logiciel de recherche développé par le groupe Interface de l'Institut de Mécanique des Fluides de Toulouse. Ce code permet l'étude des paramètres intervenant au sein d'écoulements diphasiques. Il résout les équations de Navier-Stokes en 3D dans le cadre de l'hypothèse de fluides incompressibles et pour des cas instationnaires.

Jadim utilise un maillage structuré et permet d'effectuer les simulations en DNS ou en LES

Hypothèses simplificatrices

Nous n'étudierons pas directement les phénomènes de rupture et de coalescence. En effet, l'étude de ces phénomène demande le suivi lagrangien des gouttes. Nous nous restreindrons dans le cadre de cette étude à l'étude eulerienne de l'écoulement.  Les gouttes ne sont pas ici modélisées et l'écoulement simulé est donc monophasique. Il s'agit donc ici de prendre les propriétés physique d'un fluide équivalent. Nous ne prendrons en compte aucun effet thermique. 

​Propriétés physiques

Les propriétés du fluide sont prises égales à celles de l'eau, c'est à dire : 

µf = 10-3

$\rho_{f}$= 103

​Contraintes

Il existe un certain nombre de contraintes inhérentes à l'outil numérique. Notamment, le respect des critères de stabilité entraîne une contrainte sur les dimensions des mailles et sur le pas de temps, ce qui peut entraîner des temps de calcul très important. 

Nous garderons également à l'esprit que le maillage doit permettre de prendre en compte toutes les échelles impliquées dans les phénomènes que l'on souhaite observer. Un compromis entre précision et temps de calcul doit donc être trouvé.

 

Première simulation - jet centré

Présentation

Présentation de la simulation

Cette simulation est une première approche du problème nous permettant de savoir s'il est possible de simuler directement l'écoulement dans les conditions industrielles et dans le cadre des hypothèses d'un écoulement monophasique incompressible. Ce premier cas nous a permis de poser par la suite d'autres hypothèses simplificatrices essentielles. Le calcul est effectué en 2D axisymétrique et tourne en DNS.

Une première simulation aux conditions industrielles demande cependant plusieurs étapes intermédiaires. Un premier cas pour V = 1m.s-1 nous permet de vérifier que le maillage généré permet au code de tourner. Il s'agit ensuite d'augmenter la vitesse jusqu'à la vitesse réelle d'écoulement. Cependant, l'essentiel restera d'obtenir un cas qui ne diverge pas, ce qui peut demander certains ajustements.

Géométrie

On gardera pour cette étude les dimensions du système industriel. La géométrie du système est la suivante : 

Géométrie industrielle

On ne modélisera cependant que la cavité dans laquelle s'écoule le jet au niveau de l'élargissement brusque

Domaine d'écoulement du jet

On imposera un profil de type Poiseuille en entrée du domaine plutôt que de simuler l'écoulement en amont.

Maillage

Le maillage utilisé pour cette simulation est raffiné en arctangente au niveau de la paroi d'entrée et autour du jet. La zone de jet, de même épaisseur que l'entrée, est maillée de manière uniforme de façon à obtenir 50 mailles dans les 10 micromètres de l'entrée.  Ce raffinage doit permettre de capter les forts gradients de vitesse en entrée causées par les contraintes de cisaillement et au sein du jet où la vitesse décroît rapidement et  dans lequel on observe une forte intensité de la turbulence.

 

Maillage du domaine

Maillage de l'entrée

Résultats et commentaires

Validation du maillage

Afin de valider le maillage que nous avons généré à l'aide du logiciel matlab, nous avons procédé à une première simulation très brèves pour une vitesse d'entrée de 1m.s-1. Nous avons ainsi pu vérifier le bon comportement du fluide aux limites du domaine, notamment en observant la présence d'un Poiseuille en entrée.

 

Champ de vitesse - test des conditions limites

 

   

  Vitesse en entrée du système - t=0.1 ms     Vitesse en sortie de système - t=0.1 ms

Le profil observé en entrée est une extrapolation entre le profil Poiseuille en première maille (non-observable) et la maille calculée suivantes. Ainsi on n'observe pas directement un profil de Poiseuille.

A la paroi limite nord du domaine, on observe bien une vitesse nulle à la paroi et une vitesse sortante à la sortie.

Test aux conditions industrielles

Ainsi assurés de la validité du maillage, nous avons pu lancer une seconde simulation pour une vitesse de 200m.s-1 correspondant à la vitesse d'entrée du fluide dans le système industriel. Le pas de temps est extrêmement faible, de l'ordre de 10-10s. Le calcul tourne alors 1 060 000 itérations (0.104 ms) avant de diverger.

Champ de vitesse dans le domaine - le calcul diverge

Dans le but de faire converger le calcul, nous avons donc décidé de diminuer la vitesse d'entrée du fluide. Cependant, pour conserver la physique de l'écoulement qui nous intéresse, il nous a fallu conserver le même nombre de Reynolds pour toutes les simulations.

Conservation du Reynolds

Afin de diminué la vitesse d'entrée du fluide tout en conservant le Reynolds de l'écoulement, nous avons choisi de diviser la vitesse et la viscosité du fluide par un même facteur. Ainsi nous avons choisi arbitrairement de diviser ces deux quantités par 10. Toutes les simulations réalisées par la suite l'ont été pour une vitesse d'entrée de 20m.s-1

Par ailleurs, on observe que le calcul à V=200 m.s-1 diverge lorsque les vitesses deviennent importantes en sortie. Dans cette zone, les vitesses passent d'une valeur nulle en paroi à une valeur très importante à "l'angle" de la conduite de sortie.  Dans la suite de notre travail, nous avons donc cherché à réduire l'importance de cette contrainte, lourde pour un code CFD.

Seconde simulation - jet décalé

Présentation

Démarche

Afin d'éviter d'obtenir de trop fort gradients de vitesse en sortie, nous avons décidé de décaler la position du jet de manière à éloigner la zone impactante de la sortie.

Maillage

Nous avons raffiné d'avantage le maillage au niveau des parois et notamment au niveau de la paroi opposée à l'entrée où le jet vient impacter. Le nombre de mailles selon y a été changé pour passer de 64 à 128 mailles. Ainsi nous obtenons le maillage suivant :

Maillage raffiné aux parois - jet décalé

 

 

Résultats et commentaires

Simulation du jet décalé

Le fait d'avoir décalé le jet plus loin de la sortie et d'avoir modifié le maillage nous a permis d'obtenir un écoulement beaucoup plus établi que précédement. Le calcul diverge au bout de 228 000 itérations seulement, mais le pas de temps calculé par Jadim oscille autour de 10-8 secondes au lieu de 10-10 précédement.

Champ de vitesse dans le domaine au moment où le calcul diverge

 

On peut ici observer les grands tourbillons de recirculation après impact à la paroi. A nouveau, le calcul diverge lorsque les gradients de vitesse deviennent importants au niveau de la sortie.

Cette simulation nous conforte dans l'idée que la divergence du calcul provient de variations brutales de la vitesse au niveau de la sortie. Ainsi, il s'agira dans le suite du projet d'atténuer ces brusques variations au niveau de la condition de sortie.

Troisième simulation - sortie modifiée

Présentation

Démarche

Afin d'éviter d'obtenir des gradients de vitesse en sortie trop importants, nous avons décidé de modéliser une partie de la conduite de sortie afin que le profil de vitesse ait le temps de s'homogénéiser en partie avant d'atteindre la condition de sortie.

 

Maillage

Toute la partie du domaine dont la géométrie reste inchangée demeure maillée de la même manière. La partie ajoutée en sortie apporte donc des mailles supplémentaires, raffinée en arctangente pour ne pas entraîner de "saut" dans les dimensions des mailles en entrée de conduite de sortie. On obtient ainsi le maillage suivant :

Maillage du domaine - sortie modifiée

Toute la partie du maillage qui a été ajoutée et qui n'est pas comprise dans la conduite de sortie n'est pas prises en compte par Jadim. Les équations ne sont pas résolues dans ce domaine. Ainsi le maillage "effectif" à considérer est le suivant :

Maillage du domaine de résolution des équations

 

Résultats et commentaires

La modification de la sortie a permis au calcul de converger. Nous avons cependant du l'arreter au bout d'un temps relativement court à cause des limites de temps du projet. Nous avons donc laisser le calcul tourner 764 000 itérations, ce qui correspond à 7.5ms de temps physique. Le résultat de la simulation et visible sur la vidéo suivante :

On observe bien les boucles de recirculation qui s'étendent jusqu'aux extrémités du domaine. D'autre part, nous pouvons nous intéresser aux longueurs de décroissance de la vitesse selon x et y afin d'estimer les différentes échelles de la turbulence comme expliqué au cours de l'étude théorique. 

Considérons tout d'abord la décroissance avec y de la vitesse. On se place en x égale au x du centre du jet. On trace ensuite la vitesse moyenne du fluide en fonction de y ; on obtient alors la courbe suivante : 

U moyen en fonction de Y en X = X(centre du jet)

On observe alors que la vitesse décroît de moitié sur une longueur Ydiss​ = 3.59 e-3 m. 

Pour observer la décroissance de la vitesse avec x, la question de la position en y se pose. En effet, la longueur selon x que met le fluide à perdre la moitié de sa vitesse dépend de la position en y choisie. Afin de nous placer au centre de notre domaine de dissipation, nous avons choisit de nous placer en Y=Ydiss/2. En traçant la vitesse moyenne du fluide en fonction de x, on obtient alors : 

U moyen en fonction de X en Y<Ydiss/2 (bleu), Y=Ydiss​/2 (rouge) et Y>Ydiss/2 (vert)

On trouve alors une valeur de Xdiss​=1e-5 m, c'est à dire égale à la largeur du jet. 

En reprenant les formules énoncées dans l'étude théorique, on trouve alors : 

  • Adiss=3.59 e-8 m2
  • $\epsilon$ = 5.58 e4 m2.s-3
  • $\eta$ = 2.06 e-6 m

Le maillage, en entrée du domaine, comprend des mailles de 2e-7 m selon x et de 3e-6 m selon y. Il devrait donc permettre de capter de manière satisfaisante les plus petites échelles de la turbulence. On se propose a présent d'observer le champ de vitesse du fluide en entrée du domaine : 

Champ de vitesse en entrée du domaine - visualisation des petits tourbillons (en noir)

On devine alors certaines structures tourbillonnaires, cependant la finesse ne permet pas une bonne visualisation des plus petites structures. 

D'autre part, en nous basant sur le Reynolds d'entrée et de la valeur de $\eta$ que nous venons de trouver, nous avons pu évaluer la dimension des plus grandes structures de la turbulence. On trouve alors : l0=1.1e-4 m. Visualisons à présent le champ de vitesse à l'échelle du domaine : 

Champ de vitesse dans le domaine - visualisation des grands tourbillons (en noir) - axes en mm

Les dimensions des plus grands tourbillons sont ici de quelques millimètres. Cet écart peut être du à la différence entre le Reynolds d'entrée que nous avons pris pour évaluer l0 et le Reynolds de l'écoulement, difficilement évaluable. 

 

Post-traitement

Grâce notre dernière simulation, nous avons pu établir des moyennes temporelles sur les champs de vitesse et ainsi remonter à des champs de grandeurs caractéristiques de la rupture de gouttes. Notamment, nous avons pu établir les lieux des ruptures, quantifier les grandeurs caractéristiques et les comparer à d'autres travaux effectués dans ce domaine. 

Pour observer l'influence de la turbulence sur la rupture, nous avons tracé sur matlab les contours de We : 

We dans le domaine

Pour mieux observer la répartition du We au niveau du jet, nous avons tracé le We dans cette zone uniquement : 

We au niveau du jet

On observe alors que le We est important au niveau de l'entrée et 1mm au dessus de l'entrée. Ce résultat est à mettre en relation avec les résultats obtenus par R. Maniero dans son travail sur la rupture en conduite en aval d'une restriction ; on y observe une certaine distance avant que la rupture ne se produise, distance partiellement retrouvée dans les simulations :

 

Résultats obtenus par R. Maniero

​Les ordres de grandeur des We restent cependant difficilement comparable, puisque notre valeur de la tension de surface a été prise arbitrairement.

​Nous avons opéré un travail similaire pour les effets visqueux en étudiant le nombre capillaire Ca. Celui-ci est concentré au niveau de l'entrée, où le cisaillement est très fort, et négligeable dans le reste du domaine. On se propose donc d'observer le Ca en entrée du domaine : 

Ca en entrée du domaine

Ces résultats sont cohérents avec la théorie. 

Conclusion

Conclusion

Ce travail a permis la mise en place d'une étude approfondie de la turbulence et des phénomènes de rupture au sein d'un homogénéisateur haute pression. Ce projet s'inscrit dans une démarche de plus grande ampleur et pose les bases d'un travail qui a vocation à être poursuivi avec différents jeux de paramètres. Notamment, les problèmes que nous avons rencontré au niveau de la condition de sortie ont été résolus et nous avons pu mettre en place une géométrie qui permet au code de converger. 

D'autres part, les résultats obtenus sont cohérents avec d'autres travaux effectués dans le domaine de la rupture de gouttes en écoulement diphasique turbulent et avec la théorie, ce qui valide en partie notre méthode. 

Bibliographie

[1] Modeling and simulation of inertial drop break-up in a turbulent pipe flow downstream of a restriction -  Riccardo Maniero - 2012

[2] Analysis of the flow field in a high-pressure homogenizer - Fredrik Innings, Christian Trägardh - 2007

[3] Etude expérimentale de la rupture de gouttes dans un écoulement turbulent - Sophie Galinat - 2005