Real-Time Equilibrium Reconstruction in a Tokamak

J. Blum    C. Boulbe and B. Faugeras Laboratoire J.A. Dieudonné, UMR 6621, Université de Nice Sophia Antipolis, Parc Valrose, 06108 Nice Cedex 02, France,,

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.

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 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 [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):


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 [8] and the equilibrium equation for the plasma is


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


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.

From (8) and the second relation of (1) we obtain the following expression for :


where and are the poloidal and toroidal components respectively of , and the operator is defined by


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


Relation (5) combined with the expression (9) implies that and are colinear, and therefore is likewise constant on each magnetic surface


The equilibrium relation (3) combined with the expression (8) and (9) for and implies that:


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.

From (9) it is clear that right-hand 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


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)

Figure 1: Definition of the plasma boundary (thick blue line). Left, JET (Joint European Torus) example, X-point configuration. Right, TORE SUPRA (the CEA-EURATOM Tokamak at Cadarache) example, limiter configuration (the limiter is represented by the black line). The thin blue lines represent iso-contours of .

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 .

    • , which corresponds to the component of the magnetic poloidal field, measured by the magnetic probes (see Fig. 2), which is tangent to the vacuum vessel. Indeed from Eq. (8) the tangential component of is equal to the normal component of .

  • 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 :


Figure 2: Left: the straight green lines represents the chords used for polarimetry and interferometry measurements. Right: part of the vacuum vessel. At the bottom middle an example of finite element mesh used for numerical simulations (see next Section).

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 [5].

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 [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,


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:


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:


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 .

In order to solve this problem we use an iterative algorithm based on fixed-point iterations for Eq. (21) and the normal equation of Eq. (22).

3.3.1 Algorithm

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.

Figure 3: An output of EQUINOX. The plasma is in an X-point configuration. On the left column the identified , and functions as well as the toroidal current density and the safety factor are displayed in terms of and of (in the equatorial plane).

4 Conclusion

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 [14].

5 Aknowledgements

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.



  • [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 Ill-posed problems (Winston, Washington, D.C.)
  • [10] Ciarlet P 1980 The Finite Element Method For Elliptic Problems (North-Holland)
  • [11] Bosak K 2001 Real-time numerical identification of plasma in tokamak fusion reactor Master’s thesis University of Wroclaw, Poland URL
  • [12] Blum J, Bosak K and Joffrin E 2004 12th ICPP International Congress on plasma physics (Nice) URL
  • [13] Bosak K, Blum J, Joffrin E and Sartori F 2003 EPS Conference on plasma physics (Saint-Petersbourg)
  • [14] Joffrin et al E 2003 Plasma Physics and Controlled Fusion 45 A367–A383
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description