\documentclass[XUPS,XML,SOM,Unicode,francais,NoFloatCountersInSection]{cedram}
\setcounter{tocdepth}{2}
%\XUPScorrections
\usepackage{xups15-03}
\graphicspath{{xups15-03_figures}}

\hyphenation{déri-vée déri-vées entro-pie entro-pies appli-qué varia-ble}

\newcommand{\PhG}{\relax{Ph. G}}

\begin{DefTralics}
\newcommand{\PhG}{Ph. G}
\end{DefTralics}

\begin{document}
\frontmatter

\title{Modélisation fluide des plasmas dans\nobreakspace les\nobreakspace Tokamaks}

\author[\initial{J.} \lastname{Blum}]{\firstname{Jacques} \lastname{Blum}}
\address{Laboratoire J.A. Dieudonné,
Université de Nice Sophia Antipolis,
Parc Valrose,
06108 Nice Cedex 02, France et Inria, Centre de recherche Sophia Antipolis Méditerranée, Équipe Projet CASTOR}
\email{jblum@unice.fr}

\thanks{Journées X-UPS 2015. Des problèmes à $N$ corps aux Tokamaks. Éditions de l'École polytechnique, 2015}

\begin{abstract}
Ce texte présente dans un premier temps les deux approches pour la modélisation d'un plasma (gaz ionisé) en présence d'un champ magnétique: l'approche microscopique basée sur les équations cinétiques (Vlasov, Boltzmann,...) et l'approche macroscopique basée sur les équations de la magnétohydrodynamique (MHD) et comment l'on passe de l'une à l'autre. Pour une machine du type Tokamak qui vise à réaliser la fusion de noyaux d'atomes légers par le principe du confinement magnétique, on montrera qu'il existe des échelles de temps allant de la microseconde (temps d'Alfven) à la seconde (temps de diffusion résistive) et que le modèle peut être simplifié en se plaçant à l'échelle de temps du phénomène qu'on se propose d'étudier. On établira ainsi le principe de l'évolution quasi-statique de l'équilibre du plasma dans un Tokamak et on montrera comment, avec des méthodes numériques appropriées, on peut reconstruire en temps réel l'équilibre du plasma à chaque instant de la décharge (lignes de flux, surfaces magnétiques, densité de courant, frontière libre du plasma,...). Ce problème rentre dans la catégorie des problèmes inverses, qui sont mal-posés au sens d'Hadamard et on montrera comment, à l'aide de techniques mathématiques de régularisation, on peut les rendre bien-posés. On pourra alors aborder le problème du contrôle de l'équilibre du plasma en temps réel, tel qu'il se pose dans chaque Tokamak, en utilisant des bases hilbertiennes associées à l'opérateur aux dérivées partielles qui régit le phénomène. Ce dispositif sera utilisé dans le Tokamak WEST en cours de montage au CEA à Cadarache, qui vise à tester le divertor en tungstène qui sera réalisé pour ITER.
\end{abstract}

\maketitle
\tableofcontents
\mainmatter

\section{Introduction}

La MHD (magnétohydrodynamique) est la description mono-fluide d'un plasma (gaz ionisé) placé dans un champ électromagnétique. Nous allons présenter ici la description bi-fluide (électrons et ions) connue sous le nom de modèle de Braginskii, qui est la base des processus de transport qui ont lieu dans un plasma collisionnel. Nous montrerons qu'à l'échelle de temps de la diffusion résistive, le plasma est en équilibre à tout instant.

Nous allons présenter l'équation de Grad-Shafranov dans la partie~3 pour l'équilibre dans une configuration axisymétrique par l'introduction d'un flux poloïdal qui est la principale inconnue de notre problème. Dans la partie 4, nous présentons le problème inverse qui permet de déterminer à chaque instant la configuration d'équilibre et en particulier la densité de courant toroïdal du plasma à partir de mesures expérimentales (mesures magnétiques, polarimétrie, interférométrie, pression cinétique, effet Stark...). Dans la partie 5, nous montrons comment les équations de transport peuvent être réécrites de façon consistante avec les équations d'équilibre. Le fait que la \hbox{vitesse} de diffusion tangentielle est beaucoup plus grande que la vitesse orthogonale nous permet de dire que les densités et températures électroniques et ioniques peuvent être supposées constantes sur chaque ligne de flux.

En moyennant les équations de transport sur chaque surface magné\-tique, on obtient un ensemble d'équations de diffusion paraboliques 1D pour ces densités et températures, la variable d'espace étant n'importe quel label de la ligne de flux. Cet ensemble d'équations modélise le transport entre lignes de flux et permet d'avoir un ensemble d'équations de diffusion 1D avec des coefficients géométriques calculés à partir des équations d'équilibre. Ceci permet de suivre l'évolution quasi-statique de la configuration d'équilibre de façon consistante avec les équations de transport.

\section{Les équations de la magnéto-hydrodynamique}
Un plasma est un gaz ionisé composé d'ions et d'électrons. Les équations cinétiques décrivent le plasma à l'aide d'une fonction de distribution $f_\alpha(\bmx, \bmv,t)$ (avec $\alpha=e$ pour les électrons et $\alpha=i$ pour les ions), où $\bmx$ est la position du point, et $\bmv$ la vitesse des particules. L'équation cinétique d'un plasma collisionnel est l'équation de Fokker-Planck
\begin{equation}\label{FP}
\frac{\partial f_\alpha}{\partial t}+
(\bmv\cdot\nabla_{\bmx})f_\alpha+\frac{\bF_\alpha}{m_\alpha}\cdot\nabla_{\bmv} f_\alpha=C_\alpha
\end{equation}
où $m_\alpha$ est la masse des particules et $\bF_\alpha$ la force appliquée aux particules et $C_\alpha$ est le terme dû aux collisions entre particules.
Cette approche microscopique requiert la résolution d'une équation aux dérivées partielles
6D ($6$-dimensionnelle) en espace-vitesse plus le temps. Cette équation est difficile à résoudre pour des raisons computationnelles. On dérive à partir de cette équation une approche macroscopique basée sur les équations fluides de la façon suivante. On définit la densité des particules par:
$$
n_\alpha(\bmx,t)=\int f_\alpha(\bmx,\bmv,t)d\bmv,
$$
la vitesse du fluide par:
$$
\bmu_\alpha(\bmx,t)=\frac{1}{n_\alpha}\int f_\alpha(\bmx,\bmv,t)\bmv d\bmv,
$$
et le tenseur pression par:
$$
P_\alpha (\bmx,t)=m_\alpha\int f_\alpha (\bmx,\bmv,t)(\bmv-\bmu_\alpha)(\bmv-\bmu_\alpha)d\bmv,
$$
qui dans le cas isotrope devient:
$$
p_\alpha(\bmx,t)=\frac{m_\alpha}{3}\int f_\alpha(\bmx,\bmv,t)(\bmv-\bmu_\alpha)^2d\bmv.
$$
En multipliant l'équation (\ref{FP}) par une fonction test $\phi(\bmv)$ et en intégrant sur l'espace des vitesses, on obtient les équations fluides. Le~premier moment (correspond à $\phi=1$)
nous donne l'équation de la densité des particules:
$$
\frac{\partial n_\alpha}{\partial t}+\nabla\cdot\int f_\alpha\bmv
d\bmv-\frac{1}{m_\alpha}\int \frac{\partial \bF_\alpha}{\partial \bmv}f_\alpha d\bmv=0.
$$
Comme ${\sfrac{\partial\bF_\alpha}{\partial \bmv}=0}$ pour les forces électromagnétiques et comme les collisions ne changent pas le nombre de particules, on obtient:
$$
\frac{\partial n_\alpha}{\partial t}+\nabla\cdot(n_\alpha \bmu_\alpha)=0
$$
Le second moment est obtenu en prenant $\phi=m_\alpha \bmv$, et on obtient l'équation de la quantité de mouvement:
$$
m_\alpha \frac{\partial}{\partial t} (n_\alpha\bmu_\alpha)+m_\alpha
\nabla_x\cdot\int f_\alpha \bmv\bmv d\bmv-\int \nabla_{\bmv}\cdot(\bF_\alpha\cdot\bmv)f_\alpha
d\bmv=\int m_\alpha \bmv C_\alpha d\bmv
$$
où l'on pose:
$\disp{\bmv=(\bmv-\bmu_\alpha)+\bmu_\alpha}$.
À l'aide de l'équation de conservation de densité, on obtient:
$$
m_\alpha n_\alpha \Bigl(\frac{\partial \bmu_\alpha}{\partial
t}+\bmu_\alpha\cdot\nabla \bmu_\alpha \Bigr)=-\nabla.P_\alpha +n_\alpha \overline\bF_\alpha
+{\bR_\alpha }
$$
avec $\overline\bF_\alpha = Ze (\bE+\bmu_\alpha\times \bB)$ où $Ze$ est la charge des particules et ${\bR_\alpha}$ est le taux de change de la quantité de mouvement dû aux collisions.

