Résultats

Comparaison OpenFOAM-Fluent

Dans cette partie, nous nous sommes attachés à comparer de manière succintement les différences qui pouvaient ressortir d'une simulation menée par Fluent et d'une simulation menée par OpenFOAM. La finalité étant, bien entendue, de pouvoir justifier l'utilisation de OpenFOAM pour la simulation numérique de l'hydrodynamique en cuve agitée.

 

Cas de comparaison

La géométrie utilisée ici est celle du Cas Test 3D.

Les caractéristiques de la simulation sont les suivantes :

  • l'agitateur tourne à une vitesse angulaire de 104.72 $ rad.s^{-1} $
  • la simulation se fait avec la méthode MRF
  • le fluide possède les propriétés suivantes :
    • viscosité cinématique de 1 $ Pa.s $
    • masse volumique de 1 $ kg. m^{-3} $
  • le régime laminaire est choisi
  • les schémas de discrétisation sont standards.
  • le nombre d'itérations est fixé à 500 pour la convergence en régime stationnaire

 

Résultats de la comparaison Fluent-OpenFOAM stationnaire

La comparaison a été faite sur la vitesse en norme sur trois lignes verticales définies comme dans la partie Post traitement à l'adresse suivante http://hmf.enseeiht.fr/travaux/bei/beiep/content/g22/openfoam-0 et sur les nombres adimensionnels de pompage et débit.

La comparaison des courbes s'est faite sous matlab, également comme indiquée dans la partie post-traitement.

 

Voici deux des trois graphes :

 

La comparaison des vitesses est assez surprenante puisqu'en simulant deux fois le même écoulement en cuve, on devrait s'attendre à une superposition des courbes. La comparaison des contours illustre également cette différence.

Contours de vitesse sous fluent Contours de vitesse sous OpenFOAM

Etant donnée que les données d'entrée pour chacune des deux simulations étaient identiques, nous avons regardé ce qu'il en était au niveau des schémas numériques de discrétisation des dérivées. Les schémas étaient identiques pour chacune des dérivées.
Seul l'algorithme reliant la pression à la vitesse était différent : d'un côté l'algorithme SIMPLEC pour Fluent, de l'autre, l'algorithme SIMPLE.

Le SIMPLEC n'est qu'une amélioration de modèle SIMPLE de base. Il est censé être plus consistant. Après avoir relancé la simulation avec Fluent en SIMPLE, nous obtenons les résultats suivants.

Les différences s'estompent. Dans le cas R=0.4m, le maximum de vitesse sur Fluent est même plus important. L'analyse des contours de vitesse rend également compte de ce rapprochement de comportement.

Contours de vitesse sous fluent Contours de vitesse sous OpenFOAM  

 

En poussant un peu plus les calculs pour OpenFOAM jusqu'à mille itérations, les résultats se rapprochent.

Cette partie illustre bien le fait suivant, pour l'obtention de bons résultats, il ne faut pas hésiter à faire tourner la simulation sur un nombre important d'itération, lorsque le temps vous le permet.

 

Détermination des nombres adimensionnels

        OpenFOAM                           Fluent            
Nombre de pompage 0.1 2.6
Nombre de puissance 70 74

Le taux de dissipation visqueux sous Fluent est défini par une User Defined Function que l'on intégre par la suite avec l'option report - Volume Integral.

Pour le nombre de pompage, un cylindre (rayon R =0.26m) a été créé dans la zone MRF et nous avons demandé à Fluent de sortir le flow rate de la vitesse radiale sur cette surface.

L'ordre de grandeur des valeurs trouvées pour le nombre de puissance est correct. Cependant, pour l'écoulement en question dont Re = 4, nous devrions trouver un nombre de puissance de l'ordre de 20.

 

 

Comparaison des simulations en instationnaire

Dans cette partie, nous nous sommes attachés à comparer les résultats de la simulation sur la même géométrie Cas Test 3D en instationnaire.

Les caractéristiques de la simulation sont les suivantes :

  • l'agitateur tourne à une vitesse angulaire de 104.72 $ rad.s^{-1} $
  • méthode MRF
  • le fluide possède les propriétés suivantes :
    • viscosité cinématique de 1 $ Pa.s $
    • masse volumique de 1 $ kg. m^{-3} $
  • le régime laminaire est choisi
  • les schémas de discrétisation sont standards.
  • l'algorithme de calcul pression-vitesse est le SIMPLE

 

L'élément le plus intéressant était de savoir à partir de combien de tours serait atteint le régime stationnaire.

Pour la simulation OpenFOAM, nous avons affiché sur le même graphe les valeurs de la vitesse en norme sur deux des trois lignes verticales pour plusieurs temps (comptés en tours sur le graphe).

