Real-Time Equilibrium Reconstruction in a Tokamak
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 non-linear source in the 2D Grad-Shafranov 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 real-time using a finite element method, a non-linear fixed-point algorithm and a least-square optimization procedure.
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 Grad-Shafranov equation [1, 2, 3]. The right-hand side of this equation is a non-linear source which represents the toroidal component of the plasma current density.
An important problem is the identification of this non-linearity [4, 5, 6]. The aim of this paper is to present a method for real-time 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 .
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):
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
where represents the mean velocity of particles and the mass density. At the resistive time-scale the first term can be neglected  and the equilibrium equation for the plasma is
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
Concerning the toroidal component we define by
where is the unit vector in the toroidal direction, and is the diamagnetic function. The magnetic field can be written as:
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.
where and are the poloidal and toroidal components respectively of , and the operator is defined by
In the plasma region, relation (4) implies that and are colinear, and therefore is constant on each magnetic surface. This can be denoted by
which leads to the so-called Grad-Shafranov equilibrium equation:
where is the linear elliptic operator given by (10) in which is equal to the magnetic permeability of the vacuum.
In the vacuum, the magnetic flux satisfies
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
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 cross-section of the vacuum vessel, the final equations governing the behaviour of inside the vacuum vessel, are:
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 real-time identification of the plasma current i.e. the non-linear 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 :
3.2 Statement of the inverse problem
The numerical identification problem is formulated as a least-square minimization with a Tikhonov regularization. The cost function is defined as:
and , , and are weighting parameters enabling to give more or less importance to the corresponding experimental measurements .
The inverse problem of the determination of and is ill-posed. Hence a regularization procedure can be used to transform it into a well-posed one . 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,
3.3 Numerical identification
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 …) . Let be the vector of defined by . With these notations the discretization of problem (20) can be written as follows:
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:
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 .
At the -th iteration, and are given. The non-linear 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 fixed-point 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 fixed-point 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 CEA-EURATOM Tokamak at Cadarache). Figure 3 shows a graphical output of Equinox. With all these techniques it is possible to follow the quasi-static evolution of the plasma equilibrium, either in TORE SUPRA or JET configurations, with free boundaries defined either by limiter contact or with an X-point. It is also possible to simulate ITER configurations.
We have presented an algorithm for the identification of the current
density profile in Grad-Shafranov equation from experimental measurements.
The decomposition of the unknown functions and in a reduced basis makes it possible to do the reconstruction in real-time.
The choice of this reduced basis must still be improved and optimized (robustness, precision, …).
Real-time reconstruction makes possible future real-time control of the current profile .
The authors are very thankful to K. Bosak for his enormous work in developing the real-time EQUINOX software and to E. Joffrin and S. Bremond for their contributions into the domain of real-time identification and control.
-  Grad H and Rubin H 1958 2nd U.N. Conference on the Peaceful uses of Atomic Energy vol 31 (Geneva) pp 190–197
-  Shafranov V 1958 Soviet Physics JETP 6 1013
-  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)
-  Lao L, Ferron J, Geoebner R, Howl W, St John H, Strait E and Taylor T 1990 Nuclear Fusion 30 1035
-  Blum J, Lazzaro J, O’Rourke J, Keegan B and Stefan Y 1990 Nuclear Fusion 30 1475
-  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
-  Zwingmann W and Stubberfield P 2000 27th EPS Conference on Contr. Fusion and Plasma Phys. (Budapest)
-  Grad H and Hogan J 1970 Physical review letters 24 1337–1340
-  Tikhonov A and Arsenin V 1977 Solutions of Ill-posed problems (Winston, Washington, D.C.)
-  Ciarlet P 1980 The Finite Element Method For Elliptic Problems (North-Holland)
-  Bosak K 2001 Real-time 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
-  Blum J, Bosak K and Joffrin E 2004 12th ICPP International Congress on plasma physics (Nice) URL http://fr.arxiv.org/abs/physics/0411181
-  Bosak K, Blum J, Joffrin E and Sartori F 2003 EPS Conference on plasma physics (Saint-Petersbourg)
-  Joffrin et al E 2003 Plasma Physics and Controlled Fusion 45 A367–A383