Section: New Results

Optimisation and control for magnetic fusion plasmas

Evolutive equilibrium and transport coupling and optimization of scenarii

Participants : Jacques Blum, Cédric Boulbe, Afeintou Sangam, Gael Selig, Blaise Faugeras, Holger Heumann.

Research of optimal trajectories for the monitoring of Tokamak discharges

The direct equilibrium code CEDRES++ in its static version (resp.dynamic) computes for externally applied PF currents (resp. voltages) and given plasma current density profile the (resp. evolution of the) poloidal flux and the plasma free boundary. The research of optimal trajectories is the corresponding inverse problem : find externally applied currents (resp. voltages), such that the plasma reaches a certain desired state. This desired state is mainly (resp. the evolution of) a prescribed plasma boundary. We formulate these inverse problems as so-called optimal control problems, where the PF currents (voltages) are the so-called control variables and the poloidal flux the so-called state variable. Optimal control problems are optimization problems with PDE (partial differential equations) constraints. In our case, the Grad-Shafranov equation is the constraint and the functional to be minimized is a cost-function that measures the mismatch between the computed plasma boundary and the desired plasma boundary. The Sequential Quadratic Programming (SQP) method is known to be a very efficient algorithm for solving non-linear constrained optimization problems. We implemented in CEDRES++ the SQP method for the two cases of finding either currents or voltages that corresponds to a desired boundary or a desired evolution of the boundary. These implementations are built on the orignal Newton methods for the direct non-linear problems. For optimization problems it is of great importance that the Newton methods are 'real' Newton methods in the sense that the Newton matrices are real derivatives. In the original implementation of CEDRES++ these matrices were the discretization of analytic derivatives of the non-linear operators, hence not derivatives of the discrete problem. We had to rewrite large parts of the code to eliminate this problem. Further, we added an interface to the linear solver library UMFPACK. For the current mesh resolution level, the performance of this linear solver for the stationary problems, both in the direct and in the inverse versions, is superior to iterative linear solvers. In the case of the inverse non-stationary problem, the problem of finding voltages that correspond to a desired evolution, the memory requirements forbid the use of UMFPACK. There, we used Conjugate Gradient-type iterations. In the future, we will have to investigate if other types of iterative solvers are suitable and allow a certain parallelism that will speed up the simulation time.

A new method of coupling equilibrium and resistive diffusion equations

In the framework of Gael Selig's PhD thesis, the resistive diffusion equation has been incorporated in the evolutive equilibrium system of CEDRES++. This equation has as unknown variable the derivative of the poloidal flux with respect to the averaged minor radius of the magnetic surface. This choice was made instead of the poloidal flux itself because this is the quantity directly involved in the averaged Grad-Shafranov equation used to compute the FF' term and thus this allows us not to perform a supplementary numerical differentiation which might introduce some numerical instability. An algorithm based on a successive prediction and correction method is proposed in order to ensure the consistency between the evolution of the 2D poloidal flux in the equilibrium equation and the evolution of the poloidal flux in the 1D resistive diffusion equation. The algorithm guarantees that at the end of each time step the total plasma current Ip and the mean radius of the plasma have the same values in both systems (see fig.2). The convergence of this new code (called CEDRES-DIF) has been numerically validated and the method has been successfully compared by G. Selig to the CEDRES-CRONOS coupled code which uses another coupling algorithm.

Introduction of halo currents in the equilibrium resolution

When VDE (Vertical Displacement events) instabilities occur in a Tokamak, currents flow from the plasma to the machine vessel structures, and then return to the plasma. These currents are called halo currents . In turn, these currents induce forces on the wall when crossing with Tokamak poloidal and toroidal magnetic fields. Moreover, when VDE instabilities take place, the plasma hits the wall with all its energy. Therefore, it is worth understanding the contribution of halo currents to total plasma current and other related plasma parameters, particularly the distribution, magnitude, and temporal evolution of halo currents for large scale machine such as ITER. Even if halo currents are actually 3D phenomena, it is important to take into account their effects in 2D free boundary equilibrium codes. In halo region, the pressure can be considered as negligeable so that the current follows the magnetic field lines. The magnetic field satisfies the force free equation $jxB=0,\nabla .B=0$ which can be rewritten

in an axisymmetric configuration. The function $f_H\left(\Psi \right)$ is supposed to be known. This simple model has been implemented in CEDRES++ and first tests have been done. This first simplified model has to be improved to get more realistic simulations and to be validated. The choice of the function $f_H$, the value of the total halo current, the geometry and the size of the halo region need to be enhanced with respect to experimental data.

Equilibrium reconstruction and current density profile identification

Participants : Jacques Blum, Cédric Boulbe, Blaise Faugeras.

EQUINOX is a real-time equilibrium reconstruction code. It solves the equation satisfied by the poloidal flux in a computation domain, which can be the vacuum vessel for example, using a P1 finite element method and solves the inverse problem of the identification of the current density profile by minimizing a least square cost-function. It uses as minimal input the knowledge of the flux and its normal derivative on the boundary of the computation domain. It can also use supplementary constraints to solve the inverse problem: interferometric, polarimetric and MSE measurements. Part of the work reported here has been done in the frame of a RTM-JET contract.

Direct use of the magnetic measurements

Equinox was not originally designed to take as magnetic inputs directly the magnetic measurements, as it should be the case in the ITM, but some outputs from the real-time codes Apolo at ToreSupra and Xloc at JET. These codes provide Equinox with the values of the flux and its normal derivative on a closed contour defining the boundary of the computation domain (this contour can be the limiter for example). As a consequence the main difficulty arising in the objective of integrating the code Equinox in the ITM structure was to interpolate between the magnetic measurements (flux loops and poloidal B-probes) with a machine independent method. This has already been achieved by using toroidal harmonic functions, as a basis for the decomposition of the poloidal flux in the vacuum region, in complement to the contribution of the PF coils. This method can provide an alternative tool, comparable to APOLO (for Tore Supra) and FELIX (for JET), to compute the plasma boundary in real time from the magnetic measurements. Some twin experiments for WEST (Tore Supra upgrade) have been successfully conducted. In a first step the equivalents of magnetic measurements were generated using the FBE code CEDRES++. In a second step these measurements were used by the toroidal harmonics algorithm to reconstruct the plasma boundary. The results are very promising and the work on this subject is ongoing for JET.

Boundary conditions for EQUINOX

In the present version of EQUINOX the boundary condition is a flux condition (Dirichlet boundary condition) and the tangential component of the poloidal field is incorporated in the cost-function to be minimized. This is a constant criticism which is made on EQUINOX. The idea was to inverse these two boundary conditions in order to determine if this choice is determinant in the results. We tried to use the tangential poloidal field (Neumann boundary condition for the flux) as boundary condition for the boundary value problem, and to put the flux (or its tangential derivative linked to the normal component of the poloidal field) in the cost function. However no convincing results could be obtained because the numerical resolution of the boundary value problem associated with Neumann boundary conditions proved to be unstable. This might be explained by the fact that a compatibility condition has to be satisfied between the Neumann conditions and the current density in the plasma which evolves during the mixed fixed-point and optimization iterations.

Induced currents in EQUINOX

In a disruption when the total plasma current disappears, there are very important induced currents, for example in the toroidal pumped limiter. These currents are in the domain of resolution of EQUINOX. Therefore it is necessary to take them into account in the resolution of the equilibrium reconstruction problem. This has been tested on a Tore Supra disruption case. The mesh generation has been modified in order to incorporate the real structure of the limiter. The structure of the equations being solved in the code also had to be modified in order to take into account the measured induced currents.