RealTime Equilibrium Reconstruction in a Tokamak
Abstract
This paper deals with the numerical reconstruction of the plasma current density in a Tokamak and of its equilibrium. The problem consists in the identification of a nonlinear source in the 2D GradShafranov equation, which governs the axisymmetric equilibrium of a plasma in a Tokamak. The experimental measurements that enable this identification are the magnetics on the vacuum vessel, but also polarimetric and interferometric measures on several chords, as well as motional Stark effect or pressure measurements. The reconstruction can be obtained in realtime using a finite element method, a nonlinear fixedpoint algorithm and a leastsquare optimization procedure.
1 Introduction
The problem of the equilibrium of a plasma in a Tokamak is a free boundary problem in which the plasma boundary is defined either by its contact with a limiter or as being a magnetic separatrix. Inside the plasma, the equilibrium equation in an axisymmetric configuration is called GradShafranov equation [1, 2, 3]. The righthand side of this equation is a nonlinear source which represents the toroidal component of the plasma current density.
An important problem is the identification of this nonlinearity [4, 5, 6]. The aim of this paper is to present a method for realtime identification from experimental measurements, such as magnetic measurements on the vacuum vessel, polarimetric measurements (integrals of the magnetic field over several chords), MSE (Motional Stark Effect) and pressure measurements. The pressure is supposed to be isotropic. For the anisotropic pressure case, one can refer to [7].
The next section is devoted to the mathematical modelling of the equilibrium problem in axisymmetric configurations. The inverse reconstruction problem is adressed in the last section.
2 Mathematical modelling of axisymmetric equilibrium of the plasma in a Tokamak
The equations which govern the equilibrium of a plasma in the presence of a magnetic field are on the one hand Maxwell’s equations and on the other hand the equilibrium equations for the plasma itself.
The magnetostatic Maxwell’s equations as follows are satisfied in the whole of space (including the plasma):
(1) 
where represents the magnetic field, is the magnetic permeability and is the current density. The first relation of (1) is the equation of conservation of magnetic induction and the second one is Ampere’s Theorem.
The momentum equation for a plasma is
(2) 
where represents the mean velocity of particles and the mass density. At the resistive timescale the first term can be neglected [8] and the equilibrium equation for the plasma is
(3) 
This equation (3) means that the plasma is in equilibrium when the force due the kinetic pressure is equal to the Lorentz force of the magnetic pressure . We deduce immediately from (3) that
(4) 
(5) 
Thus for a plasma in equilibrium the field lines and the current lines lie on isobaric surfaces (); these surfaces, generated by the field lines, are called magnetic surfaces. In order for them to remain within a bounded volume of space it is necessary that they have a toroidal topology. These surfaces form a family of nested tori. The innermost torus degenerates into a curve which is called the magnetic axis.
In a cylindrical coordinate system (where is the major axis of the torus) the hypothesis of axial symmetry consists in assuming that the magnetic field is independent of the toroidal angle . The magnetic field can be decomposed as , where is the poloidal component and is the toroidal component. From equation (1) one can define the poloidal flux such that
(6) 
Concerning the toroidal component we define by
(7) 
where is the unit vector in the toroidal direction, and is the diamagnetic function. The magnetic field can be written as:
(8) 
According to (8), in an axisymmetric configuration the magnetic surfaces are generated by the rotation of the flux lines around the axis of the torus.
From (8) and the second relation of (1) we obtain the following expression for :
(9) 
where and are the poloidal and toroidal components respectively of , and the operator is defined by
(10) 
Expressions (8) and (9) for and are valid in the whole of space since they involve only Maxwell’s equations and the hypothesis of axisymmetry.
In the plasma region, relation (4) implies that and are colinear, and therefore is constant on each magnetic surface. This can be denoted by
(11) 
Relation (5) combined with the expression (9) implies that and are colinear, and therefore is likewise constant on each magnetic surface
(12) 
The equilibrium relation (3) combined with the expression (8) and (9) for and implies that:
(13) 
which leads to the socalled GradShafranov equilibrium equation:
(14) 
where is the linear elliptic operator given by (10) in which is equal to the magnetic permeability of the vacuum.
From (9) it is clear that righthand side of (14) represents the toroidal component of the plasma current density. It involves functions and which are not directly measured inside the plasma.
In the vacuum, the magnetic flux satisfies
(15) 
The equilibrium of a plasma in a domain representing the vacuum region is a free boundary problem. The plasma free boundary is defined either by its contact with a limiter (outermost flux line inside the limiter) or as being a magnetic separatrix (hyperbolic line with an point, ). The region containing the plasma is defined as
(16) 
where either in the limiter configuration or in the point configuration (see Fig. 1)
Assuming Dirichlet boundary conditions, , are given on which is the poloidal crosssection of the vacuum vessel, the final equations governing the behaviour of inside the vacuum vessel, are:
(17) 
with and , in (this normalized flux is introduced so that and are defined on the fixed interval ), is the characteristic function of .
The aim of the following section of this paper is to provide a method for the realtime identification of the plasma current i.e. the nonlinear functions and in the elliptic equation (17).
3 The inverse problem
3.1 Experimental measurements
The given experimental measurements are:

the magnetic measurements

, given by the flux loops (see Fig. 2). Thanks to an interpolation between points these measurements provide the Dirichlet boundary condition .


the polarimetric measurements which give the Faraday rotation of the angle of infrared radiation crossing the section of the plasma along several chords :
where represents the electronic density which is approximately constant on each flux line, is the component of the poloidal field tangent to and represents the normal derivative of with respect to .

the interferometric measurements which give the density integrals over the chords

