Analysis of vertical displacement events in tokamaks removing the evolutionary equilibrium assumption

A key parameter for the vertical stabilisation system of future fusion devices like ITER is the maximum allowable vertical displacement which can be recovered by the magnetic control system (Gribov 2015 Nucl. Fusion 55 073021). The computation of this figure of merit is based on MHD evolutionary equilibrium models: the perturbation of the vertical position implicitly considers a perturbation of the external passive currents which guarantees MHD equilibrium. In this paper, we quantify the consequences of considering non-equilibrium perturbations of the vertical position. In this case, one has to specify separately the initial perturbation of the vertical position and of the external currents. We show in particular under which circumstances the evolution of the vertical position, as predicted by the evolutionary equilibrium models, coincides with that obtained considering also plasma inertia.


Introduction
Vertically elongated plasma cross-sections are today standard for tokamak operation, due to the great improvement of the poloidal beta and energy confinement time for such configurations [1]. The actual set up of a vertically elongated plasma cross-section requires the application of a quadrupole magnetic field, which determines an axisymmetric instability: any small perturbation of the vertical position of the plasma column grows exponentially in time [2]. The modelling and control of the plasma vertical motion is then of paramount importance for the operation of present tokamaks and for the design of future devices. * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The main developments around this subject date back to the years 1980's and 1990's, using different approaches depending on the application of interest. Some authors proposed to analyse the problem of the vertical motion based on the consistent solution of Magneto-Hydro-Dynamics (MHD) equations [3,4]. These models, although very detailed, are generally difficult to solve, and are mainly meant for aposteriori interpretation of experiments rather than for realtime control purposes. Great insight into the problem was gained approximating the tokamak plasma column with one or multiple rigid filaments carrying altogether the net toroidal plasma current [5][6][7][8][9]. In particular the seminal paper [8] evidenced that the inertia offered to the motion by the eddy currents induced in the structures surrounding the plasma, overwhelms the plasma mechanical inertia in many situations of practical interest. A synthetic figure of merit, the stability margin, was introduced to quantify the validity of the quasi-static assumption in the analysis of the tokamak plasma vertical motion. This definition, based on a single rigid filament approximation in interaction with a single external circuit, was extended both to the case where multiple degrees of freedom for the passive currents are present [10], and to the case of mass-less deformable plasmas [11][12][13]. In [14] it was observed that the 3D details of external conductors can impact significantly the actual axisymmetric vertical evolution of the plasma column, leading to the development of deformable quasi-static models of the plasma self-consistently coupled with fully volumetric 3D models of external conductors [15], today useful for the analysis of stability properties and disruptions in tokamak devices with significant three-dimensional features [16].
Besides the wide use of quasi-static models in the literature, the consequences of neglecting plasma mass were investigated in some details only in the earlier works on the subject [3,6,7], but were not addressed in subsequent works, hence taking this point for granted. Aim of the present paper is to come back to this issue, focusing on the projection of initial conditions on the quasi-static MHD model. This is particularly relevant on the one hand for the comparison of evolutionary equilibrium models with experiments and full-dynamic simulations and on the other hand because during the magnetic control design phase of tokamak devices a range of possible initial perturbations ∆z 0 is generally postulated on physical grounds, implicitly assuming that this is a quasi-static perturbation.
The paper is organized as follows. In section 2 we recall a simplified rigid-filament model for the plasma surrounded by external passive conductors. Taking the linearized version of such model, and considering only the vertically unstable eigenmode for the passive currents, in section 3 we provide a qualitative phase portrait of the dynamical system. This is accomplished via a Singular Perturbation Method, which exploits the different mechanical and electrical time-scales of the system. The analysis will reveal the relation existing between the initial conditions of the full dynamic model and the initial conditions to be used in the quasi-static model to trigger correctly the unstable mode of the dynamical system. In section 4 we apply the theory discussed in the previous section, showing the efficacy of our mapping for initial conditions. Section 5 explains how to apply the theory accounting all the possible external currents eigenmodes and draws the conclusion of our study, pointing out to possible future work.

