Prise en compte de l’hétérogénéité des propriétés radiatives

Prise en compte de l’hétérogénéité des propriétés radiatives

Les éléments essentiels de la physique du rayonnement ainsi que des Méthodes de Monte-Carlo appliquées à ce champ d’application ayant été présentés ; ce chapitre a pour but d’aborder les premiers travaux de cette thèse : la gestion des hétérogénéités des propriétés optiques du milieu participant. Ces travaux constituent la suite logique d’une dynamique collective initiée depuis plusieurs années au sein de l’équipe STARWest [Terrée, 2011, Piaud, 2010, Eymet, 2011a]. Les méthodes de Monte-Carlo en milieu semi-transparent sont aujourd’hui bien maîtrisées [Farmer et Howell, 1998, Siegel et al., 2011, Modest, 2013]. C’est en particulier le cas lorsque les propriétés optiques du milieu (coefficients d’absorption, de diffusion, etc.) sont homogènes. Mais au-delà de cas académiques simples, dès qu’il s’agit d’étudier ou de simuler du transfert radiatif en configurations réelles dans des milieux participants, en particulier dans les gaz, la prise en compte et la gestion des hétérogénéités apparaissent comme primordiales. Ces tâches deviennent cependant rapidement délicates et exigeantes lorsqu’il s’agit d’employer les méthodes de Monte-Carlo comme outil de simulation, tout en souhaitant garder le caractère exact qu’elles offrent. Si l’on se concentre sur les domaines d’application pour lesquels le transfert radiatif en milieu gazeux occupe une place importante, les hétérogénéités sont omniprésentes. En effet, que ce soit dans les systèmes de combustion, dans les atmosphères terrestre ou exoplanétaires, on rencontre généralement de fortes hétérogénéités de température, de concentrations d’espèces ou de pression, menant à d’importantes variations des propriétés optiques du milieu observé. Une prise en compte rigoureuse de ces disparités spatiales est alors nécessaire pour mener à bien l’évaluation des observables d’intérêt. Cependant les difficultés qu’elles impliquent nécessitent généralement de recourir à des hypothèses simplificatrices ou à des méthodes entraînant des erreurs non maîtrisées. Dans une volonté de préserver le caractère exact dont bénéficient les méthodes de Monte-Carlo et l’analyse statistique associée, un des principaux objets de cette thèse a été de proposer une méthode prenant en compte la complexité de ces hétérogénéités, sans faire appel à une quelconque approximation. La solution retenue : les algorithmes à collisions nulles fera l’objet de ce chapitre. Cette méthode, jusqu’alors absente de la littérature du rayonnement thermique, mais très employée dans d’autres disciplines de la physique du transport (neutronique et physique des plasmas) y sera décrite et adaptée à des problématiques radiatives. Tout au long de ce chapitre, la dimension spectrale sera ignorée. Les problèmes seront donc ramenés à des cas monochromatiques (l’intégration spectrale fera l’objet du chapitre Chap. 5). Dans un premier temps, les difficultés relatives à la gestion des hétérogénéités par les méthodes de Monte-Carlo, ainsi que les techniques couramment utilisées pour y répondre seront présentées. Un bref état de l’art de la littérature relative aux algorithmes à collisions nulles sera ensuite dressé. Puis, ces méthodes seront introduites et élargies à l’étude du transfert radiatif. Une extension du domaine de validité de ces méthodes sera également proposée. L’approche retenue, reposant sur l’introduction d’un coefficient virtuel de collision nulle kn,η, sera ensuite éprouvée face à un cas d’étude plus complexe : l’estimation d’un bilan radiatif dans un milieu tridimensionnel, absorbant, émettant et diffusant, entouré par des parois réfléchissantes. Ces simulations donneront lieu à une étude paramétrique permettant d’évaluer les influences du coefficient de collision nulle, des choix méthodologiques et de différentes propriétés optiques sur le comportement de l’algorithme de Monte-Carlo retenu. Dans la quatrième section, nous montrerons comment le caractère de solution de référence des algorithmes à collisions nulles peut être mis à profit dans la validation d’un code de calcul radiatif en géométrie complexe (la configuration retenue sera celle d’une chambre de combustion). La rédaction de ce chapitre s’appuie sur les deux publications [Galtier et al., 2013], [Eymet, 2011b] (données en Annexe D et Annexe E) qui ont fait suite aux travaux présentés dans le présent manuscrit.

