Validation de l’algorithme à collisions nulles pour différentes épaisseurs optiques

Dans l’optique de valider l’algorithme présenté à la Fig. 4.7, ce dernier a été comparé à un algorithme standard de Monte-Carlo bien maîtrisé, formulé en épaisseur optique et dans lequel l’inversion des épaisseurs optiques est réalisée en ajustant de façon très précise le champ du coefficient d’extinction, le long d’une ligne de visée, par une décomposition en splines cubiques (inversible analytiquement). Pour cette validation, le bilan radiatif monochromatique a été calculé par les deux algorithmes en deux points du milieu participant d’intérêt : le centre du cube x0 = [0, 0, 0] et le point où les épaisseurs optiques sont maximales x0 = [−D, 0, 0] (voir Fig. 4.4) 2 . Le champ de ˆkη est défini de façon uniforme comme égal à k max a,η + k max d,η (ρ = 1). Ainsi, excepté en x = [−D, 0, 0] où il est nul, le champ du coefficient de collision nulle est toujours strictement positif. En première approximation, on considère également les parois comme noires (ε = 1) de température Tw = 0K et le paramètre d’asymétrie de la fonction de phase d’Henyey-Greenstein égal à 0 (fonction de phase isotrope). Les tables 4.1 et 4.2 rassemblent les résultats obtenus respectivement pour x0 = [0, 0, 0] et x0 = [−D, 0, 0] par les deux algorithmes et pour différentes épaisseurs optiques maximales d’absorption k max a,η D et de diffusion k max d,η D. Pour chacun de ces deux algorithmes, l’estimation du bilan radiatif ainsi que son écart-type adimensionnalisés sont donnés. Ceux-ci sont respectivement définis par Sr,η(x0)/[4πka,η(x0)L eq η (x0)] et σ/[4πka,η(x0)L eq η (x0)]. De plus, pour l’algorithme à collisions nulles, les temps de calcul pour 106 réalisations indépendantes (noté t) et pour obtenir un écart-type relatif de 1% (noté t1%) sont également fournis 3 . La simulation a été effectuée avec un processeur « Intel Core i5 – 2.4GHz » sans parallélisation. Les estimations du bilan radiatif par ces deux algorithmes concordent parfaitement. Les erreurs relatives indiquent également un bon niveau de convergence de l’algorithme à collisions nulles quelles que soient les épaisseurs optiques considérées (l’erreur relative e = σ/Sr,η obtenue après 106 réalisations est inférieure à 0.2% dans tous les cas). Enfin, les temps de calcul relevés pour l’algorithme à collisions nulles sont du même ordre de grandeurs que ceux que l’on rencontrerait avec un algorithme standard de Monte-Carlo dans lequel les propriétés du milieu seraient uniformes. L’ajout de collisions nulles, sans aucun effet sur la précision du calcul, semble donc n’avoir eu qu’un effet modéré sur ces temps de calcul. 

Prise en compte de coefficients de collision nulle négatifs

 Jusqu’à présent, l’algorithme proposé ne permettait que de définir un coefficient ˆkη supérieur en tous points au coefficient d’extinction réel kη. Il est possible, comme proposé à la Sec. 4.2.4 de définir de nouvelles probabilités pour s’affranchir de cette contrainte. Dans la continuité de la probabilité d’émission/absorption proposée à la Sec. 4.2.4 pour un cas purement absorbant/émettant, nous proposons ici les probabilités suivantes pour tenir compte de la diffusion .

Comportement numérique en fonction des valeurs du coefficient de collision nulle

