1 | \chapter{Le code}
2 |
3 | \section{Introduction}
4 | \code (Parallel RadIation Solver with Spectral integration on Multicomponent mediA) % (Discrete Ordinates Method Applied with Spectral Integration on Unstructured Meshes)
5 | est un code de calcul de transfert radiatif développé en Fortran. Il permet de résoudre l'équation de transfert radiatif pour calculer la luminance dans une géométrie complexe à partir de solutions issues de calcul de CFD \footnote{Computational Fluid Dynamics} : les champs de concentration et de température. \newline
6 |
7 | Le calcul de la luminance $I$ permet d'obtenir :
8 | \begin{itemize}
9 | \item $G$, le rayonnement incident ;
10 | \item $H$, l'éclairement surfacique ;
11 | \item $S_r$, le terme source radiatif ;
12 | \item $Q$, le flux aux parois.
13 | \end{itemize}
14 |
15 | \section{Historique}
16 |
17 | Le code \code s'appuie sur une longue collaboration entre Mouna El Hafi du Centre Énergétique et Environnement de l'École des Mines d'Albi-Carmaux \footnote{\urlcolor{http://www.mines-albi.fr}} et Bénédicte Cuénot au CERFACS (Toulouse) \footnote{Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique : \urlcolor{http://www.cerfacs.fr}}, voir Fig.~\ref{Historique}. Dans sa version actuelle, le code \code a été écrit par Jorge Amaya (thèse au CERFACS 2006-2010) et développé conjointement avec Damien Poitou (thèse École des Mines d'Albi/CERFACS 2006-2009).
18 |
19 | \figScale{Historique}{Historique du code \code.}
20 |
21 | Initialement un premier code a été développé durant la thèse de David Joseph\cite{phdJoseph} (2001-2004)~: \mbox{DOMASIUM} (Discrete Ordinates Method Applied with Spectral Integration on Unstructured Meshes). Ce code a permis de mettre en place un code radiatif adapté aux maillages non structurés utilisés en CFD dans le but de réaliser des calculs de transferts radiatifs dans des configurations de combustion. Ce code est basée sur la méthode des ordonnées discrètes dans un soucis d'efficacité de temps de calcul, cette méthode offrant un bon compromis entre précision et coût de calcul. Cette méthode est utilisé conjointement avec un modèle spectral en {\it k}-distributions sur des bandes étroites du spectre infrarouge. Ce modèle permet une description précise des propriétés d'absorption du milieu composé de gaz produits de combustion (\HdeuxO, \COdeux,\CO) et éventuellement de suies. La première version du code est principalement séquentielle bien qu'une première tentative de parallélisation ait été réalisée par D. Joseph.
22 |
23 | À partir de la version originale de \mbox{DOMASIUM}, différents développements on été réalisé dans le but d'optimiser le code dans le cadre du projet CORAYL (PROJET ANR-05-CICG-012). Ce projet est une collaboration entre l'École des Mines d'Albi, le CERFACS et le laboratoire EM2C \footnote{Laboratoire d'Énergétique
24 | Moléculaire et Macroscopique, Combustion} %\urlcolor{http://www.em2c.ecp.fr/}
25 | de l'École Centrale Paris. David Joseph, Patrice Perez et Damien Poitou ont réalisé dans ce projet un étude sur la mise en place de modèles spectraux rapides pour réduire les temps de calculs radiatifs. Durant la thèse de Rogerio Dos Santos une version optimisée et parallélisée a été développée à l'EM2C~: \mbox{MPI-DOMASIUM}. Cette version a permis de réduire le temps de calcul du code et de démontrer la faisabilité d'un couplage rayonnement-combustion turbulente avec le code \mbox{AVBP}\footnote{\urlcolor{http://www.cerfacs.fr/4-26334-The-AVBP-code.php}}.
26 |
27 | Dans le but de synthétiser les différentes versions précédentes pour les intégrer dans un code intrinsèquement parallèle, un nouveau code a entièrement été réécrit~: \code. Ce code a été optimisé en mémoire afin de rechercher les meilleures performance de calcul dans un contexte de calculs parallèles. Ce travail de réécriture du code a été principalement réalisé par J. Amaya. De nouveaux modèles spectraux ont également été intégrés durant la thèse de D. Poitou\cite{th_poitou} ce qui a permis d'atteindre des rapports de temps de calcul \mbox{PRISSMA/AVBP} inférieurs à 1 montrant la faisabilité de couplage dans des configurations industrielles. Durant le développement de \code la stabilité et l'utilisation du code ont été grandement améliorés notamment par la création d'un ensemble d'outils de pré et post-processing.
28 |
29 | Nous allons rappeler dans la partie suivante les bases théoriques sur lesquelles s'appuient le code.
30 |
31 |
32 | \section{Les bases théoriques}
33 |
34 | \code \cite{IJTS,Stanford} has been designed to simulate the radiative heat transfer in coupled simulations with flow dynamics, involving unstructured grids. In the following and for sake of clarity, the intensities and radiative properties are expressed for a single wavenumber (monochromatic case) but the formulation can be easily extended to a full spectrum. \\
35 | Discrete Ordinates Method have been introduced first by \cite{Cha50} and have been widely used in radiative transfer applications. Considering an absorbing-emitting and non-scattering gray medium, the
36 | variation of the radiative intensity $I(\textbf{s})$ along a line of sight can be written as:
37 | \begin{equation}
38 | \frac{dI(\textbf{s})}{ds} = \kappa I_{b} - \kappa I(\textbf{s})
39 | \label{TheBIGONE}
40 | \end{equation}
41 | where $I(\textbf{s})$ is the radiative intensity along the directional coordinate $\textbf{s}$, $I_b$ is the blackbody radiative intensity, and $\kappa$ is the absorption coefficient. Boundary conditions for diffuse surfaces are taken from the relation giving the intensity leaving the wall $I_w$ as a function of the blackbody intensity of the wall $I_{b,w}$ and of the incident radiative intensity:
42 | \begin{equation}
43 | I_{w}(\textbf{s})=\epsilon_{w}I_{b,w}+\frac{\rho_{w}}{\pi}\int_{\textbf{n}.\textbf{s'}
44 | <0}I_{w}(\textbf{s'})|\textbf{n}.\textbf{s'}|d\Omega'
45 | \label{Boundings}
46 | \end{equation}
47 | where $\epsilon_{w}$ is the wall emissivity, $\rho_{w}$ the wall reflectivity, $\textbf{n}$ the unit vector normal to the wall and $\textbf{s'}$ the direction of propagation of the incident radiation confined within a solid angle $d\Omega'$.
48 |
49 | \subsection{Angular discretization}
50 |
51 | In the DOM, the calculation of a radiative source term at a given point is based on the discretization of the Radiative Transfer Equation (Eq.\ref{TheBIGONE}) according to a chosen number $N_{dir}$ of discrete directions $\textbf{s}_i (\mu_i, \eta_i, \xi_i)$, associated to the corresponding weights $w_i$,
52 | contained in the solid angle $4\pi$, and where ($\mu_i, \eta_i, \xi_i$) are directional cosines. Different angular discretizations may be used. In a recent study, Koch and Becker \cite{Koc03} compare the efficiency of several types of angular quadratures. They recommend the $LC_{11}$ quadrature for its better accuracy. However calculations performed with the $S_4$ quadrature satisfy a good compromise between accuracy and rapidity as shown in \cite{Stanford}, and may also be used.
53 |
54 | \subsection{Spatial discretization for hybrid grids}
55 |
56 | The RTE (Eq.\ref{TheBIGONE}) is solved for every discrete direction $\textbf{s}_i$ using a finite volume approach. The integration of the RTE over the volume $V$ of an element limited by a surface $\Sigma$, and the application of the divergence theorem yields:
57 | \begin{equation}
58 | \int_{\Sigma} I(\textbf{s}_i) . \textbf{s}_i . \textbf{n} d\Sigma = \int_{V} (\kappa I_{b} - \kappa I(\textbf{s}_i)) dV
59 | \label{EquaIntegrated1}
60 | \end{equation}
61 | The domain is discretized in three-dimensional control volumes $V$. It is assumed that $I_b$ and $I(\textbf{s}_i)$ are constant over the volume $V$ and that the intensities $I_j$ at the faces are constant over each face. Considering that $I_j$ is the averaged intensity over the $j^{th}$ face, associated with the center of the corresponding face, that $I_{b,P}$ and $I_P$ are the averaged intensities over the volume $V$, associated with the center of the cell, and assuming plane faces and vertices linked by straight lines, Eq.(\ref{EquaIntegrated1}) can be discretized as follows :
62 | \begin{equation}
63 | \sum_{j=1}^{N_{face}} I_j(\textbf{s}_i).(\textbf{s}_i.\textbf{n}_j)A_j = \kappa V (I_{b,P}-I_P(\textbf{s}_i))
64 | \label{EquaDiscrete1}
65 | \end{equation}
66 | where $\textbf{n}_j$ is the outer unit normal vector of the surface $j$.\\
67 | The scalar product of the $i^{th}$ discrete
68 | direction vector with the normal vector of the $j^{th}$ face of the considered
69 | cell is defined by $D_{ij}$:
70 | \begin{equation}
71 | D_{ij} = \textbf{s}_i . \textbf{n}_j = \mu_i n_{xj} + \eta_i n_{yj} + \xi_i n_{zj}
72 | \end{equation}
73 |
74 | The discretization of the boundary condition (Eq.(\ref{Boundings})) is straightforward:
75 | \begin{equation}
76 | I_{w} = \epsilon_w I_{b,w} + \frac{1-\epsilon_w}{\pi}
77 | \sum_{\textbf{n}.{\textbf{s}_i <0}} w_{i}I(\textbf{s}_i)\mid
78 | \textbf{n}.\textbf{s}_i\mid
79 | \label{BOUND}
80 | \end{equation}
81 | For each cell, the incident radiation $G$ is evaluated as follows:
82 | \begin{equation}
83 | G = \int_{4\pi} I(\textbf{s}) d \Omega \, \simeq
84 | \sum_{i=1}^{N_{dir}} w_i I(\textbf{s}_i)
85 | \label{Form_G}
86 | \end{equation}
87 | and the incident heat flux $H_w$ at the wall surfaces is :
88 | \begin{equation}
89 | H_w = \int_{\textbf{n}.\textbf{s} <0} I(\textbf{s}) \mid \textbf{n}.\textbf{s} \mid d \Omega \, \simeq \sum_{\textbf{n}.\textbf{s}_i <0} w_{i}I_{i} \mid \textbf{n}.\textbf{s}_i \mid
90 | \label{Form_H}
91 | \end{equation}
92 | For a gray medium, the radiative source term $S_r$ is given by:
93 | \begin{equation}
94 | S_r = \nabla.Q_r \, = \kappa (4\pi I_b - G)
95 | \label{Form_Sr}
96 | \end{equation}
97 | where $Q_r$ is the radiative heat flux, and the radiative net heat flux at the wall is:
98 | \begin{equation}
99 | Q_w = \epsilon \pi I_{b,w} - H
100 | \label{Form_Qw}
101 | \end{equation}
102 | For the evaluation of the radiative intensity $I(\textbf{s}_i)$ in Eq. (\ref{BOUND}) to (\ref{Form_Qw}) Str\"ohle et al. \cite{Str01} proposed a simple spatial differencing scheme based on the mean flux scheme that proved to be very efficient in the case of hybrid grids. This scheme relies on the following formulation:
103 | \begin{equation}
104 | I_P = \alpha \overline{I_{out}} + (1-\alpha) \overline{I_{in}}
105 | \label{Equapatialscheme}
106 | \end{equation}
107 | where $\overline{I_{in}}$ and $\overline{I_{out}}$ are respectively the intensities averaged over the entering and the exit faces of the considered cell.
108 | $\alpha$ is a weighting number between $0$ and $1$. Substituting $\overline{I_{out}}$ from Eq.(\ref{Equapatialscheme}) into Eq.(\ref{EquaDiscrete1}) yields (for more details see \cite{IJTS}):
109 | \begin{equation}
110 | I_P = \frac{\alpha V \kappa I_b - \displaystyle
111 | \sum_{\substack{j\\D_{ij}<0}} D_{ij} A_j I_j}{\alpha\kappa V +
112 | \displaystyle \sum_{\substack{j\\D_{ij}>0}}
113 | D_{ij}A_j}
114 | \label{Equaprincipal}
115 | \end{equation}
116 | The case $\alpha=1$ corresponds to the Step scheme used by Liu et al. \cite{Liu00c}. The case $\alpha=0.5$ is called Diamond Mean Flux Scheme (\textbf{DMFS}) which is formally more accurate than the Step scheme. After calculation of $I_P$ from Eq.(\ref{Equaprincipal}), the radiation intensities at cell faces such that $D_{ij}>0$ are set equal to $\overline{I_{out}}$, obtained from Eq.(\ref{Equapatialscheme}). For a given discrete direction, each face of each cell is placed either upstream or downstream of the considered cell center (a face parallel to the considered discrete direction plays no role). The control volumes are treated following a sweeping order such as the radiation intensities at upstream cell faces are known. This order depends on the discrete direction under
117 | consideration. An algorithm for the optimization of the sweeping order has been implemented \cite{IJTS}. Note that this sweeping order is stored for each
118 | discrete direction, and only depends on the chosen grid and the angular quadrature, i.e. it is independent on the physical parameters or the flow and may be calculated only once, prior to the full computation.
119 |
120 | \subsection{Spectral gas properties}
121 |
122 | The absorption coefficient $\kappa$ of the combustion products is highly dependent on the wavenumber $\nu$ as shown by line spectra of radiative gases (H$_2$O, CO$_2$ and CO). To take into account this spectral dependency, the absorption coefficient of each species is here represented by the \textit{SNB-ck} model \cite{Sou97,Liu00b,Liu01}. For the gas mixture composed of different species, the same model is used, building data according to the mixing model exposed by Liu \cite{Liu01}. The radiative solutions are obtained by computing $N_{bands} \times N_{quad}$ independent calculations where $N_{bands}=367$ is the number of narrow bands of spectral width $\Delta \eta = 25 \, cm^{-1}$, describing the spectral properties in the range $150-9300 \, cm^{-1}$, and $N_{quad}=5$ is the number of the Gauss-Legendre quadrature points used for the spectral integration over each narrow band. For non-gray media, introducing spectral dependencies in Eq. (\ref{Form_Sr}), gives for the source term :
123 | \begin{equation}
124 | S_{r,DOM}=\sum_{i=1}^{N_{band}}\sum_{j=1}^{N_{quad}} \Delta \nu_i w_{ij}
125 | \kappa_{ij}(4\pi \overline{I}_{b,ij}-G_{ij})
126 | \end{equation}
127 | where $G_{ij}$ is obtained from Eq.(\ref{Form_G}).\\
128 | The computational efficiency of such a model is strongly linked to the number of bands $N_{bands}$, that has to be optimized depending on the studied case.
129 |
130 | Other models called global models have been developed : WSGG (Weighted Sum of Gray Gas) \cite{Sou94} and FS-SNBcK (Full Spectrum SNBcK) \cite{Poi09}.
131 |
132 |
133 | \subsubsection{Modèle SNB-CK}
134 | Le modèle SNB-CK est tout d'abord un modèle SNB et il considère donc un découpage du spectre des gaz rayonnants par des bandes dans lesquelles on fait l'hypothèse que la luminance du corps noir peut \^etre considérée comme constante (hypothèse de bande étroite). Dans chaque bande étroite, les propriétés spectrales d'un gaz sont données représentées par deux paramètres $(\overline{\kappa},\phi)$ tabulés pour des températures allant de 300 K à 2900 K par pas de 200 K, et pour des longueurs d'onde allant de $150 \, cm^{-1}$ à $9300 \, cm^{-1}$. Les gaz considérés sont $H_2O$, $CO_2$ et $CO$ et les données spectrales du modèles SNB proviennent de \cite{Sou97}. Ces données sont valables pour un modèle SNB basé sur un modèle de Malkmus pour lequel les intensités de raies d'un gaz dans chaque bande étroite sont distribuées suivant une loi inverse exponentielle \cite{Mal67,Duf99}. Avec de tels modèles, le transmittivité moyenne d'une colonne de gaz de longueur $l$ homogène et isotherme pour une bande étroite de largeur $\Delta \nu$ est donnée par:
135 | \begin{equation}
136 | \overline{\tau}_{\Delta \nu}(l)=exp \left[-\frac{\phi}{\pi} \left( \left( 1+ \frac{2 \pi \overline{\kappa} l}{\phi} \right)^{1/2} -1\right) \right]
137 | \end{equation}
138 | Le calcul des grandeurs telles que le terme source radiatif ou le flux aux parois suppose une intégration sur le spectre en général et donc sur chaque bande étroite. Si on note $\kappa_{\eta}$ le coefficient d'absorption au nombre d'onde $\nu$ dans une bande étroite, le calcul d'une grandeur $H$ intégrée sur la bande revient à calculer l'intégrale suivante:
139 | \begin{equation}
140 | \overline{H}_{\Delta \nu}= \frac{1}{\Delta \nu} \int_{\Delta \nu} H(\kappa_{\nu}) d\nu
141 | \label{lintegrale}
142 | \end{equation}
143 | Par conséquent, seule la valeur intégrée sur une bande étant importante, on a toute liberté de réorganiser les raies présentes dans une bande sans altérer le résultat et les calculs dans une bande. Autrement dit, seule la fonction de distribution des coefficients d'absorption dans le bande $f(\kappa)$ est importante, et pas la fréquence à laquelle ces coefficients d'absorption sont calculés.
144 |
145 |
146 | La fonction de distribution des coefficients d'absorptions pour le modèle de Malkmus a une expression analytique \cite{Duf99} et est obtenue par~:
147 | \begin{equation}
148 | f_{\Delta \nu}(\kappa)=\sqrt{\frac{\phi \overline{\kappa}}{2 \pi \kappa^3}} exp \Bigg[-\frac{\phi}{2}\frac{(\kappa - \overline{\kappa})^2}{\kappa\overline{\kappa} } \Bigg]
149 | \label{fdistrib}
150 | \end{equation}
151 |
152 | On utilise alors l'expression de la fonction cumulée de la distribution des coefficients d'absorption comme variable de description spectrale ($g(\kappa)=\int_0^\kappa f(\kappa')d\kappa'$), qui a comme expression analytique pour un modèle de Malkmus~:
153 | \begin{equation}
154 | g_{\Delta \nu}(\kappa)=\Gamma \Bigg[ -\sqrt{\frac{\phi}{\kappa / \overline{\kappa}}} \Bigg(1-\frac{\kappa}{\overline{\kappa}} \Bigg) \Bigg] + exp \big(2 \phi \big)\Gamma \Bigg[ -\sqrt{\frac{\phi}{\kappa / \overline{\kappa}}} \Bigg(1+\frac{\kappa}{\overline{\kappa}} \Bigg) \Bigg]
155 | \label{cumulat}
156 | \end{equation}
157 |
158 | Par cette opération, l'intégrale précédente (Eq. \ref{lintegrale}) se réécrit:
159 | \begin{equation}
160 | \overline{H}_{\Delta \nu}= \frac{1}{\Delta \nu} \int_{\Delta \nu} H(\kappa_{\nu}) d\nu = \frac{1}{\Delta \nu} \int_{0}^{1} H(\kappa (g)) dg
161 | \end{equation}
162 | où $g$ est la cumulative de la fonction de distribution du coefficient d'absorption. L'intégration d'une grandeur sur une bande étroite est effectivement réalisée en utilisant une quadrature sur la fonction $g$. La quadrature utilisée par défaut dans le code \code est une quadrature de Gauss-Legendre à 5 points \cite{Liu00b}, de sorte que l'intégrale précédente devient:
163 | \begin{equation}
164 | \overline{H}_{\Delta \nu}= \sum_{j=1}^{N_{quad}} w_j H(\kappa (g_j))
165 | \end{equation}
166 | où $N_{quad}$ est le nombre de points de quadratures sur une bande étroite pour représenter la fonction g, et $w_j$ est le poids associé au $j^\textrm{ième}$ point de la quadrature.\\
167 |
168 | A ce stade, pour traiter des milieux non homogènes et non isothermes, on fait une hypothèse supplémentaire afin de prendre en compte la déformation du spectre le long d'un chemin optique. Cette hypothèse consiste à considérer que toutes les raies d'une bande étroite se déforment de la m\^eme manière en fonction de la température. Autrement dit, on peut considérer que la quantité $\frac{\kappa}{\overline{\kappa}}$ est constante quel que soit la fréquence dans la bande étroite. Cela revient également à considérer que la fonction $g$ est constante le long d'un chemin optique. C'est l'hypothèse des K corrélés ou hypothèse CK.
169 |
170 | Par ailleurs, le traitement du mélange de gaz est obtenu par le m\^eme modèle en suivant la méthode de Liu \cite{Liu01} basée sur la limite optiquement mince et qui donne les paramètres d'un mélange de gaz en fonction des paramètres de ses constituants:
171 | \begin{equation}
172 | \overline{\kappa}_{m\'elange}=\sum_{i=1}^{N_{gaz}} \overline{\kappa}_i
173 | \end{equation}
174 | et
175 | \begin{equation}
176 | \frac{\overline{\kappa}^2_{m\'elange}}{\phi_{m\'elange}}=\sum_{i=1}^{N_{gaz}} \frac{\overline{\kappa}^2_i}{\phi_i}
177 | \end{equation}
178 |
179 | Lors d'un calcul de terme source, les calculs radiatifs sont effectués $N_{bandes} \times N_{quad}$ fois de manière indépendantes, où $N_{bandes}=367$ est le nombre de bandes étroites, de largeur $\Delta \eta = 25 \, cm^{-1}$ considérées. Ces bandes permettent de considérer une région du spectre infrarouge allant de $150 \, cm^{-1}$ à $9300 \, cm^{-1}$. $N_{quad}=5$ représente le nombre de points de quadratures pour chaque bande étroite. L'équation du terme source radiatif calculé est donnée par~:
180 | \begin{equation}
181 | S_{r,MOD}=-\sum_{i=1}^{N_{bandes}}\sum_{j=1}^{N_{quad}} \Delta \nu_i w_{ij}
182 | \kappa_{ij} \Bigg(4\pi \overline{I}_{b,ij}- \sum_{k=1}^{N_{dir}}w_{k}{I}_{ij}(\textbf{s}_k) \Bigg)
183 | \end{equation}
184 |
185 | Le modèle SNB-CK étant le modèle de propriétés radiatives des gaz disponible le plus précis dans \code, les résultats produits avec ce modèle serviront de résultats de référence pour les modèles FSCK et WSGG. Le code \code utilisé avec le modèle SNB-CK a fait l'objet de benchmarks avec d'autres codes et l'on pourra trouver des points de comparaison dans \cite{Stanford}.
186 |
187 |
188 | \subsubsection{Modèle FS-SNBCK}
189 | En partant du modèle SNB-CK, Liu et al. \cite{Liu99} proposent d'opérer des regroupements de bandes étroites voisines et d'effectuer ensuite les calculs sur ces bandes élargies. Alors, pour un regroupement de $M$ bandes étroites de m\^eme largeur, la transmittivité moyenne est donnée par:
190 | \begin{equation}
191 | \overline{\tau}_{wb}=\frac{1}{M}\sum_{i=1}^{M} \overline{\tau}_i
192 | \end{equation}
193 | où $\overline{\tau}_i$ est la transmittivité moyenne de la ième bande étroite dans le regroupement des M bandes. La transmittivité moyenne sur la bande large est alors simplement la moyenne des transmittivités moyennes des bandes étroites qui la composent.
194 | Si on considère maintenant une formulation en coefficients d'absorption, la fonction de distribution de la bande large est donnée par:
195 | \begin{equation}
196 | f_{wb}(\kappa)=\frac{1}{M}\sum_{i=1}^{M} f_i(\kappa)
197 | \end{equation}
198 | où $f_i(\kappa)$ est la fonction de distribution des coefficients d'absorption dans la ième bande étroite. La fonction cumulée peut \^etre obtenue de la m\^eme façon par:
199 | \begin{equation}
200 | g_{wb}(\kappa)=\frac{1}{M}\sum_{i=1}^{M} g_i(\kappa)
201 | \end{equation}
202 | Ce principe du regroupement de bandes introduit des erreurs dans l'hypothèse de bande étroite. En effet, plus on ajoute de bandes étroits voisines, et moins on peut considérer que la luminance du corps noir dans la bande est constante. Liu et al. \cite{Liu99} soulignent cette source d'erreurs et recommandent de se limiter dans le nombre de bandes regroupées. Cependant, il est possible de tenir compte de la variation de la luminance du corps noir dans le différentes bandes regroupées par une pondération sur chaque bande étroite \cite{Liu04a}. Il est alors possible de regrouper toutes les bandes pour obtenir le modèle FSCK. Dans ce modèle, toutes les bandes étroites sont regroupés et les fonctions $f$ et $g$ obtenues de ce regroupement à une température donnée sont exprimées par:
203 | \begin{equation}
204 | f_{wb}(\kappa)=\frac{1}{\sigma T^4}\sum_{i=1}^{367} I_{b\nu,i} \, f_i(\kappa) \Delta \nu_i
205 | \end{equation}
206 | où $\sigma$ est la constante de Stefan-Boltzmann, $I_{b\nu,i}$ est la luminance du corps noir pour la bande étroite $i$ à la fréquence $\nu$. La m\^eme pondération est affectée à la fonction cumulée:\begin{equation}
207 | g_{wb}(\kappa)=\frac{1}{\sigma T^4}\sum_{i=1}^{367} I_{b\nu,i} \, g_i(\kappa) \Delta \nu_i
208 | \label{gwb}
209 | \end{equation}
210 | L'utilisation pratique de ce modèle se fait en plusieurs étapes nécessitant l'inversion de la fonction g. En effet, pour réaliser l'intégrale de l'équation \ref{lintegrale}, il faut conna\^itre la valeur de la fonction $H$ pour le coefficient d'absorption $\kappa$ correspondant au point de quadrature $g_i$. On dispose au départ de la quadrature qui donne les valeurs de $g_j$ et les poids $w_j$ associés pour le spectre regroupé. Il faut alors chercher pour chaque valeur de $g_j$ la valeur $\kappa_j$ telle que l'équation \ref{gwb} soit respectée. Pour une valeur de $g_j$, cette opération se fait par dichotomie. Pour chaque valeur de $\kappa$ testée, on calcule donc la valeur de $g(\kappa)$ pour chaque bande, puis on calcule la valeur de $g_{wb}$ en utilisant l'équation \ref{gwb}. La valeur maximale de départ prise pour $\kappa_j$ est de 10 fois la valeur maximale du paramètre $\overline{\kappa}$ pour l'ensemble des bandes étroites. Un critère d'arr\^et de la dichotomie à 50 itérations a été introduit afin d'éviter des phénomènes de boucle infinie.\\
211 | L'inversion par dichotomie peut \^etre consommatrice en temps de calcul et il sera donc souhaitable de rechercher une tabulation pour éviter cette opération . Toutefois, une fois cette opération effectuée, l'intégration spectrale se fait alors très rapidement, puisqu'il n'y a plus que $N_{quad}$ calculs spectraux à effectuer.
212 |
213 |
214 |
215 | \subsubsection{Modèle WSGG}
216 | Dans le modèle de somme pondérée de gaz gris, l'émissivité totale $\epsilon_T$ d'une colonne de gaz de longueur l à la température T est donnée par:
217 | \begin{equation}
218 | \epsilon_T(T,l,p_a)=\sum_{k=1}^n a_k(T) \left[ 1- exp \left( -\kappa_k p_a l \right) \right]
219 | \label{wsgg}
220 | \end{equation}
221 | où $\kappa_k$ est le coefficient d'absorption du $k^\textrm{ième}$ gaz gris, $p_a$ est la pression partielle des espèces absorbantes, et $a_k(T)$ est le facteur de pondération du $k^\textrm{ième}$ gaz. La pondération $a_k(T)$ peut \^etre interprétée comme la fraction du rayonnement du corps noir dont le coefficient d'absorption est proche de $\kappa_k$. Dans le cas d'un milieu homogène et isotherme, la résolution de l'équation de transfert radiatif peut \^etre remplacée par un système de n équations de transfert pour n gaz gris:
222 | \begin{eqnarray}
223 | I &=& \sum_{k=1}^n I^k \\
224 | \frac{dI^k}{ds} &=& \kappa_k p_a \left( a_k I_b - I^k \right)
225 | \end{eqnarray}
226 | Ce m\^eme système est couramment utilisé pour des milieux inhomogènes et anisothermes. Le facteur de pndération $w_k$ est alors pris à la température locale.\\
227 |
228 | Les paramètres du modèle de gaz gris qui sont le nombre de gaz gris $n$, les coefficinets d'absorption $\kappa_k$ et les pondérations $a_k(T)$ ont été calculés par de nombreux auteurs pour la vapeur d'eau seule, pour des mélanges $CO_2-H_2O$, ou m\^em en prenant en compte le rayonnement de la suie. dans le cas de mélanges $CO_2-H_2O$, la tabulation doit \^etre faire pour des valeurs données du rapport entre la pression partielle d'eau ($p_w$) et celle de $CO_2$ ($p_c$). Ces rapport de pressions $\frac{p_w}{p_c}$ sont couramment choisis égaux à 2, 1, ou 0.5. Le modèle choisi pour \code est un modèle à 3 gaz gris plus un gaz transparent, dont les paramètres ont été calculés pour des mélange $CO_2-H_2O$ à pression atmosphérique et avec un rapport $\frac{p_w}{p_c}=2$ correspondant par exemple à une situation de combustion méthane-air \cite{Sou94}. Ces paramètres ont été calculés à partir de données spectrales SNB couvrant une région allant de $150 cm^{-1}$ à $7000 cm^{-1}$. Elles ne couvrent donc pas la totalites du spectre utilisé dans le modèle SNB-CK, en particulier pour les grands nombres d'onde, mais prennent en compte les bandes les plus importantes pour $H_2O$ et $CO_2$.\\
229 |
230 | Le modèle choisi comprenant 3 gaz gris, l'équation \ref{wsgg} se reécrit:
231 | \begin{equation}
232 | \epsilon_T(T,l,p_a)=\sum_{k=1}^n a_k(T) \left[ 1- exp \left( -\kappa_k p_w l \right) \right]
233 | \label{wsgg3}
234 | \end{equation}
235 | Et comme le quatrième gaz est un gaz transparent, la relation suivante lie les coefficients $a_k$:
236 | \begin{equation}
237 | a_4=1- \sum_{k=1}^3 a_k
238 | \end{equation}
239 | LOrs du calcul de ces paramètres, les coefficients $\kappa_k$ et $a_k$ ont été ajustés pour une température $T=1700K$. Pour les autres températures couvrant une plage de $300K$ à $2500K$, les valeurs de $K_k$ sont fixées à celles obtenues pour $T=1700K$ et les coefficients $a_k$ sont ajustés pour chaque température par un polyn\^ome à 5 coefficients:
240 | \begin{equation}
241 | a_k(T)=\sum_{j=0}^5 \alpha_{kj} T_j
242 | \end{equation}
243 | La valeurs des paramètres $a_k(T)$ sont repésentées sur la figure \ref{coef_wsgg_300_2500}.
244 |
245 | \figScaleX{0.7}{coef_wsgg_300_2500}{Coefficients pour le modèle WSGG}