Problèmes liés à l’hétérogénéité des propriétés optiques du milieu participant

Incapacité d’échantillonner analytiquement des libres parcours

 Pour illustrer les difficultés rencontrées lorsque les propriétés optiques du milieu ne sont pas uniformes, reprenons la configuration de la Sec. 3.3.1 : le calcul de la luminance Lη(x0, u0) dans un milieu infini purement absorbant/émettant. Ce cas d’étude est suffisant pour aborder le problème relatif aux hétérogénéités. Le passage à un cas diffusif ou à une géométrie fermée n’entraînera aucune difficulté supplémentaire – si ce n’est de formalisme. Dans ces considérations,

Alternatives couramment proposées

 Pour répondre à cette limite, deux principales approches sont couramment employées dans la communauté du rayonnement thermique : la discrétisation des propriétés optiques du milieu et l’inversion numérique des épaisseurs optiques. Discrétisation du milieu La plus commune d’entre-elles consiste à discrétiser spatialement le volume d’intérêt et à considérer les propriétés du milieu comme uniformes à l’intérieur de chaque maille (voir Fig. 4.1). De ce fait, il devient possible d’échantillonner de façon analytique les libres parcours d’extinction, puisque les propriétés optiques sont constantes par morceaux le long du chemin optique. Une telle méthode possède cependant quelques limites. En effet, en discrétisant les propriétés optiques du milieu, le modèle physique est modifié. Les résultats de simulation dépendent alors du choix de maillage et les erreurs numériques causées par ce choix ne sont pas maîtrisées (un exemple volontairement pathologique est présenté à la Fig. 4.2). Même si les compétences développées par les spécialistes de ces approches maillées rendent généralement les erreurs causées par ce type de discrétisation faibles voire négligeables, ces dernières ne sont, en pratique, pas quantifiables, et font ainsi perdre aux méthodes de Monte-Carlo leur caractère de solution de référence. La seconde contrainte associée à ce type de résolution est d’ordre purement pratique : il est nécessaire à chaque nouveau cas d’étude, à chaque modification de géométrie ou de champs de propriétés de repenser la discrétisation du milieu et de produire un nouveau maillage. Cette étape nécessaire est généralement complexe et lourde à réaliser, en particulier lorsqu’il s’agit de valider la pertinence du maillage

Non-linéarité dans l’expression statistique de l’équation du transfert radiatif

 Les difficultés liées aux hétérogénéités des propriétés optiques du milieu d’intérêt ne se limitent toutefois pas à l’échantillonnage des positions de collision. Si l’on approche de façon purement statistique ce problème, il serait toujours possible d’insérer une nouvelle fonction densité de probabilité p˜L(l) qui elle, permettrait un échantillonnage aisé des libres parcours. Mais en pratique ce n’est pas envisageable : il serait alors nécessaire d’estimer l’épaisseur optique par une simulation de Monte-Carlo complète pour chacune des Nmc réalisations de l’algorithme permettant d’estimer Lη(x0, u0). Si l’on imagine que chaque algorithme (le premier estimant l’épaisseur optique et le second estimant la luminance) soit constitué de 106 réalisations indépendantes, il serait alors nécessaire de réaliser 1012 opérations, ce qui représenterait un temps de calcul prohibitif. S’il n’est pas possible de traiter statistiquement de façon simultanée ces deux termes intégraux : la luminance et l’épaisseur optique, c’est à cause de la fonction exponentielle qui introduit une non-linéarité dans l’expression statistique de la luminance Lη(x0, u0).

Les algorithmes à collisions nulles 

Historique des algorithmes à collisions nulles