the kinetic pressure measurements obtained from density and temperature measurements, for instance in the equatorial plane:

and MSE (Motional Stark Effect) angle measurements taken at different points :
with
3.2 Statement of the inverse problem
The numerical identification problem is formulated as a leastsquare minimization with a Tikhonov regularization. The cost function is defined as:
(18) 
with
and , , and are weighting parameters enabling to give more or less importance to the corresponding experimental measurements [5].
The inverse problem of the determination of and is illposed. Hence a regularization procedure can be used to transform it into a wellposed one [9]. The Tikhonov regularization term constrains the function , and to be smooth enough and reads:
where , and are the regularizing parameters.
It should be noticed that the electronic density does not intervene in Eq. (17). However as soon as we want to use the polarimetric measurements it is necessary to include (and hence interferometry) in the identification procedure. The inverse problem can finally be formulated as,
(19) 
3.3 Numerical identification
Problem (17) is solved using a finite element method [10]. Let and denote the usual Sobolev spaces. The finite element approximation is based on the following weak formulation:
(20) 
Classically is approximated using triangles by a polygonal domain , the space is approximated by a space of finite dimension . A finite element method is used, in which functions of are affine over each triangle and continuous on the whole domain.
Let denote the finite element stiffness matrix. Let us also (abusively) denote by the components of the magnetic flux fonction approximated in .
The unknown functions , and are approximated by a decomposition in a reduced basis
This basis can be made of different types of functions (polynomials, splines, wavelets, etc …) [6]. Let be the vector of defined by . With these notations the discretization of problem (20) can be written as follows:
(21) 
Where denotes the “plasma current matrix”, and is the stiffness matrix modified in order to impose the Dirichlet boundary condition represented by .
The discrete inverse optimization problem is:
(22) 
where is the observation operator. The quantity represents the outputs of the model corresponding to the experimental measurements, given in a previous subsection, denoted by . The matrix represents the regularization terms. The first term of in Eq. (22) corresponds to and the second to .
In order to solve this problem we use an iterative algorithm based on fixedpoint iterations for Eq. (21) and the normal equation of Eq. (22).
3.3.1 Algorithm
At the th iteration, and are given. The nonlinear mapping is approximated by the affine relation
and the cost function to be minimized by
with obvious notations. The normal equation
is solved to update to . Then a fixedpoint iteration for Eq. (21) enables the update of to
Since the algorithm is usually initialized with the equilibrium at a previous time step, two or three fixedpoint iterations are usually enough to ensure convergence.
3.3.2 Equinox software
Based on the algorithm presented above, a C++ software, called EQUINOX [11, 12, 13] has been developed in collaboration with the Fusion Department at Cadarache, and has been implemented for JET (Joint European Torus) and for TORE SUPRA (the CEAEURATOM Tokamak at Cadarache). Figure 3 shows a graphical output of Equinox. With all these techniques it is possible to follow the quasistatic evolution of the plasma equilibrium, either in TORE SUPRA or JET configurations, with free boundaries defined either by limiter contact or with an Xpoint. It is also possible to simulate ITER configurations.
4 Conclusion
We have presented an algorithm for the identification of the current
density profile in GradShafranov equation from experimental measurements.
The decomposition of the unknown functions and
in a reduced basis makes it possible to do the reconstruction in realtime.
The choice of this reduced basis must still be improved and optimized
(robustness, precision, …).
Realtime reconstruction makes possible
future realtime control of the current profile [14].
5 Aknowledgements
The authors are very thankful to K. Bosak for his enormous work in developing the realtime EQUINOX software and to E. Joffrin and S. Bremond for their contributions into the domain of realtime identification and control.
References
References
 [1] Grad H and Rubin H 1958 2nd U.N. Conference on the Peaceful uses of Atomic Energy vol 31 (Geneva) pp 190–197
 [2] Shafranov V 1958 Soviet Physics JETP 6 1013
 [3] Mercier C 1974 The MHD approach to the problem of plasma confinement in closed magnetic configurations Lectures in Plasma Physics (Luxembourg: Commission of the European Communities)
 [4] Lao L, Ferron J, Geoebner R, Howl W, St John H, Strait E and Taylor T 1990 Nuclear Fusion 30 1035
 [5] Blum J, Lazzaro J, O’Rourke J, Keegan B and Stefan Y 1990 Nuclear Fusion 30 1475
 [6] Blum J and Buvat H 1997 IMA Volumes in Mathematics and its Applications, Large Scale Optimization with applications, Part 1: Optimization in inverse problems and design vol 92 ed Biegler, Coleman, Conn and Santosa pp 17–36
 [7] Zwingmann W and Stubberfield P 2000 27th EPS Conference on Contr. Fusion and Plasma Phys. (Budapest)
 [8] Grad H and Hogan J 1970 Physical review letters 24 1337–1340
 [9] Tikhonov A and Arsenin V 1977 Solutions of Illposed problems (Winston, Washington, D.C.)
 [10] Ciarlet P 1980 The Finite Element Method For Elliptic Problems (NorthHolland)
 [11] Bosak K 2001 Realtime numerical identification of plasma in tokamak fusion reactor Master’s thesis University of Wroclaw, Poland URL http://panoramix.ift.uni.wroc.pl/~bosy/mgr/mgr.pdf
 [12] Blum J, Bosak K and Joffrin E 2004 12th ICPP International Congress on plasma physics (Nice) URL http://fr.arxiv.org/abs/physics/0411181
 [13] Bosak K, Blum J, Joffrin E and Sartori F 2003 EPS Conference on plasma physics (SaintPetersbourg)
 [14] Joffrin et al E 2003 Plasma Physics and Controlled Fusion 45 A367–A383