Filamentary rigid model
Here we briefly recall the simplified plasma model that we use in the present paper, which is in the stream of [5][6][7][8][9]. We approximate the plasma column as a rigid, axisymmetric current-carrying ring, which we indicate as plasma ring in the following. Due to these hypotheses, the radial location and the shape of this conductor are of course fixed. The plasma ring conductor can only displace vertically, constrained to be coaxial with the z-axis throughout the whole time evolution of the system. The state-space variables for our system are limited to the vertical position of the plasma column Z p , its velocity V p , plus the currents in the external conductors surrounding the plasma I w . This rigid conductor is surrounded by the magnetic field produced by itself, the equilibrium active coils currents I a and the possible induced currents flowing in passive conductors I w . We recall that, using the rigid body approximation, only external forces play a role in accelerating the ring: the magnetic field produced by the plasma ring itself does not determine any acceleration. A Lagrangian formulation of this electro-mechanical model [17] immediately provides us with the system of non-linear ordinary differential equations: Here we introduced the inductance and resistance matrix of external passive conductors, L w and R w , the mutual inductance between active and passive conductors M w,a , besides the mutual inductance terms between plasma ring and other external conductors, M w,p and M a,p . The latter inductance terms, together with their spatial derivatives, depend clearly on the vertical position of the rigid plasma ring Z p , leading the system to be non-linear. The plasma mass m p and the plasma current I p are here assumed to be constant.
In the following we will linearise model (1) around a static MHD equilibrium condition ( . Moreover, we do not allow here active currents variations, i.e. I a = I 0 a . In this section, we isolate the unstable eigen-mode of the system (i u , v u , z u ) and we assume that the unstable mode for the passive currents is the only relevant one, performing a model order reduction tout-court. The resulting linear wall-single-mode model for the system is given by (2) Here τ u = L u /R u is the electromagnetic time constant for the unstable mode, where L u = i u T · L w · i u and R u = i u T · R w · i u . The stabilizing force per external unit current was defined as while the destabilizing force per unit vertical displacement due to the magnetic field applied by the active coils was indicated by

Asymptotic solution: qualitative phase portrait
Although the third order dynamic linear system (2) can be solved analytically, we will derive here an asymptotic solution, exploiting the different time-scales of the system. The solution, obtained via a singular perturbation method, will highlight the physical properties of the system. While the electromagnetic time-scale τ u was already introduced in previous section, it is convenient to introduce as well the Alfvén time for the instability τ A , which is dominant in the no-wall limit (i.e. for i u = 0): In particular the ratio between Alfvén time and wall time τ A /τ u is usually ≪ 1, due to the extremely small plasma mass, the large plasma current, and the relative proximity of passive conductors to the plasma. The dynamic system (2) can be lumped into the third order Ordinary Differential Equation where we defined the stability margin [8]: We can now provide an asymptotic solution of the linear system (6) with typical methods used to treat multiple time-scale systems [18,19]. We expand in the small parameter ε = τ A /τ u the solution of the linear system, We introduce explicitly two distinct times, the time variable of fast dynamics t 1 = t/ε and the time variable of slow dynamics t 2 = t. The idea is essentially to think the Hamiltonian vector field in the phase space as given by a combination of commuting vector fields: In case of positive stability margin (m u > 0), assuming that the solution is given by the zero-th order approximation of expansion (8), i.e. z(t) ≃ z 0 (t) and that the separation of time variables (9) holds, we get the following asymptotic solution of the linear system: where the growth rate of the unstable mode has been defined as Details about the procedure leading to asymptotic solution (10) are reported in appendix A. The asymptotic dynamics revealed in equation (10) are correct up to first order corrections in τ A /τ u , as further confirmed by the direct computation of eigenvalues in appendix B. The asymptotic solution above reveals few important physics properties of the system. The solution is the sum of an exponentially increasing term, which represents the unstable mode, and two fast oscillatory terms, which are damped on the electromagnetic time scale [3,6,7]. The equilibrium point (z = i u = 0) is a saddle and the unstable mode solicitation is essentially independent on the initial velocity ∆v 0 . From (10) and the definition of stability margin (7), the unstable motion will be directed upwards whenever the force due to the initial current ∆i u,0 overwhelms the ideal wall reaction force for the initial displacement ∆z 0 : Simple algebra reveals that in the limit t → ∞ the system asymptotically tends to a condition of mechanical equilibrium where i 0 is the corresponding asymptotic solution for the current in passive structures. Rigorously we can draw a phase portrait only in the threedimensional (i u , v p , z p ) space. There we would find that the system gets away from the equilibrium position, spiralling around the unstable direction. However the traces of such trajectories in the (i u , z) plane, which would provide rigorously a Poincaré plot, are sketched qualitatively in figure 1, averaging out the fast oscillations. The trace in the (i u , z) plane of the unstable direction is provided by the equilibrium characteristic given by F u I i 0 (t) + F a,0 z z 0 (t) = 0. The actual upward or downward direction of the motion is determined by the initial conditions and whether these are above or below the ideal wall response characteristic, as expressed by equation (12).
The quasi-static approximation for this single-mode wall model results immediately imposing τ A /τ u → 0 in equation (6). The corresponding solution is of course In order to capture correctly the unstable mode excitation, the initial vertical displacement to impose in a quasi static model is related to the real (full-dynamic) initial conditions by as seen by direct comparison of the asymptotic solution (10) with the quasi-static solution (14). Similarly, in terms of initial condition on the current, it is possible to show:

An example of application
As an example of application of the above theory, we use a range of parameters referred to a TCV-like geometry, since the TCV (Tokamakà Configuration Variable) is particularly suited for this kind of analyses [20]. For the specific configuration under analysis, the time constant associated to the unstable current is τ u = 5.3 ms. The Alfvén time associated to the vertical motion is estimated to be more than four order of magnitude smaller, i.e. τ A /τ u ≃ 8.6 × 10 −5 . Further the single mode stability margin is m u = 0.842, making the asymptotic solution (10) robust. The single mode growth rate is calculated as γ u = 1/(m u τ u ) = 224 s −1 . First, we consider an initial perturbation of the vertical position ∆z 0 = +0.5 cm without any perturbation of the velocity and of the unstable current, i.e. ∆v 0 = 0 and ∆i u,0 = 0. Correspondingly in the quasistatic model we shift the plasma according to equation (15), so that the quasi-static shift to provide is ∆z qs,0 = 1.09 cm. The corresponding evolution in the full-dynamic and massless limit are shown in figure 2. This tuning of initial conditions reveals clearly that the full-dynamic solution oscillates around the solution of the quasi-static model (14). The dynamics neglected in the quasi-static model are oscillations in the range of the Alfvén frequency, damped on the electromagnetic time scale, as revealed by the terms in sine/cosine of the asymptotic solution (10). In figure 2, we also show the evolution of the quasi-static model corresponding to the choice ∆z qs,0 = ∆z 0 , which clearly produces a completely different solicitation of the unstable mode.
It is useful to analyse the evolution of the vertically unstable plasma-ring considering a fixed initial condition ∆z 0 for the vertical position and different initial conditions for the unstable wall current. In order to assign physically meaningful initial conditions for the current we imagine the following situation. The vertical position of the plasma ring is displaced In case of a displacement which is fast compared to the electromagnetic time constant of the wall, the initial condition for the wall current corresponds to the purely inductive reaction limit. Indeed from equation (17): On the other hand, for a very slow displacement we recover the no-wall situation, since no current is induced in wall structures, i.e. lim (∆t/τu)→∞ In figure 3, for fixed ∆z 0 , we consider various initial conditions for the unstable wall current, assigned according to equation (17). Substituting the equilibrium constraint F u I ∆i u = −F a,0 z ∆z 0 in equation (17) we can calculate the ∆t/τ u which realizes a quasi-static perturbation. Expecting ∆t/τ u ≃ 1, we may approximate e −∆t/τu ≃ (−∆t/τ u + 2)/e + · · · , hence the quasi static evolution is recovered approximately whenever For the case under exam equation (20) provides ∆t/τ u = 1.510, while the exact value is 1.377. For a motion faster than this value the wall reaction is stronger than what we get in the quasi-static limit: the unstable mode is triggered significantly less for the same ∆z 0 and much more oscillations of the vertical position appear. On the other side, above this threshold, the wall reaction is progressively reduced compared to the quasi-static perturbation, so that the unstable motion is triggered much more, as well as oscillations. For any ∆i u,0 in the range provided by equation (17) the motion has the same sign as ∆z 0 . The only possibility to get an unstable evolution in the opposite direction as compared to ∆z 0 is to have an initial current ∆i u,0 which is beyond the ideal wall reaction limit, e.g. due to fast external perturbations imposed by the vertical control circuit.

Discussion and conclusions
The correct mapping between initial conditions from full dynamic to quasi-static plasma rigid filament models, was identified in equations (15) and (16) for the single-mode description of wall structures. That relation was identified by direct comparison of the asymptotic solution (10) of the dynamic system against the solution (14) of the quasistatic model. Notice that we essentially projected the point (∆z 0 , ∆i 0 ) on the equilibrium hyperplane F I δi + F a,0 z δz = 0 along the direction of the inductive wall response to a plasma displacement L w δi + F i δz = 0. This observation allows to generalize the mapping to the multiple-mode model of passive conductors: ∆i qs,0 = ∆i 0 − L w −1 F I (∆z qs,0 − ∆z 0 ) .
Here we introduced the inductive stability margin for the multiple-mode model of external currents [10,13]: Summing up, we have identified the subspace of fast dynamics of our system with the plane identified by the velocity direction and the ideal wall response direction. It is a question left for future work whether this circumstance can be generalized to more accurate fluid models.
The above discussion highlights that the maximum allowable displacement, as defined in [21], is a quasi-static definition. This means that the displacement ∆z 0 is accompanied by wall currents variations which enforce the plasma mechanical equilibrium. This circumstance shall be stressed in order to avoid discrepancies between the predictions of full-dynamic and mass-less models. When considering a fulldynamic model, including the effect of the plasma inertia, one can in principle impose separately the initial values of plasma displacement and of wall currents. As an example, we supposed to be able to displace the plasma ring of a certain quantity ∆z 0 at constant velocity within a tuneable time interval ∆t. An instantaneous plasma displacement (∆t/τ u → 0) determines an ideal response of wall structures, i.e. ∆i u,0 = −(F u I /L u )∆z 0 . In this case, the unstable eigenmode is not triggered at all, independently on the value of ∆z 0 . On the other side of the range, for a very slow displacement (∆t/τ u → ∞) no external currents are induced in the wall. The quasi-static perturbation lies in between these two limiting cases, being approximately realized when ∆t/τ u is given by equation (20). For ∆t/τ u below this quasi-static limit, the wall reaction to the displacement is stronger, resulting in a lower solicitation of the unstable mode. Oppositely, for ∆t/τ u above the quasi static limit, the induced currents are solicited less and the unstable mode is triggered much more. In both cases, the farther we are from the limit of quasi-static perturbation, the greater is the solicitation of the damped oscillatory behaviour.
During the electromagnetic design of a tokamak, physicists and engineers shall agree not only on the actual value of maximum allowable displacement, but also on the definition of this quantity, to avoid errors while extrapolating physics considerations to engineering requirements. This paper provides some insights precisely about this point and highlights the merits and the limits of quasi-static approximations, which are largely popular in the plasma control environment.
Future work will be addressed to understanding how the simple picture illustrated in the present paper, with reference to rigid plasma motion, can be extended to more general cases, in which the plasma evolution is determined by the solution of nonlinear MHD equations including plasma inertia. In particular, it should be understood how fast transient occurring in the plasma (e.g. thermal quenches, H-to-L transitions, Edge Localized Modes etc) should be treated in evolutionary equilibrium models, in the light of the considerations reported here.

Data availability statement
The data cannot be made publicly available upon publication due to legal restrictions preventing unrestricted public distribution. The data that support the findings of this study are available upon reasonable request from the authors.
A regular perturbation method can be applied to find the rootŝ λ k of the polynomial defined in (30), assuming the expansion of the roots in the small parameter ε: λ k =λ k,0 + ελ k,1 + · · · . (31) The unperturbed rootsλ k,0 are determined immediately from equation (30) in the limit ε → 0: Given the unperturbed roots, we can calculate one by one the correction factors simply plugging (31) into (30) and claiming that the equation has to be satisfied for any ε. In particular, requiring the coefficient multiplying ε 0 to be null we get We can now wrap up the results expressed in (32)-(34), via (31), getting the approximate eigenvalues for the dynamical system: Here we have used the definition λ =λ/ε and the definition of growth rate (11). The eigenvalues expressed by (35) are the correct ones up to first order corrections in ε = τ A /τ u and confirm the asymptotic dynamics which we have found via the multiple time scale method used to derive the asymptotic solution (10). The fact that two different perturbative methods agree in the description of the dynamics of the system gives more confidence in the soundness of the results obtained.