Simulation

Présentation de la simulation

Cette partie est consacrée à la mise en place du problème sous Neptune : Géométrie, maillage, conditions aux limites et initiales.

Géométrie et maillage

Pour simuler le combusteur expérimental, on utilise une géométrie 2D.

Géométrie

On dessine une colonne de 90 cm par 12.5 cm. Les 10 cm en haut de la colonne sont occupé par un convergent pour assurer l'évacuation des fumées.

Le lit à une hauteur d'environs 36cm, mailler un peu plus de 2 fois cette hauteur permet de maintenir une zone de dégagement au dessus du lit en mouvement et de suivre l'évacuation des gaz.

schéma du lit (source : BEI 2012/2013, gazéification de la biomasse)

Maillage

Le maillage est réalisé sous Xsimail. Le logiciel Neptune ne prend en compte que les maillages 3D. Nous souhaitons réaliser une simulation en 2D pour ne pas avoir de temps de calcul trop long. Nous utilisons donc une géométrie avec une maille d'épaisseur. 

On place 33 mailles dans la largeur et 72 pour les 80 cm de hauteur droite. Le nombre total de mailles est de 2574. Les mailles font 11 mm de hauteur, 3.8 mm de largeur et 10 mm d'épaisseur.

 

conditions aux limites et initialisation

Conditions aux limites

Les tableaux suivant décrivent les conditions aux limites en haut et en bas de la colonne pour les phases et les scalaires. Les conditions pour les cotés sont de type wall.

  phase vitesse type taux de présence
BAS COLONNE 1 (0;0.49;0) inlet $1$
2 (0;0;0) wall $0$
3 (0;0;0) inlet $ 10^{-12}$
HAUT COLONNE 1 (0;0;0) outlet $1$
2 (0;0;0) outlet $ 10^{-12}$
3 (0;0;0) outlet $ 10^{-12}$
  condition
$Y_{CO_{2}}$ dirichlet $10^{-12}$
$Y_{N_{2}}$ dirichlet $0.778$
$\chi_{d}$ dirichlet $1$
$H_{gaz}$ dirichlet enthalpie à 850°C
$H_{olivine}$ dirichlet enthalpie à 850 °C
$H_{char}$ dirichlet enthalpie à 850 °C

Un dirichlet est une condition aux limites sous Neptune. Si l'écoulement est entrant, le flux est imposé à la valeur indiquée. Si le flux est sortant, le dirichlet est imposé et évacué. 

Installation du lit d'olivine

Après mise à l'échelle des conditions opératoires (5.5kg d'olivine, 2g de char, 6.9kg/h d'air) sur notre maillage, on choisit de placer 550 g d'olivine. La masse volumique de l'olivine est de 3040 $kg.m^{-3}$. Avec un taux de présence fixé à 0.4, il faut placer de l'olivine sur une hauteur de 0.3688 m pour atteindre la masse souhaitée.

 Injection du char

Dans l'installation du LGC, le char est injecté par une canne sur le coté du réacteur. Dans notre simulation, il n'y a pas de canne d'injection. Pour simuler l'injection initial du char, on place une couche de char dans la partie supérieur du réacteur. a l'instant initial cette couche tombe sur l'olivine et se mélange à celle-ci sous l'effet de la fluidisation. cette opération simule l'injection de char.

On souhaite injecter 2 g de char.Pour une masse volumique de 740 $kg.m^{3}$ et un taux de présence de 0.00128, la couche de char doit avoir une épaisseur de 0.168 m. Pour injecter d'autre masse initiale de char, on modifiera le taux de présence. Les initialisations sur le char et l'olivine sont réalisées dans la sub-routine usiniv.F

 Vitesse débitante de l'air entrant

              vitesse minimale de fluidisation

Le paramètre prépondérant pour un lit fluidisé est la vitesse minimale de fluidisation. c'est la vitesse la plus faible pour laquelle les particules du lit sont mises en mouvement. Elle se calcul à l'aide la relation empirique de Thonglimp.

\begin{equation}  U_{mf}= \left( \frac{d_{pm}~\phi}{1.75} ~ \frac{(\rho_{p}-\rho_{g})~g~\epsilon^{3}}{\rho_{g}} \right)\end{equation}