Le troisième moment nous donne l'équation d'énergie à laquelle il faut ajouter des relations de fermeture sur le flux de chaleur qui viennent d'un modèle de transport.

Les équations de la magnétohydrodynamique à un fluide s'obtiennent en définissant la densité de masse
\begin{align*}
\rho&=m_en_e+m_i n_i\\
&=m_e Z{n_i}+m_in_i \approx m_i n_i,
\end{align*}
la vitesse du fluide
$$
\bmu=\frac{m_e n_e \bmu_e + m_i n_i \bmu_i}{\rho}\approx \bmu_i,
$$
la densité de courant
\begin{align*}
\bmj&=-e n_e \bmu_e + Ze n_i\bmu_i\\
&=e n_e (\bmu_i - \bmu_e)
\end{align*}
la pression scalaire
$$
p=n_e k T_e + n_i k T_i,
$$
où $k$ est la constante de Boltzmann. Il faut y ajouter les équations de Maxwell puisqu'on est en présence d'un champ magnétique $\bB$ et d'un champ électrique $\bE$. Les équations de la MHD résistive pour un seul fluide (voir \cite{braginskii}) s'écrivent alors
\begin{equation}\label{int:mhd}
\left\{
\begin{array}{ll}
\disp{\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bmu)=s} &\hspace*{-2cm}\text{(Conservation des particules)}\\
\disp{\rho\Bigl(\frac{\partial \bmu}{\partial t}+\bmu\cdot\nabla\bmu \Bigr)+\nabla p=\bmj\times
\bB}\\
&\hspace*{-4.5cm}\text{(Conservation de la quantité de mouvement)}\\
\disp{ \frac{3}{2} \Bigl(\frac{\partial p}{\partial t}+\bmu\cdot\nabla
p \Bigr)+\frac{5}{2} p \nabla\cdot\bmu+\nabla Q = s'}\\
&\hspace*{-4.5cm}\text{(Conservation de l'énergie des particules)}\\[5pt]
\disp{\nabla\times\bE=-\frac{\partial \bB}{\partial t} }&\hspace*{-2cm}\text{(Loi de Faraday)}\\
\nabla\cdot\bB=0 &\hspace*{-2cm}\text{(Conservation de $\bB$)}\\
\bE+\bmu\times\bB=\eta \bmj &\hspace*{-2cm}\text{(Loi d'Ohm)}\\
\nabla\times \bH =\bmj &\hspace*{-2cm}\text{(Loi d'Ampère)}\\
\bB=\mu \bH &\hspace*{-2cm}\text{(Perméabilité magnétique)}\end{array}
\right.
\end{equation}
où $Q$ désigne le flux de chaleur, $\eta$ le tenseur de résistivité,
$s$ et $s'$ les termes sources. Nous donnerons dans la section 5 les équations du modèle à deux fluides comme dans \cite{braginskii}.

Pour simplifier le système (\ref{int:mhd}), nous devons définir quelques cons\-tantes de temps caractéristiques du plasma. Le temps d'\hbox{Alfvén}~$\tau_A$~est
$$
\tau_A=\frac{a(\mu_0mn)^{1/2}}{B_0}
$$
où $a$ est le petit rayon du plasma et $\bB_0$ le champ magnétique toroïdal. Son ordre de grandeur est la microseconde dans les Tokamaks actuels.

La constante de temps de diffusion de la densité des particules~$n$~est
$$
\tau_n=\frac{a^2}{D}
$$
où $D$ le coefficient de diffusion des particules. De même, la constante de temps de diffusion de chaleur des électrons et des ions sont
$$
\tau_e=\frac{n_ea^2}{K_e},\qquad
\tau_i=\frac{n_ia^2}{K_i}
$$
où $n_e$, $n_i$ sont la densité respectivement des électrons et des ions et~$K_e$, $K_i$ sont leur conductivité thermique. Ces constantes $\tau_n$, $\tau_e$, $\tau_i$ sont de l'ordre de la milliseconde dans les Tokamaks actuels.

Enfin, la constante de temps de la diffusion résistive de la densité du courant et du champ magnétique dans le plasma est donné par:
$$
\tau_r=\frac{\mu_0 a^2}{\eta}
$$
et est de l'ordre de la seconde.

Si une constante de temps globale pour la diffusion dans le plasma est définie par
$$
\tau_p=\inf(\tau_n,\tau_e,\tau_i,\tau_r)
$$
nous observons que
$$
\tau_A\ll \tau_p.
$$
À l'échelle de temps de la diffusion $\tau_p$, le terme $\rho(\sfrac{\partial
u}{\partial t}+\bmu\nabla \bmu)$ est petit comparé à $\nabla p $ (voir
\cite{maschke,nelson}) et l'équation d'équilibre
$$
\nabla p=\bmj\times\bB
$$
est satisfaite à chaque instant.

\section{Équilibre du plasma dans un Tokamak}
Les équations qui décrivent l'équilibre d'un plasma en présence d'un champ magnétique sont d'une part les équations de Maxwell et d'autre part les équations d'équilibre du plasma lui-même.
Les équations de Maxwell magnétostatiques suivantes sont satisfaites dans tout l'espace (plasma inclus):
\begin{equation}
\label{eqn:maxwell}
\left \lbrace
\begin{array}{lll}
\nabla \cdot B &=& 0\\
\nabla \times \Bigl(\dfrac{B}{\mu} \Bigr) &=& j
\end{array}
\right.
\end{equation}
où $B$ représente le champ magnétique, $\mu$ la perméabilité magnétique et $j$ la densité de courant.
La première relation de (\ref{eqn:maxwell}) est l'équation de conservation de l'induction magnétique et la seconde est le théorème d'Ampère. 

L'équation d'équilibre du plasma s'écrit:

\begin{equation}
\label{eqn:equilibrium}
\nabla p = j \times B
\end{equation}
Cette équation (\ref{eqn:equilibrium}) signifie que le plasma est en équilibre lorsque la force
$\nabla p$ est égal en tout point à la force de Lorentz de la pression magnétique $j \times B$. On déduit immédiatement de (\ref{eqn:equilibrium}) que
\begin{align}
\label{eqn:cons1}
B \cdot \nabla p &= 0\\
\label{eqn:cons2}
j \cdot \nabla p &= 0.
\end{align}
Ainsi, pour un plasma à l'équilibre, les lignes de champ et de courant reposent sur des surfaces isobares ($p=const.$); ces surfaces, engendrées par les lignes de champ, sont appelées surfaces magnétiques. Afin qu'elles restent à l'intérieur d'un volume borné de l'espace, il est nécessaire qu'elles aient une topologie toroïdale. Ces surfaces forment une famille de tores emboîtés (voir figure \ref{SM}).
Le tore le plus interne dégénère en une courbe qui est appelée axe magnétique.
\begin{figure}
\begin{center}
\includegraphics[width=.6\textwidth]{surface.jpg}
\caption{Surfaces magnétiques toroïdales pour un Tokamak}
\label{SM}
\end{center}
\end{figure}

En coordonnées cylindriques $(r,z,\phi)$ (où $r=0$ est l'axe vertical du tore), l'hypothèse de symétrie axiale consiste à supposer que le champ magnétique $B$ est indépendant de l'angle toroïdal $\phi$. Le champ magnétique peut être décomposé sous la forme $B=B_p +B_{\phi}$, où $B_p=(B_r,B_z)$ est la composante poloïdale et $B_{\phi}$ est la composante toroïdale.
À partir de l'équation (\ref{eqn:maxwell}) on peut définir le flux poloïdal $\psi(r,z)$ tel que
\begin{equation}
\left \lbrace
\begin{array}{lll}
B_r &=&-\dfrac{1}{r}\dfrac{\partial \psi}{\partial z}\\[10pt]
B_z &= & \dfrac{1}{r}\dfrac{\partial \psi}{\partial r}
\end{array}
\right.
\end{equation}
En ce qui concerne la composante toroïdale $B_{\phi}$ nous définissons $f$ par
\begin{equation}
B_{\phi}=\frac{f}{r}\,e_{\phi},
\end{equation}
où $e_{\phi}$ est le vecteur unité dans la direction toroïdale, et $f$ est la fonction diamagnétique.
Le champ magnétique s'écrit:
\begin{equation}
\label{eqn:B}
\left \lbrace
\begin{array}{lll}
B&=&B_p+B_{\phi}\\[8pt]
B_p &=&\dfrac{1}{r}[\nabla \psi \times e_{\phi}] \\[8pt]
B_{\phi} &= & \dfrac{f}{r}\, e_{\phi}
\end{array}
\right.
\end{equation}

Le champ poloïdal et le champ toroïdal sont représentés sur la figure \ref{tok} qui montre le Tokamak comme un transformateur dont le plasma est le secondaire. La figure \ref{iter} représente les Tokamaks Tore Supra (CEA Cadarache, France), JET (Joint European Torus, Culham, Angleterre), et ITER (en construction à Cadarache).
\begin{figure}[htb]
\begin{center}
\includegraphics[width=.65\textwidth]{tokamak.jpg}
\caption{Schéma de principe d'un Tokamak}
\label{tok}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=.95\textwidth]{Tokamaks-new.pdf}
\caption{Tokamaks Tore Supra (à gauche), JET (au milieu, avec l'aimable autorisation d'EUROfusion) et ITER (à~droite, crédit \copyright\ ITER Organization)}
\label{iter}
\end{center}
\end{figure}

D'après (\ref{eqn:B}), dans une configuration axisymétrique les surfaces magné\-tiques sont engendrées par la rotation des lignes de flux $\psi=const.$ autour de l'axe $r=0$ du tore.

De (\ref{eqn:B}) et de la seconde relation de (\ref{eqn:maxwell}) nous obtenons l'expression suivante pour $j$:
\begin{equation}
\label{eqn:j}
\left \lbrace
\begin{array}{lll}
j&=&j_p+j_{\phi}\\[10pt]
j_p &=&\dfrac{1}{r} \Bigl[\nabla \Bigl(\dfrac{f}{\mu} \Bigr) \times e_{\phi} \Bigr] \\[10pt]
j_{\phi} &= & (-\Delta^* \psi) e_{\phi},
\end{array}
\right.
\end{equation}
où $j_p$ et $j_{\phi}$ sont respectivement les composantes poloïdale et toroïdale de $j$, et l'opérateur $\Delta^*$ est défini par
\begin{equation}
\label{eqn:delta*}
\Delta^* .= \frac{\partial }{\partial r}\Bigl(\frac{1}{\mu r} \,\frac{ \partial\cdot}{\partial r} \Bigr)
+ \frac{\partial}{\partial z} \Bigl(\frac{1}{\mu r}\, \frac{ \partial\cdot}{\partial z} \Bigr)
\end{equation}

Les expressions (\ref{eqn:B}) et (\ref{eqn:j}) pour $B$ et $j$ sont valides dans tout l'espace car elles ne font intervenir que les équations de Maxwell et l'hypothèse d'axisymétrie.

Dans la région du plasma, la relation (\ref{eqn:cons1}) implique que $\nabla p$ et $\nabla \psi$ sont colinéaires, et donc $p$ est constant sur chaque surface magnétique. Ce qu'on note par
\begin{equation}
p=p(\psi)
\end{equation}
La relation (\ref{eqn:cons2}) combinée avec l'expression (\ref{eqn:j}) implique que $\nabla f$ et~$\nabla p$ sont colinéaires, et donc $f$ est aussi constante sur chaque surface magnétique
\begin{equation}
f=f(\psi)
\end{equation}
La relation d'équilibre (\ref{eqn:equilibrium}) combinée avec l'expression
(\ref{eqn:B}) et (\ref{eqn:j}) pour~$B$ et $j$ implique que:
\begin{equation}
\label{eqn:gradshaf0}
\nabla p = - \dfrac{\Delta^* \psi}{r}\, \nabla \psi - \frac{f}{\mu_0 r^2}\, \nabla f,
\end{equation}
ce qui conduit à l'équation d'équilibre dite de Grad-Shafranov:
\begin{equation}
\label{eqn:gradshaf}
-\Delta^* \psi = r p'(\psi) + \frac{1}{\mu_0 r}(ff')(\psi),
\end{equation}
où $\Delta^*$ est l'opérateur linéaire elliptique donné par (\ref{eqn:delta*})
où $\mu$ est égal à la perméabilité magnétique $\mu_0$ du vide, et $p'$ et $f'$ désignent respectivement les dérivées de $p$ et $f$ par rapport à $\psi$.

De (\ref{eqn:j}) il est clair que le terme de droite de (\ref{eqn:gradshaf})
représente la composante toroïdale de la densité de courant du plasma. Il fait intervenir $p(\psi)$ et $f(\psi)$, qui ne sont pas mesurées directement dans le plasma.

Dans le vide, le flux magnétique $\psi$ satisfait à
\begin{equation}
\label{eqn:gradshafvac}
-\Delta^* \psi = 0.
\end{equation}
L'équilibre d'un plasma dans un domaine $\Omega$ représentant la région du vide est un problème à frontière libre.
La frontière libre du plasma est définie soit par son contact avec un limiteur $D$
(ligne de flux la plus extérieure à l'intérieur du limiteur),
soit comme étant une séparatrice magnétique (courbe hyperbolique avec un point-$X$).
La~région $\Omega_p \subset\nobreak \Omega$ contenant le plasma est définie par
\begin{equation}
\Omega_p=\lbrace x \in \Omega \mid \psi(x) \ge \psi_b \rbrace,
\end{equation}
où soit $\psi_b=\max_{D} \psi$ dans la configuration du limiteur,
soit $\psi_b=\nobreak\psi(X)$ dans la configuration du point-$X$ (voir Figure \ref{fig:jet-tore}). La détermination numérique de la frontière libre du plasma est présentée dans \cite{Heumann-et-al}.

\begin{figure}[htb]
\centering
\includegraphics[width=.95\textwidth]{jet-tore.pdf}
\caption{\label{fig:jet-tore}Définition du bord du plasma (ligne bleue épaisse).
À gauche, l'exemple JET (Joint European Torus), configuration point-X. À droite, l'exemple TORE SUPRA (le Tokamak CEA-EURATOM à Cadarache), configuration avec limiteur (le limiteur est représenté par la ligne noire). Les lignes bleues fines représentent les iso-contours de $\psi$.}
\end{figure}

Supposons que les conditions de Dirichlet au bord, $h$, soient données sur
$\Gamma = \partial \Omega$ qui est la section poloïdale
de l'enceinte à vide,
les équations gouvernant le comportement
de $\psi(r,z)$ dans l'enceinte à vide sont finalement:
\begin{equation}
\label{eqn:final}
\left \lbrace
\begin{array}{rcl}
-\Delta^* \psi &= &\Bigl[r A(\bar{\psi}) +
\dfrac{1}{r} B(\bar{\psi}) \Bigr] \chi_{\Omega_p}\quad \text{dans } \Omega \\[10pt]
\psi&=& h\quad \text{sur } \partial\Omega
\end{array}
\right.
\end{equation}
avec
\[
A(\bar{\psi})= p'(\bar{\psi}),\quad B(\bar{\psi})=\dfrac{1}{\mu_0}(ff')(\bar{\psi}),\quad
\bar{\psi}= \dfrac{\psi - \max_\Omega \psi}{\psi_b - \max_\Omega \psi} \in [0,1]
\]
dans $\Omega_p$
(ce flux normalisé est introduit de sorte que $A$ et $B$
soient définis sur l'intervalle fixe $[0,1]$), $\chi_{\Omega_p}$
est la fonction caractéristique de $\Omega_p$.

Le but de la section suivante est de donner une méthode pour l'identification en temps réel du courant du plasma, c'est-à-dire les fonctions non linéaires $A$ et $B$ dans l'équation elliptique (\ref{eqn:final}).

\section{Le problème d'équilibre inverse}

\subsection{Mesures expérimentales}


Les données de mesures expérimentales sur JET (voir figure \ref{fig:jet-data-mesh}) sont:

\begin{itemize}

\item Les mesures magnétiques
\begin{itemize}
\item $\psi(M_i)=h_i\ \mathrm{ sur }\ \Gamma$,
donné par les boucles de flux. Grâce à une interpolation entre les points $M_i$ ces mesures produisent la condition de Dirichlet au bord~$h$.
\medskip

\item $\dfrac{1}{r} \dfrac{\partial \psi}{\partial n}(N_i)=g_i\ \mathrm{sur}\ \Gamma$,
qui correspond à la composante du champ magnétique poloïdal,
mesuré par les sondes magnétiques,
qui est tangente à l'enceinte à vide. En effet, d'après (\ref{eqn:B}) la composante tangentielle de $B_p$ est égale à la composante normale $\dfrac{1}{r}\dfrac{\partial \psi}{\partial n}$ de $\dfrac{1}{r}\nabla \psi$.
\end{itemize}

\medskip
\item Les mesures polarimétriques qui donnent la rotation de Faraday de l'angle de radiation infrarouge qui traverse la section du plasma le long de plusieurs cordes $C_i$:
$$
\int_{C_i} n_e(\bar{\psi}) B_{\|}\,dl =
\int_{C_i} \frac{n_e(\bar{\psi})}{r}\, \frac{\partial \psi}{\partial n}\,dl=\alpha_i,
$$
où $n_e$ représente la densité électronique qui est approximativement constante sur chaque ligne de champ,
$B_{\|}$ est la composante du champ poloïdal tangent à $C_i$
et $\sfrac{\partial}{\partial n}$ représente la dérivée normale
de $\psi$ par rapport à $C_i$.

\item Les mesures interférométriques qui donnent les intégrales de densité sur les cordes $C_i$
$$
\int_{C_i}n_e(\bar{\psi})\,dl=\beta_i.
$$

\item Les mesures de pression cinétique obtenues à partir des mesures de pression et de température, par exemple dans le plan équatorial:
$$
p(r,0)=p_d(r).
$$

\item Et les mesures d'angle MSE (Motional Stark Effect) prises en divers points $x_i=(r_i,z_i)$:
$$
m(B_r,B_z,B_\phi)_i=\gamma_i
$$
avec
$$
\tan (m(B_r,B_z,B_\phi)) = \dfrac{a_{1} B_r + a_{2} B_z + a_{3} B_\phi}
{a_{4} B_r + a_{5} B_z + a_{6} B_\phi}.
$$
\end{itemize}

\begin{figure}[htb]
\centering
\includegraphics[height=.9\textwidth,angle=-90]{jet-data-mesh.pdf}
\caption{\label{fig:jet-data-mesh} À gauche: les droites vertes représentent les cordes utilisées pour les mesures de polarimétrie et d'interférométrie. À droite: partie d'une enceinte à vide. Au milieu en bas, un exemple de maillage d'éléments finis utilisé pour les simulations numériques.}
\end{figure}

\subsection{Énoncé du problème inverse}

Le problème d'identification numérique est formulé comme un problème de minimisation des moindres carrés avec une régularisation de Tikhonov.
La fonction coût est définie par:
\begin{equation}
\label{eqn:costfunction}
J(A,B,n_e)=J_0 + K_1 J_1 + K_2 J_2 + K_3 J_3 + K_4 J_4 + J_\epsilon
\end{equation}
avec
$$
\begin{array}{l}
J_0= \displaystyle \sum_i
\Bigl(\dfrac{1}{r} \frac{\partial \psi}{\partial n}(N_i)
- g_i \Bigr)^2 \\[10pt]
J_1= \displaystyle \sum_i \Bigl(\displaystyle \int_{C_i}
\frac{n_e}{r} \frac{\partial \psi}{\partial n}dl - \alpha_i \Bigr)^2 \\[10pt]
J_2= \displaystyle \sum_i \Bigl(\displaystyle \int_{C_i} n_edl - \beta_i \Bigr)^2 \\[10pt]
J_3= \displaystyle \int_{R_{min}}^{R_{max}} \bigl(p(r,0)-p_d(r) \bigl)^2 dr\\[10pt]
J_4= \displaystyle \sum_i \bigl(m(B_r,B_z,B_\phi)_i - \gamma_i \bigr)^2
\end{array}
$$
et $K_1$, $K_2$, $K_3$ et $K_4$ sont des paramètres de poids permettant de donner plus ou moins d'importance aux mesures expérimentales correspondantes \cite{Blum90}.

Le problème inverse de la détermination de $A$ et $B$ est mal posé au sens d'Hadamard. Par conséquent une procédure de régularisation peut être utilisée pour le transformer en un problème bien posé
\cite{Tikhonov77}.
Le terme de régularisation de Tikhonov $J_\epsilon$ contraint les fonctions~$A$,~$B$ et~$n_e$ à être suffisamment lisses et s'écrit:
$$
J_\epsilon =
\epsilon_1 \displaystyle \int_0^1 [A''(x)]^2 dx
+
\epsilon_2 \displaystyle \int_0^1 [B''(x)]^2 dx
+
\epsilon_3 \displaystyle \int_0^1 [n_e''(x)]^2 dx
$$
où $\epsilon_1$, $\epsilon_2$ et $\epsilon_3$ sont les paramètres de régularisation.

Il faut remarquer que la densité électronique $n_e$ n'intervient pas dans (\ref{eqn:final}). Cependant, dès que l'on veut utiliser les mesures polarimétriques il est nécessaire d'inclure $n_e$ (et donc l'interférométrie) dans la procédure d'identification.
Le problème inverse peut finalement être formulé comme suit:
\begin{equation}
\left \lbrace
\begin{array}{l}
\text{Trouver $A^*$, $B^*$, $n_e^*$ tels que}\\[5pt]
J(A^*,B^*,n_{e}^*)=\inf J(A,B,n_{e})
\end{array}
\right.
\end{equation}

\subsection{Identification numérique}

%\subsubsection{Direct problem}
Le problème (\ref{eqn:final}) se résout en utilisant une méthode d'éléments finis
\cite{Ciarlet80}.
Soient $H^1(\Omega)$ et $V=H^1_0(\Omega)$ les espaces de Sobolev usuels.
L'approximation par éléments finis est basée sur la formulation faible suivante:
\begin{equation}
\label{eqn:var}
\left \lbrace
\begin{array}{l}
\text{Trouver $\psi \in H^1(\Omega)$ tel que} \\[5pt]
\displaystyle \int_\Omega \dfrac{1}{\mu_0 r}\nabla \psi \cdot \nabla v dx =
\int_{\Omega_p}[r A(\bar{\psi})+ \frac{1}{r}B(\bar{\psi})] v dx\quad \forall v\ \in V \\[10pt]
\end{array}
\right .
\end{equation}
De manière classique, $\Omega$ est approché en utilisant des triangles par un domaine polygonal $\Omega_h$,
l'espace $V$ est approché par un espace~$V_h$ de dimension finie $n$. On utilise une méthode d'éléments finis $P1$ dans laquelle les fonctions dans $V_h$ sont affines sur chaque triangle et continues sur tout le domaine.

Soit $K$ la matrice de rigidité de la méthode d'éléments finis.
Notons aussi (abusivement) $\psi \in \RR^n$ les composantes de la fonction flux magnétique approchée dans $V_h$.

Les fonctions inconnues $A$, $B$ et $n_e$ sont approchées par une décomposition
dans une base réduite $(\phi_i)_{i=1,\dots m}$
\begin{align*}
A(x)&= \sum_i a_i \phi_i(x)\\[10pt]
B(x)&= \sum_i b_i \phi_i(x)\\[10pt]
n_e(x)&= \sum_i c_i \phi_i(x)
\end{align*}
Cette base peut être constituée de divers types de fonctions
(polynômes, splines, ondelettes, etc.) \cite{Blum97}.
Soit $u$ le vecteur de $\RR^{3m}$ défini par
$u=(a_1,\dots ,a_m,b_1,\dots,b_m,c_1,\dots,c_m)$.
Avec ces notations, la discrétisation du problème (\ref{eqn:var})
s'écrit comme suit:
\begin{equation}
\label{eqn:fp}
\left \lbrace
\begin{array}{l}
\text{Étant donné $u \in \RR^{3m}$, résoudre l'équation}\\[5pt]
\tilde{K} \psi = D(\psi) u + h\quad\text{par la méthode du point fixe,}
\end{array}
\right .
\end{equation}
où $D(\psi)$ désigne la \og matrice du courant de plasma\fg $n\times 3m$, et $\tilde{K}$ est la matrice de rigidité modifiée pour imposer la condition de Dirichlet au bord représentée par $h$.

%\subsubsection{Inverse problem}

Le problème d'optimisation inverse discret est:
\begin{equation}
\label{eqn:optim}
\left \lbrace
\begin{array}{l}
\text{Trouver $u$ minimisant}:\\[5pt]
J(u)=\| C(\psi)\psi - k \|^2 + u^{T} \Lambda u\\[5pt]
\text{avec $\psi$ satisfaisant à (\ref{eqn:fp}),}
\end{array}
\right.
\end{equation}
où $C(\psi)$ est l'opérateur d'observation.
La quantité $C(\psi)\psi$ représente les sorties du modèle correspondant aux mesures expérimentales, données plus haut, notées $k$.
La matrice $\Lambda$ représente les termes de régularisation.
Le premier terme de $J$ dans (\ref{eqn:optim}) correspond à
$J_0 + K_1 J_1 + K_2 J_2 + K_3 J_3 + K_4 J_4$ et le second à $J_\epsilon$.

Pour résoudre ce problème nous utilisons un algorithme basé sur une méthode de point fixe pour (\ref{eqn:fp}) et l'équation normale pour (\ref{eqn:optim}).

\subsubsection{L'algorithme}
À la $n$\ieme itération, $\psi_n$ et $u_n$ sont donnés.
L'application non linéaire $u \mapsto \psi(u)$ est approchée par la relation affine
$$
\psi =\tilde{K}^{-1}[D(\psi_n) u +h]
$$
et la fonction coût est à minimiser par
$$
\begin{array}{lll}
J(u)&=&\| C(\psi_n)\psi - k \|^2 + u^{T} \Lambda u \\[10pt]
&=&\| C(\psi_n) \tilde{K}^{-1} D(\psi_n) u +(C(\psi_n) \tilde{K}^{-1} h - k) \|^2 + u^{T} \Lambda u \\[10pt]
&=&\| E_n u+ F_n \|^2 + u^{T} \Lambda u
\end{array}
$$
avec des notations évidentes.
L'équation normale
$$
(E^{T}_n E_n + \Lambda)u= - E^{T}_n F_n
$$
est résolue pour remplacer $u_n$ par $u_{n+1}$.
Alors une méthode du point fixe pour (\ref{eqn:fp})
permet de remplacer $\psi_n$
par $\psi_{n+1}$
$$
\psi_{n+1}= \tilde{K}^{-1}[D(\psi_n) u_{n+1} +h].
$$
Puisque l'algorithme est habituellement initialisé par l'équilibre à l'instant précédent, deux ou trois itérations sont en général suffisantes pour assurer la convergence.
Un logiciel C++ appelé EQUINOX \cite{Bosak01,BBF,JCP-Equinox}, basé sur l'algorithme présenté ci-dessus, a été développé en collaboration avec le CEA-IRFM à Cadarache, et a été implémenté pour JET (Joint European Torus)
et pour TORE SUPRA (le Tokamak CEA-EURATOM à Cadarache).
La figure \ref{fig:equinox} montre une sortie graphique d'Equinox pour JET.
Avec toutes ces techniques, il est possible de suivre l'évolution quasi statique de l'équilibre du plasma, dans les configurations TORE SUPRA ou JET,
avec des frontières libres définies par contact avec le limiteur ou par un point-$X$.
Il est aussi possible de simuler des configurations ITER.
\begin{figure}[htb]
\centering
\includegraphics[width=.95\textwidth]{equinox.pdf}
\caption{\label{fig:equinox} Une sortie d'EQUINOX. Le plasma est dans une configuration point-X.
Sur la colonne de gauche, les fonctions identifiées $p'$, $ff'$ et $n_e$ ainsi que la densité du courant toroïdal $j$ et le facteur de sécurité $q$ sont affichés en terme de $\psi$ et de $r$ (dans le plan équatorial).}
\end{figure}

\section{Le système de transport}

Plusieurs groupes ont développé des modèles de calcul pour les processus de transport dans les Tokamaks. Ils utilisent généralement une géométrie uni-dimensionnelle, supposant que les surfaces magnétiques sont des cylindres avec des sections circulaires concentriques avec l'axe magnétique. De plus, ces modèles ignorent la condition d'équilibre ($\nabla
p=\bmj\times \bB$) et utilisent une loi d'Ohm simplifiée ($\bE=\eta\bmj$). Le modèle qui est présenté ici (voir aussi le chapitre 6 de \cite{book_blum}) est compatible avec l'équation d'équilibre pour des configurations axisymétriques et prend en compte des sections non circulaires des surfaces magnétiques.

Les équations de transport restent uni-dimensionnelles en utilisant un paramètre approprié pour les lignes de flux comme variable \og \hbox{radiale}\fg comme cela a été montré dans \cite{hinton} pour le transport néoclassique.

\subsection{La technique de moyennisation}
Soit $\chi$ une coordonnée arbitraire qui paramètre la surface magnétique $S$; nous définissons la moyenne sur $S$ d'une quantité arbitraire~$A$ par
\begin{equation}\label{17}
\langle A\rangle = \frac{\partial}{\partial V}\int_V
AdV=\frac{1}{V'}\int_S \frac{AdS}{|\nabla \chi|},
\end{equation}
avec
$$
V'=\frac{\partial V}{\partial \chi}=\int_S\frac{dS}{|\nabla \chi|}.
$$
où $V$ est le volume délimité par la surface magnétique $S$. 

Une telle moyenne a les propriétés suivantes:
\begin{equation}\label{18}
\langle \div \bW\rangle=\frac{\partial}{\partial V}\langle\bW\cdot\nabla
V\rangle, \quad\forall \bW,
\end{equation}
\begin{equation}\label{19}
\frac{\partial}{\partial t}(V'\langle A
\rangle)=V'\langle\dot{A}\rangle+\frac{\partial}{\partial \chi}\langle A\bmu_\chi\cdot\nabla
V\rangle,\quad\forall A
\end{equation}
où $\dot{A}$ désigne la dérivée temporelle au point fixe $(r,z)$,
alors que ${\sfrac{\partial}{\partial t}}$ désigne la dérivée temporelle à $\chi$ fixé. Le vecteur $\bmu_\chi$ est la vitesse de la surface à $\chi$ constant, définie par
\begin{equation}\label{20}
\dot{\chi}+\bmu_\chi\cdot\nabla \chi =0.
\end{equation}
De (\ref{19}) on peut déduire que
\begin{equation}\label{21}
\frac{\partial V'}{\partial t}=\frac{\partial}{\partial \chi}\langle
u_\chi\cdot\nabla V\rangle.
\end{equation}
Nous allons maintenant moyenner les équations de transport de \hbox{Braginskii} (voir \cite{braginskii}) afin d'obtenir un ensemble d'équations de diffusion 1-D par rapport à la variable d'espace $\chi$. Comme la vitesse de diffusion parallèle à la surface magnétique est plus grande que la vitesse de diffusion perpendiculaire, nous supposerons que les densités et températures sont constantes sur chaque surface magnétique.

\subsection{L'équation de densité}
L'équation de conservation de la densité électronique $n_e$ peut s'écrire en terme de la vitesse du flot électronique $\bmu_e$,
\begin{equation}\label{22}
\dot{n_e}+\nabla\cdot(n_e\bmu_e)=S_1,
\end{equation}
où $S_1$ est un terme source.

En moyennant (\ref{22}) sur une surface magnétique, en multipliant par~$V'$ et en utilisant les relations (\ref{18}) et (\ref{19}), nous obtenons
\begin{equation}\label{23}
\frac{\partial}{\partial t}(n_e V')+\frac{\partial}{\partial
\chi}\langle n_e(\bmu_e-\bmu_\chi)\cdot\nabla V\rangle=\langle S_1\rangle V'.
\end{equation}
Si le flux de particules relativement à une surface à $\chi$ constant est défini par
\begin{equation}\label{24}
\Gamma_e=\langle n_e(\bmu_e-\bmu_\chi)\cdot\nabla \chi\rangle,
\end{equation}
l'équation (\ref{23}) peut s'écrire
\begin{equation}\label{25}
\frac{\partial}{\partial t}(n_e V')+\frac{\partial}{\partial
\chi}(\Gamma_e V')=\<S_1\>V'.
\end{equation}

\subsection{Les équations d'énergie}
Si on néglige les termes de viscosité, les équations d'équilibre d'énergie pour les électrons sont
\begin{equation}\label{26}
\frac{3}{2}\,\dot{p_e}+\nabla\cdot\Bigl(\bQ_e+\frac{5}{2}\,p_e\bmu_e \Bigr)=\bmj\cdot\bE-
Q_\Delta - \bmu_i\cdot\nabla p_i + S_2,
\end{equation}
où $p_e$ et $p_i$ sont, respectivement, la pression des électrons et des ions,
$\bQ_e$ le flux de chaleur électronique, $Q_\Delta$ l'échange de chaleur électron-ion, $\bmu_i$ la vitesse du flux ionique et $S_2$ un terme source.

Suivant \cite{braginskii}, nous avons
$$
Q_\Delta=\frac{3m_e}{m_i}\,\frac{n_e}{t_e}\,(T_e-T_i),
$$
$$
p_e=n_e k T_e,\quad  p_i=n_ikT_i, \quad  n_e=Zn_i,
$$
où $t_e$ est le temps de collision électronique et $k$ la constante de Boltzmann.

En moyennant (\ref{26}), en multipliant par $V'$ et en utilisant (\ref{18}),
(\ref{19}) et (\ref{21}), nous obtenons
\begin{multline}
\frac{3}{2(V')^{2/3}}\,\frac{\partial}{\partial
t}[p_e(V')^{5/3}]+\frac{\partial }{\partial
\chi}\Bigl[(q_e+\frac{5}{2}\,kT_e\Gamma_e)V'\Bigr]\\
= \Bigl[\<\bmj\cdot\bE\>-\<u_\chi\cdot\nabla
p_e\>-\<\bmu_i\cdot\nabla p_i\>-Q_\Delta +\<S_2\> \Bigr]V',\label{28}
\end{multline}
avec $q_e=\<\bQ_e\cdot\nabla \chi\>$.

De manière analogue, l'équation d'équilibre d'énergie ionique est
\begin{equation}\label{29}
\frac{3}{2}\,\dot{p_i}+\nabla\cdot \Bigl(\bQ_i+\frac{5}{2}\,p_i\bmu_i \Bigr)=Q_\Delta + \bmu_i\cdot\nabla p_i + S_3,
\end{equation}
où $\bQ_i$ est le flux de chaleur ionique et $S_3$ un terme source.
En utilisant la même technique que pour l'équation d'énergie électronique, on obtient
\begin{multline}\label{30}
\frac{3}{2(V')^{2/3}}\, \frac{\partial}{\partial t}(p_i
(V')^{5/3})\\+\frac{\partial}{\partial
\chi} \Bigl[\Bigl(q_i+\frac{5}{2}\,\frac{kT_i\Gamma_e}{Z} \Bigr)V'\Bigr]-\frac{\Gamma_e}{n_e}\,\frac{\partial
p_i}{\partial \chi}V'
=Q_\Delta V'+\<S_3\>V',
\end{multline}
avec $q_i=\<\bQ_i\cdot\nabla \chi\>$.

\subsection{Les équations de flux}
La vitesse $\bmu_\Psi$ d'une surface à $\Psi$ constant est définie par
\begin{equation}\label{31}
\dot{\Psi}+\bmu_\Psi\cdot\nabla \Psi =0.
\end{equation}
Nous avons
\begin{equation}\label{32}
\dot{\Psi}=\frac{\partial \Psi}{\partial t}+\frac{\partial
\Psi}{\partial \chi}\,\dot{\chi}.
\end{equation}
En utilisant les définitions (\ref{20}) et (\ref{31}) pour $\bmu_\chi$ et $\bmu_\Psi$, l'équation (\ref{32}) devient l'équation d'évolution pour le flux
poloïdal $\Psi$,
\begin{equation}\label{33}
\frac{\partial \Psi}{\partial t}=(\bmu_\chi-\bmu_\Psi)\cdot\nabla\Psi.
\end{equation}
Définissons maintenant le flux toroïdal $\Phi$ par
\begin{equation}\label{34}
\Phi=\frac{1}{2\pi}\int_V\frac{B_t}{r}dV
\end{equation}
et le facteur de sécurité $q$ par
\begin{equation}\label{35}
q=-\frac{1}{2\pi}\,\frac{\partial \Phi}{\partial \Psi}.
\end{equation}
Des définitions (\ref{33}) et (\ref{35}) on peut déduire
\begin{equation}\label{36}
q=-\frac{1}{4\pi^2}\,\frac{\partial V}{\partial \Psi}\,\<r^{-2}\>f.
\end{equation}
De plus, on a
\begin{equation}\label{37}
\frac{\partial \Phi}{\partial t}=\frac{\partial \Phi}{\partial
t}\Big|_\Psi-2\pi q\, \frac{\partial \Psi}{\partial t}.
\end{equation}
Suivant \cite{hinton} et \cite{pereversev}, la dérivée temporelle de
$\Phi$ à $\Psi$ fixé est égale à
\begin{equation}\label{38}
\frac{\partial \Phi}{\partial t}\Big|_\Psi=\frac{1}{2\pi}\,\frac{\partial
V}{\partial \Psi}\,\<\bE\cdot\bB\>.
\end{equation}
On peut alors déduire de (\ref{33}), (\ref{36})-(\ref{38}) l'équation d'évolution pour le flux toroïdal $\Phi$:
\begin{equation}\label{39}
\frac{\partial \Phi}{\partial t}=\frac{1}{2\pi}\,\frac{\partial
V}{\partial \Psi} \Bigl[f\<r^{-2}\>(\bmu_\chi-\bmu_\Psi)\cdot\nabla \Psi+\<\bE\cdot\bB\> \Bigr].
\end{equation}

\Subsection{Le choix du paramètre \texorpdfstring{$\chi$}{chi} pour la surface de flux}
Définissons $N_e$, $\sigma_e$ et $\sigma_i$ par
\begin{equation}\label{40}
N_e=n_eV',\quad \sigma_e=p_e(V')^{5/3},\quad \sigma_i=p_i(V')^{5/3}.
\end{equation}

La variable $N_e$ est le taux d'électrons entre deux surfaces magnétiques adjacentes, tandis que $\sigma_e$ et $\sigma_i$ sont respectivement
les entropies électronique et ionique. Il suit clairement de (\ref{25}),
(\ref{28}) et (\ref{30}) que~$N_e$,~$\sigma_e$ et~$\sigma_i$ restent constants pendant les évolutions adiabatiques; on les appelle des variables adiabatiques.

Les équations (\ref{25}), (\ref{28}), (\ref{30}), (\ref{33}) et (\ref{39}) constituent un ensemble d'équations $1D$ par rapport à $\chi$,
pour les cinq variables $N_e$, $\sigma_e$, $\sigma_i$,~$\Psi$ et $\Phi$. Si on prend une de ces quantités adiabatiques comme variable~$\chi$, nous verrons que le système se réduit à un ensemble de quatre équations de diffusion.

L'idée la plus simple est de prendre $\Psi$ comme variable indépendante. C'est ce qu'ont fait Byrne et Klein dans la première version de leur code G2M
(voir \cite{byrne}) et Holmes et \textit{al.} dans leur étude des FCT (Flux Conserving Tokamaks)
(voir \cite{holmes}). Mais, dans les
tokamaks conventionnels, le flux poloïdal varie avec le courant du plasma et les tensions d'induction externes. Il semble donc mieux de prendre pour
$\chi$ le flux poloïdal normalisé $\overline{\Psi}$, qui reste dans l'intervalle $[0,1]$. C'est ce qu'ont fait Shumaker et \textit{al.}
(voir \cite{shumaker}) et aussi Turnbull et Storer (voir
\cite{turnbull}).

Mais dans les tokamaks actuels, le flux toroïdal $\Phi$ est moins mobile que le flux poloïdal $\Psi$. Le champ toroïdal externe appliqué est si grand que, grossièrement parlant, le flux toroïdal est presque invariant
pendant la décharge d'un tokamak. C'est pourquoi nous pouvons choisir pour $\chi$ le rayon moyen de section de la surface magnétique, défini à partir du flux toroïdal par
\begin{equation}\label{42}
\chi=\Bigl(\frac{\Phi}{\pi B_0} \Bigr)^{1/2},
\end{equation}
où $B_0$ est une constante, par exemple le champ toroïdal externe appliqué. C'est le choix fait par Hinton et Hazeltine (voir
\cite{hinton}) et utilisé par Miller (voir \cite{miller}) pour le contrôle de forme des Doublets et par Blum et \textit{al.} (voir
\cite{blum1981}).

Il faut aussi indiquer que Helton et {\textit{al}.} (voir~\cite{helton}) et Hogan (voir~\cite{hogan78}) utilisent le paramètre de flux de surface
\begin{equation}\label{43}
\chi= \Bigl(\frac{V}{2\pi^2R_0} \Bigr)^{1/2},
\end{equation}
où $R_0$ est le grand rayon du Tokamak. Le volume $V$ n'est pas une variable adiabatique, et c'est pourquoi chaque phase de l'évolution doit être scindée en deux étapes: une phase d'évolution où la géométrie de la surface magnétique est fixée et où $\chi$, donné par l'équation (\ref{43}), est pris comme paramètre de flux de surface, et une
phase adiabatique où un nouvel équilibre est calculé, en maintenant les quantités adiabatiques $N_e$, $\sigma_e$ et $\sigma_i$ constantes.

Finalement, indiquons qu'Hirshman et Jardin (voir
\cite{hirshman,jardin}) laissent $\chi$ représenter un paramètre arbitraire de flux de surface, et la vitesse des surfaces à $\Phi$ constant est calculée en dérivant par rapport au temps l'équation de Grad-Shafranov.

Le paramètre de flux de surface $\chi$ étant défini par (\ref{42}), les
équations (\ref{25}),
(\ref{28}), (\ref{30}), (\ref{33}) et (\ref{39}) peuvent se récrire comme un système de quatre équations de diffusion 1-D par rapport à $\chi$. En~effet, on déduit de (\ref{39}), puisque $\sfrac{\partial
\Phi}{\partial t}=0$, que
\begin{equation}\label{44}
(\bmu_\chi-\bmu_\Psi)\cdot\nabla\Psi=-\frac{\<\bE\cdot\bB\>}{f\<r^{-2}\>}.
\end{equation}
Par conséquent, (\ref{33}) donne
\begin{equation}\label{45}
\frac{\partial\Psi}{\partial t}=-\frac{\<\bE\cdot\bB\>}{f\<r^{-2}\>}.
\end{equation}
À partir de la loi d'Ohm
\begin{equation}\label{46}
\bE+\bmu\times\bB=\eta \bmj
\end{equation}
on peut déduire que
\begin{equation}\label{47}
\<\bE\cdot\bB\>=\eta\<\bmj\cdot\bB\>,
\end{equation}
où $\eta$ est la résistivité de Spitzer. De plus, on peut déduire
\begin{equation}\label{48}
\<\bmj\cdot\bB\>=-\frac{f^2}{\mu_0V'}\,\frac{\partial}{\partial\chi} \Bigl(\frac{C_2}{f}\,\frac{\partial\Psi}{\partial\chi} \Bigr),
\end{equation}
avec
$$
C_2=V'\<r^{-2}\nabla^2\chi\>.
$$
Ainsi (\ref{45}) devient
\begin{equation}\label{49}
\frac{\partial\Psi}{\partial t}-\frac{\eta f }{\mu_0
C_3}\,\frac{\partial}{\partial
\chi} \Bigl(\frac{C_2}{f}\,\frac{\partial\Psi}{\partial \chi} \Bigr)=0,
\end{equation}
avec $C_3=V'\<r^{-2}\>$. En différenciant (\ref{49}), on obtient l'équation de diffusion suivante pour $\Psi'=\sfrac{\partial\Psi}{\partial \chi}$:
\begin{equation}\label{50}
\frac{\partial\Psi'}{\partial t} -\frac{1}{\mu_0}\,\frac{\partial}{\partial \chi} \Bigl[\frac{\eta f}{C_3}\,\frac{\partial}{\partial \chi} \Bigl(\frac{C_2 \Psi'}{f} \Bigr) \Bigr]=0.
\end{equation}
On déduit des définitions (\ref{34}) et (\ref{42}) de $\Phi$ et $\chi$,
la relation suivante pour $f$:
\begin{equation}\label{51}
f=\frac{4\pi^2B_0\chi}{C_3}.
\end{equation}
L'équation (\ref{50}) peut s'écrire
\begin{equation}\label{52}
\frac{\partial\Psi'}{\partial t}-\frac{1}{\mu_0}\,\frac{\partial}{\partial
\chi} \Bigl[\frac{\eta\chi}{C_3^2}\,\frac{\partial}{\partial \chi} \Bigl(\frac{C_2C_3}{\chi}\,\Psi'\Bigr) \Bigr]=0.
\end{equation}
La quantité adiabatique $\Psi'$ est reliée au facteur de sécurité $q$ par
\begin{equation}\label{53}
\Psi'=-\frac{B_0 \chi}{q}.
\end{equation}
L'équation (\ref{52}) est une équation parabolique 1-D qui caractérise
la diffusion résistive du flux poloïdal par rapport au flux toroïdal.

Considérons maintenant l'équation d'énergie électronique (\ref{28}). Comme il a été montré dans \cite{hinton}, on a
\begin{equation}\label{54}
\<\bmj\cdot\bE\>=\frac{\<j_\phi/r\>\<\bE\cdot\bB\>}{f\<r^{-2}\>}+\<\bmu_\chi\cdot\nabla p\>,
\end{equation}
avec
\[
{\<j_\phi/r\>=-\frac{1}{\mu_0V'}\,\frac{\partial(C_2\Psi')}{\partial \chi}}.
\]
Mais, à l'aide des relations (\ref{47}), (\ref{48}), (\ref{51}) et (\ref{54}), le terme de droite de (\ref{28}) se simplifie, de sorte que l'équation d'énergie électronique n'est autre que\vspace*{-3pt}
\begin{multline}
\frac{3}{2(V')^{2/3}}\,\frac{\partial\sigma_e}{\partial
t}+\frac{\partial}{\partial
\chi} \Bigl[\Bigl(q_e+\frac{5}{2}\,kT_e\Gamma_e \Bigr)V'\Bigr]\\[5pt]
=-\frac{\Gamma_eV'}{n_e}\,\frac{\partial p_i}{\partial \chi}+\frac{\eta\chi}{\mu_0^2C_3^2}\,\frac{\partial}{\partial
\chi}(C_2\Psi')\,\frac{\partial}{\partial\chi} \Bigl(\frac{C_2C_3}{\chi}\,\Psi'\Bigr)-Q_\Delta V'+\<S_2\>V'\cdot\label{55}
\end{multline}
Les équations de densité électronique et d'énergie ionique (\ref{25}) et (\ref{30}) sont inchangées, de sorte que le système d'équations
(\ref{25}), (\ref{55}), (\ref{30}) et (\ref{52}) constitue un ensemble d'équations 1-D pour $N_e$, $\sigma_e$, $\sigma_i$ et $\Psi'$ en terme de $\chi$. Pour clore ce système d'équations il est nécessaire d'exprimer les variables adiabatiques $N_e$, $\sigma_e$, $\sigma_i$ et
$\Psi'$. C'est ce qui est fait dans \cite{hinton}, où les flux classique,
néoclassique et anormal sont exprimés comme des combinaisons quasi-linéaires de $\partial n_e/\partial \chi$, $\partial p_e/\partial
\chi$, $\partial p_i/\partial \chi$, et $\<\bE\cdot\bB\>$. Hirshman et Jardin (voir \cite{hirshman}) donnent une description complète de ces flux dans les régimes Pfirsch-Schluter et banane.

Dans la section suivante, nous allons utiliser des expressions simpli\-fiées pour $\Gamma_e$, $q_e$ et $q_i$ qui correspondent au modèle appelé
\og \hbox{diagonal}\fg.

\subsection{Le modèle diagonal}
Dans le modèle \og diagonal\fg nous supposons que
\begin{equation}\label{56}
n_e(\bmu_e-\bmu_\chi)=-D\nabla n_e,\quad \bQ_e=-K_e\nabla
T_e,\quad \bQ_i=-K_i\nabla T_i,
\end{equation}
où $D$ est le coefficient de diffusion électronique et $K_e$ et $K_i$
sont les conductivités thermiques électronique et ionique respectivement.

Nous avons donc
\begin{equation}\label{57}
\begin{split}
\Gamma_e&=-D\<\nabla^2\chi\>\,\frac{\partial n_e}{\partial
\chi},\\
q_e&=-K_e\<\nabla^2\chi\>\,\frac{\partial T_e}{\partial
\chi},\\
q_i&=-K_i\<\nabla^2\chi\>\,\frac{\partial T_i}{\partial\chi},
\end{split}
\end{equation}
et les équations (\ref{25}), (\ref{55}) et (\ref{30}) deviennent
\begin{equation}\label{58}
\frac{\partial N_e}{\partial t}-\frac{\partial}{\partial
\chi} \Bigl(DC_1\,\frac{\partial n_e}{\partial \chi} \Bigr)=\<S_1\>V',
\end{equation}
\begin{multline}
\frac{3}{2(V')^{2/3}}\,\frac{\partial\sigma_e}{\partial t}-\frac{\partial}{\partial\chi} \Bigl(K_eC_1\,\frac{\partial T_e}{\partial\chi} \Bigr)-\frac{5k}{2}\,\frac{\partial}{\partial\chi} \Bigl(DC_1 T_e\,\frac{\partial n_e}{\partial\chi} \Bigr)\\
=\frac{\eta\chi}{\mu_0^2C_3^2}\,\frac{\partial}{\partial\chi}(C_2\Psi')\frac{\partial}{\partial\chi} \Bigl(\frac{C_2C_3}{\chi}\Psi'\Bigr)+\frac{DC_1}{n_e}\,\frac{\partial p_i}{\partial\chi}\,\frac{\partial n_e}{\partial \chi} -Q_\Delta V'+ \<S_2\>V',\label{59}
\end{multline}
\begin{multline}\label{60}
\frac{3}{2(V')^{2/3}}\,\frac{\partial\sigma_i}{\partial
t}-\frac{\partial}{\partial \chi} \Bigl(K_iC_1\,\frac{\partial T_i}{\partial
\chi} \Bigr)-\frac{5k}{2}\,\frac{\partial}{\partial\chi} \Bigl(\frac{DC_1T_i}{Z}\,\frac{\partial n_e}{\partial \chi} \Bigr)\\
+\frac{DC_1}{n_e}\,\frac{\partial n_e}{\partial\chi}\,\frac{\partial p_i}{\partial \chi}
=Q_\Delta V'+\<S_3\> V'
\end{multline}
avec $C_1=V'\<\nabla^2\chi\>$.

Le modèle de diffusion diagonal est alors constitué par les équations (\ref{58}),
(\ref{59}), (\ref{60}) et (\ref{52}).

\subsection{L'équation d'équilibre moyennée}
L'équation d'équilibre de Grad-Shafranov peut s'écrire
\begin{equation}\label{61}
-\nabla\cdot(\frac{1}{\mu_0r^2}\nabla\Psi)=\frac{\partial p}{\partial
\Psi}+\frac{1}{2\mu_0r^2}\frac{\partial f^2}{\partial \Psi}
\end{equation}
La moyenne de l'équation (\ref{61}) sur une surface magnétique et l'utilisation de la relation (\ref{18}) conduisent à
\begin{equation}\label{62}
-\Psi'\frac{\partial}{\partial \chi}(C_2\Psi')=\mu_0 V'\frac{\partial
p }{\partial \chi}+\frac{1}{2}C_3\frac{\partial f^2}{\partial \chi}.
\end{equation}
Puisque ${p=\sfrac{(\sigma_e+\sigma_i)}{(V')^{5/3}}}$, les équations
(\ref{51}) et (\ref{62}) nous permettent de déterminer les fonctions
$f(\chi)$ et $V'(\chi)$ si les profils $\Psi'(\chi)$,
$\sigma_e(\chi)$, $\sigma_i(\chi)$ et les coefficients géométriques
$\<r^{-2}\nabla^2\chi\>$ et $\<r^{-2}\>$ ont été calculés auparavant.

L'intégration de (\ref{62}) sur le plasma donne
\begin{equation}\label{63}
I_p=-\frac{C_2\Psi'(\chi_{\max})}{2\pi\mu_0},
\end{equation}
où $I_p$ est le courant de plasma total et $\chi_{\max}$ la valeur de
$\chi$ au bord du plasma.

\subsection{Couplage des systèmes d'équilibre et de diffusion. Conditions au bord}

L'équation de Grad-Shafranov requiert la connaissance des fonctions
$p(\Psi)$ et $\bmf(\Psi)$. Les équations de diffusion (\ref{25}),
(\ref{55}), (\ref{30}) et~(\ref{52}) (ou~(\ref{58}), (\ref{59}),
(\ref{60}) et (\ref{52}) pour le modèle diagonal) permettent de calculer $N_e(\chi,t)$, $\sigma_e(\chi,t)$, $\sigma_i(\chi,t)$,
et $\Psi'(\chi,t)$. Les équations (\ref{51}) et~(\ref{62}) permettent de calculer $f(\chi)$ et $V'(\chi)$ à chaque instant. Les
dérivées ${\sfrac{\partial p}{\partial \Psi}}$ et ${\sfrac{\partial f^2}{\partial \Psi}}$
peuvent alors se déduire de l'ensemble des équations de diffusion 1-D.

Mais d'autre part, pour résoudre ces équations de diffusion, il est nécessaire de connaître les coefficients géométriques $\<\nabla^2 \chi\>$,
$\<r^{-2}\nabla^2 \chi\>$ et $\<r^{-2}\>$, qui sont déduits des calculs d'équilibre.

Les équations de diffusion ne nécessitent aucune condition au bord sur l'axe magnétique,
puisque $C_1$, $C_2$ et $C_3$ s'annulent avec $V'$ pour $\chi=0$. Sur le bord du plasma ($\chi=\chi_{\max}$), les valeurs de la densité électronique $n_e$ et des températures électronique et ionique $T_e$ et~$T_i$ sont spécifiées. L'équation (\ref{62}) est une équation différentielle du premier ordre pour $f$ par rapport à $\chi$; sa condition au bord est donnée par la valeur de la fonction $f$ sur le bord du plasma, qui est égale au champ toroïdal externe appliqué $B_0$ multiplié par $R_0$.

La condition (\ref{63}) est une condition de Dirichlet au bord pour $\Psi'$.
En résumé, nous avons
\begin{equation}\label{64}
\begin{aligned}
n_e(\chi_{\max})&=n_e^0,& T_e(\chi_{\max})&=T_e^0,\\ T_i(\chi_{\max})&=T_i^0,& f(\chi_{\max})&=R_0B_0,
\end{aligned}
\end{equation}
et
\begin{equation}
\Psi'(\chi_{\max})=-2\pi \mu_0\frac{I_p}{C_2}.
\end{equation}

Le couplage entre l'équation d'équilibre 2D de Grad-Shafranov et les équations de transport1D, qui dépend de la géométrie des lignes de flux est appelé modèle de transport \hbox{1 1/2 D}. Il permet de suivre l'évolution quasi-statique de l'équilibre à l'échelle de temps du transport (voir \cite{Blum-Lefoll}). Il est bien sûr nécessaire d'ajouter dans les équations de transport les termes sources dus à un chauffage additionnel (injec\-tion de neutres, chauffage haute fréquence HF) ainsi que, dans la loi d'Ohm, le courant engendré par le système HF, comme réalisé par exemple dans le logiciel CRONOS (voir \cite{Artaud}).

\subsubsection*{Remerciements}
L'auteur remercie chaleureusement Cédric Boulbe et Claude Sabbah pour leur aide précieuse dans la préparation de ce texte.

\backmatter
\bibliographystyle{jepalpha+eid}
\bibliography{xups15-03}
\end{document}