Pour la simulation avec OpenFOAM, il semblerait que l'état stationnaire soit atteint pour mille tours d'agitateur, soit environ soixante secondes.

 

Les résultats pour la simulation sur Fluent sont complètement différents.
 

Les maximums de vitesse en norme sont quasiment les mêmes, quelque soit le logiciel support de simulation. Par contre le temps  mis pour atteindre le régime stationnaire diffère : 84 tours sont nécessaires sur Fluent alors qu'il en faut quasiment 1000 sur OpenFOAM.

Cas de la géométrie réelle

Dans cette partie, nous nous attachons à réaliser le post-traitement défini dans la rubrique post-traitement ainsi que quelques comparaisons entre le cas stationnaire et le cas instationnaire.

 

Données d'entrée du calcul

Les simulations ont été lancées en régime laminaire avec un fluide industriel dans les conditions réelles . Ces caractéristiques sont les suivantes :

  • Fluide XX :
    • Densité 1500 $ kg. m^{-3} $
    • Viscosité dynamique 40000 Cp = 40 Pa.s
    • Viscosité cinématique $ 2,66.10^{-2} m^{2}.s^{-1} $
  • Paramétrage pour les dictionnaries
    • Méthode MRF
    • Rotation à - 4,39 rad/s (signe négatif pour produire un pompage de fluide vers le fond de la cuve
    • Régime laminaire : Re = 105
  • Maillage : 124000 cellules tetraédriques

 

Cependant afin de pouvoir illustrer nos post-traitements, nous avons besoin de réaliser des simulations en régime turbulent. Les données sont les suivantes :

  • Paramétrage pour les dictionnaries
    • Méthode MRF
    • Rotation à - 4,39 rad/s (signe négatif pour produire un pompage de fluide vers le fond de la cuve
  • Fluide en régime turbulent :
    • Densité 1500 $ kg. m^{-3} $
    • viscosité cinématique de $ 2,66.10^{-5} m^{2}.s^{-1} $
    • Re = 76000
  • Maillage : 330000 cellules tétraédriques

 

 

Post-traitement des résultats

Cas stationnaire

L'axe Y est l'axe de rotation de l'agitateur. Les zones rouges représentent les zones dans lesquelles le fluide développe un mouvement ascendant fort. Les zones bleues représentent les zones dans lesquelles le fluide développe un mouvement descendant fort.

On constate bien qu'il y a un pompage du fluide dans le sens des Y descendants le long de l'arbre d'agitation et un refoulement dans le sens des Y ascendants près des parois.

Les flèches représentent le vecteur vitesse en norme ; elles retranscrivent bien le mouvement dans le sens horaire du fluide.

 

Post traitement

Par l'utilisation des de OpenFOAM et Matlab comme décrit dans la partie post-traitement, nous avons cherché à visualiser quelques grandeurs pertinentes.

On rappelle que le tracé s'est fait sur les lignes dessinées sur l'image suivante.

 
Vitesse en norme à plusieurs altitudes

Les pâles ont un rayon d'environ 0.85m. Les mouvements de fluide ont lieu essentiellement dans les régions proches des extrémités de pâles.

 

Vitesse selon Y

On constate ici que le fluide remonte près des parois pour redescendre le long de l'arbre de rotation.

 

Taux de dissipation turbulente

La dissipation d'énergie turbulente est importante dans la région ou les fluctuations de vitesse sont importantes, à savoir près de l'extremité des pâles et entre la paroi et les pâles.

 

 

Les nombres caractéristiques issus du post traitement

Le cas du nombre de Reynolds est trivial. Il est donné en haut de page, le régime est turbulent, il légitime l'utilisation d'un modèle de turbulence.

 

Le nombre de puissance

\begin{equation}
N_p= \frac {P} {\rho N^3 D^5}
\end{equation}

La méthode de post-traitement donne une valeur de \begin{equation}  \iiint\limits_V  \mu \phi_v \, dv = 1138 W
\end{equation}

Attention : $ \phi_v $ le taux de dissipation visqueux est calculé à partir des composantes du tenseur des contraintes qui sont normalisées par la densité du fluide sous OpenFOAM. Il faut donc diviser non pas par $ \mu^2 $ mais par $ \nu^2 $ dans l'expression de $ \phi_v $ . Pour plus de précision, voir le détail dans la rubrique Post traitement.

Nous obtenons donc la valeur suivante : \begin{equation}  N_p= \frac {\iiint\limits_V  \mu \phi_v \, dv} {\rho N^3 D^5} = 0.16 \end{equation}

Par rapport aux valeurs de la littérature, le nombre de puissance trouvé se trouve dans la gamme des valeurs usuelles. La méthode de post traitement semble fonctionner correctement.

 

 

 

Le nombre de pompage

La détermination du nombre de pompage transférée au fluide selon paraview donne :

Le slice est sur l'image réalisé à une hauteur Y =0.75m, soit juste au dessus des pâles. Pour ne prendre en compte que le débit de pompage, le slice doit être fait juste en dessous de l'agitateur.

L'intégration sur le disque à une hauteur Y=0.15m donne un débit de pompage de $ Q_{P}  = 0.024 m^{3}.s^{-1} $ , soit un nombre de pompage de :

\begin{equation}
N_{Q_P} = \frac {Q_P} {N D^3} = \frac {0.024} {0.7* 1.7^3} = 0.007
\end{equation}

Avec :

$ N=0.7 tr .s^{-1} $
$ D= 1.7 m $

 

D'après les valeurs trouvées et répertoriées de manière non exhaustive dans la nomenclature de la partie Ecoulement (http://hmf.enseeiht.fr/travaux/bei/beiep/content/g22/ecoulement), le nombre de pompage est du bon ordre de grandeur bien plus faible. Cependant, on pourra noter que l'agitateur est à une hauteur de 24% de celle de la cuve ce qui est en dessous des hauteurs classiques de 33% de la hauteur de cuve pour lesquels les systèmes sont optimisés. De plus, le fond de cuve étant plat, la circulation du fluide ne se fait pas correctement comme indiqué sur cet URL http://hmf.enseeiht.fr/travaux/bei/beiep/content/g22/inventaire-mobiles-dagitation .

De plus, le débit de circulation, calculé sur un slice à une hauteur de 0.60m donne un débit de $ Q_{P}  = 0.38 m^{3}.s^{-1} $ soit \begin{equation}  N_{Q_C} = \frac {Q_C} {N D^3} = \frac {0.38} {0.7* 1.7^3} = 0.11 \end{equation}.
Selon les valeurs de la littérature on devrait être aux alentours de  1,2 à 1,5 $ N_{Q_p} $. Or ici, nous sommes à 15.7 $ N_{Q_p} $

Par conséquent, en se fiant à ces deux constats, nous pouvons en déduire que

  • soit le système est ici non optimisé.
  • soit la méthode de post-traitement n'est pas correcte. Au vu des résultats correct du nombre de puissance déterminé grâce aux mêmes outils sous paraview, on peut privilégier un défaut dans la conception de la cuve.

Il ne faudra pas oublier que les paramètres de simulation ne sont pas les meilleurs non plus. Entre autres, le modèle de turbulence k-Epsilon qui a été choisi pour des raisons liées à la vitesse de convergence des calculs et non pas pour sa capacité à modéliser les phénomènes en cuve agitée.

 

 

Cas instationnaire

 

Le nombre de puissance

Le graphe ci dessous représente l'évolution du nombre de puissance en fonction du nombre de tours de l'agitateur.

Deux éléments sont à noter :

  • l'évolution du nombre de puissance subit des évolutions notables dans la phase instationnaire de la simulation
  • l'état instationnaire n'est pas atteint avant au moins 120 tours, ce qui correspond à un temps de stabilisation d'au moins une minute et 25 secondes

 

Le nombre de pompage

Le graphe ci dessous représente l'évolution du nombre de pompage sur plusieurs surfaces d'intégration en fonction du nombre de tours de l'agitateur.

Selon la définition du nombre de pompage donnée dans la rubrique "Les mobiles d'agitation", il faut regarder la courbe bleue.

Tout comme pour le nombre de puissance, l'évolution de nombre de pompage subit d'importantes variations avant d'atteindre l'état stationnaire.

Discussion des résultats

Les résultats obtenues ne sont pas à prendre pour argent comptant.

En effet, les paramètres de simulation ont été choisis afin de pouvoir développer un post-traitement basé sur des simulations relativement rapides en de temps de calcul. Il n'est donc pas concevable de considérer les chiffres obtenus comme physiquement acceptables.

Plusieurs points peuvent aider à améliorer la qualité des résultats :

  • les fluides choisis ne sont pas représentatifs de la majorité des fluides rencontrés dans les procédés industriels ( viscosité cinématique de 1 $ Pa.s $ et masse volumique de 1 $ kg. m^{-3} $ ).
  • les maillages utilisés ne sont pas assez raffinés : 330000 cellules pour mailler une cuve de plusieurs mètres cubes et pour capter le taux de dissipation turbulent des petites échelles de la turbulence est très clairement trop peu.
  • la définition du nombre de puissance doit prendre en compte le taux de dissipation turbulent pour les très grands Reynolds.
  • dans le cas d'une prise en compte de la turbulence, les modèles de fermeture tels que le modèle k-Epsilon ne sont absolument pas adaptés aux écoulements tournant en cuve agitées. Un modèle k-Omega SST serait plus approprié voir peut être même un modèle k-Epsilon $ V^{2}F $. A voir