$\phi$ est la facteur de forme(supposé égal à 1 pour des sphères). On détermine $\epsilon$ (taux de vide) à l'aide d'une autre corrélation :

\begin{equation}  \epsilon = 1.57~Re_{p}^{0.29}~Ga^{-0.19} \end{equation}

La vitesse minimale de fluidisation est celle de l'olivine, média porteur. La quantité d'olivine est très supérieur à celle de char, celle-ci adopté le comportement hydrodynamique du média porteur.

$U_{mf}=0.075 m.s^{-1}$

             vitesse débitante

On injecte l'air à une vitesse de 0.49 $m.s^{-1}$ (soit un débit de $6.9~kg.h^{-1}$). Cette vitesse est égale à $6~U_{mf}$. Cette vitesse permet d'assurer une bonne fluidisation du lit. De plus, elle permet d'opérer dans maintien conditions où le char est réactif limitant. Le fluide entrant est soit du dioxygène pur, soit de l'air.

Etude de la cinétique de réaction

Dans cette série de simulations, nous avons cherché à tester la différence entre une constante de vitesse constante et une constante de vitesse dépendant des paramètres physiques et géométrique de la réaction. Nous avons codé l'expression de cette constante variable dans les subroutines du logiciel neptune selon le modèle décrit à la page "mise en place du problème : réaction".

Résultats

Les test ont été réalisés avec un gaz entrant composé uniquement de dioxygène pur. Cinq cas ont été testés : constante de vitesse constante à $10^{-3}$, constante de vitesse calculée à partir d'un nombre de Sherwood égal à 2 et trois constante dépendant d'un Sherwood calculé avec des corrélations empiriques. Les corrélations utilisées sont [2]:

  • Analogie Ranz-Marshall pour les tranferts thermiques : $Sh=2+0.6~Re^{1/2}~Sc^{1/3}$,
  • Frössling (1938) : $Sh=2~\epsilon_{mf}+0.61~Rep^{0.48}~Sc^{1/3}$, avec $\epsilon_{mf}$ fixé à 0.4
  • Chakraborty & Howard (1981) : $Sh=2~\epsilon+0.69~Re^{1/2}~Sc^{1/3}$, avec $\epsilon= \alpha_{g}$

Les simulations pour kc en $10^{-3}$ ont été réalisée par le groupe de l'année dernière. Nous allons pouvoir comparer ce cas avec les autres que nous avons codés cette année.

La constante de réaction chimique sera établit avec la loi d'Arrhénius pour une température constante de 850 °C (température raisonnable pour une combustion).

                      Test pour une alimentation en dioxygène pur

Dans cette série de simulation, on alimente le lit avec du dioxygène pur.

Un écart marqué est visible entre le cas $k_{c}$ à  $10^{-3}$ et les cas $k_{c}$ dépendant des caractéristiques du lit. On peut penser qu'un kc tenant compte des paramètres de la réaction représente mieux celle-ci qu'un $k_{c}$ fixé arbitrairement. Cette supposition sera validée dans le paragraphe suivant avec les résultats expérimentaux. Pour les $k_{c}$ fonction de Sherwood, la conversion totale est obtenue pour 180 s de réaction. sur la même durée, la conversion atteint pour le kc en $10^{-3}$ reste inférieure à 0.1. La prise en compte des paramètres du réacteur dan sla cinétique amèliore les résultats antérieures.

Par contre, l'utilisation de différentes corrélations donne des profils de conversion semblables. La corrélation Ranz-Marshall "classique" sera utilisé dans la suite des simulations. De léger écart apparaissent au milieu de la réaction, mais la conversion totale est toujours obtenue pour la même durée de simulation.

                      Test pour une alimentation en air ($O_{2}/N_{2}$)

Dans cette seconde série de simulation, nous avons ajouté au code Neptune un scalaire $N_{2}$ afin de simuler le fluide entrant comme de l'air ambiant (20 % dioxygène et 80 % d'azote). Cela permet d'affiner les résultats par rapport aux résultats obtenue jusqu'ici par les différents groupes de BEI.

Ce gaz entrant est similaire à celui utilisé dans les expérience au LGC. Les simulations sont effectuées pour différentes masse de char injectées (2,4 et 10 g).

