Méthode MN

La méthode M_N est une méthode utilisant les moments de l'équation de Boltzmann introduite par Harold Grad en théorie cinétique des gaz. Elle est appliquée à la résolution du transfert radiatif. Dans cette méthode la fermeture du système utilise la maximisation de l'entropie de Bose-Einstein dans le cas des photons, celle de Boltzmann dans d'autres cas, par exemple en irradiation médicale. Cette méthode a été introduite par Gerald Minerbo en 1978. Elle est essentiellement utilisée dans les problèmes de transfert énergétique mais donne une solution trop imprécise angulairement pour les problèmes optiques.

Températures dans un jet de disque protoplanétaire. Calcul M₁ couplé à l'hydrodynamique.

Moments de la luminance

Valeurs spectrales

L'équation de transfert radiatif porte sur la luminance spectrale L (x, t, θ, ν) correspondante à l'intervalle de fréquence (ν, ν + dν). On se place en une dimension d'espace et la distribution angulaire est de révolution, décrite par le seul angle de colatitude θ pour des coordonnées sphériques. On supposera de plus que les coefficients d'absorption et d'extinction sont indépendants de la fréquence et que le terme d'émission en volume correspond à un corps noir

$B(\nu ,T)={\frac {2h\nu ^{3}}{c^{2}}}n_{eq}\,,\qquad n_{eq}={\frac {1}{e^{\frac {h\nu }{kT}}-1}}$

où n_eq est le nombre d'occupation à l'équibre thermodynamique, h est la constante de Planck, k la constante de Stefan-Boltzmann et c la vitesse de la lumière. La température T (x, t) est supposée connue.

La diffusion est élastique (sans changement de fréquence) et caractérisée par la fonction de phase P (α) fonction de la seule déviation α.

L'équation de Boltzmann instationnaire qui décrit la conservation de l'énergie dans la direction correspondante à μ s'écrit ici

${\frac {1}{c}}{\frac {\partial L}{\partial t}}+\mu {\frac {\mathrm {d} L}{\mathrm {d} x}}+\kappa _{t}L=\kappa _{a}B+\kappa _{d}\int _{-1}^{1}P(\alpha )L\mathrm {d} \mu$

