Modélisation

Le but du projet va être de reproduire de façon simplifiée la maquette du Laboratoire de Génie Chimique de Toulouse sur un logiciel de simulation. Pour cela, il faut réaliser un maillage, définir des modèles et imposer des conditions aux limites.

L'ensemble de cette partie aura pour but de présenter ces travaux, de détailler ce qui a été récupéré de précédents projets et ce qui a été développé au cours de ce BEI.

La quasi totalité des simulations ont été réalisées à distance en utilisant le supercalculateur du calmip - groupement de 17 laboratoires de Toulouse et de la région Midi-Pyrénées qui a pour but d'encourager le développement de nouvelles techniques de calcul au sein de la communauté scientifique -, HYPERION.

HYPERION nous aura permis de lancer un calcul employant 128 processeurs durant près de 48 heures. Au vue de la complexité du maillage et du problème, ce fut la seule façon d'arriver à un résultat. Voici le supercalculateur :

NEPTUNE_CFD

La simulation est effectuée avec le logiciel NEPTUNE_CFD V1.08@Toulouse.

Cette branche de NEPTUNE_CFD developpée à l'IMFT est destinée à la simulation des écoulements gaz-particules, tout particuliérement les lits fluidisés.

NEPTUNE_CFD est un code de calcul développé par le consortium EDF (Electricité de France) - CEA (Commissariat à l'Energie Atomique) - AREVA - IRSN (Institut de Radioprotection et de Sûreté Nucléaire) dédié aux écoulements multiphasiques réactifs 3D turbulents autant à l'échelle locale que dans des dimensions industrielles.

Le codage se réalise en Fortran 77 pour la majorité du code et le portage se fait par Linux principalement.

 

NEPTUNE_CFD s'associe par défaut au mailleur Simail 7.0.4. et travaille en pseudo-2D, c'est-à-dire en 3D sur une épaisseur d'une maille.

C'est un logiciel intéractif graphique de génération et de manipulation des maillages linéiques, surfaciques et volumiques. Ses algorithmes performants le rendent adapté au traitement des géométries les plus complexes.

 

Puis le logiciel Edamox permet de définir les paramètres de calcul telles que les conditions au limites, les propriétés de chaque phases, les paramètres de calculs, les modèles physiques adoptés.

Récupération des travaux précédents

Dans le cadre de l'étude de l'hydrodynamique de l'écoulement, une maquette de simulation a déjà été mise en place par Hervé Neau de l'IMFT - Institut de Mécanique des Fluides de Toulouse - et Hadrien Benoit, au cours d'un précédent BEI (année scolaire 2011 - 2012). Les différents maillages ont été réalisés par Hervé Neau et la première simulation par Hadrien Benoit.

Hadrien Benoit a également mis en place un système de codage de la régulation en pression du dispositif via un fichier utilisateur uskpdc.F.

L'ensemble de ce travail va être repris et adapté. Nous allons détailler dans les deux sous-rubriques les caractéristiques du travail repris et ce que nous avons modifié.

 

Types de maillages

Le maillage très complexe que l'on voit ci-dessous à était réalisé sous Simail, logiciel qui a permis le recollement, la translation, etc des différents maillages.

On peut voir ici le caisson avec une alimentation en solide latérale et une vanne de régulation au dessus. Le tube est plongé dans le caisson, il a une hauteur de 2 m et il est placé à 10 cm de la base du bac.

 

Les simulations se font avec deux types de maillage :

  • "maillage grossier" de 105 128 cellules qui permet la prise en main rapide des calculs (moins d'une heure) afin de déterminer les bons paramètres de la simulation
  • "maillage fin" de 1 649 044 cellules (très raffiné) qui prend donc beaucoup plus de temps (72 h) mais qui permet d'obtenir des résultats précis

Ci-dessous, nous pouvons observer de plus près le maillage fin au niveau du tube et l'emplacement de l'aération. Chaque cellule est de type hexahédrique, avec une hauteur de 1.5 mm et un largeur de 1.2 mm.

 

Contrôle de la pression dans le dispositif

Le contrôle de la pression dans le dispositif est quelque chose de primordial. En effet, ce sont les forces motrices de pression qui vont permettre aux particules du lit fluidisé de monter dans le tube et ainsi de faire fonctionner le dispositif.

Le système doit donc imposer une pression dans le ciel du lit - entre le haut des particules dans le caisson et le haut du caisson - pour créer une perte de charge suffisante entre le ciel du lit et le sommet du tube qui est à pression atmosphérique.

Cette pression imposée est contrôlée par un dispositif vanne - capteur. L'image suivante représente un schéma très simplifié du dispositif avec la position de la vanne et les deux extrémités où va s'appliquer la perte de charge :

Lorsque la vanne est fermée, tout l'air injecté par la base du caisson va passer par le tube et va entraîner avec lui le solide. Si l'on veut ralentir le débit, il suffira d'ouvrir la vanne pour réduire la pression dans le ciel du lit et ainsi diminuer le débit d'air passant par le tube et donc la vitesse à laquelle le solide en sort.

Il va donc falloir reproduire ce système vanne - capteur dans NEPTUNE_CFD. Pour ce faire, nous avons repris les travaux d'Hadrien Benoit qui consistent en un fichier utilisateur fortran uskpdc.F dans lequel est codée cette perte de charge. Le fichier code les pertes de charge pour chaque phase de la façon suivante :

$\dfrac{ \mathrm{d} U_{k,i}}{ \mathrm{d} t} = -K_{k,ij}U_{k,j}$

où $U_{k,i}$ est la vitesse de la phase k dans la direction i et les $K_{k,ij}$ sont les coefficients de la matrice $K_k$ des pertes de charge de la phase k. Les matrices $K_k$ sont définies pour chaque cellules du maillage et doivent être diagonales avec des coefficients positifs ou nuls. La différence de pression dans la direction i, résultant de l'application de ces pertes de charge est :

$\Delta P_i = - K_{k,ij} \rho_k L_i U_{k,i}$

où $L_i$ est la longueur des cellules avec perte de charge dans la direction i et $\rho_k$ la masse volumique de la phase k.

Le fichier nécessite plusieurs données pratiques telles que la surface de la vanne $S$ déterminée sur le maillage via simail, le débit d'air total entrant dans le caisson $Q_{tot}$ ainsi que la consigne de perte de charge $\Delta P_c$ qui eux ont été déterminés grâce à l'expérimentation au LGC. Nous imposerons :

$\Delta P_c = 290 \, \mathrm{mbar}$

Le fichier utilisateur écrit par Hadrien Benoit a été adapté avec ces nouvelles valeurs. Il a également été adapté au maillage fin au niveau du capteur de pression, c'est-à-dire au niveau de la sélection d'une cellule dans le ciel du lit et dans un voisinnage proche de la vanne de façon à en mesurer la pression. Cette cellule, unique pendant le calcul et toujours la même fait office de capteur de pression pour le dispositif de contrôle.

On se référera au document d'Hadrien Benoit "Codage d'une vanne sous NEPTUNE_CFD" pour de plus amples détails sur le code en lui même.

Adaptation au pilote expérimental

Géométrie

L'équipe s'est rendu au LGC de Toulouse afin de voir et de bien comprendre le fonctionnement du pilote expérimental qui fait l'objet de la thèse de M.Benjamin Boissière.

Sur place, il est apparu que la géométrie de la précédente simulation ne correspondait pas exactement à celle du pilote. Il a donc été nécessaire de modifier le maillage afin de se rapprocher de la réalité expérimentale.

  • le bac a été agrandi, il a à présent pour dimensions : 20 cm de largeur, 20 cm de profondeur et 40 cm de hauteur
  • la diamètre interne du tube à été diminué à 34 mm

 

 

Ces modifications au niveau de la géométrie amènent à une dilatation du maillage ce qui nous amène à changer quelques valeurs dans le fichier utilisateur uskpdc.F :

  • la section de la vanne qui est plus grande et que l'on détermine sous Paraview
  • la vitesse de l'air dans la vanne doit diminuer d'environ 50% car le tube est plus petit et la vanne est deux fois plus grande
  • la recherche de la cellule dans le ciel du lit et au voisinage de la vanne doit être recodée puisque la vanne n'est plus au même endroit
  • la consigne de pression doit être diminuée pour atteindre le régime permanent

 

Aération

Nous rappelons que l'aération permet de remplacer la chauffe dans le pilote. En effet, les deux phénomènes permettent l'expension du gaz et l'allégement du lit. Cela évite également la formation de paquets dans le lit fluidisé pouvant boucher les tubes et donc stopper le processus.

Les travaux précédents n'avaient pas encore pris en compte l'aération réalisée sur le pilote expérimental du LGC.

Nous avons ajouté dans le maillage cette entrée d'air à une hauteur à partir du tube de 56 cm.

Nous avons testé trois façon d'intégrer l'aération à la simulation :

  • Aération dès le début mais cela peut entraîner des court-circuitages avec les autres sorties
  • Rampe d'aération qui l'on code dans le fichier utilisateur usclim.F afin que le débit d'air augmente petit à petit jusqu'à atteindre la valeur demandée après 5 secondes
  • Aération ne démarrant que après 5 secondes de simulation, ainsi l'air entre à un endroit non vide du tube et permet l'expension du gaz sans redescendre dans le tube ou sans faire d'autres choses inattendues.

 

 

 

Une autre façon de placer l'aération pourrait être d'installer une rampe d'aération mais cette fois-ci à partir de 5 secondes. Ainsi, l'air arrive à un endroit du tube où celui-ci est rempli de lit fluidisé, et l'aération n'est pas trop brusque et ne peut donc pas pertuber le processus.

Modèles utilisés

NEPTUNE_CFD par l'utilisation de EDAMOX permet de choisir les modèles adéquats qui correspondent aux écoulements en jeu.

 

Force de traînée

La force de traînée s'écrie :

$$ F_D=-\frac{\alpha_p\rho_p}{\tau_{fp}^F}V_{r,i} $$

où le temps de relaxation est :

$$ \frac{1}{\tau_{fp}^F}=\frac{3}{4}\frac{\rho_g}{\rho_p}\frac{<|v_r|>}{d_p}C_d $$

Le code NEPTUNE_CFD permet de tenir compte du caractère variable du coefficient de traînée $C_d$.

Le coefficient de traînée de Wen and Yu Ergun, typique pour des lits fluidisés gaz-particules, est choisi pour les particules.

$$
C_d= \left\{
    \begin{array}{ll}
        C_{d,WY} & \mbox{si } \alpha \le 0.3 \\
       min [C_{d,WY}, C_{d,Erg}] & \mbox{si } \alpha > 0.3
    \end{array}
\right.
$$

 

Les coefficients de traînée d'Ergun et de Wen & Yu sont  :

$$ C_{d,Erg}=200\frac{\alpha_p}{Re_p}+\frac{7}{3} $$

et

$$
C_{d,WY}= \left\{
    \begin{array}{ll}
        \frac{24}{Re_p}(1+0.15Re_p^{0.687})\alpha_p^{-17} & \mbox{si } Re_p<1000 \\
       0.44\alpha_g^{-1.7} & \mbox{si } Re_p \geq 1000
    \end{array}
\right.
$$

 

 

avec

$$ Re_p=\alpha_g\frac{\rho_g<|v_r|>d_p}{\mu_g} $$ 

 

Écoulement compressible

Le fichier utilisateur usphyv.F qui permet le codage des propriétés physiques variables est configuré pour tenir compte de la variabilité de la masse volumique.

Elle s'écrit au moyen de la loi des gaz parfait. La simulation se fait à froid, il n'y a donc pas d'échange de chaleur et la température reste constante donc la masse volumique est une fonction linéaire de la pression.

L'influence de la variation de la masse volumique sur l'écoulement le rend compressible.

 

 

Modèles de turbulence

Pour la phase gaz, le modèle de turbulence isotrope $ k-\epsilon $ est adopté.

Ce modèle fait intervenir des termes dans les équations de transport, de nouvelles constantes sont donc nécessaires afin de fermer le calcul :

$ C_\mu$ $ \sigma_k$ $\sigma_{\epsilon}$ $C_{\epsilon1}$ $C_{\epsilon2}$
0.09 1 1.3 1.44 1.92

Ce modèle est limité dans la mesure où l'énergie cinétique turbulente est surestimée dans les régions d'impact et de ré-attachement, le développement des jets ou sillages ainsi que le rattachement après un décollement sont en général mal prédits.

Cependant c'est un modèle simple et robuste très utilisé car son comportement est bien appréhendé notamment grâce à son ancienneté.

 

Pour la phase particules, le modèle à deux équations de transport q²-q12 est adopté car nous sommes dans le cas d'un écoulement gaz-particules où $ \rho_p \gg \rho_f $.

L'équation de l'agitation des particules est prise en compte avec :

$$ q²=\frac{1}{2} <u_{p,i}^{'} u_{p,i}^{'}>_p $$

Et l'équation de transport de la covariance avec :

$$ q_{fp}= <u_{g,i}^{'} u_{p,j}^{'}> $$

Les collisions sont prises en compte au travers d'un temps caractéristique $ \tau_p^c $.

Et la modélisation du tenseur de contraintes particulaires $ \Sigma_{p, ij} $ se fait avec la viscosité particulaire qui possède deux contributions :

$$ \mu_p=\alpha_p\rho_p(\nu_p^{kin}+\nu_p^{col}) $$

 

Modèle intéractions particules-particules

Les modèles choisi côté particules sont :

  • modèle frictionnel (écoulements très denses)
  • modèle granulaire (écoulements denses, modèle collisionnel)
  • modèle cinétique (écoulements dilués)

Ainsi : $$ \mu_p=\mu^{coll}+\mu^{cin}+\mu^{fric} $$

Paramètres de la simulation

 

Les paramètres du cas de référence se déterminent grâce au pilote expérimental. Le schéma ci-dessous présente les conditions aux limites données par Benjamin Boissière.

 

 

Conditions aux parois

Nous adoptons la condition de non glissement (no-slip) en paroi pour les particules, elles auront donc une vitesse nulle à cette surface.

Cela se traduit par :

$$ (U_{p,\tau})_{wall}=0 $$

$$ K_p(\frac{\partial q_p^2}{\partial n})_{wall}=0 $$

Et du côté de l'air, la condition de friction est choisie.

Conditions aux limites

Nous initialisons le niveau de particules à $z = 0,3 \, \mathrm{m}$ donc avant   la fluidisation le bac est rempli de particules avec une fraction massique de   0.52.

L'aération est imposée à 150 NL/h d'entrée d'air soit $Q_{aération}=5,017 \cdot 10^{-5} \, \mathrm{kg/s}$. Nous choisissons d'activer l'aération au bout de 5 secondes. Pour cela, nous codons la condition dans le fichier fortran usclim.F.

L'entrée d'air qui permet de fluidiser le lit est fixée à 17,5 mm/s donc, pour une section de caisson de 0,04 m. Cela correspond à $ Q_{e,air} = 8,43 \cdot 10^{-4} \, \mathrm{kg/s}$ .

 

L'entrée des particules de SiC est de 50 kg/h de solide. On désire avoir 50% de volume de particules, cela nous donne les débits en air et en particule d'entrée :

$ Q_{air\_entrée\_particules} = 5,21 \cdot 10^{-6} \, \mathrm{kg/s}$

$ Q_{particules\_entrée\_particules} = 1,36 \cdot 10^{-2} \, \mathrm{kg/s}$

 

Le tableau ci-dessous donne les conditions aux limites en entrée imposées dans EDAMOX.

  Entrée d'air Entrée de particules Aération
Débit d'air (kg/s) $ 8.43.10^4 $ $ 5.21.10^6 $ $ 5.02.10^5 $
Débit de particules (kg/s) $ 0 $ $ 1.36.10^2 $ $ 0 $