Mise en place du problème

Le logiciel NEPTUNE_CFD

Présentation du logiciel NEPTUNE_CFD

Le logiciel utilisé pour réaliser les simulations numériques est le logiciel NEPTUNE_CFD développé par un consortium comprenant EDF, le CEA, AREVA et l'IRSN. La version utilisée ici est la branche NEPTUNE_CFD V1.08@Tlse, développée par l'IMFT et dédiée à la simulation des écoulements gaz-particules et en particulier des lits fluidisés.

NEPTUNE_CFD est basé sur une méthode RANS (Reynolds Average Navier Stokes Equations) et utilise une méthode des volumes finis de type cell-centered qui co-localise toutes les variables au centre de chaque cellule du maillage.

Les équations sont des moyennes d'ensemble de phase, pondérées par les fractions volumiques de gaz pour la phase gazeuse et complétée par la théorie des écoulements granulaires et les effets de la turbulence pour la phase dispersée.

NEPTUNE_CFD est un code de calcul non structuré et parallèle qui permet :

  • de modéliser des écoulements multiphasiques (de 1 à 20 phases), réactifs, 3D et turbulents
  • de modéliser des géométries simples ou complexes (échelles industrielles)
  • d'utiliser de nombreux modèles de turbulence, de masse, de quantité de mouvement, d'enthalpie totale mais aussi de modéliser des scalaires passifs

 

Définition des paramètres

Dans le cadre de notre étude, 3 phases entrent en jeu :

  • Le char, réactif de la réaction de combustion
  • L'olivine, media fluidisé permettant de transmettre la cinétique du lit fluidisé aux particules du char de quelques millimètres, difficilement fluidisables.
  • La phase gazeuse, qui comprend un mélange de dioxyde de carbone et de dioxygène, pondéré par les fractions massiques de chacune des espèces.

L'air est le gaz de fluidisation : il est injecté à pression atmosphérique et à teméprature ambiante (environ 20°C)

Paramètres :

La température de référence est fixée à 293,14 °C et la capacité thermique massique de l'air est prise égale à 1400 J/kg.K. On choisira les propriétés de l'air comme propriétés de la phase gaz, le dioxygène étant largement majoritaire dans le mélange $O_{2}/CO_{2}$

Le tableau ci-dessous regroupe les principales caractéristiques des phases présentes avant la réaction (gaz, olivine et char) :

Rappelons qu'alpha correspond au taux de présence de l'espèce k.

Attention : la densité donnée dans le tableau ci-dessus n'est valable qu'à une température ambiante d'environ 20°C. La réaction de combustion se fait quant à elle à une température de l'ordre de 855°C : il faut donc recalculer la masse volumique de l'air et sa viscosité.

On trouve à 855°C :

  • $\rho=0.313 kg/m^{3}$
  • $\mu=4,16. 10^{-5} Pa.s$

 

Hypothèses et modèles choisis

Hypothèses :

Plusieurs hypothèses doivent être faites afin de faciliter les calculs. Elles sont énumérées ci-après :

  • Le système est adiabiatique
  • L'échauffement et l'exothermicité de la combustion sont dans un premier temps négligés
  • Le coefficient de transfert de matière autour de la particule p est dans un premier temps considéré comme constant et égal à 0.001 S.I
  • Les gaz sont considérés comme parfaits
  • L'air est constitué uniquement de dioxygène (on ne prend pas en compte la présence du diazote)
  • Le système est constitué de 3 phases : char, olivine et gaz (on regroupe sous une seule phase les gaz)
  • La diffusion est un phénomène limitant
  • Les réactions en chaîne de la combustion ne sont pas prises en compte

 

Choix des modèles :

 

 

Description des modèles :

  • Modèles de turbulence :

Modèle k-$\epsilon$ : il modélise l'énergie cinétique turbulente ainsi que la dissipation. C'est un modèle simple,  à deux équations. Il présente l'avantage d'être connu et prévisible. C'est ce modèle que l'on choisit pour la phase continue (gaz).

