REAL-TIME EQUILIBRIUM RECONSTRUCTION IN A TOKAMAK

# 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 jblum@unice.fr, boulbe@unice.fr, faugeras@unice.fr
###### 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 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):

 ⎧⎪⎨⎪⎩∇⋅B=0∇×(Bμ)=j (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

 ρdudt=j×B−∇p (2)

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

 ∇p=j×B (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

 B⋅∇p=0 (4)
 j⋅∇p=0 (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

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Br=−1r∂ψ∂zBz=1r∂ψ∂r (6)

Concerning the toroidal component we define by

 Bϕ=freϕ (7)

where is the unit vector in the toroidal direction, and is the diamagnetic function. The magnetic field can be written as:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩B=Bp+BϕBp=1r[∇ψ×eϕ]Bϕ=freϕ (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 :

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩j=jp+jϕjp=1r[∇(fμ)×eϕ]jϕ=(−Δ∗ψ)eϕ (9)

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

 Δ∗.=∂∂r(1μr∂.∂r)+∂∂z(1μr∂.∂z) (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

 p=p(ψ) (11)

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

 f=f(ψ) (12)

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

 ∇p=−Δ∗ψr∇ψ−fμ0r2∇f (13)

 −Δ∗ψ=rp′(ψ)+1μ0r(ff′)(ψ) (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 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

 −Δ∗ψ=0 (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

 Ωp={x∈Ω, ψ(x)≥ψb} (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 cross-section of the vacuum vessel, the final equations governing the behaviour of inside the vacuum vessel, are:

 ⎧⎨⎩−Δ∗ψ=[rA(¯ψ)+1rB(¯ψ)]χΩpin Ωψ=hon Γ (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 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 :

 ∫Cine(¯ψ)B∥dl=∫Cine(¯ψ)r∂ψ∂ndl=αi

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

 ∫Cine(¯ψ)dl=βi
• the kinetic pressure measurements obtained from density and temperature measurements, for instance in the equatorial plane:

 p(r,0)=pd(r)
• and MSE (Motional Stark Effect) angle measurements taken at different points :

 m(Br,Bz,Bϕ)i=γi

with

 tan(m(Br,Bz,Bϕ))=a1Br+a2Bz+a3Bϕa4Br+a5Bz+a6Bϕ

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

 J(A,B,ne)=J0+K1J1+K2J2+K3J3+K4J4+Jϵ (18)

with

 J0=∑i(1r∂ψ∂n(Ni)−gi)2J1=∑i(∫Ciner∂ψ∂ndl−αi)2J2=∑i(∫Cinedl−βi)2J3=∫RmaxRmin(p(r,0)−pd(r))2drJ4=∑i(m(Br,Bz,Bϕ)i−γi)2

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:

 Jϵ=ϵ1∫10[A′′(x)]2dx+ϵ2∫10[B′′(x)]2dx+ϵ3∫10[n′′e(x)]2dx

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,

 {Find A∗, B∗, n∗e such that:J(A∗,B∗,n∗e)=infJ(A,B,ne) (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:

 ⎧⎪⎨⎪⎩Find ψ∈H1(Ω), such that ψ=h on Γ, and∫Ω1μ0r∇ψ⋅∇vdx=∫Ωp[rA(¯ψ)+1rB(¯ψ)]vdx∀v ∈V (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

 A(x)=∑iaiϕi(x)B(x)=∑ibiϕi(x)ne(x)=∑iciϕi(x)

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:

 {Given u∈R3m, solve the fixed−point equation~Kψ=D(ψ)u+h (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:

 ⎧⎪⎨⎪⎩Find u minimizing:J(u)=∥C(ψ)ψ−k∥2+uTΛuwith ψ satisfying (???) (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 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

 ψ=~K−1[D(ψn)u+h]

and the cost function to be minimized by

 J(u)=∥C(ψn)ψ−k∥2+uTΛu=∥C(ψn)~K−1D(ψn)u+(C(ψn)~K−1h−k)∥2+uTΛu=∥Enu+Fn∥2+uTΛu

with obvious notations. The normal equation

 (ETnEn+Λ)u=−ETnFn

is solved to update to . Then a fixed-point iteration for Eq. (21) enables the update of to

 ψn+1=~K−1[D(ψn)un+1+h].

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.

## 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.

## 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 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 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 (Saint-Petersbourg)
• [14] Joffrin et al E 2003 Plasma Physics and Controlled Fusion 45 A367–A383
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