$\mu =\cos(\theta )$	donne la direction de propagation pour une fonction de phase supposée de révolution,
$\alpha =\mu \mu '+{\sqrt {1-\mu ^{2}}}{\sqrt {1-\mu '^{2}}}$	déviation (μ, μ'),
$\kappa _{t}(x,t)=\kappa _{a}(x,t)+\kappa _{d}(x,t)$	coefficient d'extinction totale,
$\kappa _{a}(x,t)$	coefficient d'absorption,
$\kappa _{d}(x,t)$	coefficient d'extinction,

Les premiers moments sont

• l'énergie volumique

$E(x,t,\nu )={\frac {2\pi }{c}}\int _{-1}^{1}L\mathrm {d} \mu$

• l'exitance (le flux) M et sa valeur normalisée f

$M(x,t,\nu )=2\pi \int _{-1}^{1}L\mu \mathrm {d} \mu \,,\qquad f={\frac {M}{cE}}$

• la pression et sa valeur normalisée (coefficient d'Eddington)

$P(x,t,\nu )={\frac {2\pi }{c}}\int _{-1}^{1}L\mu ^{2}\mathrm {d} \mu \,,\qquad \chi ={\frac {P}{E}}$

À l'équilibre thermodynamique

$L_{eq}=B\,,\quad E_{eq}={\frac {4\pi B}{c}}\,,\quad M_{eq}=0\,,\quad P_{eq}^{R}={\frac {4\pi B}{3c}}={\frac {E_{eq}}{3}}$

Valeurs intégrées en fréquence

Compte tenu du choix de l'indépendance des coefficients vis-à-vis de la fréquence toutes les équations ci-dessus sont valables pour les quantités intégrées en fréquence

$L^{R}=\int _{0}^{\infty }L\mathrm {d} \nu$

$E^{R}={\frac {2\pi }{c}}\int _{0}^{\infty }\int _{-1}^{1}L\mathrm {d} \mu \mathrm {d} \nu ={\frac {1}{c}}\langle L\rangle =aT_{R}^{4}$

$M^{R}=\int _{0}^{\infty }\int _{-1}^{1}L\mu \mathrm {d} \mu \mathrm {d} \nu =\langle L\mu \rangle$

$P^{R}={\frac {1}{c}}\int _{0}^{\infty }\int _{-1}^{1}L\mu ^{2}\mathrm {d} \mu \mathrm {d} \nu ={\frac {1}{c}}\langle L\mu ^{2}\rangle$

On a défini la notation ⟨ . ⟩ pour l'intégrale en fréquence et en angle, ainsi que la température radiative T_R.

À l'équilibre thermodynamique

$L_{eq}^{R}={\frac {\sigma T^{4}}{\pi }}\,,\quad E_{eq}^{R}=aT^{4}\,,\quad M_{eq}^{R}=0\,,\quad P_{eq}^{R}={\frac {aT^{4}}{3}}={\frac {E_{eq}^{R}}{3}}\,,\quad a={\frac {4\sigma }{c}}$

où σ est la constante de Stefan-Boltzmann.

Système aux moments

En prémultipliant l'équation de Boltzmann par 1, μ, ... et en intégrant sur μ on obtient un système sur les moments de la luminance

${\begin{array}{lcl}{\frac {\partial E}{\partial t}}+\mu {\frac {\partial M}{\partial x}}&=&\kappa _{a}(4\pi B-cE)\\[0.6em]{\frac {\partial M}{\partial t}}+c^{2}\mu {\frac {\partial P}{\partial x}}&=&-c\kappa _{f}(x,t)M\\[0.6em]...\end{array}}$

où apparaît le premier moment de la fonction de phase g₁

$\kappa _{f}(x,t)=\kappa _{a}(x,t)+[1-g_{1}\kappa _{d}(x,t)]\,,\qquad g_{1}=\int _{-1}^{1}P(\alpha )\alpha \mathrm {d} \alpha$

g₁ = 0 pour toute distribution symétrique : distribution isotrope ou diffusion Thomson par exemple.

Le système sur les valeurs intégrées en fréquence est le même (en remplaçant 4 π B par a T⁴) puisqu'on a supposé les coefficients indépendants de celle-ci.

Propriétés

Par substitution dans le système ci-dessus on voit que celui-ci est entièrement décrit par E, f et χ, ces deux dernières valeurs jouant le rôle de nombres adimensionnels. Ils sont appelés premier et second facteur d'Eddington dans la théorie à N moments.

Ces moments sont des moments de Hausdorff et possèdent un certain nombre de propriétés

• les symétries du problème entraînent    ${\frac {1}{3}}\leq \chi \leq 1$

χ = 1 / 3 correspond à l'isotropie et χ = 1 à un rayon.

• l'inégalité de Cauchy-Schwarz impose  $\chi \geq f^{2}$

Par ailleurs :

• $|f|\leq 1$ : cette inégalité découle des expressions précédentes. Elle exprime que le flux maximal possible correspond au transfert de l'énergie à la vitesse c

$|M_{max}|=cE$

• le système décrit ci-dessus est hyperbolique. En supposant que les valeurs propres du système croissent avec | f | on montre que χ est une fonction convexe de f et que la dérivée χ' (f) vérifie | χ' (±1) | = 2 afin que la vitesse caractéristique en χ = ±1 (correspondante à un rayon lumineux) soit ±c.

On voit que f et χ sont sévèrement limités. Ceci a incité les physiciens à proposer un certain nombre de fermetures du type χ = g (f) appelées méthodes à facteur d'Eddington variable.

Fermeture entropique

Distribution angulaire pour diverses valeurs de f (0, 0.3, 0.7, 0.9) dans la méthode M₁ (la direction du flux est ici π / 4 ).

Le système ci-dessus est, comme tout système aux moments, incomplet : chaque équation sur un moment fait apparaître un moment d'ordre supérieur inconnu. Au-delà du travail initial de Gerald Minerbo une méthode générale a été introduite par David Levermore. Elle utilise l'entropie volumique spectrale d'un gaz de photons

${\mathcal {S}}(L)=[(n+1)\ln(n+1)-n\ln n]\,{\frac {2k\nu ^{2}}{c^{3}}}\,,\qquad n={\frac {c^{2}}{2h\nu ^{3}}}\,L$

où n est le nombre d'occupation.

La fermeture consiste à choisir la luminance spectrale qui maximise l'entropie sous la contrainte de respecter les moments pris en compte. Ceci est toujours possible, l'entropie étant une fonction convexe de L. Un exemple est celui de la recherche de la fonction d'équilibre thermodynamique B supposée inconnue

${\mathcal {S}}(B)=\min _{L}\{{\mathcal {S}}(L):{\frac {1}{c}}\langle L\rangle =aT^{4}\}$

La solution de cette équation est la distribution de Planck définie plus haut.

Démonstration

Le problème de la recherche du minimum est équivalent à la celui du point-selle du lagrangien

${\mathcal {L}}(L,\alpha )={\mathcal {S}}(L)+\alpha \left({\frac {1}{c}}\langle L\rangle -aT^{4}\right)$

où α est le multiplicateur de Lagrange. La contrainte étant réalisée par définition le minimum est défini par

${\frac {\partial {\mathcal {L}}}{\partial L}}=0$

soit

$\left\langle \left[-{\frac {k}{hc\nu }}\log {\left(1+{\frac {1}{n}}\right)}+\alpha \right]\right\rangle =0$

dont la solution est

$\alpha ={\frac {k}{hc\nu }}\log {\left(1+{\frac {1}{n}}\right)}$

D'où la solution

$B={\frac {2h\nu ^{3}}{c^{2}}}{\frac {1}{e^{\frac {hc\alpha \nu }{k}}-1}}$

La contrainte E^R = a T⁴ permet de calculer la constante

$\alpha ={\frac {1}{cT}}$

Dans le cas général on note m = (1, μ, μ²,...). Le problème de minimisation s'écrit

${\mathcal {S}}(B)=\min _{L}\{{\mathcal {S}}(L):\langle \mathbf {m} L\rangle ={\mathcal {M}}_{i}\}$

où ℳ₀, ℳ₁, ℳ₂,...(ces quantités sont l'énergie, le flux, la pression,... et sont ici des paramètres du problème) sont les moments à respecter. On leur substitue les moments adimensionnés f, χ,.... La solution de ce problème est la fonction duale au sens de Legendre de l'entropie

$B={\frac {2h\nu ^{3}}{c^{2}}}{\frac {1}{e^{{\frac {hc\nu }{k}}(\alpha _{0}+\alpha _{1}\mu +\alpha _{2}\mu ^{2}+...)}-1}}$

La démonstration est identique à celle de l'encadré précédent. La solution obtenue correspond à un spectre de corps noir et à une distribution angulaire d'autant plus complexe que l'ordre est élevé.

Dans le cas où l'on prend les deux premiers moments la solution appelée M₁ est analytique et conduit à la relation de fermeture

$\chi ={\frac {3+4f^{2}}{5+2{\sqrt {4-3f^{2}}}}}$

Il n'est pas nécessaire de connaître la distribution angulaire qui est sous-jacente dans le problème résolu. Celle-ci est régulière (voir courbe jointe) pouvant aller de l'isotropie au rayon. La méthode M₂ permet d'obtenir des distributions plus complexes.

Propriétés de la solution

Le caractère hyperbolique permet d'utiliser des solveurs numériques efficaces développés pour la mécanique des fluides. Le couplage avec l'hydrodynamique en est grandement facilité. La solution du système couplé est entropique (l'entropie du système croît au cours du temps) ce qui assure un bon comportement numérique.

Une pathologie de cette méthode est la possibilité de discontinuités de la solution comme en mécanique des fluides : des « chocs radiatifs » sans existence physique rencontrés lors du croisement de deux faisceaux. Ce problème n'existe que pour les développements de moments à l'ordre impair. Ils peuvent être surmontés par une méthode de partition de l'espace analogue à la méthode à deux flux, laquelle n'est guère utilisée, le problème ne se posant généralement pas.

La méthode a été étendue au problème multidimensionnel en espace où M est une quantité vectorielle et P et les moments suivants des quantités tensorielles et à un spectre quelconque en utilisant le découpage en bandes de fréquences. Il existe également des travaux sur la méthode M₂ utilisant l'entropie de Maxwell-Boltzmann.

Un exemple

Distributions angulaires dans une lame isotherme en x / e = 0.95 et κ e = 3 (rouge : solution exacte, vert : méthode M₁).

Prenons l'exemple très simple du rayonnement dans et hors d'une lame d'épaisseur e, de température T, émissive et absorbante mais non diffusante, en régime stationnaire. Le repère utilisé a son origine sur l'une des faces et pointe vers l'intérieur du domaine. L'équation de Boltzmann est ici

$\mu {\frac {\mathrm {d} L^{R}}{\mathrm {d} x}}=\kappa _{a}(aT^{4}-L^{R})$

La solution de ce problème différentiel linéaire du premier ordre est

$4\pi L^{R}=\left\{{\begin{array}{lc}acT^{4}\left(1-e^{-{\frac {\kappa x}{\mu }}}\right)&\mu \geq 0\\[0.5em]acT^{4}\left[1-e^{-{\frac {\kappa (e-x)}{\mu }}}\right]&\mu <0\end{array}}\right.$

$M^{R}=2\sigma T^{4}\left\{E_{3}[\kappa _{a}(x-e)]-E_{3}(\kappa _{a}x)\right\}$

$E^{R}={\frac {aT^{4}}{2}}\left\{2-E_{2}[\kappa _{a}(x-e)]-E_{2}(\kappa _{a}x)\right\}$

où E_i est l'exponentielle intégrale d'ordre i.

Cet exemple a le mérite d'être simple et donc de permettre de tester le cas plus général d'un spectre quelconque en permettant un calcul de référence raie par raie. La comparaison avec la solution M₁ (voir courbe ci-contre) montre que la distribution angulaire n'est que médiocrement approchée, ce qui ne fait pas de la méthode une bonne candidate pour les problèmes d'optique. Pourtant l'erreur sur le flux sortant n'excède pas 3%. Ce chiffre se compare à la méthode SN ou à la méthode PN avec N ≈ 10, pour un coût calcul notablement moindre. Une telle erreur doit être rapprochée de la connaissance sur les coefficients d'absorption qui sont de l'ordre de 10% dans les meilleurs cas.

Voir aussi

• Méthode PN