Quelque soit la masse injectée les profils de conversion du char simulés correspondent à ceux établit expérimentalement. Ces résultats confirme que le modèle utilisé pour représenter la cinétique de réaction est adapté. L'allure et l'évolution des profils est proche quelque soit la masse injecté. Cette observation est normale tant que $O_{2}$ reste limitant. La prise en compte des transferts thermiques devrait permettre améliorer encore la précision la proximité simulation/réaction.

La conversion totale est atteinte au bout d'un peu plus de 200s pour les masses de 2 et 4g. Il faut atteindre près de 250 s pour atteindre une conversion totale de 10g de char. On note que le profil simulé pour 10 g de char présente un écart significatif avec l'expérience. Dans les simulations, le taux de conversion change avec la masse injectée, alors que les expériences ont montré le contraire. La prise en compte des transferts thermiques pourra corriger ce problème. En effet, pour une masse de char plus importante, l'énergie libérée par la combustion est plus forte. La cinétique est liée à la température du milieu par la loi d'Arrhenius. La prise en compte de la température doit augmenter la vitesse de réaction et permettre d'éviter le retard de conversion quand la masse augmente.

                     Taux de $CO_{2}$ en sortie

Les profils expérimentaux et simulés ont la même allure. Les profils simulés  atteignent le même niveau de valeur (on prend en considération les plaquettes dont la géométrie est la plus proche de notre cas de simulation). Les rejets de $CO_{2}$ sont élevés en début de réaction puis diminuent.

Evolution du diamètre de la particule

L'évolution du diamètre des particules est une propriété importante. Au dela d'un certain diamètre, les particules peuvent être élutrier, c'est à dire projeter hors du réacteur.

Le diamètre est calculée en fonction de la masse de char à chaque instant.

\begin{equation} d=d_{0}~(1-X)^{1/3}  \end{equation}

Où X est le taux de conversion du char.

Le graphique suivant présente l'évolution de la taille des particules (en noir) et le nombre de particule sortant de la colonne (en rouge). Pour obtenir cette seconde donnée, on a placé une sonde sur chaque cellule de la dernière ligne de cellules. Dans chaque cellule on relève le nombre de particules, ensuite on somme pour avoir le nombre total de particules éjecter à chaque instant.

La taille de particule décroît au cours du temps au fur et à mesure que le char est consommé. L'élutrition commence seulement quand les particules deviennent très petit. Pendant L'essentiel de la réaction, aucune particule n'est éjectée hors du réacteur. Toutes les particules sont entièrement consommée pendant la réaction.

On peut utiliser la corrélation de Heider & Levenspiel pour estimer la vitesse d'envollement d'une particule en fonction du diamètre. Si cette vitesse dépasse la vitesse du gaz entrant, les particules seront entraînées hors du lit.

\begin{equation}U_{t}= \left( \frac{18}{{d_{p}^*}^{2}}+\frac{2.335-1.744~\phi}{{d_{p}^*}^{0.5}} \right)^{-1} \left( \frac{\rho_{g}^{2}}{\mu_{g} \left( \rho_{p}-\rho_{g} \right) g} \right)^{-1/3} \end{equation}

Avec $d_{p}^*=\left(d_{p}^{3}~\frac{\rho_{g}~(\rho_{p}-\rho_{g})~g}{\mu_{g}^{2}} \right)^{1/3}$

Avec cette corrélation, on trouve une vitesse d'entrainement inférieure à la vitesse du gaz entrante pour des particules de char d'environ 200 µm. Ce qui est confirmé par la lecture graphique. On trouve par le calcul que les particules sont expulsées à partir de 230 µm. L'élutriation des particules débutent à 220 s.

Comparaison des résistances aux transferts

Une simulation a été réalisée avec une constante de réaction chimique nulle. Ainsi les particules de char sont sonsommées sous le seul effetdu transfert de matière. Le but de cette simulation est de comparé l'influence respectivement des deux phénomènes gouvernant la réaction.

La conversion du char sous le seul controle du transfert de matière est plus rapide que celle prennnant en compte les deux résistance au tranfert. Ce résultat n'est pas surprenant. Comme on peut le voir sur le graphique suivant, quand le diamètre des particules diminue, $k_{t}$ devient très important tandis que $k_{c}$ reste stable. Ainsi sous le seul controle du transfert de matière, la réaction s'intensifierais à mesure de la diminution de $k_{t}$, mais la résistance chimique vient limiter cette effet.

 

