JOURNAL INTERNATIONAL DE TECHNOLOGIE, DE L'INNOVATION,
DE LA PHYSIQUE, DE L'ENERGIE ET DE L'ENVIRONNEMENT

# Étude théorique de la modélisation du rayonnement d'un plasma thermique : application à l'arc libre

A. Truilhé, C. Gouze, Y. Cressault, F. Valensi

<figure><img src="image_1.png" /></figure>

ISSN : 2428-8500

DOI : 10.52497/jitipee.v9i1.383

Le sujet de cet article a été présenté lors du colloque sur les arcs électriques (CAE XVII) les 17 et 18 mars 2025 à Limoges.

L'article publié a fait l'objet d'une expertise indépendante par deux spécialistes du domaine.
# Étude théorique de la modélisation du rayonnement d'un plasma thermique : application à l'arc libre

Antoine Truilhé (1), Camille Gouze (1)(2), Yann Cressault(1), Flavien Valensi(1)

antoine.truilhe@laplace.univ-tlse.fr

**Résumé** – Le rayonnement est un aspect important de la physique des plasmas thermiques, notamment dans sa modélisation. Dans ces travaux seront présentées les méthodes couramment utilisées pour la modélisation du rayonnement des plasmas thermiques. On rappelle et décrit assez précisément les étapes de calculs qui aboutissent aux équations utilisées pour simuler le transfert radiatif dans les logiciels de modélisation COMSOL Multiphysics et ANSYS Fluent. A titre d'exemples et donc sans détailler l'implémentation de ces méthodes dans les deux logiciels, nous présentons, en fin d'étude, des résultats sur un arc libre d'argon à 200A afin de comparer les différentes méthodes entre elles et ainsi montrer leur impact sur la physique de l'arc.

## Nomenclature

<table><tr><td>B</td><td>[T]</td><td>champ magnétique</td></tr><tr><td>C<sub>p</sub></td><td>[J.K<sup>-1</sup>.kg<sup>-1</sup>]</td><td>chaleur spécifique à pression constante</td></tr><tr><td>e</td><td>[C]</td><td>charge élémentaire de l'électron (1.6021892.10<sup>-19</sup> C)</td></tr><tr><td>E</td><td>[V.m<sup>-1</sup>]</td><td>champ électrique</td></tr><tr><td>g</td><td>[m.s<sup>-2</sup>]</td><td>accélération de la pesanteur</td></tr><tr><td>h</td><td>[J.kg<sup>-1</sup>.m<sup>-3</sup>]</td><td>enthalpie massique par unité de volume</td></tr><tr><td>j</td><td>[A.m<sup>-2</sup>]</td><td>vecteur densité de courant</td></tr><tr><td>k</td><td>[W.m<sup>-1</sup>.K<sup>-1</sup>]</td><td>conductivité thermique</td></tr><tr><td>k<sub>B</sub></td><td>[J.K<sup>-1</sup>]</td><td>constante de Boltzmann (1.380662.10<sup>-23</sup>)</td></tr><tr><td>K'<sub>λ</sub></td><td>[m<sup>-1</sup>]</td><td>coefficient d'absorption corrigé de l'émission induite</td></tr><tr><td>L<sub>λ</sub></td><td>[W.m<sup>-3</sup>.sr<sup>-1</sup>]</td><td>luminance spectrale</td></tr><tr><td>m</td><td>[kg]</td><td>masse</td></tr><tr><td>P</td><td>[Pa]</td><td>pression</td></tr><tr><td>Q</td><td>[J.m<sup>-3</sup>]</td><td>énergie thermique par unité de volume</td></tr><tr><td>S<sub>rad</sub></td><td>[J.m<sup>-3</sup>]</td><td>pertes radiatives par unité de volume</td></tr><tr><td>T</td><td>[K]</td><td>température</td></tr><tr><td>v</td><td>[m.s<sup>-1</sup>]</td><td>vecteur vitesse</td></tr><tr><td>ε<sub>λ</sub></td><td>[W.m<sup>-3</sup>.sr<sup>-1</sup>.m<sup>-1</sup>]</td><td>émission propre d'un milieu</td></tr><tr><td>λ</td><td>[m]</td><td>longueur d'onde</td></tr><tr><td>μ</td><td>[kg.m<sup>-1</sup>.s<sup>-1</sup>]</td><td>viscosité</td></tr><tr><td>ρ</td><td>[kg.m<sup>-3</sup>]</td><td>masse volumique</td></tr><tr><td>σ</td><td>[S.m<sup>-1</sup>]</td><td>conductivité électrique</td></tr><tr><td>τ</td><td>[N.m<sup>-2</sup>]</td><td>tenseur des forces visqueuses</td></tr></table>

**Mots clés** : modélisation magnétohydrodynamique, Coefficient d'Émission Nette, harmoniques sphériques, méthode P-1, Méthode des Ordonnées Discrètes, arc libre, COMSOL Multiphysics, ANSYS Fluent.

DOI: 10.52497/jitipee.v9i1.383

<footnote>(1) Université de Toulouse, Toulouse INP, CNRS, LAPLACE, Toulouse, France</footnote>

<footnote>(2) IRT Saint Exupéry, Toulouse, France</footnote>
## Introduction

Ces travaux ont pour objectif l'étude des méthodes de calcul du rayonnement d'un plasma thermique en s'appuyant sur la modélisation d'un arc libre d'argon à 200A entre deux électrodes. L'enjeu et l'importance des différents logiciels de simulation, notamment dans ce domaine d'étude où la mise en place d'expériences s'avère très complexe du fait des grandes valeurs de températures au sein du milieu (plusieurs dizaines de milliers de Kelvin), mais également des grandes valeurs de tension ou de courant nécessaires afin de réaliser le claquage de l'arc. En effet, les manipulations expérimentales restent faisables et intéressantes mais le modèle numérique permet d'accéder à de nombreux paramètres, dans différents régimes, et ceci assez rapidement en manipulant la géométrie (faire des coupes, ...) ou même de pénétrer au cœur du plasma sans le perturber, ce qui s'avère parfois infaisable avec les méthodes expérimentales.

## 1. Rappels théoriques : équations de la magnétohydrodynamique

Les plasmas thermiques sont des plasmas pour lesquels les densités électroniques varient de $10^{20}m^{-3}$ à $10^{26}m^{-3}$ à 10 kK. Ce type de plasmas, constitués de molécules, d'atomes, d'ions et d'électrons mais globalement électriquement neutre, se retrouve aussi bien dans les appareils de coupure que dans les torches à plasma. Cependant, dans les plasmas thermiques, il existe plusieurs types d'équilibre thermodynamique comme les plasmas à l'équilibre thermodynamique complet (ETC), les plasmas à l'équilibre thermodynamique local (ETL) ou encore les plasmas hors d'équilibre. Le plasma est dit à l'ETC lorsqu'il est entièrement caractérisé par quatre lois avec une seule température : la loi de distribution de Maxwell-Boltzmann, la loi de distribution de Boltzmann, la loi d'action de masse et la loi de Planck. On parle d'ETL lorsque l'ETC est vérifié localement. Afin de modéliser un tel plasma thermique, on utilise les équations de la magnétohydrodynamique (MHD). Tout d'abord, afin d'obtenir un modèle mathématique cohérent, un certain nombre d'hypothèses doivent être prises en compte pour les équations de la MHD :

- le plasma est en équilibre thermodynamique local (ETL) ;

- le plasma d'arc est considéré comme un milieu continu ;

- le plasma d'arc est traité comme un fluide neutre conducteur ;

- le fluide a une description eulérienne.

Ainsi, la modélisation la plus simple d'un arc le considère comme stable, axisymétrique et optiquement fin ; sans prise en compte de la gravité et de la dissipation de chaleur due à la viscosité [1].

Les équations de la magnétohydrodynamique sont constituées d'une première famille d'équations : les équations de conservations.

* Conservation de la masse :

La conservation de la masse impose qu'il n'y a ni apparition ni disparition de matière dans la particule de fluide considérée.