Modèle $q_{2}-q_{12}$ : Ce modèle est utilisé pour modéliser la turbulence dans la phase dispersée, ie dans la phase solide soit olivine et char. C'est un modèle à 2 équations de transport : l'une sur l'agitation des particules $q^{2}_{2}$;autre sur le transport de la covariance $q_{12}=\overline{u'_{1,i}u'_{2,j}}$. Les collisions sont prises en compte et le tenseur des contraintes est modélisé au travers de la viscosité particulaire.

  • Modèle de trainée : Wen&Yu Ergun

Le choix d'un modèle de trainée adapté est primordial. Il faut trouver la bonne expression du coefficient de trainée. Le modèle choisi pour les phases dispersées est le modèle Wen&Yu Ergun : la loi de trainée est celle de Wen & Yu, limitée par celle d'Ergun pour les écoulements denses. Le coefficient de trainée est alors défini comme suit :

$$
C_{d} = \left\{
    \begin{array}{ll}
        C_{d,WY} & \mbox{si } \alpha_{p} \le 0,3 \\
        min[C_{d,WY},C_{d,Erg}] & \mbox{si} \alpha_{p} > 0,3
    \end{array}
\right.
$$

avec $$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_{g}^{-1,7} Re_{p} < 1000 \\
        0,44 \alpha_{g}^{-1,7} Re_{p} \ge 1000 \\
    \end{array}
\right.
$$

On rappelle que :

$ Re_{p} = \alpha_{g} \frac{\rho_{g} \langle | v_{r} |\rangle d_{p}}{\mu_{g}} $

p désigne la pième phase solide et g la phase gazeuse.

 

 

 

 

 

Modélisation des équations

Équation de la combustion : 

L'équation de la combustion s'écrit :

$C+O_2\longrightarrow CO_2$

k=1 gaz g (air+$CO_2$)

k= 2 (olivine)

k= 3 (char C)

Conservation de la masse :

Le bilan de masse pour chaque phase k s’écrit :

$\frac{\partial \alpha_k \rho_k}{\partial t} + \frac{\partial \alpha_k \rho_k U_{k,i}}{\partial x_i}=\Gamma_k$ avec  $\sum \Gamma_k​=0$

  • Pour la phase gaz ($O_{2}+CO_{2}$) :

$\Gamma_g=\sum \Gamma_{p \rightarrow g}=-\Gamma_C$

  • Pour l'olivine (k=2) :

$\Gamma_{olivine}=0 $ car l'olivine est inerte

Expression du taux de transfert de masse d'une particule réactive de char :

Soit $n_k$ le nombre de moles de l'espèce k

On choisit un modèle à coeur rétrécissant pour modéliser la variation du nombre de moles de char en fonction du temps. Le choix de ce modèle sera justifié par la suite :

$\frac{dn_c}{dt}=k_{c}S(C_{{O_2}_{s}}-C_{{O_2}_{\infty}})$

avec $n_c$ en $mol$, kc en $m/s$, C en $mol/m^3$ et $S$ en $m^2$

Cette expression sera validée par les expériences réalisées en laboratoire.

La concentration en dioxygène à la surface de la particule solide étant nulle on obtient :

$\frac{dn_c}{dt}=-k_{c}SC_{{O_2}_{\infty}}$ (1)

avec kc le coefficient de transfert de matière autour de la particule p et nc le nombre de moles de carbone.

On prendra dans un premier temps la constante cinétique $k_c$ constante est égale à $10^{-3} m^³ $de gaz$/m^²$ de solide.s

(1) peut s'écrire : $\frac{dm_c}{dt}=-k_{c}M_{c}SC_{{O_2}_{\infty}}$

afin d'obtenir la variation de la masse de char en fonction du temps.

Expression du taux de transfert de masse $\Gamma_c$ :

Le taux de transfert de masse $\Gamma_{C}$ en $kg/m³/s$ d'une particule réactive de char s'écrit :

$\Gamma_c=n(\frac{dm_c}{dt})$ avec n le nombre de particules de char par unité de volume en $m^{-3}$.

n s'exprime de la façon suivante :

$n=\frac{\rho_{c} \alpha_c}{m_c}$ avec m la masse d'une particule de char, $\alpha \rho$ la masse apparente du char.

Donc, $n=\frac{\rho_{c} \alpha_c}{V_{c} \rho_c}=\frac{\alpha_c}{V_c}$ avec $V_c$ le volume d'une particule de char.

On obtient donc :

$$\Gamma_c=-\frac{\alpha_{c}k_{c}M_{c}SC_{{O_2}_{\infty}}}{V_c}$$

 

Avec $C_{O_2}=\frac{Y_{O_2} \rho_{mélange} \alpha_{g}}{M_{O_2}}$ soit $C_{O_2}=\frac{(1-Y_{CO_2}) \rho_{mélange} \alpha_{g}}{M_{O_2}}$ avec $C_{0_2}$ en $mol/m^3$

$Y_{CO_2}$ est le scalaire transporté par la phase gazeuse

$\rho_{mélange}$ est calculé à l'aide de la loi des gaz parfaits dans le scripts usphyv.F

 

 

Bilan de quantité de mouvement :

Pour chaque phase k, le bilan de quantité de mouvement s'écrit :

$\alpha_{k}\rho_{k}\frac{\partial​ U_{k,i}}{\partial t}+\alpha_k \rho_k U_{k,i}\frac{\partial U_{k,i}}{\partial x_j} $

$=$

$ -\alpha_k \frac{\partial​ P_g}{\partial x_i}+\alpha_k \rho_k g_i+\frac{\partial}{\partial x_i}[-\alpha_k \rho_k < u'_{k,i} u'_{k,j}>+\Theta_{k,ij}]+I'_{k,i}+[U_{\sigma,i}-U_{k,i}]\Gamma_k+\delta_{kp}\sum_{n \neq  k}S_{kn,i}$

 où $I'_{k,i}$ est tel que:

$I'_{g,i}​=I'_{k,i}​=\alpha_{k}\rho_{k}\frac{V_{(r,k),i}}{\tau _{gk}^{F}}​$ avec $V_{(r,k),i}$ la ième

composante de la vitesse relative locale entre la kième phase particulaire et l'écoulement gazeux localement non perturbé, i.e $V_{(r,k),i}=U_{k,i}-U_{g,i}-V_{(d,k),i}$.

$\delta_{kp}$ est le symbole de Kronecker (où l’indice p définit une phase particulaire donnée) et $S_{kn,i}$ modélise les effets de collisions entre les différentes phases particulaires. Le terme $U_{\sigma,i} \Gamma_k$ représente le transfert de quantité de mouvement dû aux échanges de masses entre la phase gazeuse et les particules réactives.

On suppose que la vitesse de la phase gazeuse après la réaction (ici $C_{O_2}$) $U_{\sigma,i}$ est égale à la vitesse de la particule (ie le char C): $U_{\sigma,i}=U_{C,i}$.

Les scalaires

Les scalaires résolus ici sont $Y_{C_{O_2}}$ et $\chi_d$ avec $Y_{O_2}=1-Y_{C_{O_2}}$.

Sous l'hypothèse d'un mélange de gaz idéaux dans le réacteur, la densité du mélange $\rho_g$ et la fraction massique de $ Y_{CO{_2}}$ sont contrôlés :

$\rho_g=\frac{P_ref W_g}{R_g T_g}=\frac{P_ref}{R_g T_g}[\frac{Y_{O_2}}{W_{O_2}}+\frac{Y_{C{O_2}}}{W_{Y_{C{O_2}}}}]^{-1}$

L'équation de transport de $Y_{C_{O_2}}$

$\alpha_{g} \rho_{g} \frac{\partial Y_{CO_2}}{\partial t} + \alpha_{g} \rho_{g} U_{g,j} \frac{\partial Y_{CO_2}}{\partial x_{j}} =  \frac{\partial}{\partial x_{j}} (\alpha_{g} \rho_{g} \frac{\nu_{CO_{2}}}{\sigma_{g}} \frac{\partial Y_{CO_{2}}}{\partial x_{j}}) - Y_{CO_{2}} \Gamma_{g} + \psi_{Y_{CO_{2}}}$

$\nu_g^{t}$ et $\sigma_g$ sont respectivement la viscosité cinématique et le nombre de Schmidt turbulent. $W_g$ et $R_g$ sont respectivement la masse molaire du gaz et la constante  de masse du gaz, alors que $P_{ref}$ est la pression de référence au sein du réacteur.

Conservation du nombre de particules par unité de masse ($\chi_d$)

$\chi_d$ est le scalaire transporté pour modéliser la taille des particules. 

$\chi_d=\frac{\rho_{p_0}}{\rho_p}(\frac{d_{p_0}}{dp})^3$

$\rho_{p,0}$ et $dp_{p,0}$ représentent respectivement la masse volumique et le diamètre de référence du char. Ici, la masse volumique du char est constante donc :

$\chi_{d}=(\frac{d_{p_0}}{d_{p}})^3$

L'équation de transport du scalaire $\chi_d$ s'écrit :

$\alpha_{c}\rho_{c}\frac{\partial \chi_{d}}{\partial t}+\alpha_{c} \rho_{c} U_{c}\frac{\partial \chi_{d}}{\partial x_i} = - \frac{\partial}{\partial x_i}(\alpha_c \rho_c {<u'_{c} \chi_{d} >}_{c})+\Gamma_d -\chi_d \Gamma_c$

Or dans notre cas l'attrition n'est pas prise en compte donc $\Gamma_c$ qui représente la variation de $chi_d$ due à l'attrition et/ou l'agglomération des particules de char est nul. L'équation devient alors :

$\alpha_{c}\rho_{c}\frac{\partial \chi_{d}}{\partial t}+\alpha_{c} \rho_{c} U_{c}\frac{\partial \chi_{d}}{\partial x_i} = - \frac{\partial}{\partial x_i}(\alpha_c \rho_c {<u'_{c} \chi_{d} >}_{c}) -\chi_d \Gamma_c$

 

 

Implémentation des fichiers sources

Equation de combustion : $C_s + {O_2}_{(g)} \rightarrow {CO_2}_{(g)}$

Scalaires : $\chi_{d}$ et $Y_{CO_{2}}$

Conservation du nombre de particules par unité de masse $\chi_{d}$ :

$\alpha_{p} \rho_{p} \chi_{d}+\frac{\partial}{\partial x_i}\alpha_{p} \rho_{p} U_{p,i} \chi_{d}=\frac{\partial}{\partial x_i}(\alpha_{p} \rho_p D_{p,\chi_{d}} \frac{\partial \chi_d}{\partial x_{i}}) +\psi_{p}$

$\psi_{p,\chi_c}=0$ car il n'y a pas d'attrition ou d'agglomération

$\alpha_{p}\rho_{p}\frac{\partial \chi_d}{\partial t} +\alpha_{p}\rho_{p}U_{p,i}\frac{\partial \chi_d}{\partial x_i}= \frac{\partial}{\partial x_i}(\alpha_{d} \rho_{p} D_{p,\chi_d} \frac{\partial \chi_o}{\partial x_i}) - \chi_d \Gamma_c$

 

d'où $ \left\{
\begin{array}{ll}
TSA =\frac{-\Gamma_{c}\chi_d}{\alpha_c} \\
TSB = \frac{\partial TSA}{\partial \chi_d}=-\frac{\Gamma_c}{\alpha_c}\\
\end{array}
\right.$

On prendra $\chi_{d}$ égal à 1 à l'instant initial dans l'injecteur : en effet, les particules de char n'ont pas encore subi de réaction de combustion donc $d=d_0$.

Conservation de l'espèce gazeuse $CO_{2}$ :

$\Psi_{CO_2}=-\frac{W_{CO_2}}{W_c}\Gamma_c$

$\alpha_g \rho_{g}Y_{CO_2}+\frac{\partial}{\partial x_i} \alpha_{g} \rho_{g} U_{g,i} Y_{CO_2} = \frac {\partial}{\partial x_i}(\alpha_{g}\rho_{g}D_{g,Y_{CO_2}}\frac{\partial Y_{CO_2}}{\partial x_i}) +\Psi_{g,Y_{CO_2}}$

$\alpha_{g} \rho_{g} \frac{\partial Y_{CO_2}}{\partial t} + Y_{CO_2} \frac{\partial (\alpha_{g} \rho_{g})}{\partial t} + Y_{CO_2} \frac{\partial}{\partial x_i}(\alpha_{g}\rho_{g}U_{g,i})+\alpha_{g}\rho_{g}U_{g,i}\frac{\partial Y_{CO_2}}{\partial x_i} = \frac{\partial}{\partial x_i}(\alpha_{g}\rho_{g}D_{g,Y_{CO_2}}\frac{\partial Y_{CO_2}}{\partial x_i}) - \frac{W_{CO_2}}{W_c}\Gamma_c$

Or on a : $\frac{\partial}{\partial t} (\alpha_{g}\rho_{g})+\frac{\partial}{\partial x_i}(\alpha_{g} \rho_{g} U_{g,i})=\Gamma_g=-\Gamma_c$

d'où $\alpha_{g} \rho_{g} \frac{\partial Y_{CO_2}}{\partial t}alpha_{g} \rho_{g} U_{g,i} \frac{\partial}{\partial x_i} Y_{CO_2} - \Gamma_{c}Y_{CO_2}=\frac{\partial}{\partial x_i} (\alpha{g} \rho_{g} D_{g,Y_ {CO_2}} \frac{\partial}{\partial x_i} Y_{CO_2}) _\frac{W_{CO_2}}{W_c}\Gamma_c$

d'où $\alpha_{g} \rho_{g} \frac{\partial Y_{CO_2}}{\partial t}alpha_{g} \rho_{g} U_{g,i} \frac{\partial}{\partial x_i} Y_{CO_2}=\frac{\partial}{\partial x_i} (\alpha{g} \rho_{g} D_{g,Y_ {CO_2}} \frac{\partial}{\partial x_i} Y_{CO_2})+\Gamma_{c}(Y_{CO_2}-\frac{W_{CO_2}}{W_c})$

 

d'où $ \left\{
\begin{array}{ll}
TSA =\frac{-\Gamma_{c}}{\alpha_g}(-\frac{W_{CO_2}}{W_c}+Y_{CO_2}) \\
TSB = \frac{1}{\alpha_g}(\frac{W_{CO_2}}{W_c}\frac{\partial \Gamma_c}{\partial Y_{CO_2}}+\Gamma_c)=\frac{\Gamma_c}{\alpha_g}\\
\end{array}
\right.$

On prendra $Y_{CO_2}$ nul à l'instant initial.

Conservation de l'espèce gazeuse $CO_{2}$ :

Dans les fichiers sources de Neptune il nous faut coder la densité volumique de transfert de masse $\Gamma=-\frac{k_{c}SC_{O_2}}{V}$ avec $S=4\Pi R^2$ et $V=\frac{4}{3}\Pi R^3$ respectivement la surface externe et le volume de la particule sphérique de char.

$C_{02}$ est calculée dans la partie $\href{http://hmf.enseeiht.fr/travaux/bei/beiep/content/g12/modelisation-equations}{\textbf{Implémentation des équations}}$, k vaut $10^-{3} m/s$ et R=d/2 avec d=0,00344 m le diamètre initial d'une particule de char.