Mise en place du problème

Cette partie est consacrée à l'explication de la mise en équation du problème.

La modélisation du comportement chimique et physique du lit est présentée. Les bilans utilisés et leurs implémentations sous Neptune sont décrits.

Hypothèses et modèles

Hypothèses

Un certain nombre d'hypothèses ont été faites pour simplifier le problème :

  • gaz parfait
  • gaz entrant composé exclusivement de dioxygène (une étude prenant en compte une injection d'air a aussi été réalisée par la suite)
  • température en paroi constante
  • réaction de combustion unique et complète (pas de sous-produit)

Modèles

Un certains nombre de modèles codés dans Neptune sont utilisés pour représenter le comportement des phases. Ces différents modèles sont résumés dans le tableau suivant.

 

gas olivine char
turbulence model k- epsilon q2-q12 q2-q12
turbulence reverse coupling   small inclusions small inclusions
particles interactions : friction   fluxes fluxes
particles interactions   granulometrie   fluxes fluxes
particles interactions : kinetic   uncorellated collisions uncorellated collisions

polydisperse model

  activated activated
wall boundary condition   no slip no slip
drag law   Wen&Yu Ergun Wen&Yu Ergun

Le modèle $k-\epsilon$

Au sein du lit fluidisé se développe des écoulements pleinement turbulents. C'est pourquoi le modèle de turbulence $k-\epsilon$ est utilisé. Ce modèle est basé sur l'écriture des vitesses dans les équations de Navier-Stokes comme la somme d'un écoulement moyen et d'un écoulement instantané: $u=\bar{U}+u'$ ainsi que sur l'hypothèse de Boussinesq: $R_{ij}=\mu_{t}\left( \frac{\partial U_{i}}{\partial x_{j}} + \frac{\partial U_{j}}{\partial x_{i}}\right)$.

La réécriture des équations permet de faire apparaître les équations suivantes:

\begin{equation} \frac{\partial k}{\partial t} + U_{k} \frac{ \partial k}{ \partial x_{k}} = \frac{\partial}{\partial x_{k}} \left[ \left( \nu+ \frac{C_{\mu}k^{2}}{\sigma_{k} \epsilon}\right) \frac{\partial k}{\partial x_{k}} \right] + \frac{C_{\mu}k^{2}}{\epsilon} \left( \frac{\partial U_{i}}{\partial x_{k}} + \frac{\partial U_{k}}{\partial x_{i}} \right) \frac{\partial{U_{i}}}{{\partial x_{k}}} - \epsilon\end{equation}

\begin{equation} \frac{\partial \epsilon}{\partial t} + U_{k} \frac{ \partial \epsilon}{ \partial x_{k}} = \frac{\partial}{\partial x_{k}} \left[ \left( \nu+ \frac{C_{\mu}k^{2}}{\sigma_{\epsilon} \epsilon}\right) \frac{\partial \epsilon}{\partial x_{k}} \right] + C_{\epsilon_{1}} C_{_\mu} k \left( \frac{\partial U_{i}}{\partial x_{k}} + \frac{\partial U_{k}}{\partial x_{i}} \right) \frac{\partial{U_{i}}}{{\partial x_{k}}} - C_{\epsilon_{2}} \frac{\epsilon^2}{k}\end{equation}

Les constantes de ce modèle sont répertoriées dans le tableau suivant:

$C_{\mu}$ $\sigma_{\epsilon}$ $\sigma_{k}$ $C_{\epsilon_{1}}$ $C_{\epsilon_{2}}$
$0.09$ $1.22$ $1$ $1.44$ $1.92$

Le modèle q2-q12

Le modèle $q_{2}-q_{12}$ est un modèle à 2 équations, l'une sur l'agitation des particules $q_{2}^2=1/2 <u_{2,i}'u_{2,i}'>_{2}$ et l'autre sur le transport de la covariance $q_{12}=\overline{u_{1,i}'u_{2,j}'}$

Ce modèle prend en compte les transfert interfaciaux via la loi de traînée (cf § suivant) ainsi que les collisions à travers un temps de collision $\tau_{p}^{c}$. Le tenseur des contraintes particulaire $\Sigma_{p,ij}$ est modélisé à travers d'une viscosité particulaire comprenant 2 contributions (cinétique et collisions) $\mu_{p} = \alpha_{p} \rho_{p} \left( \nu_{p}^{kin}+\nu_{p}^{col} \right) $

Le modèle est adapté pour des écoulements gaz-particules avec $\rho_{p}\gg \rho_{g}$

Les lois de trainées Wen & Yu Ergun

La force de trainée s'oppose au mouvement de la particule: $F_{D}=-\frac{\alpha_{p} \rho_{p}}{\tau_{fp}^{F}} V_{r,i}~\mathrm{avec}~\frac{1}{\tau_{fp}^{F}}=\frac{3 \rho_{g} <V_{r}>}{4 \rho_{p} d_{p}} C_{D}$

Dans le logiciel, la loi de traînée est modélisée à l'aide de la relation de Wen & Yu - Ergun:

\begin{equation} C_{D}=\left\{ \begin{array}{ll} C_{D,wy} & \mathrm{si}~\alpha_{p}\le 0,3 \\ \mathrm{min}[C_{D,wy};C_{D,erg}]& \mathrm{sinon} \end{array} \right. \end{equation}

Dans cette équation:

\begin{equation} C_{D,erg}=200\frac{\alpha_{p}}{Re_{p}}+\frac{7}{3} \end{equation}

\begin{equation} C_{D,wy}=\left\{ \begin{array}{ll} C_{D,wy}=\frac{24}{Re_{p}} \left( 1+ 0.15 Re_{p}^{0.687}\right) \alpha_{p}^{-1.7}& \mathrm{si}~Re_{p} < 1000 \\ C_{D,wy}=0.44\alpha_{p}^{-1.7}& \mathrm{sinon} \end{array} \right. \end{equation}

 

 

 

Réaction

La réaction de combustion dans le lit est celle du char (phase solide) par du dioxygène (gazeux).

Réaction de combustion

La réaction considérée dans le lit est celle mettant en jeu le carbone et le dioxygène d'ordre 1 par rapport au dioxygène telle que :

\begin{equation*} C + O_{2} \to CO_{2} \end{equation*}

Dans la suite les phases seront désignées par la lettre k avec :

  • k=1, le gaz
  • k=2, l'olivine
  • k=3 le char

Expression de la constante de vitesse

L'évolution d'une particule de char pendant la combustion est décrite par le modèle du "coeur rétrécissant". La réaction se déroule de manière homogène, en surface de la particule. La constante de vitesse tient compte de deux phénomènes: la réaction et transfert de matière.

 

Cinétique de réaction : kr

La constante de réaction caractérise l'équilibre de la réaction. Elle ne dépend que de la réaction considérée et de la température. La valeur de kr est relié à ces paramètres par la loi d'Arrhenius.

\begin{equation} k_{r}=k_{0}T \mathrm{exp}\left(-\frac{E_{a}}{R~T}\right) \end{equation}

Avec $E_{a}$= énergie d'activation, $R$ = constante d'état des gaz parfait et $k_{0}$ = facteur pré-exponentiel

dans notre cas $Ea=75189~J/mol$ et $k_{0}=1.699$

Transfert de matière kt

Le second phénomène limitant est le transport de l'oxygène dans la couche gazeuse entourant la particule jusqu'a la surface de cette particules pour réagir. Ce phénomène est modélisé à l'aide d'une corrélation de type Ranz-Marshall.

On commence par calculer les nombres de Schmidt et Reynolds.

\begin{equation*} Sc=\frac{\mu}{\rho D_{O_{2}}}  \end{equation*}

\begin{equation*} Re=\frac{\rho~d~V_{d}}{\mu} \end{equation*}

$V_{d}$= vitesse relative des phase gaz et dipsersée.

On utilise ensuite la corrélation pour calculer le nombre de Sherwood :

\begin{equation} Sh=2+0.6 Re^{1/2} Sc^{1/3} \end{equation}

Pour Re<200 et Sc<250

On peut alors calculer kr

\begin{equation} k_{t}=\frac{D_{O_{2}}~Sh}{d} \end{equation}

L'expression du calcul de $D_{O_{2}}$ est développée ici.

Expression du taux de réaction

La constante globale de réaction est la combinaison de ces deux résistances :

\begin{equation} k_{c}=\left(\frac{1}{k_{r}}+\frac{1}{k_{t}}\right)^{-1} \end{equation}

On calcul alors le taux de réaction de la phase char :

\begin{equation} \Gamma_{c}=-\frac{\alpha_{c}~k_{c}~M_{c}~C_{O_{2}}^{\infty}~S}{V_{p}} \end{equation}

Avec $C_{O_{2}}^{\infty}=$ concentration en dioxygène loin de la particule, $S$:surface de la particule et $V_{p}$: volume de la particule.

calcul de $CO_{2}^{\infty}$

On écrit avec loi des gaz parfaits : $CO_{2}^{\infty}=\frac{n_{O_{2}}}{V_{g}}=\frac{m_{O_{2}}}{M_{O_{2}}~V_{g}}=\frac{(1-(Y_{O_{2}}+Y_{N_{2}}))~\rho_{g}~V_{g}}{M_{O_{2}}~V_{g}}=\frac{(1-(Y_{O_{2}}+Y_{N_{2}}))~\rho_{g}}{M_{O_{2}}}$.

Avec :

  • $n_{O_{2}}$ : quantité dioxygène  (mol)
  • $M_{O_{2}}$ : masse molaire de l'oxygène (kg)
  • $V_{g}$ : volume de gaz ($m^{3}$)

Le rapport $\frac{S}{V_{p}}$ pour une sphère vaut $\frac{6}{d}$, avec $d$ diamètre d'une particule.

On obtient l'expression suivante pour le taux de réaction :

\begin{equation} \Gamma_{c}=-\frac{6~\alpha_{c}~k_{c}~M_{c}~\rho_{g}~(1-(Y_{CO_{2}}+Y_{N_{2}}))}{d~M_{o_{2}}} \end{equation}

On implémente les expressions de $k_{c}$ et de $\Gamma_{c}$ dans le fichier user usphyv.F du logiciel. On définit $\Gamma$ via une fonction user qui est ensuite pris en compte dans les bilans de masse, quantité de mouvement et enthalpie. C'est une fonction globale définie dans un fichier et utilisé dans l'ensemble des subroutines.

Coefficient de diffusion

Cette partie est dédiée au calcul des coefficients de diffusion. La méthode de Lennard-Jones est utilisée pour ces calculs [3].

Expression

\begin{equation}  D_{AB}=1,858 \times 10^{-7}~\frac{[T^{3}~(M_{A}+M_{B})/(M_{A}~M_{B})]^{1/2}}{P~\sigma_{AB}^{2}~\Omega_{AB}^{(1.1)*}} \end{equation}

Avec :

  • $D_{AB}$ le coefficient de diffusion (m2/s)
  • $M_{A}$ et $M_{B}$ les masses molaires des constituants (g/mol)
  • $T$ la température absolue (K)
  • $P$ la pression (bar)

$  \sigma_{AB}=1/2~(\sigma_{A}+\sigma_{B}) $ le rayon moyen des molécules (Å)

$ \Omega_{AB}^{(1,1)*}=\Omega^{(1,1)*}$ les intégrales de collisions du modèle de potentiel 12-6 de Lennard-Jones (hypothèse de sphères dures).

On évalue $\Omega^{(1,1)*}$ avec la relation suivante :

\begin{equation*} \Omega^{(1,1)*}= 1,06036~(T^{*})^{-0,1561}+0,19300~\mathrm{exp}[-0,47635~T^{*}] \end{equation*} \begin{equation} +1,03587~\mathrm{exp}[-1,52996T^{*}]+1,76474~\mathrm{exp}[-3,89411T^{*}] \end{equation}

$T^{*}$ est défini par : \begin{equation}T^{*}=\frac{k~T}{\epsilon_{AB}} \end{equation} et \begin{equation} \frac{\epsilon_{AB}}{k}=[(\epsilon_{A}/k)~(\epsilon_{B}/k)]^{1/2}\end{equation}

Résultats

composé M (g/mol)  $\sigma$ (Å) $\epsilon / k$ (K)
$O_{2}$ 32 3,323 137
$CO_{2}$ 44,01 3,703 266,1
$N_{2}$ 28 3,568 113

Pour une pression de 1 bar et à une température de 850°C, on obtient pour de l'air composé de dioxygène pur:

$D_{O_{2}/O_{2}}=2,19.10^{-4}$ m2/s

$D_{CO_{2}/O_{2}}=1,59.10^{-4}$ m2/s

Afin d'éviter trop de calculs sous NEPTUNE_CFD, des régressions polynomiales sur une gamme de température [800°C;2000°C] ont été réalisées pour un gaz entrant composé de dioxygène pur ou d'air ($O_{2}+N_{2}$)

  • Coefficients de diffusion pour du dioxygène pur

  • Coefficients de diffusion pour de l'air

Afin de prendre en compte la diffusion de l'oxygène dans chacune des deux autres espèce $CO_{2} ~ \mathrm{et} ~ N_{2} $, on utilise la relation de Wilke:

\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}

Bilan de masse et de quantité de mouvement

Bilan de masse

Sur chaque phase, le bilan de masse s'écrit :

\begin{equation} \frac{\partial \alpha_{k} \rho_{k}}{\partial t}+ \frac{\partial \alpha_{k} \rho_{k} U_{k}}{\partial x_{i}} = \Gamma_{k} \end{equation}

De plus, $\Gamma_{1}=-\Gamma_{3}$=$\Gamma_{c}$ (car une mole de char réagit avec une mole de dioxygène) et $\Gamma_{2}=0$ car l'olivine est inerte.

Bilan de quantité de mouvement

Le bilan de quantité de mouvement pour chaque phase s'écrit de la même façon :

\begin{equation*} \alpha_{k}~\rho_{k}~\frac{\partial U_{k}}{\partial t} +  \alpha_{k}~\rho{k}~U_{k,i}~\frac{\partial U_{k,i}}{\partial x_{j}} = \end{equation*} \begin{equation}-\alpha_{k}~\frac{\partial P_{g}}{\partial x_{i}}+\alpha_{k}~\rho{k}~g_{k} + \frac{\partial[\alpha_{k}~\rho_{k}~<u_{k,i}'~u_{k,j}'+\Theta>]}{\partial x_{i}} + I_{k,i}'+ [U_{\sigma,i}-U_{k,i}]~\Gamma_{k} + \delta_{k,i}~\sum{S_{k,j}}  \end{equation}

$I_{k,i}=\frac{\alpha_{k}~\rho{k}~V_{k,i}}{\tau_{k}}$ avec$V_{k,i}$ la ième vitesse de la phase k. Ce bilan modélise le transfert de quantité de mouvement entre le fluide et les particules. $\delta{k,i}$ est l symbole de Kronecker (où i et k désignent les différentes phases). $S_{k,j}$ est le terme modélisant l'effets des collisions entre phases. le termes en $\Gamma_{k}$ décrit les transferts de quantité de mouvement liés à la réaction. La vitesse de la phase gazeuse après réaction est supposée égale à la vitesse de la particule.

Scalaires

Deux scalaires résolus sont définis YCO2 et $\chi_{d}$ (un scalaire $N_{2}$ sera ajouté pour l'étude en injection d'air). On écrit pour chacun d'entre eux une équation de transport. Ces deux scalaires représentent le transport de deux paramètres clés de la simulation : le taux de CO2 rejeté et la conversion du char.

\begin{equation} \alpha_{1}~\rho_{1}~\frac{\partial Y_{CO_{2}}}{\partial t} + \alpha_{1}~\rho_{1}~U_{1,j}~\frac{\partial Y_{CO_{2}}}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} \left ( \alpha_{1}~\rho_{1}~\frac{\mu_{1}}{\sigma_{1}}~\frac{\partial Y_{CO_{2}}}{\partial x_{j}} \right) - Y_{CO_{2}}~\Gamma_{1}+\phi_{1} \end{equation}

On définit $\rho_{g}=\frac{P_{ref}}{R~T_{g}}~\left(\frac{Y_{O_{2}}}{W_{O_{2}}}+\frac{Y_{CO_{2}}}{W_{CO_{2}}}\right)^{-1}$

\begin{equation} \alpha_{3}~\rho_{3}~\frac{\partial \chi_{d}}{\partial t} + \alpha_{3}~\rho_{3}~U_{3}~\frac{\partial \chi_{d}}{\partial x_{i}}= -\frac{\partial \alpha_{3}~\rho_{3}~<u_{3}'~\chi_{d}>}{\partial x_{i}} - \chi_{d}~\Gamma_{c} \end{equation}

$\chi_{d}$ étant défini par $\chi_{d}=\left(\frac{d}{d_{0}}\right)^{3}$, où $d$ est le diamètre d'une particule. Le calcul de $\chi_{d}$ permet d'évaluer l'évolution du diamètre de particule avec la réaction. Il est aussi possible de remonter à un $\chi_{d}$ moyen peut aussi être obtenu via le bilan de masse sur la phase char (hypothèse : toutes les particules réagissent identiquement et qu'il n'y a pas d'élutriation).

Par la suite, on ajoutera aussi un troisième scalaire : $Y_{N_{2}}$. Ce scalaire, codé pour correspondre à la fraction d'azote dans le milieu, permet de faire évoluer le gaz de fluidisation de dioxygène pur à air ambiant. Cette évolution permet de faire coller le gaz de fluidisation à celui utilisé dans l'expérience du LGC. On prendra une composition 80%$N_{2}$ / 20%$O_{2}$. Ce scalaire suit l'équation de transport suivante :

\begin{equation} \alpha_{1}~\rho_{1}~\frac{\partial Y_{N_{2}}}{\partial t} + \alpha_{1}~\rho_{1}~U_{1,j}~\frac{\partial Y_{N_{2}}}{\partial x_{j}} = \frac{\partial}{\partial x_{j}} \left ( \alpha_{1}~\rho_{1}~\frac{\mu_{1}}{\sigma_{1}}~\frac{\partial Y_{N_{2}}}{\partial x_{j}} \right) - Y_{N_{2}}~\Gamma_{1}+\phi_{1} \end{equation}

Les codes correspondant à ces deux scalaires sont dans le user ustssp.F de NEPTUNE_CFD. L'expression du diamètre est rentré via une fonction user dans usphyv.F.

         Implémentation du fichier source ustssp.F.

Le terme sources pour le bilan de quantité de mouvement sur un scalaire s'écrit sous la forme TS = TSA + TSBdX

TSA correspond aux termes source explicite et TSB à sa dérivé par rapport aux scalaires.

Pour $\chi_{d}$

\begin{equation*} TSA =  \frac{\Gamma_{c}}{\alpha_{3}}~\chi_{d} \end{equation*}

\begin{equation*}  TSB = \frac{\Gamma_{c}}{\alpha_{c}} \end{equation*}

Pour $Y_{CO_{2}}$

\begin{equation*} TSA =  \frac{\Gamma_{c}}{\alpha_{1}}~\left( \frac{M_{CO_{2}}}{M_{c}} + Y_{CO_{2}} \right) \end{equation*}

\begin{equation*}  TSB = \frac{\Gamma_{c}}{\alpha_{1}} \end{equation*}

Pour $Y_{N_{2}}$

La quantité de ce scalaire dans le lit ou le courant entrant ne change pas au cours de la simulation. Mais la proportion d'azote varie avec la production de $CO_{2}$ et la consommation $O_{2}$. Les termes sources entrés dans le code modélisant ce comportement sont :

\begin{equation*} TSA =  \frac{\Gamma_{c}~Y_{N_{2}}}{\alpha_{1}}  \end{equation*}

\begin{equation*}  TSB = \frac{\Gamma_{c}}{\alpha_{1}} \end{equation*}

 

Bilan enthalpique

Bilan

Sur chaque phase, le bilan enthalpique s'écrit omme suit. Les échanges thermiques se décomposent en trois phénomènes : convection-diffusion, rayonnement et transfert de chaleur dû à la réaction. Le transfert du aux choc particules-particules sont négligés. En l'absence de source de chaleur, le bilan s'écrit :

\begin{equation*} \alpha_{k}~\rho_{k}~\frac{\partial H_{k}}{\partial t}+\alpha_{k}~\rho_{k}~U_{k}~\frac{\partial H_{k}}{\partial x_{j}} = \end{equation*}

\begin{equation} \frac{\partial}{\partial x_{j}} (\alpha_{k}~\rho_{k}~U_{k}~\frac{\partial H_{k}}{\partial x_{j}})+[H_{\sigma}-H_{k}]~\Gamma_{k}  \end{equation}

Il faut ajouter à ce bilan un terme de radiation. Ce terme est codé dans les sub-routines de Neptune avec un modèle de Rosseland. Les émissivités du char et de l'olivine sont respectivement 0.9 et 0.82 [1].

$H_{\sigma}$ est l'enthalpie de transfert qui prend en compte les échanges de chaleur entre les phases réactives. 

Calcul de l'enthalpie de transfert.

Cette enthalpie modélise l'échauffement du char et du gaz du à la réaction entre le char et l'oxygène. Elle s'écrit:

\begin{equation} H_{\sigma}=H_{CO_{2}}^{0}+(T_{3}-T_{ref})~(Cp_{c}-Cp_{CO_{2}})-(T_{1}-T_{ref})(Cp_{O_{2}}) \end{equation}

Avec :

$H_{\sigma}$=enthalpie formation du $CO_{2}$ = -393500 J/mol

$Cp_{c}$= capacité calorifique du char = 1700 J/kg

$Cp_{CO_{2}}$= capacité calorifique du dioxyde de carbone = 2000 J/kg

$Cp_{O_{2}}$ capacité calorifique du dioxygène

T désigne les températures de chaque phases avec $T_{ref}$ = 273.15 K

L'ensemble des capacités calorifiques varient avec la température. Cependant, ces variations ne seront pas prises en compte lors des simulations. Intégrer ces variations pourrait faire l'objet d'une amélioration future. Dans cette optique, le tableau ci-dessous récapitule l'évolution de la capacité calorifique des espèces gazeuses avec la température selon la formule [6]:

\begin{equation} Cp_{i} = A+B~t+C~t^{2}+D~{t^3}+\frac{E}{t^2}\end{equation}

Avec $t=T/1000$ et $T$ en $K$. Les $Cp$ sont données en $\mathrm{J/mol/K}$

Coefficient $N_{2}$ $O_{2}$ $CO_{2}$
Gamme T $[500;2000]$ $[700;2000]$ $[1200;6000]$
A $19.50583$ $30.03235$ $58.16639$
B $19.88705$ $8.772972$ $2.720074$
C $-8.598535$ $-3.988133$ $-0.492289$
D $1.369784$ $0.788313$ $0.038844$
E $0.527601$ $-0.741599$ $-6.447293$

Pour des températures comprises entre $273 \mathrm{K}$ et $1273 \mathrm{K}$, la capacité calorifique du char suit la corrélation [7] ($Cp_{c} \mathrm{kJ/kgK}$:

\begin{equation} Cp_{c}=0.42+2.09.10^{-3}T+6.85.10^{-3}T^{2} \end{equation}

Implémentation des fichiers sources

Il s'agit de coder dans l'user ushig.F l'expression du terme source $H_{\sigma}$. Neptune considère cette source sous la forme \begin{equation} HSHIG=BHGIG + \Sigma~ AHSIG(i)~dH(i) \end{equation}

           Pour la phase gaz

BHGIG = $-H_{\sigma}$

AHSIG(k) = $ \delta_{1,k}$

          Pour l'olivine

BHGIG = 0 , car l'olivine est inerte

AHSIG(1) = 0

        Pour le char

BHGIG = $H_{char}$

AHSIG(k) = $ \delta_{3,k}$

Où $\delta_{i,k}$ est le symbole de Kronecker.