$$ \frac{\partial \rho}{\partial t} + \text{div} \rho v = 0 \qquad (1) $$

avec $\rho$ est la masse volumique et $v$ la vitesse de la particule.
* Conservation de la quantité de mouvement :

La conservation de la quantité de mouvement est l'application du Principe Fondamental de la Dynamique.

$$ \frac{\partial \rho v}{\partial t} + \text{div}(\rho v \cdot v) = \rho g - \text{grad}P - \text{div}\tau + j \times B \quad (2) $$

* Conservation de l'énergie :

L'équation de conservation de l'énergie découle du premier principe de la thermodynamique qui définit la variation de l'énergie totale d'une particule fluide comme étant égale à la somme de la quantité de chaleur et du travail reçus par la particule fluide (variation de l'énergie interne nulle) :

$$ \frac{\partial \rho h}{\partial t} + \text{div}(\rho v h) = -\text{div}Q + \frac{\partial P}{\partial t} + v \cdot \text{grad}P + j \cdot E + \frac{5 \cdot k_p}{2 \cdot e} j \cdot \text{grad} \frac{h}{c_p} - S_{rad} \quad (3) $$

Dans cette dernière équation figure le terme $S_{rad}$ correspondant à l'énergie transportée par rayonnement dans le plasma et qui fait l'objet de cette présente étude.

## 2. Modèles de Rayonnement

L'énergie transportée par rayonnement dans le plasma s'obtient par la résolution de l'Équation du Transfert Radiatif (ETR) au sein de celui-ci.

### 2.1. Équation du Transfert Radiatif

<figure><img src="image_2.png" /><figcaption>Figure 1 : Flux de rayonnement</figcaption></figure>

Considérons un élément de surface $dS$, perpendiculaire à la direction de propagation $\vec{s}$, défini par un point $P(\vec{r})$ et d'un angle solide $d\Omega$ (Fig. 1.) [2].

La luminance spectrale $L_\lambda$ au point P suivant la direction $\vec{s}$ est la puissance rayonnée par angle solide $d\Omega$, par unité de surface S, par unité de longueur d'onde $\lambda$ (en $W \cdot m^{-2} \cdot sr^{-1} \cdot m^{-1}$). Ainsi, en l'intégrant sur tout le domaine spectral, nous obtenons la luminance totale $L$ (en $W \cdot m^{-2} \cdot sr^{-1}$) :

$$ L(\vec{r}, \vec{s}) = \int_{\lambda} L_{\lambda}(\vec{r}, \vec{s}) d\lambda \quad (4) $$

On peut alors définir le flux de rayonnement $\vec{F}_R$ dans toutes les directions par intégration sur l'angle solide :

$$ \vec{F}_R(\vec{r}, t) = \int_{\Omega} L(\vec{r}, \vec{s}) \vec{s} \, d\Omega = \int_{\lambda} \int_{\Omega} L_{\lambda}(\vec{r}, \vec{s}) \vec{s} \, d\lambda \, d\Omega \quad (5) $$
L'Équation du Transfert Radiatif (ETR) ci-dessous, décrit l'évolution de la luminance spectrale à travers un élément de volume :

$$ \frac{dL_{\lambda}(\vec{r}, \vec{s})}{ds} = K'(\lambda, T(\vec{r}, \vec{s})) (L_{\lambda}^{0}(T(\vec{r}, \vec{s})) - L_{\lambda}(T(\vec{r}, \vec{s}))) \quad (6) $$

Qui s'écrit le plus souvent :

$$ \frac{dL_{\lambda}}{ds} = K'_{\lambda} L_{\lambda}^{0} - K'_{\lambda} L_{\lambda} \quad (7) $$

$$ \text{où } \varepsilon_{\lambda} = K'_{\lambda} L_{\lambda}^{0} $$

avec

* $L_\lambda$ la luminance spectrale (en W.m⁻².sr⁻¹.m⁻¹) variant en fonction du gaz, de la température, de la pression, de la longueur d'onde et de la direction ;

* $L_\lambda^0$ la fonction de Planck en fonction de la température (en W.m⁻².sr⁻¹.m⁻¹) ;

* $K'_\lambda$ le coefficient d'absorption corrigé de l'émission induite à la longueur d'onde $\lambda$ (en m⁻¹) variant en fonction du gaz, de la température, de la pression et de la longueur d'onde ;

* $\varepsilon_\lambda$ l'émission propre du milieu à la longueur d'onde $\lambda$ et la position $(\vec{r})$

Cette équation traduit l'évolution du rayonnement $L_\lambda$ à travers un milieu. Elle peut se décomposer en deux termes : le premier, noté $\varepsilon_\lambda$, correspondant à l'émission propre du milieu et relié au coefficient d'absorption $K'_\lambda$ par la loi de Kirchhoff, le second correspondant à l'absorption de celui-ci.

<figure><img src="image_3.png" /><figcaption>Figure 2 : Intégration de l'ETR</figcaption></figure>

On résout l'ETR suivant une direction $\vec{s}$ (ici $\overline{OX}$). Après avoir intégré les deux termes entre 0 et X (Fig. 2), la luminance $L_\lambda(X)$ se met alors sous la forme suivante [2] :

$$ L_{\lambda}(X) = L_{\lambda}(0) \exp\left(-\int_{0}^{X} K'_{\lambda}(\xi) d\xi\right) + \int_{0}^{X} L_{\lambda}^{0}(\xi) K'_{\lambda}(\xi) \exp\left(-\int_{\xi}^{X} K'_{\lambda}(\eta) d\eta\right) d\xi \quad (8) $$

**Terme I :** solution homogène de l'ETR, représentant la luminance initiale en X = 0, c'est-à-dire l'énergie radiative provenant de l'extérieur, dans la direction considérée, et qui est absorbée dans le milieu d'épaisseur X.

**Terme II :** solution particulière de l'ETR, correspondant à l'émission propre du milieu, représentant le rayonnement émis en $\xi$ et absorbé par le milieu d'épaisseur $(X - \xi)$.
On suppose alors qu'il n'y a aucune énergie radiative venant de l'extérieur vers le milieu ($L_λ(0) = 0$), on ne s'intéresse qu'à l'émission propre, que l'on peut écrire en intégrant sur l'ensemble des longueurs d'onde :

$$L(X) = \int_{0}^{\infty} \int_{0}^{X} L_{\lambda}^{0}(\xi) K'_{\lambda}(\xi) \exp \left( -\int_{\xi}^{X} K'_{\lambda}(\eta) d\eta \right) d\xi d\lambda \quad (9)$$

Compte tenu du nombre de paramètres pouvant influer la résolution de cette équation, il est inimaginable de la résoudre numériquement avec exactitude à cause du nombre de dépendances : spatiale (correspondant à la géométrie du problème), directionnelle (correspondant au nombre de directions possibles) et spectrale (correspondant au nombre de longueurs d'onde à considérer) auxquelles s'ajoutent la composition du gaz, la température et la pression. Pour cela, plusieurs approches simplifiées sont utilisées dans les modèles afin de résoudre cette équation avec le plus de précision possible en fonction des cas, tout en conservant un temps de calcul raisonnable.

## 2.2. Méthode du Coefficient d'Emission Nette

<figure><img src="image_4.png" /><figcaption>**Figure 3 :** Simplification géométrique du plasma pour la méthode du CEN</figcaption></figure>

La méthode du Coefficient d'Emission Nette (CEN) est une simplification géométrique de l'ETR. Le plasma est alors une sphère de rayon $R_p$ supposée isotherme, isobare et homogène (Fig. 3.) [3]. On peut définir ce coefficient comme étant la part de puissance rayonnée par unité de volume et par stéradian d'un élément sur l'axe de l'arc électrique et qui s'échappe de la sphère après avoir traversé une épaisseur $R_p$ du plasma. Il ne dépend donc que de la température, du rayon du plasma, de la pression et de la nature du gaz.

Par cette méthode, il est possible d'établir une relation simple entre $ε_N$ et la Divergence du Flux Radiatif (DFR), permettant en tout point du plasma de décrire le bilan énergétique radiatif.

On part de l'équation du transfert radiatif (6). Le premier terme de gauche représente la divergence de la luminance spectrale. En l'intégrant sur tout l'angle solide $Ω$, nous obtenons la divergence du flux radiatif spectral [4] :

$$4\pi \frac{dL_{\lambda}(\vec{s})}{ds} = 4\pi K'(\lambda, T) L_{\lambda}^{0}(T) - K'(\lambda, T) \int_{\Omega} L_{\lambda}(T) d\Omega = \vec{\nabla} \cdot \overrightarrow{q_{\lambda}} \quad (10)$$

Avec $\overrightarrow{q_{\lambda}}$ le flux radiatif spectral.
On peut alors poser le CEN dépendant de la longueur d'onde $\varepsilon_{N\lambda}$ sous la forme :

$$ \varepsilon_{N\lambda}(T) = K'(\lambda, T) L_{\lambda}^{0}(T) - K'(\lambda, T) J_{\lambda} \quad (11) $$

$$ \text{avec } J_{\lambda} = \frac{\int_{\Omega} L_{\lambda}(T) d\Omega}{4\pi} $$

$J_{\lambda}$ représente la luminance moyenne sur l'angle solide $d\Omega$ pour une longueur d'onde donnée.

Dans le cas d'une sphère isotherme, on a :

$$ J_{\lambda} = \frac{1}{2\pi} \int_{\rho'=0}^{\rho} \int_{\theta=0}^{\frac{\pi}{2}} \int_{\varphi=0}^{2\pi} K'(\lambda, T) L_{\lambda}^{0}(T) \exp\left(-\int_{0}^{\rho'} K'(\lambda, T) d\rho''\right) \sin\theta \, d\theta \, d\varphi \, d\rho' \quad (12) $$

en posant que $\rho' = 0$ correspondant au centre de la sphère et $\rho' = \rho$ le bord, $d\rho'' = \frac{dr''}{\sin\theta}$ et $d\rho' = \frac{dr}{\sin\theta}$, on obtient :

$$ J_{\lambda} = \frac{1}{2\pi} \int_{r=0}^{R_p} \int_{\theta=0}^{\frac{\pi}{2}} \int_{\varphi=0}^{2\pi} K' L_{\lambda}^{0} \exp\left(-\int_{0}^{r} \frac{K'}{\sin\theta} dr''\right) d\theta \, d\varphi \, dr \quad (13) $$

Etant une sphère, on peut simplifier en :

$$ J_{\lambda} = \int_{r=0}^{R_p} \int_{\theta=0}^{\frac{\pi}{2}} K' L_{\lambda}^{0} \exp\left(-\frac{K'r}{\sin\theta}\right) dr \, d\theta \quad (14) $$

On pose alors $x = K'r$ et $dx = K'dr$, pour trouver :

$$ J_{\lambda} = \int_{x=0}^{K'R_p} \int_{\theta=0}^{\frac{\pi}{2}} K' L_{\lambda}^{0} \exp\left(-\frac{x}{\sin\theta}\right) \frac{dx}{K'} \, d\theta = \int_{x=0}^{K'R_p} \int_{\theta=0}^{\frac{\pi}{2}} L_{\lambda}^{0} \exp\left(-\frac{x}{\sin\theta}\right) \, d\theta \, dx \quad (15) $$

On pose alors $G_n(x)$ tel que :

$$ G_n(x) = \int_{\theta=0}^{\frac{\pi}{2}} \sin^n \theta \exp\left(-\frac{x}{\sin\theta}\right) d\theta \quad (16) $$

J.J. Lowke [5] a montré que :

$$ \int_{0}^{x} G_{0}(x') dx' = 1 - G_{1}(x) \quad (17) $$

Ainsi, on peut simplifier l'expression (14) par :

$$ J_{\lambda} = L_{\lambda}^{0} \left(1 - G_{1}(K'R_{p})\right) \text{ avec } G_{1}(x) = \int_{\theta=0}^{\frac{\pi}{2}} \sin\theta \exp\left(-\frac{x}{\sin\theta}\right) d\theta \quad (18) $$

Cependant, Liebermann et Lowke [6] ont montré que cette fonction $G_1(x)$ se comporte de la même manière que $\exp(-x)$ dans le cas d'une sphère isotherme.

Donc, dans le cas d'une sphère isotherme de rayon $R_p$, nous obtenons :

$$ J_{\lambda} = L_{\lambda}^{0} \left(1 - \exp(-K'R_{p})\right) \quad (19) $$
En reprenant alors la formule du CEN par unité de longueur d'onde (11), et en prenant en compte l'expression de $J_λ$ on obtient :

$$ \epsilon_{N\lambda} = K'(\lambda, T) \cdot L_{\lambda}^{0}(T) - K'(\lambda, T) \cdot J_{\lambda} = K' L_{\lambda}^{0} - K' L_{\lambda}^{0} (1 - \exp(-K' R_{p})) \quad (20) $$

On peut alors intégrer cette expression sur l'ensemble des longueurs d'ondes :

$$ \epsilon_N = \int_0^\infty [K'_\lambda L_\lambda^0 - K'_\lambda L_\lambda^0 (1 - \exp(-K'_\lambda R_p))] d\lambda = \int_0^\infty [K'_\lambda L_\lambda^0 (1 - (1 - \exp(-K'_\lambda R_p)))] d\lambda \quad (21) $$

On trouve alors l'expression du coefficient d'émission nette en fonction de la température et du rayon du plasma :

$$ \epsilon_N(R_p, T) = \int_0^\infty L_\lambda^0(T) K'_\lambda(T) \exp(-K'_\lambda(T) \cdot R_p) d\lambda \quad (22) $$

Avec cette méthode, on peut ainsi considérer la divergence du flux radiatif dans l'équation de l'énergie grâce au terme source linéaire $S_{rad}$ défini par :

$$ S_{rad} = \vec{\nabla} \cdot \overline{q_{\lambda}} = 4\pi\epsilon_N \quad (23) $$

Cette méthode est bien adaptée au calcul de la puissance perdue par rayonnement par les régions chaudes du plasma, lorsque les gradients de température ne sont pas trop importants. Cependant, le rayonnement absorbé dans les régions tièdes et froides n'est pas pris en compte mais peut jouer un rôle important ou non dans le bilan d'énergie selon l'importance des phénomènes conducto-convectifs et turbulents dans ces zones. Concernant le paramètre $R_p$, le rayonnement étant absorbé en grande partie dans les deux premiers millimètres, sa valeur a peu d'incidence sur le résultat lorsqu'elle est choisie supérieure à 2 mm [7]. Cependant, une méthode simple de détermination du $R_p$ adéquat en fonction de son étude est de le choisir tel que 80% de la température maximale atteinte par le plasma soit contenue à l'intérieur de ce rayon : cela permet de prendre en compte la région où l'émission radiative est la plus intense. Diminuer ce rayon revient à considérer des pertes radiatives plus importantes et donc à un plasma plus froid.

## 2.3. Harmoniques Sphériques : la Méthode P-N

Reprenons l'Équation du Transfert Radiatif (6) :

$$ \frac{dL_{\lambda}(\vec{r}, \vec{s})}{ds} = K'(\lambda, T(\vec{r}, \vec{s})). (L_{\lambda}^{0}(T(\vec{r}, \vec{s})) - L_{\lambda}(T(\vec{r}, \vec{s}))) \quad (24) $$

Une des difficultés essentielles du transfert radiatif est de prendre en compte le champ directionnel de la luminance $L(\vec{r}, \vec{s})$ en tout point $P(r)$, sur tout le spectre et pour toute direction caractérisée par $\vec{s}(\theta, \varphi)$. L'approche du modèle P-1 (P-N à l'ordre 1), ou des harmoniques sphériques, consiste à décomposer $L(\vec{r}, \vec{s})$ sous la forme d'une double somme constituée d'une partie spatiale et d'une partie directionnelle. Cette nouvelle base complète de fonctions dépendant séparément de $\vec{r}$ et de $\vec{s}$ correspond à exprimer la luminance sous forme d'une série d'harmoniques sphériques sous la forme suivante [8] :

$$ L(\vec{r}, \vec{s}) = L(\vec{r}, \theta, \varphi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{m=l} A_l^m(\vec{r}) Y_l^m(\theta, \varphi) \quad (25) $$
Les fonctions $Y_l^m(\Omega)$ dépendant uniquement de l'angle solide $\Omega$ et donc des variables
angulaires $\theta$ et $\varphi$, s'expriment de manière connue sous la forme de polynômes de Legendre de
première espèce $P_l^m(\mu = \cos(\theta))$ :

$$
Y_l^m(\theta, \varphi) = \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} e^{im\varphi} P_l^m(\cos \theta) \quad (26)
$$

avec :

$$
P_l^m(\mu = \cos\theta) = (-1)^m \frac{(1-\mu^2)^{\frac{m}{2}}}{2^{l+m} l!} \frac{d^{l+m}}{d\mu^{l+m}} (\mu^2 - 1)^l \quad (27)
$$

De plus, $P_l^m$ et $P_l^{-m}$ sont liés par la relation suivante :

$$
P_l^{-m}(\mu) = (-1)^m \frac{(l-m)!}{(l+m)!} P_l^m(\mu) \quad (28)
$$

Les expressions des polynômes de Legendre $P_l^m(\mu = \cos \theta)$ pour $l$ allant de 0 à 3 sont
présentées dans le tableau 1.

<table><caption>**Tableau 1 :** Expression des polynômes de Legendre pour *l* allant de 0 à 3</caption>
  <thead>
    <tr>
      <th><i>l</i></th>
      <th><b>m=0</b></th>
      <th><b>m=1</b></th>
      <th><b>m=2</b></th>
      <th><b>m=3</b></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th><b>0</b></th>
      <td>1</td>
      <td>-</td>
      <td>-</td>
      <td>-</td>
    </tr>
    <tr>
      <th><b>1</b></th>
      <td>cos &theta;</td>
      <td>- sin &theta;</td>
      <td>-</td>
      <td>-</td>
    </tr>
    <tr>
      <th><b>2</b></th>
      <td>(3 cos<sup>2</sup> &theta; &minus; 1)/2</td>
      <td>3cos &theta; sin &theta;</td>
      <td>3sin<sup>2</sup> &theta;</td>
      <td>-</td>
    </tr>
    <tr>
      <th><b>3</b></th>
      <td>1/2(5 cos<sup>2</sup> &theta; &minus; 3) cos &theta;</td>
      <td>3/2(&minus;5 cos<sup>2</sup> &theta; + 1) sin &theta;</td>
      <td>15 cos &theta; sin<sup>2</sup> &theta;</td>
      <td>15 sin<sup>3/2</sup> &theta;</td>
    </tr>
  </tbody>
</table>

<caption>**Tableau 1 :** Expression des polynômes de Legendre pour *l* allant de 0 à 3</caption>

Dans la littérature, on note qu'il existe plusieurs façons de décrire ces expressions en fonction
de la manière dont on traite le cas de *m* négatif : l'équation (28) est parfois remplacée par de
simples valeurs absolues [9][10].

Les $A_l^m(\vec{r})$ sont des fonctions inconnues qui ne dépendent que de la position. Le principe de la méthode P-N consiste à identifier ces différentes fonctions à partir d'un système d'équations aux moments. Avec la méthode P-1, on arrête le développement à l'ordre 1. On reprend alors la double somme de l'équation (25) que l'on écrit avec $l=1$ :

$$
L'(\vec{r}, \theta, \varphi) = L'(\vec{r}, \theta, \varphi) = \sum_{l=0}^{1} \sum_{m=-1}^{m=1} A_l^m(\vec{r}) Y_l^m(\theta, \varphi) \quad (29)
$$

Donc :

$$
L'(\vec{r}, \theta, \varphi) = A_0^0 Y_0^0 + A_1^{-1} Y_1^{-1} + A_1^0 Y_1^0 + A_1^{-1} Y_1^{-1} \quad (30)
$$

Or, en utilisant les équations (27) et (28), ainsi que le tableau 1, nous déterminons les $Y_l^m$ :

$$
Y_0^0 = \sqrt{\frac{1}{4\pi}} ; \quad Y_1^{-1} = \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta ; \quad Y_1^0 = \sqrt{\frac{3}{4\pi}} \cos\theta ; \quad Y_1^1 = -\sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta
$$
$L'(\vec{r}, \theta, \varphi)$ devient donc :

$$L'(\vec{r}, \theta, \varphi) = A_0^0 \sqrt{\frac{1}{4\pi}} + A_1^{-1} \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta + A_1^0 \sqrt{\frac{3}{4\pi}} \cos\theta - A_1^1 \sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta \quad (31)$$

Il nous faut déterminer les $A_l^m(\vec{r})$ pour trouver $L'(\vec{r}, \theta, \varphi)$. Ces coefficients $A_l^m(\vec{r})$ s'expriment en fonction des moments de la luminance, eux même obtenus en multipliant l'ETR par des puissances de cosinus directeurs $l_i$ et en intégrant sur l'angle solide $\Omega$. En s'appuyant sur la figure 4, les cosinus directeurs de $\vec{s}$ sont : $l_1 = \cos\theta$, $l_2 = \sin\theta \cos\varphi$ et $l_3 = \sin\theta \sin\varphi$.

En reprenant les différentes expressions des moments de la luminance, on aboutit alors à un système dont les $A_l^m(\vec{r})$ sont les inconnues. Nous donnons la démonstration dans ce qui suit.

### a) <u>Moment d'ordre 0 de la luminance</u>

Le moment d'ordre 0 de la luminance $L_0$, correspond à la luminance moyenne :

$$L_0(\vec{r}) = \int_{\Omega=0}^{4\pi} L'(\vec{r}, \Omega) d\Omega = \int_{\theta=0}^{\pi} \int_{\varphi=0}^{2\pi} L'(\vec{r}, \theta, \varphi) \sin\theta \, d\theta \, d\varphi \quad (32)$$

$$L_0(\vec{r}) = \int_{\theta=0}^{\pi} \int_{\varphi=0}^{2\pi} \left( A_0^0 \sqrt{\frac{1}{4\pi}} + A_1^{-1} \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta + A_1^0 \sqrt{\frac{3}{4\pi}} \cos\theta - A_1^1 \sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta \right) \sin\theta \, d\theta \, d\varphi$$

$$L_0 = 4\pi A_0^0 \sqrt{\frac{1}{4\pi}} \qquad (33)$$

$$A_0^0 = L_0 \frac{\sqrt{\pi}}{2\pi} \qquad (34)$$

### b) <u>Moments d'ordres 1 de la luminance :</u>

Les moments d'ordre 1 ($L_1, L_2, L_3$), représentent les composantes du flux radiatif selon les différentes directions $x_i$ du système de coordonnées ($x_1, x_2$ et $x_3$).

$$L_i(\vec{r}) = \int_{\Omega=0}^{4\pi} l_i L'(\vec{r}, \Omega) d\Omega \quad \text{avec } (i = 1, 2, 3) \qquad (35)$$

On calcule les 3 différents moments d'ordre 1 :

$$L_1(\vec{r}) = \int_{\Omega=0}^{4\pi} \cos\theta \left( A_0^0 \sqrt{\frac{1}{4\pi}} + A_1^{-1} \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta + A_1^0 \sqrt{\frac{3}{4\pi}} \cos\theta - A_1^1 \sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta \right) d\Omega$$
D'où

$$L_1 = \frac{2\sqrt{3\pi}}{3} A_1^0 \Leftrightarrow A_1^0 = \frac{\sqrt{3\pi}}{2\pi} L_1 \qquad (36)$$

$$L_2(\vec{r}) = \int_{\Omega=0}^{4\pi} \sin\theta \cos\varphi \left( A_0^0 \sqrt{\frac{1}{4\pi}} + A_1^{-1} \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta + A_1^0 \sqrt{\frac{3}{4\pi}} \cos\theta - A_1^1 \sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta \right) d\Omega$$

D'où

$$L_2 = \frac{\sqrt{6\pi}}{3} A_1^{-1} - \frac{\sqrt{6\pi}}{3} A_1^1 \qquad (37)$$

$$L_3(\vec{r}) = \int_{\Omega=0}^{4\pi} \sin\theta \sin\varphi \left( A_0^0 \sqrt{\frac{1}{4\pi}} + A_1^{-1} \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta + A_1^0 \sqrt{\frac{3}{4\pi}} \cos\theta - A_1^1 \sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta \right) d\Omega$$

D'où

$$L_3 = \frac{-i\sqrt{6\pi}}{3} A_1^{-1} - \frac{i\sqrt{6\pi}}{3} A_1^1$$

Pour déterminer $A_1^{-1}$ et $A_1^1$, nous devons résoudre le système constitué de $L_2$ et $L_3$:

$$\left\{ \begin{array}{l} L_2 = \frac{\sqrt{6\pi}}{3} A_1^{-1} - \frac{\sqrt{6\pi}}{3} A_1^1 \\ \\ L_3 = \frac{-i\sqrt{6\pi}}{3} A_1^{-1} - \frac{i\sqrt{6\pi}}{3} A_1^1 \end{array} \right.$$

Ce système linéaire de deux équations à deux inconnues peut être résolu par la méthode de Cramer en passant par une formulation matricielle :

$$\begin{pmatrix} \frac{\sqrt{6\pi}}{3} & -\frac{\sqrt{6\pi}}{3} \\ -\frac{i\sqrt{6\pi}}{3} & -\frac{i\sqrt{6\pi}}{3} \end{pmatrix} \begin{pmatrix} A_1^{-1} \\ A_1^1 \end{pmatrix} = \begin{pmatrix} L_2 \\ L_3 \end{pmatrix}$$

Nous obtenons :

$$\left\{ \begin{array}{l} A_1^{-1} = \frac{3}{2\sqrt{6\pi}} (L_2 + iL_3) \\ \\ A_1^1 = \frac{-3}{2\sqrt{6\pi}} (L_2 - iL_3) \end{array} \right. \qquad (38)$$
Les expressions des coefficients $A_l^m(\vec{r})$ avec $l = 1$ (29) en fonction des moments de la luminance, en posant $\beta = \frac{1}{2\sqrt{\pi}}$, s'écrivent :

$$ egin{cases} A_0^0 = \beta L_0 \\ A_1^0 = \sqrt{3\beta} L_1 \\ A_1^{-1} = \sqrt{\frac{3}{2}} \beta (L_2 + iL_3) \\ A_1^1 = -\sqrt{\frac{3}{2}} \beta (L_2 - iL_3) \end{cases} \qquad (39) $$

Nous avons obtenu les expressions des $Y_l^m$ et des $A_l^m$, l'expression de la luminance (30) s'écrit :

$$ L'(\vec{r}, \theta, \varphi) = A_0^0 Y_0^0 + A_1^{-1} Y_1^{-1} + A_1^0 Y_1^0 + A_1^1 Y_1^1 \qquad (40) $$

$$ L'(\vec{r}, \theta, \varphi) = \beta L_0 \sqrt{\frac{1}{4\pi}} + \sqrt{\frac{3}{2}} \beta (L_2 + iL_3) \times \frac{1}{4} \sqrt{\frac{6}{\pi}} e^{-i\varphi} \sin\theta + \sqrt{3\beta} L_1 \sqrt{\frac{3}{4\pi}} \cos\theta \\ + \sqrt{\frac{3}{2}} \beta (L_2 - iL_3) \times \sqrt{\frac{3}{8\pi}} e^{i\varphi} \sin\theta \qquad (41) $$

Pour la méthode P-1, nous obtenons l'expression de la luminance suivante :

$$ L'(\vec{r}, \theta, \varphi) = \frac{1}{4\pi} L_0 + \frac{1}{4\pi} (3L_1 \cos\theta + \frac{1}{4\pi} (3L_2 \cos\varphi \sin\theta + 3L_3 \sin\varphi \sin\theta)) \\ \Leftrightarrow \\ L'(\vec{r}, \theta, \varphi) = \frac{1}{4\pi} (L_0 + 3L_1 \cos\theta + 3L_2 \cos\varphi \sin\theta + 3L_3 \sin\varphi \sin\theta) \qquad (42) $$

Appliquons nos résultats à la résolution de l'ETR par la méthode P-1 dans le cas 1D d'un plan infini dont la configuration est présentée sur la figure 4 [8].

<figure><img src="image_5.png" /><figcaption>**Figure 4** : Description de la luminance dans le système de coordonnées</figcaption></figure>
Dans ce cas, l'ETR selon $x_1$ s'écrit :

$$ \cos \theta \frac{dL'}{dx} = K'L^0 - K'L' \quad (43) $$

Or, d'après la démonstration précédente, la luminance $L'(x, \theta)$ peut se développer en série d'harmoniques sphériques en simplifiant l'équation (42) sous une seule dimension.

On obtient alors :

$$ L'(x, \theta) = \frac{1}{4\pi} (L_0 + 3L_1 \cos \theta) \quad (44) $$

avec :

$$ L_0 = \int_{\Omega=0}^{4\pi} L'(x, \theta) d\Omega \quad (45) $$

Représente le rayonnement incident $\vec{G}_\lambda$ suivant x (27)

$$ L_1 = \int_{\Omega=0}^{4\pi} \cos \theta L'(x, \theta) d\Omega \quad (46) $$

Représente le flux radiatif $\vec{q}_r$ suivant x (28) avec $l_i = \cos \theta$)

Pour obtenir la luminance $L'(x, \theta)$, nous devons déterminer les deux moments inconnus $L_0$ et $L_1$ à partir de l'ETR. En la multipliant par 1 et par le cosinus directeur de la luminance $\cos \theta$, puis en intégrant sur l'angle solide $\Omega$, on obtient le système d'équations suivant :

$$ (\times 1) \quad \int_{\Omega=0}^{4\pi} \cos \theta \frac{dL'}{dx} d\Omega + \int_{\Omega=0}^{4\pi} K'L' d\Omega = \int_{\Omega=0}^{4\pi} K'L^0 d\Omega \quad (47) $$

$$ (\times \cos \theta) \quad \underbrace{\int_{\Omega=0}^{4\pi} \cos \theta \frac{dL'}{dx} \cos \theta \, d\Omega}_{\frac{dL_{11}}{dx}} + \int_{\Omega=0}^{4\pi} K'L' \cos \theta \, d\Omega = \int_{\Omega=0}^{4\pi} K'L^0 \cos \theta \, d\Omega \quad (48) $$

Nous remarquons que le premier terme de la deuxième équation (48) du système correspond à la dérivée du moment d'ordre 2 de la luminance, d'après sa définition :

$$ L_{11} = L_{ij}(\vec{r}) = \int_{\Omega=0}^{4\pi} l_i l_j L'(\vec{r}, \Omega) d\Omega = \int_{\Omega=0}^{4\pi} \cos^2 \theta L'(x, \theta) d\Omega \quad \text{avec } (i, j = 1, 2, 3) \quad (49) $$

Finalement, le système d'équations des moments à résoudre s'écrit :

$$ \left\{ \begin{array}{l} \frac{dL_1}{dx} + K'L_0 = 4\pi K'L^0 \\ \frac{dL_{11}}{dx} + K'L_1 = 0 \end{array} \right. \quad (50) $$
En partant de la définition du moment d'ordre 2 (49) et en prenant en compte l'expression développée en harmonique sphérique de la luminance (44), nous obtenons :

$$L_{11} = \frac{1}{4\pi} \left( \int_{\Omega=0}^{4\pi} L_0 \cos^2 \theta \, d\Omega + \int_{\Omega=0}^{4\pi} 3L_1 \cos^3 \theta \, d\Omega \right) \quad (51)$$

En remplaçant l'angle solide élémentaire $d\Omega$ par son expression $d\Omega = \sin \theta \, d\theta \, d\varphi$, on a :

$$L_{11} = \frac{1}{4\pi} \left( \int_{\varphi=0}^{2\pi} \int_{\theta=0}^{\pi} L_0 \cos^2 \theta \sin \theta \, d\theta \, d\varphi + \int_{\varphi=0}^{2\pi} \int_{\theta=0}^{\pi} 3L_1 \cos^3 \theta \sin \theta \, d\theta \, d\varphi \right) \quad (52)$$

Nous obtenons :

$$L_{11} = \frac{1}{3} L_0 \quad (53)$$

Enfin, en substituant cette expression de $L_{11}$ dans le système des équations au moment, on résout le système d'équation suivant :

$$\left\{ \begin{aligned} \frac{dL_1}{dx} + K'L_0 &= 4\pi K' L^0 \\ \frac{1}{3} \frac{dL_0}{dx} + K'L_1 &= 0 \end{aligned} \right. \qquad (54)$$

⇔

$$\left\{ \begin{aligned} \frac{1}{3K'^2} \frac{d^2 L_0}{dx^2} - L_0 &= -4\pi L^0 \\ \frac{1}{K'^2} \frac{d^2 L_1}{dx^2} - 3L_1 &= \frac{4\pi}{K'} \frac{dL^0}{dx} \end{aligned} \right. \quad (55)$$

Nous remarquons que ce nouveau système d'équation et un système d'équations différentielles du second ordre avec termes sources, beaucoup plus simple à résoudre.

De plus, en pratique, seule l'équation différentielle (55) portant sur la variable $L_0$ devra être résolue. En effet, une fois le rayonnement incident $L_0$ déterminé, il est facile d'en déduire le flux radiatif $L_1$ en l'isolant dans la relation (54) :

$$\frac{1}{3} \frac{dL_0}{dx} + K' L_1 = 0 \quad (56)$$

⇔

$$L_1 = -\frac{1}{3K'} \frac{dL_0}{dx} \quad (57)$$

Nous généralisons ces résultats (55), en notant $L_0 = \vec{G}_\lambda$ la radiation incidente et $L_1 = \vec{q}_\lambda$ le flux radiatif pour une longueur d'onde donnée :

$$\vec{q}_\lambda(\vec{r}) = \frac{-1}{3K'(\lambda, T(\vec{r}))} \vec{\nabla} \cdot \vec{G}_\lambda(\vec{r}) \quad (58)$$

De plus, en généralisant également l'équation (54), nous obtenons :

$$\vec{\nabla} \cdot \vec{q}_{\lambda}(\vec{r}) = K'(\lambda, T(\vec{r})) (4\pi L_{\lambda}^{0}(\vec{r}) - G_{\lambda}(\vec{r})) \quad (59)$$
Finalement, nous obtenons :

$$ \vec{\nabla} \frac{-1}{3K'(\lambda, T(\vec{r}))} \vec{\nabla} G_{\lambda}(\vec{r}) = K'(\lambda, T(\vec{r})) (4\pi L_{\lambda}^{0}(\vec{r}) - G_{\lambda}(\vec{r})) \quad (60) $$

Cette dernière équation (60) peut se mettre sous la forme généralisée de Patankar [4] :

$$ \frac{\partial}{\partial x_i} \left( -\Gamma_k \frac{\partial \varphi_k}{\partial x_i} \right) = S_{\varphi_k} \quad (61) $$

$$ \text{En posant : } \varphi_k = G_\lambda, \quad \Gamma_k = \frac{1}{3K'(\lambda)} \text{ et } S_{\varphi_k} = K'(\lambda) (4\pi L_\lambda^0 - G_\lambda(\vec{r})) $$

L'intégration sur toutes les longueurs d'onde demande de longs temps de calcul : nous pouvons alors découper le spectre en groupes de longueurs d'ondes et calculer des Coefficients Moyens d'Absorption (CMAs) définis sur différents intervalles spectraux, limitant ainsi le nombre d'équations à résoudre au nombre d'intervalles spectraux choisi.

Les pertes radiatives se mettent alors sous la forme suivante [4] :

$$ S_{\text{rad}} = \vec{\nabla} \cdot \vec{q}_{\lambda} = \sum_{l=0}^{N\text{bandes}} \vec{\nabla} \cdot \vec{q}_{l} = \sum_{l=0}^{N\text{bandes}} \bar{K}_{l} (4\pi L_{l}^{0} - G_{l}(\vec{r})) \quad (62) $$

Où $L_l^0$ est l'intégration de la fonction de Planck sur l'intervalle spectral $l$ et $\bar{K}_l$ le CMA de l'intervalle $l$. Ainsi, sous cette forme, l'équation peut directement être injectée dans l'équation de l'énergie en tant que terme source. On aura cependant à résoudre autant d'équation (61) que d'intervalles définis.

## 2.4. Méthode des Ordonnées Discrètes

La Méthode des Ordonnées Discrètes (DOM) est basée sur la discrétisation de l'Equation du Transfert Radiatif (ETR) sur un certain nombre $N_{dir}$ de directions discrètes $s_i(\mu_i, \eta_i, \xi_i)$ contenu dans l'angle solide $\Omega = 4\pi$ et chacune associée à son poids $\omega_i$, s'interprétant comme la contribution de chaque sous-espace angulaire à la divergence totale. On peut alors utiliser différents types de quadratures comme les harmoniques sphériques (méthode $S_N$), en choisissant l'ordre afin d'indiquer le nombre de directions angulaires utilisées pour l'intégration : plus l'ordre est élevé, plus fine sera la résolution angulaire. En effet, la précision du calcul sera linéairement impactée par une augmentation de cet ordre, contrairement à la méthode P-N où le gain en véracité est non-linéaire. Cependant, un compromis doit être trouvé entre la précision du calcul et la durée de celui-ci [11]. Il est de coutume de choisir un ensemble de directions et de poids symétriques (invariant suivant une rotation de 90°) et satisfaisant les équations aux moments d'ordre 0, 1 et 2 [12] :

$$ \int_{4\pi} d\Omega = 4\pi = \sum_{i=1}^{n} \omega_i \quad (63) $$

$$ \int_{4\pi} \vec{s}_l d\Omega = \vec{0} = \sum_{i=1}^{n} \omega_i \vec{s}_l \quad (64) $$

$$ \int_{4\pi} \vec{s}_l \cdot \vec{s}_l d\Omega = \frac{4\pi}{3} \vec{\delta} = \sum_{i=1}^{n} \omega_i \vec{s}_l \cdot \vec{s}_l \quad (65) $$
Où $\vec{\delta}$ est le tenseur unitaire. Des ensembles de directions et de poids satisfaisant ces critères ont été tabulé, par exemple, par Lathrop et Carlson [13].

La quadrature $S_4$ est souvent utilisé car maximisant dans la plupart des cas le compromis entre précision et temps de calcul, équivalente à $N_{dir} = N(N + 2) = 24$ directions. L'ETR est alors résolue pour chacune des directions discrètes $\vec{s}_i$ par une approche de volumes finis. L'intégration de l'ETR sur le volume $V$ d'un élément limité par sa surface $\Sigma$, de vecteur unitaire normal $\vec{n}$ donne [11] :

$$ \int_{\Sigma} L_{\lambda}(\vec{s}_i) \cdot \vec{n} \, d\Sigma = \int_{V} (K'_{\lambda} L_{\lambda}^{0} - K'_{\lambda} L_{\lambda}(\vec{s}_i)) dV \quad (66) $$

Le domaine est discrétisé en volumes de contrôle, pouvant être réguliers ou non.

En prenant $L_j$ l'intensité moyenne de la face $j$ associée à son centre et à sa surface $A_j$, et $L_p$ l'intensité moyenne sur le volume $V$ associé au centre de la cellule représenté par un point P, on obtient alors [14] :

$$ \sum_{j=1}^{N_{faces}} L_j(\vec{s}_i \cdot \vec{n}_j) A_j = K'_{\lambda} V (L_{\lambda}^{0} - L_p) \quad (67) $$

Où $L_j$ est considéré constant sur chaque face et le produit scalaire entre le $i$-ème vecteur directeur discret et la $j$-ème face de la cellule considérée, est défini par :

$$ D_{ij} = \vec{s}_i \cdot \vec{n}_j = \mu_i n_{xj} + \eta_i n_{yj} + \xi_i n_{zj} \quad (68) $$

Pour résoudre cette équation (67), un schéma de discrétisation spatiale basé sur le flux moyen (DMFS) proposé par Ströhle et al. [15] peut être par exemple utilisé. Ce schéma repose sur la décomposition en deux termes de l'intensité moyenne $L_p$ sous la forme suivante :

$$ L_p = \alpha \overline{L_{out}} + (1 - \alpha) \overline{L_{in}} \quad (69) $$

avec $\alpha \in [0,1]$ et $\overline{L_{in}}$ est la moyenne des intensités aux faces entrantes de la cellule, $\overline{L_{out}}$ celle aux faces sortantes, soit :

$$ \overline{L_{in}} = \frac{\sum_{j, D_{ij}<0} D_{ij} A_j L_j}{\sum_{j, D_{ij}<0} D_{ij} A_j} \quad (70) $$

$$ \overline{L_{out}} = \frac{\sum_{j, D_{ij}>0} D_{ij} A_j L_j}{\sum_{j, D_{ij}>0} D_{ij} A_j} \quad (71) $$

Ainsi, en substituant $\overline{L_{in}}$ de l'équation (69) dans l'équation (67), nous obtenons, après quelques opérations algébriques, l'expression de $L_p$ suivante [14] :

$$ L_p = \frac{\alpha V K'_{\lambda} L^{0} - \theta \sum_{j, D_{ij}<0} D_{ij} A_j L_j}{\alpha V K'_{\lambda} + \sum_{j, D_{ij}>0} D_{ij} A_j} \quad (72) $$

Avec :

$$ \theta = \alpha - (1 - \alpha) \frac{\sum_{j, D_{ij}>0} D_{ij} A_j}{\sum_{j, D_{ij}<0} D_{ij} A_j} \quad (73) $$
Si on considère que le produit scalaire $D_{ij}$ est le terme de transformation de la projection de la surface $A_j$, suivant le vecteur $\vec{s}_i$, sur le plan orthogonal à ce dernier, on peut écrire :

$$ \sum_{j, D_{ij}>0} D_{ij} A_j = - \sum_{j, D_{ij}<0} D_{ij} A_j \qquad (74) $$

Le terme $\Theta$ devient alors égal à 1 et on obtient :

$$ L_p = \frac{\alpha V K_{\lambda}' L^0 - \sum_{j, D_{ij}<0} D_{ij} A_j L_j}{\alpha V K_{\lambda}' + \sum_{j, D_{ij}>0} D_{ij} A_j} \qquad (75) $$

La radiation incidente $G_{\lambda}$ est évaluée au centre de la cellule P par :

$$ G_{\lambda} = \sum_{i=1}^{N_{dir}} \omega_i L_p(\vec{s}_i) \qquad (76) $$

On en déduit le terme source $S_{rad}$ tel que :

$$ S_{rad} = \vec{\nabla} \cdot \overline{q_{\lambda}} = K_{\lambda}' (4\pi L_{\lambda}^0 - \sum_{i=1}^{N_{dir}} \omega_i L_p(\vec{s}_i)) \qquad (77) $$

Pour exemple, considérons la colonne cylindrique de gaz décrite sur la figure 5 [16]. Cette figure illustre en 3D les pertes par rayonnement dans une direction donnée AW (définie par les angles $\theta$ et $\varphi$), pour un point A situé dans la coupe transversale de la colonne de gaz, ayant pour centre C. On calcule la divergence du flux radiatif pour des points le long du segment CG, tel que A, c'est-à-dire du centre vers le bord de la colonne de gaz. W est un point à la frontière du cylindre (représenté à $x_W$), $x_j$ dans l'équation (78) représente un point sur le segment AW, tel que D, et $L_{\lambda}^0(0)$ la fonction de Planck au point cible A ($x_j = 0$). $\xi_j$ n'est qu'une variable intermédiaire utilisée pour l'intégration.

<figure><img src="image_6.png" /><figcaption>**Figure 5** : Coupe transversale de la colonne de plasma cylindrique</figcaption></figure>
En utilisant la méthode DOM, l'angle solide de $4\pi$ en un point cible A est ainsi divisé en directions ordonnées, et une équation différentielle partielle est générée pour décrire le transfert radiatif dans chaque direction. Ainsi, l'intégration sur l'ensemble de l'angle solide est convertie en une sommation pondérée basée sur la quadrature numérique choisie. Nous pouvons écrire [16] :

$$S_{rad} = K'_{\lambda} \sum_{i=1}^{n} \omega_{i} \left[ L_{\lambda}^{0}(0) \exp \left( -\int_{0}^{x_{w,j}} K'_{\lambda}(x_{j}) dx_{j} \right) - \int_{0}^{x_{w,j}} \left[ L_{\lambda}^{0}(x_{j}) - L_{\lambda}^{0}(0) \right] K'_{\lambda}(x_{j}) \exp \left( -\int_{0}^{x_{j}} K'_{\lambda}(\xi_{j}) d\xi_{j} \right) dx_{j} \right] \quad (78)$$

Où $i$ représente une direction parmi les $n$, et $\omega_i$ le poids associé à cette direction $i$.

### 3. Application au cas d'un Arc Libre à dans l'Argon à la Pression Atmosphérique.

Nous avons étudié théoriquement trois différentes méthodes de calcul des pertes de chaleur par rayonnement utilisables dans les modèles. Afin de visualiser leurs différences nous allons appliquer deux d'entre elles, en utilisant deux logiciels de modélisation différents (COMSOL Multiphysics et ANSYS Fluent) sur un cas classique de modélisation du plasma thermique : l'arc libre.

#### 3.1. Description du Modèle

Pour comparer et valider notre modèle, nous avons repris les hypothèses et conditions aux limites utilisées par Fulbert Baudoin [1] dans sa thèse. Un arc libre dans l'argon à 200 A est modélisé selon les conditions suivantes :

- Le plasma est en équilibre thermodynamique local ;

- Le modèle est stationnaire et en 2D axisymétrique ;

- Le plasma est un fluide visqueux, laminaire, newtonien et compressible ($Ma<0.3$)

- La gravité est négligée devant les forces de pression ;

- Les phénomènes physiques dans les zones anodique et cathodique ne sont pas pris en compte, la modélisation concerne uniquement la colonne de l'arc.

Après validation avec notre modèle de référence de F. Baudoin, ainsi qu'avec l'expérience que lui-même cherchait à modéliser de K.C. Hsu [17], on compare alors les résultats en température obtenus en fonction de la méthode de calcul du rayonnement (CEN et P-1) et des logiciels utilisés (COMSOL Multiphysics et ANSYS Fluent).

#### 3.2. Résultats sous COMSOL Multiphysics

L'argon étant un gaz non moléculaire dans nos conditions, l'absorption du rayonnement sur les bords du plasma se fait majoritairement par le continuum atomique, ce qui explique le peu de différences observées entre le CEN et le P-1 pour la température radiale (Figure 7.a) et la divergence du flux radiatif (Figure 7.b) sur ce rayon. On observe toutefois des températures qui diminuent plus rapidement en s'éloignant de l'arc lorsque les pertes radiatives sont modélisées par la méthode P-1.
<figure><img src="image_7.png" /><figcaption>**Figure 6** : Profil de température (K) de l'arc libre à 200 A dans l'argon : comparaison entre les méthodes du CEN et du P-1 avec R<sub>p</sub> = 1 mm</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7** : Divergence du flux radiatif (W·m<sup>-3</sup>)
(a) et températures radiales (K) (b) selon trois coupes à 2, 5 et 9 mm sous la cathode</figcaption></figure>

3.3. Résultats sous ANSYS Fluent

Avec le logiciel ANSYS Fluent et pour des conditions d'études similaires, les différences entre les deux méthodes sont encore plus faibles, que ce soit sur le profil de température radial (Figure 9.a) ou celui de la divergence du flux radiatif (Figure 9.b). Les valeurs obtenues pour des coupes radiales à 2 mm, 5 mm et 9 mm au-dessous de la cathode sont quasiment superposées sur l'intégralité du rayon.
<figure><img src="image_9.png" /><figcaption>**Figure 8 :** Profil de température (K) de l'arc libre à 200 A dans l'argon. Comparaison entre les méthodes du CEN et du P-1 avec $R_p = 1 \text{ mm}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9 :** Divergence du flux radiatif ($W \cdot m^{-3}$) (a) et températures radiales (K) (b) selon trois coupes à 2, 5 et 9 mm sous la cathode</figcaption></figure>
Conclusion

Le rayonnement est un aspect important de la physique des plasmas thermiques, notamment
dans sa modélisation. Cependant, il existe différentes méthodes afin de calculer les pertes
radiatives, chacune possédant ses avantages et inconvénients. Choisir une méthode plutôt
qu'une autre dépend alors du cas étudié ou d'un compromis entre la précision souhaitée et temps
de calcul.

Dans le cas d'un arc libre dans l'argon, pour des calculs souhaités rapides et acceptables, la méthode du CEN semble donc suffisante. Les méthodes P-1 et DOM quant à elles demandent plus de pré- traitements et de ressources : la connaissance des coefficients moyens d'absorption qui requiert un découpage optimal en intervalles spectraux du spectre radiatif, l'utilisation de fonctions moyennes adaptées en température, pression et taille du plasma, un nombre d'équations à résoudre qui augmente avec le nombre d'intervalle spectraux et tend par conséquent à augmenter les temps de calcul et étudier de façon plus méticuleuse les solveurs proposés.

L'utilisation de logiciels différents appuie cette observation. En effet, les deux logiciels aboutissent à la même conclusion, avec des temps de calcul comparables, malgré une différence fondamentale dans le principe de résolution des équations : COMSOL Multiphysics utilise la méthode des éléments finis, tandis que ANSYS Fluent repose sur celle des volumes finis.

De façon générale, COMSOL Multiphysics est particulièrement adapté aux géométries complexes nécessitant des maillages non structurés. Il présente l'avantage d'une interface intuitive et d'un environnement multiphysique intégré, qui facilite le couplage automatique entre différentes physiques, ce qui en fait un outil bien adapté pour débuter. En revanche, l'utilisateur a peu de contrôle direct sur les équations résolues, et le solveur peut parfois rencontrer des problèmes de convergence, notamment dans le cas de systèmes fortement non linéaires ou couplés.

À l'inverse, ANSYS Fluent repose sur une méthode de résolution conservant les flux
localement dans chaque cellule, ce qui le rend mieux adapté à la résolution des équations de
Navier-Stokes. Le logiciel est en revanche moins intuitif à prendre en main et nécessite une
maîtrise du codage pour intégrer des développements numériques personnalisés (UDF ou
UDS). Le couplage multiphysique y est moins automatisé, mais l'utilisateur garde un contrôle
fin sur les équations, ce qui permet, à terme, une meilleure robustesse du solveur, à condition
de bien en maîtriser les paramètres.

Références

[1] F. Baudoin, « Contribution à la modélisation d'un arc électrique dans les appareils de coupure basse tension », Thèse de doctorat de l'Université Blaise Pascal, Clermont-Ferrand II, 2004.
https://theses.hal.science/tel-00011591/file/2004_Baudoin.pdf

[2] H. Riad, « Calcul du transfert radiatif dans des arcs et des plasmas thermiques : application à l'hydrogène et au méthane », Thèse de doctorat de l'Université Paul Sabatier - Toulouse III, 1996.
https://theses.fr/1996TOU30178

[3] Y. Cressault, “Propriétés des plasmas thermiques dans des mélanges argon-hydrogène-cuivre”, Thèse de doctorat de l’Université Paul Sabatier - Toulouse III, 2001.
<https://theses.hal.science/tel-00278927>

[4] N. Kabbaj, Etude du transfert radiatif d'un plasma thermique d'air : influence des propriétés radiatives dans la modélisation d'un arc libre, Thèse de doctorat de l'Université Paul Sabatier - Toulouse III, 2019.
<https://theses.hal.science/tel-02936544v1>

[5] J.J. Lowke, (1974), "Predictions of arc temperature profiles usinf approximate emission coefficients for radiation losses", J. Quant. Spectrosc. Radiat. Transfer. Vol. 14, pp. 111-122. Pergamon Press.
<https://doi.org/10.1016/0022-4073(74)90004-1>

[6] R.W. Liebermann, J.J. Lowke, (1976), "Radiation emission coefficients for sulfur hexafluoride arc plasmas", J. Quant. Spectrosc. Radiat. Transfer, Vol. 16, pp. 253-264. Pergamon Press
<https://doi.org/10.1016/0022-4073(76)90067-4>

[7] J.-M. Bauchire, « Modélisation numérique d'une torche de projection à plasma : influence de la géométrie et de la turbulence sur les propriétés du plasma », Thèse de doctorat de l'Université Paul Sabatier - Toulouse III, 1997.

[8] H.-Q Nguyen, « Modélisation et simulation du remplissage de moules verriers : Prise en compte du transfert radiatif », Thèse de doctorat de l'Institut National Polytechnique de Lorraine, Laboratoire d'Energétique et de Mécanique Théorique et Appliquée Nancy-Université, CNRS UMR, 2009
<https://theses.hal.science/tel-01748913>

[9] R. Siegel, J. Howell, "Thermal radiation heat transfer – 4th edition", Taylor & Francis, 667-672, 2002

[10] M.F. Modest, "Radiative heat transfer", McGraw-Hill International Editions - Mechanical Engineering Series, pp.503-519, 1993.

[11] K.A. Jensen et al., (2007) "On various modeling approaches to radiative heat transfer in pool fires", Combustion and Flame, pp.263-279.
<https://doi.org/10.1016/j.combustflame.2006.09.008>

[12] M.F. Modest, Radiative heat transfer, McGraw-Hill International Editions - Mechanical Engineering Series, pp.541-546, 1993.

[13] K.D. Lathrop, B. G. Carlson, "Discrete-ordinates angular quadrature of the neutron transport equation, Technical information series report", Los Alamos Scientific Laboratory, 1965.

[14] D. Joseph et al., (2005); "Comparison of three spatial differencing schemes in discrete ordinates method using three-dimensional unstructured meshes", International Journal of Thermal Sciences, pp.851-864.
<https://doi.org/10.1016/j.ijthermalsci.2005.02.010>

[15] J. Ströhle, U. Schnell, K.R.G. Hein, (2001), "A mean flux discrete ordinates interpolation scheme for general co-ordinates", Int. Conference on Heat Transfer, Antalya, vol. 3, 2001.
<https://doi.org/10.1615/ICHMT.2001.RadiationSymp.80>

[16] S. Esmaeili, P. Kloc, J. D. Yan, (2022), "Spectral and Spatial Characteristic Behaviour of Three-dimensional Radiation Transfer in SF6 Switching Arcs", J. Phys. D : Appl. Phys, Vol 55, p 185502.
<https://doi.org/10.1088/1361-6463/ac4f94>

[17] K.C. Hsu, K. Etemadi, E. Pfender, (1983), "Study of the Free-Burning High-Intensity Argon Arc", Journal of Applied Physics, Volume 54, Issue 3, pp.1293-1301.
<https://doi.org/10.1063/1.332195>