Ce graphique à été établit pour une température constante de 850°C, avec une vitesse relative gaz-particule et un coefficient de diffusion moyens.

Echauffement des phases

Avant de relier échange thermique et réaction, il faut étudier les transferts de chaleur entre phases sous le seul effet des mécanisme de convection-diffusion. La prise en compte de l'évolution des propriétés physiques du gaz avec la température est réalisée selon les relations suivantes :

Masse volumique de l'air:

\begin{equation}  \rho_{g}=\rho_{g_{0}}~\frac{T_{0}~M_{g}}{T~M_{g_{0}}}\end{equation}

Dans cette équation $T_{0}$ correspond à la température de référence pour la densité de l'air. Il a été choisi de prendre cette référence à 855°C ($\rho_{g_{0}}=0.313$).

Viscosité laminaire de l'air:

Loi de Sutherland

\begin{equation}  \mu_{g}=\mu_{g_{0}}~\left(\frac{T}{T_{0}} \right)^{3/2}~\frac{T_{0}+S}{T+S} \end{equation}

Avec $S$ : facteur spécifique au gaz ($110,4 K$ pour l'air sur la plage $[170K;1900K]$ et $T_{0}=273,15K$)

Coefficient de diffusion:

Les relations utilisées pour la variation des coefficients de diffusion avec la température sont explicitées ici.

Afin de prendre en compte la diffusion dans chacune des espèce, la relation de Wilke est utilisée:

\begin{equation} D_{O_{2}}=\left( \frac{X_{CO_{2}}}{D_{O_{2}/CO_{2}}} + \frac{X_{N_{2}}}{D_{O_{2}/N_{2}}} + \frac{1-X_{CO_{2}}-X_{N_{2}}}{D_{O_{2}/O_{2}}} \right)^{-1} \end{equation}

Résultats

Une simulation ne prenant en compte que l'échauffement du lit (pas de phénomène de radiation) a été confronté à un modèle 0D codé sous Matlab simulant ce même échauffement. Les conditions initiales sont telles que $T_{lit}=1148.15K$ et que la température de l'air en entrée $T_{air,e}=1168.15K$. L'échauffement sans terme radiatif de l'ensemble du lit par l'air prenant beaucoup trop de temps, il a été décidé de n'observer que le début de l'échauffement et d'observer les pentes des courbes. Le résultat est présenté dans le graphique suivant:

Il apparaît que les résultats NEPTUNE_CFD collent bien au modèle 0D de Matlab. Ceci semble confirmé que les transfert thermiques codés dans le user usthgp. sont bien codés. Il semble y avoir un léger décrochage à la toute fin de simulation, ceci pouvant être dû au fait que sous NEPTUNE, la densité, la viscosité laminaire (et plus tard les $Cp$) varient avec la température. L'écart relatif calculé entre les deux coefficient des droites est de 1,5%. Le cas de réchauffement est donc validé.

Etude du couplage réaction/tranferts thermiques

Cette partie est un combinaison des précédentes : On regarde l'évolution de la cinétique et de la température dans le réacteur. On utilise l'analogie de Ranz-Marshall pour le nombre de Sherwood. Le gaz de fluidisation est de l'air à $0,49~\mathrm{m.s^{-1}}$. Les propriétés physiques (densité, viscosité, coefficient de diffusion) évoluent avec la température.

Resultats

A l'heure de la rédaction de ce site, l'implémentation et le codage du user ushsig.F ne sont pas encore finalisés. Cependant, des premiers résultats encourageants ont été observés comme présentés dans le graphe ci-dessous:

Le graphique présente l'évolution de la conversion ainsi que de la température du char en fonction du temps. Il apparaît que l'introduction de la température semble réhausser la courbe "sans température" des données expérimentales et semble même la sur-estimer pour des temps longs. De son coté la température du char semble augmenter dans des proportions raisonnables. Cependant, un tracé de la température de l'air en fonction de la température montre une élévation très limitée (1,5°C maximum) tout au long de la simulation, bien loin des 25°C observés lors des expériences. Des améliorations dans le codage du fichier user ushsig.F sont à faire afin d'améliorer ce point.