Les algorithmes de Monte-Carlo à collisions nulles sont apparus au début des années soixante dans deux champs disciplinaires : la physique des Plasmas et la Neutronique. Il est intéressant de constater que cette méthode a vu le jour de manière totalement indépendante dans chacune de ces deux communautés, menant ainsi à deux ensembles distincts de travaux qui ont semblé s’ignorer jusqu’à aujourd’hui (à l’exception de la publication [Boeuf et Marode, 1982] faisant un lien entre ces deux communautés). Aussi, cette sous-section a pour objectif de présenter succinctement cette littérature plutôt complexe de par sa duplicité et de par la variété des termes employés pour décrire une même technique. Les algorithmes à collisions nulles ont été développés à la fin des années soixante dans le domaine de la physique des plasmas. Très utilisés dans ce champ d’application, ils permettent notamment de tenir compte des sections efficaces d’interaction dépendant de la vitesse des particules. On les rencontre dans cette communauté sous les dénominations : Null-Collisions, Fictitious-Collisions, PseudoCollisions, Null-Events ou encore Fictitious-Events. H.R. Skullerud est le premier à aborder dans [Skullerud, 1968] un Algorithme à Collisions Nulles, sans encore le dénommer ainsi, dans le but de pouvoir tirer statistiquement des « temps libres » entre deux collisions ion/molécule produites dans un gaz soumis à un champ électrique. De nombreux travaux vont alors s’ensuivre [Lin et Bardsley, 1977, Lin et Bardsley, 1978, Boeuf et Marode, 1982, Heifetz et al., 1982, Andreucci, 1985, Brennan, 1991, Longo, 2002, Longo et Diomede, 2004] visant pour la plupart à simuler les interactions entre particules chargées et molécules neutres sous l’influence d’un champ électrique. Les travaux de Skullerud ont également mené la communauté étudiant la dynamique des gaz raréfiés à s’intéresser aux Algorithmes à Collisions Nulles [Koura, 1986, Khisamutdinov et Sidorenko, 1995, Rjasanow et Wagner, 1998]. E. Woodcock a été, de son côté, à l’origine des algorithmes à collisions nulles dans le domaine de la neutronique [Woodcock et al., 1965]. Cette technique, étendue d’un point de vue théorique par Coleman [Coleman, 1968], sera alors intensivement utilisée dans ce champ applicatif. Parmi les principaux travaux, on peut citer [MacMillan, 1967, Spanier, 1970, Androsenko et al., 1991, Martin et Brown, 2001, Brown et Martin, 2003]. La place qu’occupera cette méthode dans cette communauté sera telle qu’elle sera implémentée nativement dans plusieurs codes de simulation de transport particulaire tels que SERPENT [Leppänen, 2007b, Leppänen, 2007a, Leppänen, 2010] ou encore MORET [Miss et al., 2007, Forestier et al., 2008]. Ces travaux conduiront des spécialistes d’autres domaines applicatifs tels que ceux de la synthèse d’image [Szirmay-Kalos et Tóth, 2010, Szirmay-Kalos et al., 2011], de la radiothérapie [Wang et al., 1997] et de la tomographie [Kawrakow et Fippel, 2000, Rehfeld et Stute, 2008, Kawrakow et al., 2008, Rehfeld et al., 2009, Badal et Badano, 2009, Tóth et Magdics, 2010] à s’en inspirer. On rencontre les algorithmes à collisions nulles dans la littérature associée à ces champs d’étude sous différentes dénominations : Woodcock-Tracking, Delta-Tracking, Hole-Tracking, Woodcock-Scattering, Delta-Scattering, Pseudo-Scattering ou encore Fictitious-Scattering. Toutefois, bien que très usités dans de nombreux domaines d’application de la physique du transport corpusculaire, les algorithmes à collisions nulles semblent, à notre connaissance, absents de la littérature propre à l’étude du rayonnement thermique. Aussi, ces travaux de thèse proposent, en partie, d’étendre des algorithmes à collisions nulles à ce champ applicatif.

Principe des algorithmes à collisions nulles

 Le principe des algorithmes à collisions nulles repose sur l’addition arbitraire d’un champ positif de coefficient de collision nulle kn,η dans le champ d’extinction réel : ˆkη = ka,η + kd,η + kn,η (4.10) Ce champ fictif de collisions nulles doit être défini de façon à rendre le champ du nouveau coefficient d’extinction ˆkη suffisamment simple pour permettre un échantillonnage aisé des libres parcours selon la fonction densité probabilité de Beer-Lambert : pˆL(l) = ˆkη(x0 − lu0)exp  − Z l 0 ˆkη(x0 − l 0u0)dl0  (4.11) Le champ du coefficient de collision nulle kn,η peut par exemple être défini de sorte à rendre celui de ˆkη uniforme (voir Fig. 4.3). En pratique, c’est le nouveau champ de coefficient d’extinction ˆkη qui est défini arbitrairement, le champ de kn,η n’étant jamais explicité mais seulement défini comme kn,η = ˆkη − ka,η.