Optimisation d’une méthode de calcul pour les cuves agitées avec le logiciel OpenFOAM

 

 

Ligne de courant dans une cuve agitée de type axial

 

FOREST Matthieu - TANG Ing Yii

3e Année Elèves Ingénieurs ENSEEIHT (Mécanique des fluides)

Encadré par : LEGENDRE Dominique et TEKAM Sydney

(INP et TECH-AM INGENIERIE)


OpenFOAM - Paraview - Cuves agitées - MRFSimpleFOAM

Hélices - Moving Mesh - Matlab - SRFSimpleFOAM - Turbines

Fluent - Surfaces non conformes
 


Objectifs :

Recherches bibliographiques pour la caractérisation des systèmes agités
existants dans le monde industriel

Développement de méthodes de post traitement des données avec OpenFOAM, Parview et Matlab

Comparaison de différentes stratégies de simulation de systèmes agités avec OpenFoam et Fluent (SRF-MRF (in)stationnaire, moving mesh)


           

Présentation du sujet de BEI

Le sujet de BEI que nous avons choisi a été proposé par M. Sydney Tekam de l'entreprise TECH-AM Ingénierie.

http://techam-ingenierie.fr/

 

La société TECH-AM ingénierie propose des services en simulations numériques de mécanique des fluides auprès de clients dans de multiples domaines technologiques et industriels.
Les projets de développement sont de plus en plus complexes et requièrent une attention importante. L'expertise CFD de l'entreprise permet de réaliser des projets ambitieux de diverse nature.
De la conception de géométrie (design Modeler, T-grid,...) à la simulation numérique (Ansys Fluent, OpenFoam, Star CD) en passant par le maillage (blockMesh, Meshing), TECH-AM Ingénierie propose des solutions pouvant répondre à des problématiques complexes.
Que ce soit sous windows ou sous linux, TECH-AM Ingénierie développe également des codes en C, fortran, Scheme ainsi qu'au format UDF Fluent.
L'entreprise dispose de structures leur permettant de faire des essais en soufflerie et sur banc d'essai moteur et consacre également une partie de son activité à la veille technologique.

 

Contexte du BEI :

Les systèmes agités sont omniprésents dans les procédés industriels. Du fait de leurs dimensions et de leur construction, il est difficile de mener des campagnes d'essais pour optimiser leur conception. Le dimensionnement empirique et/ou la simulation numérique sont fréquemment utilisés pour la conception de ces agitateurs. Cependant la simulation numérique peut s'avèver très coûteuse en temps et en argent. L'utilisation de certains logiciels commerciaux fonctionnent sur le principe de licences et jetons. L'achat de licences permet de disposer de jetons qui, eux-mêmes permettent l'utilisation des logiciels. Dans le cas du logiciel commercial Fluent, un jeton est équivalent à une simulation sur un processeur. Pour des géométries complexes nécessitant un maillage raffiné  et une parallélisation des calculs, on comprend facilement que le coût d'une simulation puisse devenir très important.

Des alternatives se développent afin de proposer des solutions compétitives gratuites. C'est le cas du logiciel OpenFOAM(Open Field Operating and Manipulation), dit open source puisqu'il est gratuit dans son utilisation. L'intérêt de ce logiciel réside dans la possibilité de pouvoir implémenter facilement des solveurs, des modèles de turbulence ou des modèles thermodynamiques correspondant à une étude spécifique. Toutefois, avant de se lancer dans de la programmation en C++, il est conseillé de s'appuyer sur l'un des multiples solveurs fournis avec OpenFOAM. En effet, la bibliothèque de solveurs de OpenFOAM est très complète et permet de réaliser des simulations dans une large gamme de domaines comme les écoulements incompressibles, les échanges thermiques et les milieux poreux.

TECH-AM Ingéniérie souhaite développer des outils de post-traitement et réaliser des études sur la faisabilité des simulations numériques de cuve agitées sur OpenFOAM.

 

Sujet du BEI :

Le but de ce BEI est multiple. Une recherche bibliographique a été menée afin d'établir une nomenclature des mobiles à écoulement sur laquelle puisse s'appuyer TECH-AM Ingénierie. Le développement des méthodes de post-traitement sur des cas test et leur application à des cas réel est envisagée. Et enfin, une comparaison de différentes méthodes de simulation sous OpenFOAM vise à comprendre les limites de certaines de celles-ci faces aux problématiques des systèmes agités.

 

La prochaine partie a pour but de renseigner les principales informations sur les mobiles à écoulement que l'on retrouve dans les procédés industriels.

 

Les mobiles d'agitation

L'agitation :

L'opération de mélange est requise à de très nombreux stades des procédés. Sans cette agitation, un mélange placé en cuve serait gouverné par une loi de diffusion de Fick. Pour accélérer les processus de transfert au sein du mélange il faut conférer un mouvement au fluide afin de renouveler les contacts entre phases. C'est l'agitateur qui va entretenir ce mouvement.

Les mobiles d'agitations peuvent être de géométrie différente. Les caractéristiques des mobiles sont déterminantes dans le type d'utilisation qui veut en être fait. Certains agitateurs sont utilisés pour homogénéiser des mélanges avec ou sans différence de viscosité, de concentration ou de température. D'autres sont conçus de sorte à pouvoir mélanger des fluides très peu miscibles afin de permettre des réactions chimiques en phase hétérogène ou pour créer des émulsions. Il est aussi possible de rencontrer des agitateurs dont le rôle sera de permettre une mise en suspension de particules solides dans un liquide. Par exemple, la dissolution, la cristallisation ou encore le raclage du fond de cuve pour éviter l'amassement de particules. Enfin, la dispersion de gaz dans un liquide peut aussi être assurée par un agitateur pour des opérations d'absorption.

Pour plus d'information, voir Agitation et mélange - Aspect fondamentaux et applications industrielles de Catherine Xuereb, DUNOD

Inventaire des mobiles d'agitation

Les mobiles d'agitation peuvent être classés en fonction de l'écoulement qu'ils créent.

On distingue trois grandes catégories :

  • les mobiles à écoulement axial
  • les mobiles à écoulement radial
  • les mobiles à écoulement tangentiel
       Axial                Radial     Tangentiel

 

 

Les mobiles à écoulement axial :