Il devient ainsi possible de réaliser en toute généralité une étude de l’effet du ˆkη sur le comportement de l’algorithme de Monte-Carlo. Pour mener à bien cette analyse, plusieurs calculs du bilan radiatif monochromatique ont été effectués, à partir de l’algorithme de la Fig. 4.8, pour plusieurs valeurs de ρ = ˆkη/kmax η allant de ρ = 0.5 (où ˆkη ne majore que localement le coefficient d’extinction réel kη) à ρ = 5 (où ˆkη majore en tout point et très largement kη). Les Fig. 4.9, Fig. 4.10 et Fig. 4.11 décrivent respectivement les évolutions de l’erreur relative, du temps de calcul pour 106 réalisations et du temps de calcul pour obtenir une erreur relative de 1% en fonction de ρ, pour différentes épaisseurs optiques et deux points d’intérêts : x0 = [0, 0, 0] et x0 = [−D, 0, 0]. Pour ρ > 1, l’écart-type de l’estimation du bilan radiatif monochromatique est indépendant du coefficient de collision nulle (voir Fig. 4.9). En effet, les algorithmes à collisions nulles ne constituent qu’un artefact statistique et numérique permettant un échantillonnage plus aisé des libres parcours d’extinction. L’ajout de collisions nulles a seulement une incidence sur le temps de calcul : plus il y aura de collisions nulles (sans effet sur la convergence du calcul) plus le temps de calcul sera long (voir Fig. 4.10 et Fig. 4.11). Au contraire, plus la valeur de ρ est petite devant 1 , plus l’écart-type associé à l’estimation du bilan radiatif est important et croît de façon rapide, voir Fig. 4.8. Ce comportement était en effet attendu. L’introduction de nouvelles probabilités permettant de prendre en compte des occurrences négatives de kn,η engendre dans l’expression des poids de Monte-Carlo l’apparition d’un produit correctif. La valeur de ce dernier croît et est susceptible de changer de signe à chaque fois qu’une collision a lieu dans une région où les coefficients de collision nulle sont négatifs (cf. Sec. 4.2.4). Si un grand nombre d’événements de diffusion ou de collisions nulles se produisent le long du chemin optique dans une région où kn < 0, les poids de Monte-Carlo peuvent alors avoir des valeurs absolues très importantes et d’une très grande variance, expliquant ainsi l’accroissement conséquent de l’écart-type associé à l’estimation de Sr,η(x0). Cet effet est naturellement plus prononcé lorsque le point sonde x0 appartient à la zone où kn,η < 0 (voir Fig. 4.9b) que lorsque les chemins optiques partent d’une zone où kn,η > 0 (voir Fig. 4.9a). Toutefois, cette augmentation brutale de l’écart-type doit être relativisée : la proposition de la Sec. 4.2.4 permettant d’autoriser des occurrences de kn,η > 0 est faite pour éviter un biais des résultats de simulations dans le cas où le champ de ˆkη choisi ne majorerait pas parfaitement le champ du coefficient d’absorption réel (elle n’entraîne aucune modification si ˆkη(x) > kη(x)). Ainsi, on remarque que si le choix du champ de ˆkη est suffisamment bien pensé (ρ > 0.9), l’augmentation d’écart-type due aux coefficients négatifs de collisions nulles reste mesurée (elle est multipliée par 4 dans le cas le plus défavorable). Enfin, ces trois jeux de graphiques nous permettent de constater que les temps de calcul décroissent avec la valeur de ρ pour un nombre donné Nmc de réalisations indépendantes, (voir Fig. 4.10). Cela vient simplement du fait que plus faible est la valeur de ρ, moins il y a d’événements de diffusion et de collisions nulles. En effet, plus ρ est faible, plus la valeur de ˆk l’est aussi, les libres parcours échantillonnés selon pˆL(l) sont alors beaucoup plus longs, favorisant une absorption rapide aux parois. Toutefois, excepté pour le cas particulier d’un milieu très mince, cette décroissance des temps de calcul ne compense pas l’augmentation de l’erreur relative due aux coefficients négatifs de collision nulle (voir Fig. 4.11). Pour une erreur relative désirée, le temps de calcul est donc à la fois conditionné par l’effet de ρ sur la variance (voir Fig. 4.11 pour ρ < 1) et sur la quantité de collisions nulles qui augmente mécaniquement, mais dans une moindre mesure, ce temps de calcul (voir Fig. 4.11 pour ρ > 1). À la vue de l’évolution des temps de calcul nécessaires à l’obtention d’une erreur relative d’1%, il semble préférable, en cas de doute, de définir un champ de ˆk majorant largement celui du coefficient d’extinction que de risquer une explosion de variance causée par un trop grand nombre de collisions nulles caractérisées par des coefficients négatifs.

Influence de l’émissivité et du paramètre d’asymétrie de la fonction de phase

Plusieurs simulations ont également été réalisées pour des paramètres d’asymétrie de la fonction de phase différents de 0 et des émissivités de paroi inférieures à 1. Le paramètre d’asymétrie n’a que très peu d’effet sur le comportement numérique de l’algorithme à collisions nulles. Inévitablement, l’estimation de Sr,η(x0) est sensible à ces changements (puisque le modèle physique est modifié), mais les effets sur les erreurs relatives et les temps de calcul semblent négligeables. La réflexion multiple aux parois a, quant à elle, des conséquences sur le comportement de l’algorithme. Elle agit de la même manière que la diffusion : elle a une tendance à accroître la longueur des chemins parcourus par les photons avant d’être absorbés. Ainsi, un plus grand nombre de collisions ont lieu, augmentant alors les temps de calcul. Dans les zones où les coefficients de collision nulle sont négatifs, la réflexion multiple a tendance, au même titre que la diffusion (cf. Fig. 4.11) mais dans une moindre mesure, à augmenter davantage l’écart-type estimé.

Traitement déterministe des événements d’émission – approche par « Energy-partitioning » 

Jusqu’à présent, à chaque fois qu’un point de collision était identifié, un test statistique, de type roulette russe, avait lieu pour savoir s’il s’agissait d’une émission (par le milieu ou par la paroi), d’une réflexion, d’une diffusion ou d’une collision nulle. S’il s’agissait d’une émission, la réalisation serait stoppée et un poids wi serait calculé. Il pourrait cependant être intéressant de traiter de façon déterministe ces événements d’émission. Ainsi, après une émission/absorption les photons continueraient leur chemin, mais une information locale relative à l’émission serait prise compte. Les poids de Monte-Carlo wi seraient alors définis comme une somme de contributions d’émission rencontrées le long d’un chemin optique. Une telle approche, présentée à la Sec. 3.4.2, est connue sous le nom d’energy-partitioning ou de pathlengthmethod. L’application de ce traitement déterministe à l’estimation du bilan radiatif monochromatique Sr,η(x0) et ses conséquences sur le comportement numérique des algorithmes à collisions sont discutés dans cette section. 

 Influence de ζ sur le comportement numérique de l’algorithme à collisions nulles

Aussi, il est intéressant d’étudier l’influence que joue ce seuil ζ sur le comportement de l’algorithme de Monte-Carlo. Des bilans radiatifs monochromatiques Sr,η(x0) aux points x0 = [0, 0, 0] et x0 = [−D, 0, 0] ont alors été estimés par l’algorithme présenté à la Fig. 4.12, pour différentes valeurs de ζ, de ρ et pour différentes épaisseurs optiques maximales k max a,η D et k max d,η D. La valeur du champ d’extinction arbitraire ˆkη est définie comme égale à k max η (ρ = 1, les coefficients de collision nulle sont positifs en tout point).