Ces mobiles plus communément appelés hélices créent un écoulement axial dans la cuve. Les hélices sont généralement utilisées pour homogénéiser des mélanges qui soient à faibles voire moyenne viscosité (de l'ordre de 1 Pa.s). On les trouve dans les industries agroalimentaire et pétrolière pour de la mise en suspension de solides, pour de la cristallisation ou encore pour de la création de dispersion liquide-liquide.

Hélice classique     Hélice marine    Hélice à pale mince Hélice à double flux

Sources : http://www.mixel.fr
http://www.techniques-ingenieur.fr/res/media/docbase/table/sl2981801-web/SL2981801TBL-web.xml

 

Les mobiles à écoulement radial :

Ces mobiles sont plus simplement appelés turbines. Elles sont utilisées dans des procédés tels que la fermentation, les réactions chimiques gaz-liquide ou encore dans la fabrication d'émulsion. Ces processus nécessitent un apport d'énergie lequel est assuré par la rotation très rapide d'agitateur au design permettant la création d'écoulements fortement turbulents.

    Turbines à disque (de Rushton) 

    Turbines sans disque à pâles    
droites ou incurvées

 

 

Les mobiles à écoulement tangentiel :

Ces mobiles sont plus communément appelés ancre ou barrière. Leur structure particulière permet de racler les fluides très visqueux aux parois de cuves afin de s'assurer de la non formation de zones mortes et d'assurer les échanges thermiques pour certains procédés où l'influence de la température est importante.

Ancres Barrières

 

On pourra être amené à rencontrer d'autres mobiles à écoulement dans l'industrie. Ce ne sont pas les plus communs mais ils existent. Entre autres, les vis et rubans hélicoïdaux conçus pour assurer le mélange de fluide très visqueux par un pompage axial ascendant et les agitateurs avec rotor/stator conçus pour mélanger par la création de cisaillement très fort.

Rubans hélicoïdaux et Vis

 

La partie suivante donne quelques éléments de caractérisation des mobiles d'agitation.

 

Géométrie de cuve

En général, les géométries de cuve sont structurées et leurs dimensions sont standardisées.

Sachant que :

$ D $ est le diamètre de l'agitateur
$ T $ est le diamètre intérieur de la cuve
$ C $ est la hauteur d'implantation du mobile
 

Type d'agitateur le ratio $ D/T $
Agitateur axial $ 0,2<D/T<0,7 $
Agitateur radial $ 0,2<D/T<0,5 $
Agitateur raclant 0,95

Tableau ci-dessous montre le ratio standard entre le diamètre de l'agitateur
et le diamètre de la cuve pour des différents types de cuves.

 

Il est important de bien choisir la hauteur d'implantation du mobile dans la cuve pour avoir un écoulement bien mélangé. La hauteur de fluide est égale au diamètre de la cuve dans la plupart des cas.

Pour les mobiles axial ou radial, le ratio entre la hauteur d'implantation du mobile et le diamètre de la cuve, $ C/T $ est toujours $ 1/3 $.

Si la hauteur de liquide est supérieure à $ 1,2 T $, il faut opter pour un système multi-étagé comportant deux mobiles ou plus sur l'arbre.

Caractérisation des mobiles

L'agitation de fluides en cuves est une opération qui peut être relativement simple à réaliser mais toujours complexe à caractériser à cause de la nature des écoulements et de la complexité des systèmes cuve-agitateur. La caractérisation par des grandeurs globales moyennées dans le temps devient nécessaire afin de pouvoir quantifier les performances des systèmes et de les dimensionner ensuite.

Afin de déterminer les nombres adimensionnels caractérisant une cuve agitée, une analyse dimensionnelle est nécessaire.

 

La mise en mouvement d'un liquide en vue de son mélange avec un autre liquide, miscible ou non, ou bien en vue de la création d'une suspension solide, ou encore pour favoriser la dispersion d'une phase gazeuse demande de l'énergie.

 

 

Analyse adimensionelle

Les paramètres physiques qui entrent en jeu :

  1. N : la vitesse de rotation de l'agitateur en tours/seconde ; $[T^{-1}]$
  2. D : le diamètre de l'agitateur ; $[L]$
  3. $ \rho $ : la masse volumique du fluide ; $[ML^{-3}]$
  4. P : la puissance dissipée dans le fluide ; $[ML^2 T^{-3}]$
  5. g : l'accélération de la gravité ; $[LT^{-2}]$
  6. $ \mu $ : la viscosité dynamique du fluide ; $[ML^{-1} T^{-1}]$

On a donc : n = 6 avec n est le nombre de paramètres physiques.

Les unités fondamentales sont la masse [M], la longueur [L] et le temps [T]. On a donc p = 3 avec p est le nombre des unités fondamentales.

Selon le théorème de Vaschy-Buckingham, une équation impliquant six variables, exprimées à l'aide de trois dimensions, peut être réduite à une relation impliquant $ n-p = 6-3 = 3 $ termes sans dimensions et indépendants.

Pour un phénomène représenté par trois groupes adimensionnels, l'équation exprimant la relation entre les variables possède alors une solution de la forme :

\begin{equation}
f(\pi_1 , \pi_2, \pi_3) = 0
\end{equation}

On a trois groupes adimensionnels qui s'expriment comme ci-dessous :

\begin{equation}
\pi_1 = N^{x_1} D^{y_1} \rho^{z_1} P
\end{equation}

\begin{equation}
\pi_2 = N^{x_2} D^{y_2} \rho^{z_2} g
\end{equation}

\begin{equation}
\pi_3 = N^{x_3} D^{y_3} \rho^{z_3} \mu
\end{equation}

Pour que $ \pi $ reste adimensionnel, il faut que la somme des exposants de chaque grandeur fondamentale soit nulle. Pour l'équation (2):

\begin{equation}
0 = T^{x_1}  L^{y_1}  M^{z_1} L^{-3z_1}  ML^2 T^{-3}
\end{equation}

On trouve donc :

\begin{equation}
x_1=-3  ;  y_1=-5 ;  z_1=-1
\end{equation}

D'où l'expression du nombre adimensionnel :

$$ \pi_1= \frac {P} {\rho N^3 D^5} = N_p $$

$ N_p $ est appelé nombre de puissance.

Avec la même méthode, on trouve pour $ \pi_2 $ :

$$ \pi_2=\frac {g} {N^2 D}= \frac {1} {Fr}$$

$ Fr $ est appelé nombre de Froude, Fr.

Pour $ \pi_3 $, on a :

$$ \pi_3  =\frac {\mu} {\rho N D^2}= \frac {1} {Re} $$

$ Re $ est appelé nombre de Reynolds.

 

 

Nombre de Reynolds

Le nombre de Reynolds exprime le rapport des forces d'inertie sur les forces de viscosité. Si on a un nombre de Reynolds faible, c'est-à-dire que les forces de viscosité prédominent, l'écoulement est appelé laminaire. Un nombre de Reynolds élevé indique que les forces d'inertie prédominent dans l'écoulement, le régime est dit turbulent.

$$ R e  =\frac {\rho N D^2} {\mu} $$

 

 

Nombre de puissance

$ N_p $ est appelé nombre de puissance. Il caractérise l'énergie transmise au fluide par le système d'agitation, c'est-à-dire l'énergie adimensionnalisée dissipée dans le fluide.

$$ N_p = \frac {P} {\rho N^3 D^5} $$

Pour être un peu plus précis sur la définition de $ P $, $ P $ s'exprime en fonction de  $ P_a $ la puissance fourni à l'arbre pour tourner et $ P_0 $ la puissance à vide selon :

$$ P =  P_a - P_0 $$

Allure des courbes de puissance pour des différents mobiles ( Agitation. Mélange - Caractéristiques des mobiles d'agitation publié par Michael Rushton le 10/03/2005 )

A partir de données expérimentales, les courbes de puissance peuvent être établies. La courbe exprime la variation du nombre de puissance en fonction du nombre de Reynolds et tous les deux axes sont en coordonnées logarithmiques.

  1. Pour $ Re < 10 $, $ N_p $ dépend linéairement du nombre de Reynolds, Re. Dans cette zone, le régimeq est laminaire, le fluide de viscosité très élevée et le produit de $ N_p Re $ est égal à une constante qui est souvent notée A. Cette constante A est définie pour chaque ensemble agitateur-cuve.
  2. Pour $ Re > 10^4 $, le nombre de puissance est invariant avec Re. On est dans le régime tubulent, le fluide est peu visqueux. On peut aussi déduire que, dans un système données, le nombre de puissance est indépendant de la vitesse de rotation. Lorsque le système comporte des chicanes, le nombre de puissance est plus élevé parce qu'en effet, le rôle principal des chicanes est de casser l'écoulement tangentiel induit par la rotation. Elles transforment le mouvement tangentiel en un mouvement axial et c'est pour cela que l'on a une légère augmentation de la puissance consommée.
  3. Pour $ 10<Re<10^4 $, on est dans une zone de transition entre le régime laminaire et le régime turbulent.

 

Régime d'écoulement turbulent

Quelques valeurs du nombre de puissance en régime turbulent pour des différents types de mobiles.

Type de mobile Nombre de puissance $ N_p $
Ancre 0,2
Hélice 0,3 à 1
Turbines à pales inclinées 1,2 à 2
Turbine à disque type Rushton
Turbine à pale concave
4 à 5

Deux facteurs majeurs sont à l'origine d'une consommation énergétique élevé : la surface entre les pales et le fluide et le profil des pales. Par exemple, la turbine à disque type Rushton consomme beaucoup d'énergies parce qu'elle génère un mouvement intense qui crée un cisaillement très fort.

 

Régime d'écoulement laminaire

Comme annoncé précédemment, le nombre de puissance pour le régime laminaire s'écrit comme ci-dessous :

\begin{equation}
N_p Re=A =cste
\end{equation}

D'après l'analyse adimensionnelle, on peut aussi écrire la formule du nombre de puissance comme :

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

En reportant l'équation (7) dans l'expression (8), on obtient donc :

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

En remplaçant le nombre de Reynolds, l'expression se réduit à :

\begin{equation}
P=A \mu N^2 D^3
\end{equation}

 

 

Détermination du nombre de puissance

Numériquement, on peut déterminer le nombre de puissance par deux méthodes :

 

La première méthode :

L'hypothèse : On considère que toute la puissance du fluide est transmise au fluide. On peut déterminer la puissance comme :

\begin{equation}
P= \iiint\limits_V  \mu \phi_v \, dv
\end{equation}

avec $ \phi_v $ est le taux de dissipation visqueux et $v$ est le volume de la cuve. On peut donc définir le nombre de puissance comme ci-dessous :

\begin{equation}
 N_p= \frac {\iiint\limits_V  \mu \phi_v \, dv} {\rho N^3 D^5}
\end{equation}

pour un fluide newtonien avec

\begin{equation}
 \phi_v = [2 {\tau_{rr}} ^2+ 2 {\tau_{\theta \theta}} ^2 + 2 {\tau_{zz}} ^2 + {\tau_{rz}} ^2 + {\tau_{r\theta}} ^2 + {\tau_{z\theta}} ^2] / {\mu ^2}
\end{equation}

Les contraintes de cizaillement, $ \tau_{rr}, \tau_{r \theta}, \tau_{rz} $ sont données par :

\begin{equation}
\tau_{rr}= -2 \mu \frac {\partial U_r} {\partial r}
\end{equation}

\begin{equation}
\tau_{r \theta}= - \mu [ r \frac {\partial (U_{\theta}/r)} {\partial r} + {\frac {1} {r}} {\frac {\partial U_r} {\partial \theta}}]
\end{equation}

\begin{equation}
\tau_rz = - \mu [ \frac {\partial U_r} {\partial z}+\frac {\partial U_z} {\partial r}]
\end{equation}

 

La deuxième méthode :

La puissance peut aussi être déterminé par le calcul des efforts du fluide sur les pâles :

\begin{equation}
 P= \iint\limits_S  p \vec {n} . \vec {v} \, dS
\end{equation}

avec $ \vec {v} $ la vitesse de la pâle et $ \vec {n} $ est le normal à la paroi.

 

 

Nombre de Froude

Le nombre de Froude défini le rapport des forces d'inerties aux forces de pesanteur.

$$ Fr =\frac {N^2 D} {g}$$

 

 

Ecoulement

Dans cette partie, on va définir les principales grandeurs caractéristiques des écoulements pour compléter l'analyse énergétique du système. On va introduire deux nombres adimensionnels : le nombre de circulation et le nombre de pompage.

 

 

Pompage

Le débit de pompage $ Q_p $ est défini comme le débit de fluide qui traverse le volume balayé par l'agitateur en rotation. Selon le type de pompage, le débit de pompage diffère.

Hélices Turbines Ancres
$ Q_p = 2 \pi \int_{r_s}^{R}  W(r) r  dr $ $ Q_p = \pi d \int_{- \frac {h}{2}}^{ \frac {h} {2}} U(z) dz $ $ Q_p = h \int_{{R} -a}^{R} V(r) dr $

Avec :

$ r_s $ est le rayon d'arbre de rotation
$ R $ est le rayon de l'agitateur
$ h $ est le hauteur de vitesse de la pâle d'agitation
$ T $ est le diamètre intérieur de la cuve
$ a $ est la largeur de pâle d'ancre
 

 

Le nombre de pompage $ N_{Q_P}$ est l'expression adimensionnelle du débit de pompage :

\begin{equation}
N_{Q_P} = \frac {Q_P} {N D^3}
\end{equation}

 

Ce nombre est constant en régime turbulent.

 

 

Le débit de pompage et le débit de circulation pour une hélice. (Source : Cours de Cathérine XUEREB - INP-ENSEEIHT)

Le débit de pompage et le débit de circulation pour une turbine. (Source : Cours de Cathérine XUEREB - INP-ENSEEIHT)

 

 

Circulation

Le débit de pompage est défini comme le débit qui traverse strictement le volume balayé par l'agitateur en rotation. Or, le débit total de fluide mis en mouvement est supérieur à ce débit de pompage. Ce débit est appelé le débit de circulation.

Pour tenir compte de la totalité du débit de fluide induit par l'agitation, il faut d'abord récupérer le ou les nœuds de circulation. Ce sont les points dans la cuve où les composantes axiale et radiale de la vitesse sont nulles.

Selon le type de pompage, l'expression pour calculer le débit de circulation diffère, voici les expressions dans le tableau suivant :

Hélices Turbines Ancres
$ Q_c = 2 \pi \int_{r_s}^{r_c} W(r) r  dr $ $ Q_c = \pi d \int_{- \frac {h}{2}-e }^{ \frac {h} {2}+ e} U(z) dz $ $ Q_c = 2 \pi \int_{r_s}^{r_c} W(r) r dr $

Avec :

$ D $ est le diamètre intérieur de la cuve
$ h $ est le hauteur de vitesse de la pâle d'agitation
$ r_s $ est le rayon de l'arbre de rotation
$ r_c $ est le rayon critique où la vitesse axiale (tangentielle pour le cas d'un ancre) est nulle.
$ e $ est la distance de l'extrémité de la pâle à la hauteur à laquelle le noeuds de circulation se trouve

 

Le nombre de circulation $ N_{Q_c}$ est l'expression adimensionnelle du débit de circulation :

\begin{equation}
N_{Q_c} = \frac {Q_c} {N D^3}
\end{equation}

 

Les débits sont liés l'un à l'autre. Le tableau suivant donne une idée des différences entre nombre de pompage et nombre de circulation pour des différents types de mobiles.

Mobiles  $ N_{Q_p} $  $ N_{Q_c} $
Hélice 0,3 à 1  1,2 à 1,5 $ N_{Q_p} $
Turbine 0,5 à 1,2 (turbine à pales inclinées)
0,45 à 0,9 (turbines à pales concaves)
 1,5 à 2,0 $ N_{Q_p} $

 

OpenFOAM

OpenFOAM (Open Field Operational and Manipulation) est un logiciel de simulation multi-physiques qui résout principalement les équations différentielles de la mécanique des fluides en utilisant la méthode des volumes finis.

OpenFOAM est essentiellement constitué d'une bibliothèque logicielle en langage C++ libre. Le code de OpenFOAM est très intéressant lorsqu'il s'agit d'implémenter de nouveaux modèles. Aucune connaissance approfondie du C++ n'est nécessaire pour écrire son propre modèle sur OpenFOAM.

 

 

Génération ce maillage

OpenFOAM propose deux outils pour générer des maillages : blockMesh et snappyHexMesh. Notre BEI n'a pas pour but de se focaliser sur la manipulation de maillage. Par conséquent nous ne développons pas cette partie dans ce projet.
Pour plus d'information sur l'utilisation de ces deux outils, rendez-vous sur l'URL suivante.

http://hmf.enseeiht.fr/travaux/bei/beiep/content/g21/2-creation-maillage
http://hmf.enseeiht.fr/travaux/projnum/content/utilisation-basique-de-snappyhexmesh
 

Il est important de noter que OpenFOAM est capable de convertir un certain nombre de maillages développé par d'autres logiciels comme DesignModeler, Salome ou CFX. OpenFOAM dispose de plusieurs d'outils spécifiques à chaque mailleur pour convertir les maillages en extension lisible par les solvers OpenFOAM. On peut citer entre autre l'outil fluent3DMeshToFoam pour convertir les .msh et l'outil star4ToFoam pour convertir les maillages Star-CD. La liste des convertisseur est disponible sur l'URL suivante.

http://www.openfoam.org/docs/user/standard-utilities.php

 

 

Calcul avec des solveurs

Comme annonce dans la présentation du BEI, il existe de nombreux solveurs sur OpenFOAM qui couvrent une large gamme de domaines de la physique (combustion, écoulements incompressibles, compressibles, turbulents, multiphasiques, avec ou sans réaction chimique, transferts thermiques etc). Pour les écoulements turbulents, il existe de différents modèles implantés dans le code OpenFOAM comme RANS (k-epsilon, k-omega), LES et DNS. Pour une liste complète, nous vous envoyons une nouvelle fois sur le site OpenFOAM qui regorge d'information utiles.

http://www.openfoam.org/features/standard-solvers.php
http://www.openfoam.org/features/turbulence.php

 

 

Post-traitement

OpenFOAM est distribué avec Paraview, un logiciel de post-traitement open-source. Pour les utilisateurs qui préfèrent d'autres outils de visualisation, il existe des modules d'export pour Fluent, Ensight, FieldView. Là encore, on trouve la liste de ces modules sur le site officiel de OpenFOAM dans la rubrique standard utilities.

http://www.openfoam.org/docs/user/standard-utilities.php

 

Pour ce qui concerne le fonctionnement des répertoires et dictionnaires de OpenFOAM, nous ne détaillerons pas leur principe d'utilisation. En revanche, pour ce qui concerne notre projet, nous expliquerons les choix effectués dans nos démarches. Pour comprendre comment OpenFOAM fonctionne, un excellent tutoriel ce trouve sur l'URL suivante.

https://www.imft.fr/IMG/pdf/INITIATION_OPENFOAM_TOME_1v3-05.pdf

 

 

Intérêt de OpenFOAM

Le code OpenFOAM est un logiciel open-source et il peut être téléchargé gratuitement depuis le site officiel. Comme c'est un logiciel open-source, l'utilisation de ce logiciel est beaucoup plus souple par rapport aux autres logiciels commerciaux et tout le monde a l'accès directement au code pour pouvoir créer son propre solveur.

Voici un example d'équation différentielle :

$\frac{\partial {\rho U}}{\partial t}+\nabla . {\phi U} - \nabla . {\mu \nabla U}=-\nabla p $

Cette équation est écrit comme ci-dessous dans le code OpenFOAM :

solve

(

        fvm::ddt(rho, U)

     + fvm::div(phi, U)

     -  fvm::laplacian(mu, U)

        ==

     -  fvc::grad(p)

);

On peut remarquer que OpenFOAM est logiciel intéressant dans son utilisation pour plusieurs points :

 

 

 

Présentation du problème sur OpenFOAM

Après avoir établi une nomenclature des cuves agitées et développé les méthodes de post traitement, l'objectif est de comparer les différentes stratégies de simulation d'un système agité (cuve agitée) pour un régime turbulent avec le code OpenFOAM. Souvent, on utilise l'une des trois stratégies suivantes. L'approche à référentiel Simple (Simple Reference Frame/SRF), l'approche à référentiels multiples (Multiple Reference Frame/MRF) et l'approche maillage glissant (Moving Mesh).

 

 

Simple Reference Frame

On associe le système agité à un seul référentiel tournant avec la vitesse de rotation du rotor. Dans ce cas, les équations dans le référentiel tournant sont résolues pour tout le domaine avec la force de Coriolis.

 

 

Multiple Reference Frame

Dans cette approche, le domaine est divisé en deux parties. Une partie autour de l'agitateur est simulée dans un repère tournant avec l'agitateur (un référentiel tournant) et les équations sont résolues avec un terme en plus : le terme de Coriolis. Pour le reste du domaine (la deuxième partie), les équations sont résolues en repère fixe.

image

L'approche MRF ne tient pas compte du mouvement relatif  d'une zone mobiles par rapport à des zones adjacentes (qui peuvent être mobiles ou stationnaires), le maillage reste immobile pour le calcul. Cela est analogue à la congélation du mouvement de la partie mobile dans une position spécifique et on observe le champ de l'écoulement instantané autour le rotor à cette position.  On appelle cette aussi cette méthode la méthode du frozen rotor ou agitateur gelé.

Cette approche est approximative mais elle a l'avantage de mettre en jeu des calculs stationnaires (le maillage est immobile) et donc plus rapide.

 

 

Moving Mesh

La méthode de simulation la plus exacte consiste à utiliser le moving mesh. Pour le moving mesh, le maillage est construit de la même façon que pour MRF, mais les équations sont résolues en instationnaire. A chaque pas de temps, le maillage de la partie centrale tourne avec le rotor et l'interface entre les deux zones est recalculé. Pour obtenir des résultats fiables, il ne faut pas que la rotation de l'agitateur(et donc du maillage de la partie centrale) dépasse 5 degrés par pas de temps.

Image : OpenFOAM Project : Different ways to treat rotating geometries par Olivier Petit

Ici, on voit bien que le maillage du référentiel tournant tourne en fonction de temps.

 

 

Les géométries de la cuve

Pour ce projet, TECH-AM Ing met à disposition ses maillages de cuves agitées. Comme annoncé précédemment, les objectifs de se BEI sont de comparer les approches numériques avec le post traitement développé et non pas de concevoir les géométries et maillages.

Voici les géométries fournies par TECH-AM Ing pour le projet.

 

 

Cas Test 3D

Voici la géométrie (coupe intérieure) de la cuve agitée sur laquelle nous allons mettre en place les post traitements. Sur paraview, on peut aller dans la rubrique properties et cocher Patch Names pour afficher les noms de différents patch.

La cuve est constituée de quatre parties (que l'on appelle patch sur OpenFOAM) :

  • la pale qui tourne au centre de la cuve avec une certaine vitesse de rotation (définie par le nom helice dans le maillage)
  • le contour extérieur limitant la cuve (défini par le nom pourtour_externe dans le maillage)
  • le disque haut de la cuve ; on ne le voit pas sur cette image (défini par le nom disque_haut dans le maillage)
  • le disque bas de la cuve (défini par le nom disque_bas dans le maillage).

Il faut bien connaître les différents patches prédéfinis dans le maillage pour qu'on puisse bien définir les conditions aux limites sur OpenFOAM.

 

 

Cas de la géométrie réel :

Voici les géométries avec les patches en configuration MRF.

Voici la géométrie (coupe intérieure) de la cuve agitée.

La cuve est constituée de deux parties principales :

  • La zone_mrf (le cylindre de la taille plus petite que la cuve) contient Wall_fond_zone_mrf, agitateur, agitateur:011 et Wall_arbre_zone_mrf. Wall_fond_zone_mrf, agitateur, agitateur:011 et Wall_arbre_zone_mrf font partie de la pale qui tourne au centre de la cuve avec une certaine vitesse de rotation.
  • La zone_fluide est le reste de la cuve sauf la zone_mrf qui contient quatre composantes de la géométrie : Wall_fond_cuve, wall_arbre_cuve, Wall_virole, et surface_libre (qu'on ne le voit pas sur la géométrie ci-dessus). Wall_fond_cuve, wall_arbre_cuve et Wall_virole font partie de la cuve agitée et surface_libre représente la surface libre du fluide dans la cuve agitée.

Solveurs Choisis

Il existe différents solveurs prêts à être utilisé sur OpenFOAM.  Listons quelques solveurs pour un écoulement incompressible :

  • icoFoam : les écoulements incompressibles laminaires (instationnaires) pour des fluides Newtoniens
  • nonNewtonianIcoFoam : les écoulements incompressibles laminaires (instationnaires) pour des fluides non-Newtoniens/Newtoniens
  • simpleFoam : les écoulements incompressibles turbulents (stationnaires)
  • SFRSimpleFoam : les écoulements incompressibles turbulents (stationnaires) pour des fluides non-Newtoniens/Newtoniens dans un seul référentiel tournant
  • MRFSimpleFoam : les écoulements incompressibles turbulents (stationnaires) pour des fluides non-Newtoniens/Newtoniens en référentiels multiples (souvent un référentiel tournant : le rotor et le fluide environnant, un référentiel immobile : le fluide près de la paroi)
  • pimpleDyMFoam : écoulements incompressibles (instationnaire) pour des fluides Newtoniens en maillage glissant (moving mesh).

( cf  :  http://www.openfoam.org/features/standard-solvers.php pour les solveurs standards sur OpenFOAM.

Pour notre BEI, on utilise trois solveurs : SRFSimpleFoam, MRFSimpleFoam et pimpleDyMFoam. Le MRFSimpleFoam sera le plus utilisé pour les démarches à suivre.

Prise en main sur OpenFOAM

Dans cette partie, nous avons détailler comment sont définis les conditions aux limites et les conditions initiales pour pouvoir faire le calcul numérique à partir du maillage fourni par TECH-AM Ing.

Avant de passer dans un cas de simulation, voici quelques explications sur le fonctionnement du code OpenFOAM. Comme introduit dans les parties précédentes, il existe une multitude de solveurs dans OpenFOAM. On choisi le cas qui nous intéresse  et on le copie sur notre propre espace de travail personnel afin de pouvoir le modifier. Ce système de fichier contient initialement au moins trois dossiers : 0, constant et system.

Image : http://hmf.enseeiht.fr/travaux/bei/beiep/content/g18/structure-generale-...

Le dossier 0 contient toutes les informations sur les grandeurs à l'état initial. Il faut définir les conditions aux limites et initiale pour les différents patches en se basant sur le tutoriel. Pour prendre connaissance des différents types de conditions aux limites, beaucoup d'information se trouvent sur l'URL suivante.

http://www.openfoam.org/docs/user/boundaries.php

 

Le dossier constant contient au moins trois fichiers (transportProperties, turbulenceProperties et RASProperties) dans lesquelles on définit les propriétés du fluide, le modèle de turbulence utilisé ainsi que les paramètres du modèle de turbulence et un sous-dossier qui, lui-même, contient les éléments nécessaires à la description du maillage.

Le dossier system contient au moins trois fichiers que sont controlDict, fvSchemes et fvSolution. Dans le fichier controlDict, on a accès aux paramètres temporels comme le pas de temps et la fréquence d'enregistrement de calcul. Le fichier fvSchemes nous permet de choisir les schémas numériques de discrétisation des équations aux dérivées partielles. Le fichier fvSolution nous permet de préciser les critères de convergences pour les différentes grandeurs. L'URL suivante propose une description plus en détail de l'utilité de ces fichiers.

http://hmf.enseeiht.fr/travaux/bei/beiep/content/g18/description-fichiers

 

Les sous parties suivantes exhibent les configurations que nous avons utilisé pour les deux géométrie dans le cas MRF.

Cas Test 3D

Pour le Cas Test 3D, le modèle k-epsilon est choisi comme modèle de fermeture pour résoudre l'équation de Navier-Stokes moyennée. A noter que le modèle k-Epsilon n'est pas adapté à notre situation puisqu'il est utilisé pour des écoulements parallèles pleinement turbulent (typiquement un écoulement en conduite). Néanmoins, pour disposer de vitesses de convergence rapide il sera choisi. Préférez le Reynolds Stress Model ou le k-omega SST pour l'obtention de résultats plus cohérents.

La suite est destinée à l'explicitation des divers fichiers à paramètrer pour lancer la simulation.

 

 

Turbulent-Stationnaire

0

Puisque nous choisissons de simuler un écoulement turbulent, il est nécessaire de définir les conditions aux limites pour trois grandeurs en plus de la pression et de la vitesse. Il s'agit de la viscosité turbulente, du taux de dissipation visqueux et de l'énergie cinétique turbulente.

 

 

P

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object       p;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /
dimensions             [0 2 -2 0 0 0 0]        // attention, P est normalisé par la densité du fluide

internalField           uniform 0;        // le champ de pression relatif intérieur est nul

boundaryField
{
     disque_haut
     {
          type                      zeroGradient;        //pas de gradient normal pour toutes les grandeurs
      }
     pourtour_externe
     {
          type                      zeroGradient;    
      }
     disque_externe
     {
          type                      zeroGradient;    
      }
     helice                              
     {                                     
          type                      zeroGradient;
      }
}
 

/ / ******************************************************************************** / /

 

 

U

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object       U;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /
dimensions                  [0 1 -1 0 0 0 0]

internalField                uniform (0 0 0);       // le champ de vitesse intérieur est nul

boundaryField
{
     disque_haut
     {
          type                      fixedValue;    
          value                    uniform (0 0 0);
      }
     pourtour_externe
     {
          type                      fixedValue;    
          value                    uniform (0 0 0);
      }
     disque_externe
     {
          type                      fixedValue;    
          value                    uniform (0 0 0);
      }
     helice                        // le vecteur vitesse est nul sur disque_haut, pourtour_externe,
     {                               // disque_bas et helice, on considère la condition de non glissement
          type                      fixedValue;    
          value                    uniform (0 0 0);
      }
}
 

/ / ******************************************************************************** / /

 

 

Les trois grandeurs supplémentaires à configurer.

 

nut

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object       nut;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /
dimensions                   [0 2 -1 0 0 0 0]

internalField                 uniform 0;      

boundaryField
{
     disque_haut
     {
          type                      nutkWallFunction;       //traitement de la viscosité à la paroi
          value                    uniform 0;
      }
     pourtour_externe
     {
          type                      nutkWallFunction;    
          value                    uniform 0;
      }
     disque_externe
     {
          type                      nutkWallFunction;    
          value                    uniform 0;
      }
     helice             
     {                              
          type                      nutkWallFunction;    
          value                    uniform 0;
      }
}
/ / ******************************************************************************** / /

 

 

k

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object       k;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /
dimensions               [0 2 -2 0 0 0 0]

internalField             uniform 1;      

boundaryField
{
     disque_haut
     {
          type                      kqWallFunction;    
          value                    uniform 0;
      }
     pourtour_externe
     {
          type                      kqWallFunction;    
          value                    uniform 0;
      }
     disque_externe
     {
          type                      kqWallFunction;    
          value                    uniform 0;
      }
     helice                  
     {                                   
          type                      kqWallFunction;    
          value                    uniform 0;
      }
}
 

/ / ******************************************************************************** / /

 

 

epsilon

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object       epsilon;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /
dimensions                     [0 2 -3 0 0 0 0]

internalField                   uniform 20;      

boundaryField
{
     disque_haut
     {
          type                      epsilonWallFunction;    
          value                    uniform 0;
      }
     pourtour_externe
     {
          type                      epsilonWallFunction;    
          value                    uniform 0;
      }
     disque_externe
     {
          type                      epsilonWallFunction;    
          value                    uniform 0;
      }
     helice                         
     {                                 
          type                      epsilonWallFunction;    
          value                    uniform 0;
      }
}
 

/ / ******************************************************************************** / /

 

 

Constant

 

RASProperties

Comme expliqué dans la partie précédente : on modifie le fichier RASProperties pour activer le modèle de k-epsilon comme ci-dessous :

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      RASProperties;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /

RASModel         kEpsilon;      // Le modèle de fermeture utilisé est k-epsilon

turbulence          on;                // L'option turbulent est activée

printCoeffs          on;              // Affichage des coefficients de turbulence dans le shell
 

/ / ******************************************************************************** / /

Pour plus d'information sur les modèles de turbulence disponibles dans les librairies OpenFOAM, voir le site à l'URL suivante : http://www.openfoam.org/features/RAS.php

 

 

transportProperties

C'est dans ce fichier qu'est paramétré les propriétés du fluide présent dans la cuve, pour nous il sera newtonien avec une viscosité cinémétique de 1e-5 m²/s.

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version       2.0;
    format         ascii;
    class          dictionary;
    location     "constant";
    object         transportProperties;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /

transportModel     Newtonian;

nu                           nu [ 0 2 -1 0 0 0 0  ] 1e-5;  // La viscosité cinématique de fluide ( $ m^2 s^{-1} $)

 

/ / ******************************************************************************** / /

 

 

MRFZones

C'est dans ce fichier qu'est spécifiée la zone tournante MRF. La zone tournante est, comme son nom l'indique une zone. Elle doit être définie dans /constant/polyMesh/cellZones.

 

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version       2.0;
    format         ascii;
    class          dictionary;
    location     "constant";
    object         MRFZones;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /

1
(
     zone_mrf
     {
          
/ / Fixed patches (by default they 'move' with the MRF zone)

          nonRotatingPatches ();
          origin      origin      [0 1 0 0 0 0 0 ]    (0 0 0);     // Le centre de la cuve tournante
          axis         axis        [0 0 0 0 0 0 0 ]    (0 0 1);     // L'axe de rotation : ici, la rotation se fait selon l'axe  Z
          omega    omega   [0 0 -1 0 0 0 0 ]   40;           // la vitesse de rotation de la pale en radians par seconde
!!
     }

)
/ / ******************************************************************************** / /

 

 

 

System

 

fvSchemes

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version       2.0;
    format         ascii;
    class          dictionary;
    location     "system";
    object         fvSchemes;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /

ddtSchemes
{
    default         steadyState;        // pour le régime stationnaire
}

gradSchemes                             // les schémas de discrétisation des différents termes dans les équations
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss limitedLinearV 1;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,epsilon) Gauss limitedLinear 1;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         none;
    laplacian(nuEff,U) Gauss linear corrected;
    laplacian((1|A(U)),p) Gauss linear corrected;
    laplacian(DkEff,k) Gauss linear corrected;
    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
    interpolate(U)  linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p               ;
}

/ / ******************************************************************************** / /

 

 

fvSolution

/*--------------------------------*- C++ -*--------------------------------------------------------*\
| ======                        |                                                                                       
| \\         /  F ield             | OpenFOAM:   The Open Source CFD Toolbox     
|  \\      /   O peration      | Version:          2.1.0                                                   
|   \\   /    A nd                 | Web:                www.OpenFOAM.org                       
|    \\/     M anipulation  |                                                                                       
\*----------------------------------------------------------------------------------------------------*/
FoamFile
{
    version       2.0;
    format         ascii;
    class          dictionary;
    location     "system";
    object         fvSolution;

/ / * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * / /

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.001;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         5;
        tolerance       1e-07;                      //valeur sous laquelle le solveur doit descendre dans la convergence
                                                                //
des calculs   res < tolerance
        relTol          0.001;                        //ou alors,  relTol < residu_initial
    }                                                         // la première condition atteinte stoppe la simulation

    k
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         10;
        tolerance       1e-07;
        relTol          0.001;
    }

    epsilon
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         10;
        tolerance       1e-07;
        relTol          0.001;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 5;     // Quand les centres de cellules sont pas alignés, cela permet de corriger les
    pRefCell        0;                              // gradients
    pRefValue       0;
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.3;
        k               0.1;
        epsilon         0.1;
    }
}

/ / ******************************************************************************** / /

 

 

Turbulent-Instationnaire

Pour faire la même simulation mais instationnaire, il faut modifer le schéma de discrétisation du terme dt dans le fichier fvSchemes :

ddtSchemes
{
    default         Euler;
}

 

 

Laminaire-Stationnaire

Pour faire une simulation pour un très petit nombre de Reynolds (on est dans un régime laminaire), il faut modifier le fichier RASProperties comme ci-dessous :

// RASModel         kEpsilon;  
RASModel           laminar     

turbulence          off;                // L'option turbulent est déactivée car on est dans le régime laminaire

printCoeffs          off;             

 

Pour les conditions aux limites et conditions initiales, il suffit juste de définir les deux paramètres de $ P $ et $ U $ .

 

 

 

Cas réel 3D (une hélice)

Pour le cas Réel 3D, les paramètres sont exactement les mêmes que pour le cas Test 3D. Seul change la condition sur la partie supérieure.

  • agitateur
  • wall_virole
  • wall_fond_cuve
  • wall_fond_mrf
  • wall_arbre_zone_mrf
  • wall_arbre_cuve
  • agitateur:011
  • surface_libre

Tous les patches dans le cas Test 3D (disque_haut, pourtour_externe, disque_bas et helice) sont les murs. Reprenons notre cas Réel 3D : tous les patches sauf surface_libre sont les murs et on peut faire pareil que la partie précédente pour définir les conditions aux limites et conditions initiales dans les fichier p , U, epsilon, k et nut. Par contre, on définit le patch de surface_libre dans les fichiers p , U, epsilon, k et nut comme suivant :

    surface_libre
    {
        type            symmetryPlane;
    }

 

Pour pouvoir lancer le cas Réel 3D, on reprend les démarches expliquées dans le cas Test 3D.

Post-traitement des données

OpenFOAM

Cette rubrique a pour but d'illustrer les post-traitements que l'on a utilisé à l'aide des outils OpenFOAM. En effet, OpenFOAM dispose d'un certain nombre d'outils, majoritairement regroupé sous le nom de "Post-processing - Standard Utilities", qui permet d'analyser l'écoulement simulé.

 

foamCalc

foamCalc est un outil simple qui permet de calculer des grandeurs nécessaires au post-traitement. Pour la détermination de la vitesse en norme, on pourra utiliser la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation/system$    foamCalc mag U

Pour l'examination des composantes du tenseur des contraintes, on pourra utiliser la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation/system$    foamCalc sigmaxx

Pour connaître toutes les possibilités de l'utility foamCalc, on peut utiliser la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation/system$    foamCalc xxx

 

 

Sample

L'utility "sample" permet d'enregistrer des données de calcul le long de lignes pour tracer des graphes ou sur des surfaces pour l'affichage  d'isograndeurs en deux dimensions.
Pour bénéficier de cet outil, il suffit de spécifier les éléments nécessaires dans un fichier nommé "sampleDict" à placer dans le dossier "system".
Comme pour la plupart des fichiers OpenFOAM, il est intéressant d'aller chercher un exemplaire "sampleDict" dans les tutoriels proposés par OpenFOAM. Pour éviter de perdre du temps à chercher ce "dictionary", la commande suivante permet de recenser tous les fichiers du même nom, c'est une commande de recherche.

rmartin@tanche : /mnt/hmf/OpenFOAM/OpenFOAM-2.1.0$    find -name sampleDict

 

Une fois le fichier placé dans le répertoire /work/MaSimulation/system, ouvrir le fichier avec gedit.

rmartin@tanche : /work/rmartin/MaSimulation/system$    gedit sampleDict &

 

Le "SampeDict" ci-dessous permet de créer un répertoire /MaSimulation/sets contenant la grandeur spécifiée, ici Umag, sur les 3 lignes line_R=0.*m.

 

/*--------------------------------*- C++ -*---------------------------------------------------*\
| =====                        |                                               
| \\      /  F ield               | OpenFOAM: The Open Source CFD Toolbox          
|  \\    /   O peration       | Version:  2.1.0                                
|   \\  /    A nd                | Web:      www.OpenFOAM.org                    
|    \\/     M anipulation    |                                               
\*----------------------------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      sampleDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

interpolationScheme       cellPoint;              //choix d'interpolation des grandeurs demandées
                                                                // cellPoint : interpolation linéaire basée sur la taille des cellules

setFormat                     raw;                      //choix du format d'écriture 1D : raw - format ASCII en colonne                                                                           //exploitable sur matab. surfaceFormat pour une écriture 2D
sets                                                         //nom du répertoire dans lequel sera écrit les données
(

    line_R=0.2m                                          //nom d'écriture des fichiers
    {
        type         uniform;                             //répartition uniforme des 146 points sur la ligne
        axis          z;                                      //coordonnées de référence
        start         ( 0 0.2 -0.25 );                   //coordonnées du premier et du dernier point de la ligne
        end          ( 0 0.2 0.75 );
        nPoints     146;                                  //nombre de points sur la ligne
    }

    line_R=0.3m
    {
        type         uniform;
        axis          x;
        start         ( 0 0.3 -0.25 );
        end          ( 0 0.3 0.75 );
        nPoints     134;
    }

    line_R=0.4m
    {
        type         uniform;
        axis          x;
        start         ( 0 0.4 -0.25 );
        end          ( 0 0.4 0.75 );
        nPoints     135;
    }
);

fields          (magU);                               //liste des grandeurs à enregistrer sur les lignes

// ******************************************************************************* //

 

 

L'image ci-dessus illustre bien la position de ces trois lignes dans la cuve agitée.

 

Dans le cas de l'hélice en cuve réelle, nous tracerons les données selon quatre lignes horizontales à différentes hauteurs selon l'image suivante.

 

 

 

stressComponents

Comme nous avons expliqué dans la partie Caractérisation des Mobiles, nous pouvons calculer la puissance dissipée par la formule ci-dessous :

\begin{equation}
 P=\iiint\limits_V  \mu \phi_v \, dv
\end{equation}

Dans les coordonées cartésiennes, nous pouvons impliciter le terme de $ \phi_v $ comme suivant :

\begin{equation}
 \phi_v = [2 {\tau_{xx}} ^2+ 2 {\tau_{yy}} ^2 + 2 {\tau_{zz}} ^2 + {\tau_{xz}} ^2 + {\tau_{yz}} ^2 + {\tau_{xy}} ^2] / {\mu ^2}
\end{equation}

Avec :

\begin{equation}
{\tau_{xx}}= - 2 {\mu} \frac { \partial U} {\partial x}
\end{equation}

\begin{equation}
{\tau_{yy}}= - 2 {\mu} \frac { \partial V} {\partial y}
\end{equation}

\begin{equation}
{\tau_{zz}}= - 2 {\mu} \frac { \partial W} {\partial z}
\end{equation}

\begin{equation}
{\tau_{xy}}= - {\mu} (\frac { \partial U} {\partial y}+\frac { \partial V} {\partial x})
\end{equation}

\begin{equation}
{\tau_{xz}}= - {\mu} (\frac { \partial W} {\partial x}+\frac { \partial U} {\partial z})
\end{equation}

\begin{equation}
{\tau_{yz}}= - {\mu} (\frac { \partial V} {\partial z}+\frac { \partial W} {\partial y})
\end{equation}

Pour pouvoir calculer le puissance, il faut demander à OpenFOAM à sortir les valeurs des contraintes. Il existe une fonction auxiliaire sur l'OpenFOAM : stressComponents qui permet de sortir toutes les valeurs de $ \sigma_{ij} $.

L'équation qui relie $ \sigma_{ij} $ et $ \tau_{ij} $ est comme suit :

\begin{equation}
\sigma_{ij} = -p I_{ij}+\lambda (div U_i) I_{ij} +2 \mu D_{ij}
\end{equation}

avec :

$ I_{ij} $ est la matrice d'indentité

$ div (U_i) = 0 $ car c'est incompressible

$ D_{ij} = \frac {1} {2} (grad U_i + grad^t U_i) $

Pour sortir les composantes du tenseur sigma, lancez la commande suivante :

rmartin@tanche : /work/rmartin/MaSimulation$    stressComponents

 

Afin de mieux comprendre comment OpenFoam calcule $ \sigma_{ij} $ à partir des gradients de vitesse, nous allons sur la librairie d'OpenFOAM :

/mnt/hmf/OpenFOAM/OpenFOAM-2.1.0/applications/utilities/postProcessing/stressField/stressComponents/

 

Essayons de comprendre le script de stressComponents.C :

volSymmTensorField sigma
            (
                IOobject
                (
                    "sigma",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                laminarTransport.nu()*2*dev(symm(fvc::grad(U)))   // sigma est le produit de $ \nu $ et des gradients de U. OpenFoam
            );                                                                               

L'unité de la puissance est le Watt et sa dimension est $ [M{L^2} {T^{-3}}] $. La dimension de la viscosité dynamique est $ [M L^{-1} T^{-1}] $ et la dimension de volume est $ [L^3] $. A l'aide de l'analyse dimensionnelle, on peut déduire que la dimension de $ \phi_v $ est $ [T^{-2}] $. 

Nous pouvons donc déduire qu'il suffit de normaliser le terme de $ \sigma_{ij} $ obtenu avec OpenFOAM par $ \nu $ pour avoir le terme de $ \tau_{ij} $.

Sur paraview, afin d'obtenir le terme de $ \phi_v $, il faut normaliser le terme de $ \sigma_{ij} $ calculé par OpenFOAM comme suit :

$$ \phi_v= [2 \sigma_{xx}^2 +2 \sigma_{yy}^2 + 2 \sigma_{zz}^2 + \sigma_{xy}^2 +\sigma_{xz}^2 + \sigma_{yz}^2 ]/ {\nu^2}  $$

 

 

pyFoamPlotWatcher

pyFoamPlotWatcher.py est un script python très utile puisqu'il peut afficher les résidus et la continuité d'une simulation en cours avec gnuplot.

Pour une visualisation des résidus, il faut rentrer les deux commandes suivantes.

rmartin@tanche : /work/rmartin/MaSimulation/system$ MRFSimpleFoam > log &
rmartin@tanche : /work/rmartin/MaSimulation/system$ pyFoamPlotWatcher.py log &

 

Voici ce que l'on peut obtenir :

Matlab

L'ensemble des données extraites par les utilities OpenFOAM ont été traité sur Matlab. Nous allons détailler en quelques lignes le script Matlab utilisé.

Nous avons précisé une extension en .xy pour toutes les données que le sampleDict sortirait grâce à l'option setformat : raw. L'exploitation sous Matlbab se fait sans soucis.

L'exploitation des données se fait en deux étapes :

  • l'importation sous Matlab des données
  • le tracé des graphes correctement renseignés

 

Importation des données

L'importation des données se fait avec l'outil "importdata

F02_SIMPLE = importdata('500_Umag_0.2m_Fluent_SIMPLE.xy');
OF02 = importdata('line_R=0.2m_magU.xy')

 

 

Tracé des graphes

Afin de disposer de graphes lisibles sur la majorité des supports, optez pour les commandes set, fontsize et linewidth pour grossir et passer un gras les caractères.

figure(4);
hold on;
grid on;
plot(F02_SIMPLE(:,1),F02_SIMPLE(:,2), 'g','linewidth',3);
plot(F02_SIMPLE(:,1),OF02(:,2), 'b','linewidth',3);
xlabel('Height Z  m ', 'fontsize',22);
ylabel('Velocity Mag  m/s', 'fontsize',22);
set(gca, 'fontsize',22);
grid on;
title('Velocity Mag for R=0,2m at it=500', 'fontsize',22);
legend('Fluent 13.0 - SIMPLE', 'OpenFoam 2.1.0 - SIMPLE');
set(legend, 'fontsize',22);

 

Le graphe associé :

 

Paraview

Une partie du post-traitement des données calculées par OpenFOAM peut se faire à l'aide de paraview, notamment pour la détermination des nombres adimensionnels caractéristiques.

En effet, d'après la définition des nombres adimensionnels de débit, l'utilisation de paraview peut se révéler utile pour la création de surfaces à travers desquelles peuvent être déterminés des débits.

Cette rubrique a pour but d'illustrer la méthode employée pour la détermination grandeurs globales.

 

Méthode à suivre :

Commencer par ouvrir un terminal et aller dans le répertoire ou se trouve l'ensemble des fichiers de la simulation OpenFOAM :

rmartin@tanche :   cd /work/rmartin/MaSimulation
rmartin@tanche : /work/rmartin/MaSimulation$   ls
0   10   20   30   40   constant   system   log

Pour pouvoir lancer les commandes OpenFOAM depuis n'importe quel répertoire :

rmartin@tanche : /work/rmartin/MaSimulation$    source /mnt/hmf/OpenFoam/OF21.sh

Pour bénéficier de tous les outils de Paraview, on préférera lancer Paraview avec les deux commandes suivantes plutôt que paraFoam directement :

rmartin@tanche : /work/rmartin/MaSimulation$   touch .OpenFoam
rmartin@tanche : /work/rmartin/MaSimulation$   /usr/local/bin/paraview

Pour pouvoir exploiter les données OpenFOAM avec paraview, il faut convertir les données dans un format que peut lire Paraview. Le format .vtk est adapté. Pour convertir les données en .vtk il suffit de lancer la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation$   foamToVTK

Ouverture de paraview avec la commande suivante.

rmartin@tanche : /work/rmartin/MaSimulation$   paraview

L'ouverture des fichiers .vtk se fait manuellement dans un premier temps : => Files => Open => *.vtk

Pour la suite des commandes, on divisera les cas en fonction du type de pompage de l'agitateur : axial ou radial.

 

 

Calcul du débit de pompage

 

Cas des mobiles d'agitation à pompage axial :

Afin de déterminer le débit de pompage d'un agitateur à pompage axial, il faut créer une surface spécifique qui servira de support à l'intégration de la vitesse axiale. La détermination de ce débit de pompage se divise en 3 parties.

Voici le cas d'une cuve à agitation axiale (une hélice qui tourne selon l'axe-y) qui servira d'illustration :

 

Création d'un cylindre avec le Programmable Filter (attention, cette option n'est pas disponible dans paraFoam)

L'outil Programmable Filter est un outil puissant puisqu'il permet de créer un filtre coder en python correspondant à la fonction recherchée. Nous allons créer un Clip Cylinder, filtre qui n'existe pas dans la version utilisée, à savoir la 12.6.0 :

Pour bénéficier de l'outil Programmable Filter : => Filters => Alphabetical => Programmable Filter

Une fenêtre permet de rentrer l'ensemble des lignes de code nécessaire. Pour ce qui concerne la création d'un cylindre par un Clip, nous avons rentré l'ensemble des lignes de code suivante :

input = self.GetInputDataObject(0, 0)
inp_copy = input.NewInstance()
inp_copy.ShallowCopy(input)
inp_copy.UnRegister(None)

transf = vtk.vtkTransform()
transf.RotateY(90)

cyl = vtk.vtkCylinder()
cyl.SetCenter(0,0,0)
cyl.SetRadius(1)
cyl.SetTransform(transf)

clipper = vtk.vtkClipDataSet()
clipper.SetClipFunction(cyl)
clipper.SetInput(inp_copy)
clipper.InsideOutOn()
clipper.Update()

self.GetOutputDataObject(0).ShallowCopy(clipper.GetOutputDataObject(0))

Ce programme permet de créer un cylindre  de dimension infini selon Y grâce à une fonction implicite.

Le paramètrage des dimensions du cylindre dans les deux autres directions se fait via les deux commandes suivantes :

cyl.SetCenter(0,0,0) permet de fixer le centre du cylindre
cyl.SetRadius(1) permet de choisir le rayon du cylindre

On obtient la visualisation suivante :

C'est bien un cylindre de rayon de 0.95 m ayant un diamètre plus petit que celui de la cuve.

 

Création du slice disque :

Le Clip_Cylindre précédent a permis de délimiter la zone dans les plans normaux à Y. La création de la surface en forme de disque se fait à l'aide du filtre Slice :

=> Filters => Alphabetical => Slice

On choisit de créer un Slice_plane qui, délimité par le Clip_cylindre, revient à la création d'un disque.

Il faut bien choisir la bonne normale qui est selon l'axe Y pour créer un disque.

On obtient la visualisation suivante dans le cas d'un slice fait trop prêt de l'hélice :

 

Calcul de l'intégrale de la vitesse verticale sur le disque avec Integrate Variables

On peut utiliser la fonction IntegrateVariables sur paraview pour calculer le débit volumique de fluide qui passe à travers le disque défini dans la partie précédent.

=> Filters => Alphabetical => Integrate Variables

La fonction Integrate Variables est capable de calculer la grandeur souhaitée sur une surface ou un volume prédéfini. Attention à bien appliqué le filtre IntegrateVariables sur le Slice : ceci se fait en cochant le Slice défini dans la colonne de droite avant d'utiliser le filtre.

Le calcul par IntegrateVariables montre que le débit de fluide qui passe par le disque selon l'axe Y est d'environ  0,21 $ m^3 /s $.

 

 

Cas des mobiles d'agitation à pompage radial :

 

Afin de déterminer le débit de pompage d'un agitateur à pompage radial qui tourne selon l'axe-Z, il faut créer une surface d'un cylindre de la même hauteur de l'agitateur qui servira de support à l'intégration de la vitesse radiale.

La géométrie du cas test 3D sera utilisée pour l'illustration de la procédure à suivre dans le cas des agitateurs à pompage radial.

 

Calcul de vitesse radiale avec calculator

Le logiciel de l'OpenFOAM travaille dans un système de coordonnées cartésiennes. Afin de calculer le débit de pompage à travers la surface d'un cylindre, il nous faut tout d'abord la vitesse radiale à l'aide de l'outil calculator sur paraview.

=> Filters => Alphabetical => Calculator

On a des équations comme ci_dessous qui nous permettent de passer du système de coordonées cartesiennes au système de coordonées cylindriques :

\begin{equation}
U_r=U_x cos (\theta)+U_y sin(\theta)
\end{equation}

\begin{equation}
U_\theta=- U_x sin (\theta) +U_y cos(\theta)
\end{equation}

avec :

$$ \theta= arctg (\frac {coordY} {coordX} ) $$

Attention à la fonction non définie pour coord-x $ =0 $ car la valeur de $ \theta $ tend vers l'infini. Pour éviter ce problème, on peut définir les équations de $ U_r $ et $ U_{\theta} $ autrement :

\begin{equation}
U_r=U_x \frac {coordX} {\sqrt {{coordX}^2+{coordY}^2}} +U_y \frac {coordY} {\sqrt {{coordX}^2+{coordY}^2}}
\end{equation}

\begin{equation}
U_\theta=- U_x \frac {coordY} {\sqrt {{coordX}^2+{coordY}^2}} +U_y \frac {coordX} {\sqrt {{coordX}^2+{coordY}^2}}
\end{equation}

On peut modifier la nouvelle variable qu'on définit par la fonction calculatrice comme U_r (la vitesse radiale).

 

Création d'un cylindre avec le Programmable Filter (attention, cette option n'est pas disponible dans paraFoam)

input = self.GetInputDataObject(0,0)

inp_copy = input.NewInstance()

inp_copy.ShallowCopy(input)

inp_copy.UnRegister(None)

cutter = vtk.vtkCutter()

transf = vtk.vtkTransform()

transf.RotateX(90)

cyl = vtk.vtkCylinder()

cyl.SetCenter(0,0,0)

cyl.SetRadius(0.26)

cyl.SetTransform(transf)

cutter.SetCutFunction(cyl)

cutter.SetInput(inp_copy)

cutter.Update()

self.GetOutputDataObject(0).ShallowCopy(cutter.GetOutputDataObject(0))

Ce script de python ne marche que sur Paraview 6.14.

On obtient la visualistion suivante :

 

Création d'une boîte d'une hauteur 0.5m contenant l'hélice

On se rappelle qu'on veut faire l'intégral de la vitesse radiale sur la surface transversale d'un cylindre ayant la même hauteur que les pâles d'agitateur. Il faut donc enlever la partie haute et la partie basse de la cuve où il n'y pas d'agitateur à l'aide de la fonction clip sur paraview:

=> Filters => Alphabetical => Clip

On choisit le clip type : box et il faut absolument cocher Inside Out pour garder la partie contenant l'hélice.

On obtient donc la visualisation suivante :

 

Calcul de l'intégrale de la vitesse radiale sur la surface transversale du cylindre avec Integrate Variables

On peut utiliser la fonction Integrate Variables sur paraview pour calculer le débit volumique de fluide qui passe la surface transversale du cylindre.

=> Filters => Alphabetical => Integrate Variables

 

 

Calcul de la puissance pour des mobiles à pompage axial ou radial

 

Comme nous avons expliqué dans la partie OpenFOAM, il faut intégrer le champ de sigma (les contraintes de cisaillement) ou $ \phi_v $ avec la viscosité dynamique sur tout le volume de la cuve.

 

 

Calcul de $ \mu \phi_v $ à l'aide de calculator

Afin de pouvoir utiliser la fonction Integrate Variables pour calculer l'intégral sur tout le volume, il faut calculer abord le produit du terme de $ \phi_v $ normalisé avec la viscosité dynamique à l'aide de la calculator. Il faut rentrer la formule ci-dessous sur la fonction de calculator :

$$ \mu [2 \sigma_{xx}^2 +2 \sigma_{yy}^2 + 2 \sigma_{zz}^2 + \sigma_{xy}^2 +\sigma_{xz}^2 + \sigma_{yz}^2 ]/ {\nu^2}  $$

Avec : $ \mu=0,04 SI $ et $ \nu=2,67.10^{-5} SI $

L'image illustre la formule de $ \mu \phi_v $ entrée sur la fonction calculator sur Paraview.

 

Calcul de l'intégrale de $ \mu \ phi_v$ sur tout le volume avec Integrate Variables

Une fois le produit de $ \mu \phi_v $ calculé, nous utilisons la fonction Integrate Variables pour calculer l'intégrale de $ \mu \phi_v $ sur tout le volume de la cuve.

=> Filters => Alphabetical => Integrate Variables

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

Tutorial - Moving Mesh

Ce petit tutorial a pour but de convertir les maillages du fluent sous un format d'OpenFOAM et aussi de les fusionner afin de pouvoir faire Moving Mesh :

interface_zone_mrf
    {
        type            patch;
        nFaces          10749;
        startFace       631265;
    }
    interface_zone_cuve
    {
        type            patch;
        nFaces          10749;
        startFace       663905;
    }

   

PAR :
   

AMI1
    {
        type            cyclicAMI;
        nFaces          10749;
        startFace       631265;
        matchTolerance  0.0001;
        neighbourPatch  AMI2;
        transform       noOrdering;
    }
    AMI2
    {
        type            cyclicAMI;
        nFaces          10749;
        startFace       663905;
        matchTolerance  0.0001;
        neighbourPatch  AMI1;
        transform       noOrdering;
    }

 

 

Conclusion

Le sujet sur les cuves agitées que nous avons choisi a été proposé par TECH-AM Ing, une entreprise spécialisée dans les études CFD.  Le sujet portait sur la simulation numérique des cuves agitées avec OpenFOAM. Les cuves agitées sont omniprésentes dans les procédés. Cependant, du fait de leurs dimensions et de leur construction il est difficile d'optimiser leur conception. Le dimensionnement empirique et/ou la simulation numérique sont fréquemment utilisés pour la conception de ces agitateurs. Le choix du logiciel OpenFOAM est légitime pour tous les avantages cités dans les parties "Présentation du sujet de BEI" et "OpenFOAM".

Le sujet sur lequel nous avons travaillé depuis début Novembre 2013 était assez ambitieux dans son contenu puisqu'il y avait beaucoup de tâches à remplir. Voici les points principaux :

 

Les quatre premiers points ont été globalement bien abordés.

Malgré le peu de temps à disposition, la nomenclature que nous avons créée a été enrichie par une recherche bibliographique assez longue du fait du peu de sources disponibles et récentes. La caractérisation des systèmes d'agitation à l'aide de nombres adimensionnels pertinents a découlé de la recherche bibliographique. Nous nous sommes attachés à bien renseigner les plages de valeur des nombres adimensionnels correspondant aux domaines de fonctionnement de chaque type de mobile d'agitation. Ce travail était essentiel pour la suite des événements.

La phase suivante du projet a été la prise en main de OpenFOAM puisque nous n'avions jamais utilisé ce logiciel. Il faut avoué que son utilisation n'a rien d'intuitif au premier abord puisqu'il n'y a pas d'interface graphique. Fort heureusement, avec le temps, nous avons réussi à comprendre le fonctionnement et sa logique qui en font un logiciel à la mode.

Le post-traitement a été fait à l'aide de trois logiciels :

Le développement des méthodes de post-traitement nous a pris énormément de temps pour plusieurs raisons. La prise en main de OpenFOAM n'a pas été des plus rapides. Les étapes de post-traitement ont été longues à développer car nous avons dû faire face à plusieurs problèmes. Premièrement, la stabilité de paraview qui n'était pas bonne pour l'utilisation qu'on en faisait (crash à répétition, messages d'erreur au lancement mais pas sur toutes les sessions utilisateurs...). Deuxièmement, une des étapes sous paraview consistait en la mise en place d'un code en python ; code que n'ont n'avions jamais utilisé, beaucoup de temps a été consacré à cette étape qu'était la création d'un slice_cylinder et la création d'un clip_cylinder. Troisièmement, la recherche de toutes les commandes nécessaires sous OpenFOAM, paraview et ScriptPython pour mener à bien le projet. En prime, des problèmes de session utilisateurs nous ont empêché d'aborder dans les temps le dernier point cité précédemment et malheureusement l'un des plus importants.

Finalement, nous avons réussi à surmonter un nombre important de difficultés ce qui nous a permis de tout de même mettre en place correctement les méthodes de post-traitement pour les deux types d'agitateurs les plus utilisées dans l'industrie, à savoir les agitateurs à écoulement axial et les agitateur à écoulement radial.

Bibliographie

 [1] Ecoulement pour les procédés - Mathieu MORY

 [2] Agitation et mélange - Cathérine Xuereb - Edition DUNOD

 [3] Thèse - Agitation mécanique des fluides visqueux, contribution à l'étude des éoulements tridimensionnels - 1989 - Mohamed Salah ABID

 [4] Thèse - Approche numérique tridimensionnelle de l'agitation mécanique en régime turbulent - 1992 - Jean MOUREH

 [5] Thèse - Hydrodynamique d'écoulement en cuve mécaniquement agitée - 2009 - Hamid LAKHDARIL

 [6] Thèse - Agitation de fluides visqueux - cas de mobile à pales, d'ancres et de barrières -  Joël BERTRAND

 [7] Thèse - Contribution à l'étude de l'agitation d'un fluide pseudo-plastique dans une cuve avec un agitateur à ancre -1980 - Chairit SATAYAPRASERT

 [8] Thèse - Agitation de fluides visqueux pseudoplastiques par un doubles ruban héicoïdal - 1985 - Maher Gamal SOLIMAN - Jean Pierre COUDERC

 [9] http://www.cfd-online.com/Forums/openfoam-paraview/95729-cut-clip-cylinder-instead-sphere-plane.html - création d'un cylindre avec le programmable filter à l'aide de clipper

 [10] http://www.cfd-online.com/Forums/openfoam-paraview/117815-paraview-slice-type-cylinder.html - création d'une surface transversale du cylindre avec le programmable filter à l'aide de slice/cutter

 [11] OpenFOAM - The Open Source CFD Toolbox Programmer's Guide - Version 2.1.0 - 15th December 2011

 [12] Paraview User's Guide (v3.12)

 [13] http://www.openfoam.org/docs/user/index.php - The OpenFOAM Foundation

 [14] http://paraview.org/OnlineHelpCurrent/

 [15] ANSYS, Inc. Novembre 2010. ANSYS FLUENT User's Guide

 [16] ModéLisation numérique des écoulements gênérés dans une cuve mécaniquement agitée par une turbine de Rushton - 20ème Congrès Français de Mécanique 29 août au 2 septembre 2011 - H. Ameur, M. Bouzit, M. Helmaoui, F. Bouzit