Impact of poloidal field coil current disturbance on tokamak plasma shape

Abstract

This paper studies the impact of Poloidal Field (PF) coil current disturbance on plasma shape in Experimental Advanced Superconductive Tokamak (EAST). First, the concept of the free-boundary equilibrium problem is studied. Based on the non-rigid plasma model, a free-boundary equilibrium solver developed in MATLAB software is presented to calculate the flux distribution and plasma profile of EAST plasma equilibrium, and the results are validated against the experiments. Second, given a typical double-null plasma equilibrium, the responses of plasma shape to the lower PF coil currents are computed using the orthogonal test method. The dependences of the elongation and triangularity of plasma shape on PF coil currents are obtained, showing that the currents in PF6 and PF8 have the greatest impact on triangularity and elongation for lower single-null (LSN) divertor configuration, respectively. And for double-null (DN) divertor plasma configuration, PF5 and PF7 play equally important roles due to the up-down symmetry. The conclusions will be useful for the design of EAST plasma shape feedback control.

Keywords

EAST tokamak plasma shape response free-boundary equilibrium calculation orthogonal test method

1. Introduction

A tokamak reactor is one of the most promising approaches to implementing nuclear fusion [1]. EAST is a D-shaped cross-section tokamak with 12 independent full superconductive PF coils that are symmetrically located outside the reactor vessel [2], as shown in Fig. 1. The plasma is feedback controlled by adjusting the currents in PF coils [3]. Since the PF coils are highly coupled and each is used to provide the magnetic field necessary to maintain and shape the plasma, it is essential to decide the contribution of each PF coil to the control of plasma parameters [4]. Given this problem, in the ongoing design of an optional plasma-shape feedback controller, it is necessary to calculate the responses of the plasma shape parameters (namely, elongation, triangularity, squareness and X-point position) to each PF coils, which can be obtained by solving the free boundary equilibrium problem [5]. More specifically, based on a given plasma equilibrium, recalculating the plasma equilibrium after an artificial disturbance is added to the PF coil current.

The present plasma control system of EAST tokamak is carefully introduced in [6]. As a matter of fact, the decoupling of the plasma shape controller from the vertical stabilization has been achieved by a frequency separation approach [2]. And to decouple the control for plasma parameters, the optimized control vectors are used to distribute on the PF circuits the contribution of each control loop [7]. In order to achieve a high performance of plasma decoupled control for tokamak, many researches have been conducted on response of plasma to PF coils. R. Albanese et al. studied the dynamic responses of plasma boundary and flux expansion to PF coil currents in EAST tokamak [7]. Furthermore, they proposed an integrated control scheme to control both plasma boundary and flux expansion. S.L. Chen et al. compared the different linearized plasma response models while studying the electromagnetic interaction between the plasma and surrounding conductors [8], and the CarMa0 code was used to simulate the tokamak discharges. Q.L. QIU et al. employed the Tokamak Simulation Code (TSC) to compute the responses of vertical position and radial position to the PF coil currents [9]. However, the calculation was complex and took a long time. In this paper, calculation of responses of plasma shape to the PF coil current disturbance is presented by solving free-boundary equilibrium problem. First, A free boundary equilibrium solver based on GAQ plasma model (B.B. Brown, 1977) is developed to calculate the plasma shape and flux distribution in the vacuum vessel by solving the GradShafranov equation [10], the calculation result of EAST shoot 73625 is compared to the experiment to verify the solver. Afterwards, the orthogonal test method is used to choose the combination of PF current disturbances applied in the tests.

2. East plasma equilibrium calculation

2.1. Free boundary equilibrium calculation for EAST

The quasi-static plasma equilibrium is described by the Grad–Shafranov (GS) equation derived from the steady state ideal Magneto-hydrodynamic(MHD) equations, and is given by $\begin{eqnarray}\displaystyle r\displaystyle \frac{\partial }{\partial r}\displaystyle \frac{1}{r}\frac{\partial {\psi}}{\partial r}+\displaystyle \frac{\partial ^{2}{\psi}}{\partial ^{2}z}=-{\mu}_{0}r^{2}P^{\prime }({\psi})-F({\psi})F^{\prime }({\psi})=-{\mu}_{0}rj_{{\varphi}} & & \displaystyle\end{eqnarray}$ (1) where j _𝜑 is toroidal component of the plasma current density distribution in cylindrical coordinates (R, 𝛷, Z), 𝜇₀ is magnetic permeability of free space, 𝜓 is the normalized magnetic flux, P (𝜓) and F (𝜓) are the pressure and toroidal flux profile functions, respectively.

The GS equation is a nonlinear second order partial differential equation, in order to solve it, a free-boundary equilibrium solver is developed in MATLAB software (MathWorks, USA) based on the non-rigid plasma model. The solver has the same kernel algorithm as in the EFIT code. It uses the reconstructed equilibrium of EFIT as input, and for free boundary equilibrium calculation, it can change the initial values of PF coil currents directly in the MATLAB command window, which compared to EFIT is more convenient and user-friendly. The solver divides the computational domain into a 64 × 64 grid to calculate the flux distribution. The plasma is divided into several current filaments and the GS equation is linearized according to the GAQ current profile assumption, which is described as $\begin{eqnarray}\displaystyle j_{{\varphi}}(\overline{{\psi}})=\frac{{\lambda}{\beta}_{0}r}{r_{0}}(1-\overline{{\psi}}^{{\alpha}_{\text{m}}})^{{\alpha}_{\text{n}}}+\frac{{\lambda}(1-{\beta}_{0})r_{0}}{r}(1-\overline{{\psi}}^{{\gamma}_{\text{m}}})^{{\gamma}_{\text{n}}}=rP^{\prime }+\frac{1}{{\mu}_{0}r}FF^{\prime }. & & \displaystyle\end{eqnarray}$ (2) Here, $\bar{{\psi}}=({\psi}-{\psi}_{\text{a}})/({\psi}_{\text{b}}-{\psi}_{\text{a}})$ , where 𝜓_a is the value of 𝜓 on the magnetic axis and 𝜓_b is the value on the plasma boundary. The parameter 𝜆 is a normalization factor specified in terms of the total current I _p. The quantity r ₀ is major radius of the plasma. The quantities 𝛼_m, 𝛼_n, 𝛾_m, 𝛾_n are related to the plasma internal inductance l _i, which describes the peakedness of the current profile. The parameter 𝛽₀ is a scalar quantity related to the poloidal beta.

Fig. 1.

Cross-section of EAST and plasma.

The adjustment of PF coil currents will cause a flux disturbance in the entire region: $\begin{eqnarray}\displaystyle {\delta}\bar{{\psi}}(r,z)=\displaystyle \frac{\partial \bar{{\psi}}^{0}}{\partial r}{\delta}r+{\displaystyle \frac{\partial \bar{{\psi}}^{0}}{\partial z}}{\delta}z. & & \displaystyle\end{eqnarray}$ (3) Here, (δr, δz) is the displacement of the plasma current. The variation of the normalized flux will cause perturbance of plasma current distribution, which is expressed as: $\begin{eqnarray}\displaystyle j_{{\varphi}}^{1}=j_{{\varphi}}^{0}+{\delta}j_{{\varphi}} & & \displaystyle\end{eqnarray}$ (4) $\begin{eqnarray}\displaystyle {\delta}j_{{\varphi}}=\displaystyle \frac{\partial j_{{\varphi}}^{0}}{\partial r}{\delta}r+\displaystyle \frac{\partial j_{{\varphi}}^{0}}{\partial z}{\delta}z & & \displaystyle\end{eqnarray}$ (5) where superscript 0 identifies the result of the previous iteration and superscript 1 identifies the current iteration. According to (2), $\frac{\partial j_{{\varphi}}^{0}}{\partial r}$ and $\frac{\partial j_{{\varphi}}^{0}}{\partial z}$ can be calculated by: $\begin{eqnarray}\displaystyle \frac{\partial j_{{\varphi}}^{0}}{\partial r} & = & \displaystyle \frac{{\lambda}{\beta}_{0}}{r_{0}}(1-\bar{{\psi}}^{0{\alpha}_{\text{m}}})^{{\alpha}_{\text{n}}}-\frac{{\lambda}r_{0}(1-{\beta}_{0})}{r^{2}}(1-\bar{{\psi}}^{0{\gamma}_{\text{m}}})^{{\gamma}_{\text{n}}}+\left(\frac{{\lambda}{\beta}_{\text{0}}r{\alpha}_{\text{n}}}{r_{0}}(1-\bar{{\psi}}^{0{\alpha}_{\text{m}}})^{{\alpha}_{\text{n}}-1}(-{\alpha}_{\text{m}}\bar{{\psi}}^{0{\alpha}_{\text{m}}-1})\right.\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \left.+∼\frac{{\lambda}r_{0}(1-{\beta}_{0}){\gamma}_{\text{n}}}{r}(1-\bar{{\psi}}^{0{\gamma}_{\text{m}}})^{{\gamma}_{\text{n}}-1}(-{{\gamma}_{\text{m}}\bar{{\psi}}}^{0{\gamma}_{\text{m}}-1})\right)\frac{\partial \bar{{\psi}}^{0}}{\partial r}\end{eqnarray}$ (6) $\begin{eqnarray}\displaystyle {\displaystyle \frac{\partial j_{{\varphi}}^{0}}{\partial z}} & = & \displaystyle \left({\displaystyle \frac{{\lambda}{\beta}_{0}r{\alpha}_{\text{n}}}{r_{0}}}(1-\bar{{\psi}}^{0{\alpha}_{\text{m}}})^{{\alpha}_{\text{n}}-1}(-{\alpha}_{\text{m}}\bar{{\psi}}^{0{\alpha}_{\text{m}}-1})+{\displaystyle \frac{{\lambda}r_{0}(1-{\beta}_{0}){\gamma}_{\text{n}}}{r}}(1-\bar{{\psi}}^{0{\gamma}_{\text{m}}})^{{\gamma}_{\text{n}}-1}(-{{\gamma}_{\text{m}}\bar{{\psi}}}^{0{\gamma}_{\text{m}}-1})\right){\displaystyle \frac{\partial \bar{{\psi}}^{0}}{\partial z}}\end{eqnarray}$ (7) Eventually, the plasma current distribution on computational grids can be updated by solving the (4)–(7). The flow chart of the free boundary equilibrium calculation is shown in Fig. 2:

Step-1: Initialize the plasma equilibrium before PF coil current disturbances are added. We load from gfile the data of plasma boundary, magnetic flux distribution, P’ profile (pprime_0) and FF’ profile (ffprime_0). Then calculating the normalized flux distribution. Next, initializing the coefficients of GAQ current model: $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}pprime\_0=P^{^{\prime }}={\displaystyle \frac{{\lambda}{\beta}_{0}}{r_{0}}}(1-\overline{{\psi}}^{{\alpha}_{\text{m}}})^{{\alpha}_{\text{n}}}\\ ffprime\_0=FF^{^{\prime }}={\mu}_{0}r_{0}{\lambda}(1-{\beta}_{\text{0}})(1-\overline{{\psi}}^{{\gamma}_{\text{m}}})^{{\gamma}_{\text{n}}}\end{array}\right. & & \displaystyle\end{eqnarray}$ (8) According to (8), 𝜆 and 𝛽₀ can be computed by: $\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}{\lambda}=r_{0}\cdot pprime\_0(1)+\displaystyle \frac{1}{{\mu}_{0}r_{0}}ffprime\_0(1)\\ {\beta}_{0}=\displaystyle \frac{r_{0}\cdot pprime\_0(1)}{{\lambda}}\end{array}\right. & & \displaystyle\end{eqnarray}$ (9) where pprime_0 (1) and ffprime_0 (1) are the values on magnetic axis where $\bar{{\psi}}=∼0$ . After that, we use MATLAB fitting function “nlinfit” to calculate 𝛼_m, 𝛼_n, 𝛾_m, 𝛾_n. Finally, the plasma current distribution $j_{{\varphi}}(\bar{{\psi}})$ can be calculated according to (2).

Step-2: Adjust the PF coil currents.

Step-3: Calculate the plasma current profile: Firstly, calculating the displacement of the center of plasma (δrδz) and variation of normalized flux between two iterations according to (3). Secondly, calculating the distribution of plasma current using (4)-(7). After that, coefficient 𝜆 is recalculated to keep the total plasma current constant, also the center of the plasma current is recalculated, and the P’ and FF’ profiles are updated by (8).

Step-4: Calculate the flux produced by PF coil currents and plasma current.

Step-5: Find out the plasma boundary and magnetic axis, calculating the normalized flux.

Step-6: Compute the flux variations on the grids and judge the convergence of equilibrium calculation by checking the equation: $\begin{eqnarray}\displaystyle \left|\frac{({\psi}^{m+1}-{\psi}^{m})}{{\psi}_{\text{b}}-{\psi}_{\text{a}}}\right|<{\varepsilon} & & \displaystyle\end{eqnarray}$ (10) where m is the iteration count, ϵ is the predefined error limit. If the result is convergent, proceed to step-7, otherwise proceed to step-3.

Step-7: Calculate the equilibrium parameters and output the result.

2.2. Validation of the simulation environment

To verify the developed free boundary equilibrium solver, the simulation of EAST pulse 73625 at 1.8 s was carried out, and the result was compared to the experiment (EFIT reconstruction). In particular, starting from the reconstructed plasma equilibrium at 1.2 s by off-line EFIT code, the plasma equilibrium at 1.8 s was predicted with the PF coil currents and total plasma current set as target values, while the PF currents were inputs. A comparison of the key plasma parameters is shown in Table 1 and the plasma shape comparison is shown in Fig. 3. It can be seen that the key parameters of plasma profile calculated by simulation are almost consistent with the experiment, and the biggest error of plasma shape in simulation is 0.4 cm as compared with the experiment, which proves the precision as well as reliability of the simulation tool.

Fig. 2.

Flow chart of the free boundary equilibrium calculation.

Table 1

Comparison of equilibrium parameters

Parameters	Experiment 1.2 s	Experiment 1.8 s	Simulation 1.8 s
Plasma current I _p (kA)	342	380	380
Toroidal field B _t (T)	1.8	1.9	1.9
Poloidal beta 𝛽_p	0.58	0.6	0.56
Internal inductance l _i	1.54	1.61	1.64
Major radius R ₀ (m)	1.88	1.9	1.89
Minor radius a ₀ (m)	0.46	0.45	0.45
Elongation 𝜅	1.31	1.35	1.35
Triangularity δ	0.41	0.46	0.46

Fig. 3.

Comparison of plasma shape between the simulation and experiment. The calculated plasma boundary (red line) is superimposed on the EFIT boundary (black dashed line) at 1.8 s. The black dotted line represents the plasma boundary at 1.2 s.

3. Response of plasma shape to PF coils

EAST plasma control system (PCS) controls the plasma shape by controlling the applied magnetic field generated by PF coil currents. Currently, there are two control algorithms: “RZIP” control, which controls the plasma radial position (R), vertical position (Z) and current (Ip); and “isoflux” control, which controls the magnetic flux at a set of control points. Both two control algorithms adopt the proportional-integral-derivative (PID) controller, and for each controlled variable, optimized control vectors (M_matrix column in EAST PCS) which links the variations of the plasma parameters to the variations of PF coil currents is used to determine the contribution of 12 PF coil currents to the variation control. For other controlled shape variables (elongation and triangularity in this paper), it is possible to develop an integrated control algorithm that simultaneously controls multiple parameters. The M_matrix mentioned above can be decided through calculating the responses of shape parameters to PF coil currents exploiting the simulation tool presented in Section 2.

Taking shot 61782 for example, it is a typical double-null discharge with elongation 1.8, and plasma current up to 400 kA. The plasma equilibrium at 5 s is set as the initial condition. Since the PF coils are symmetrically powered, only the lower coils (PF2, PF4, PF6, PF8, and PF12) are studied. Recalculation of plasma equilibrium is conducted to obtain the plasma response after current disturbances are added to PF2, PF4, PF6, PF8, and PF12. To be specific, each of the PF coil currents is supposed to set to $I_{\text{pfn}}^0 − 200∼\text{A}, I_{\text{pfn}}^0 − 100∼\text{A}, I_{\text{pfn}}^0, I_{\text{pfn}}^0 + 100∼\text{A}, I_{\text{pfn}}^0 + 200∼\text{A}$ in turn, where superscript 0 identifies the initial condition, and subscript n indicates the PF coil number. However, the total number of possible combinations of PF coil currents is large, which brings a great computational burden. Thus, the orthogonal method is applied to reduce the number of simulations.

The orthogonal test method is an efficient way to study pair-wise interactions [11]. It reduces the number of tests by choosing representative variable combinations based on the combinatorial theory. With this method, the parameters that affect the output are defined as “factors”, and the different values of factors are called “levels”. According to the numbers of the factors and levels, the orthogonal array L _m(p ^q) can be constructed. Here, L denotes the orthogonal test method, m is total number of simulations, p and q are numbers of levels and factors, respectively.

To calculate the plasma response, the shape parameters elongation 𝜅 and triangularity δ are set as target functions, the six PF coil currents are considered factors, five levels are selected corresponding to the current values. If we test all the possible combinations of these PF coil currents, a set of 5⁶ = 15625 simulation cases is needed. But instead of testing all the combination of parameters, a L ₂₅(5⁶) orthogonal table as demonstrated in Table 2 is constructed to maximize the test coverage while minimizing the number of test cases to 25. To analyze the relationship between target functions and factors, the average value of each level for six factors and the range of each factor are computed by: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle k_{i}=\frac{1}{5}(J_{i1}+J_{i2}+J_{i3}+J_{i4}+J_{i5})\quad (i=1,2,3,4,5)\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle R=\max (k_{i})-\text{min}(k_{i})\end{eqnarray}$ (12) where i is the level number, J _in (n = 1,2,3,4,5) is the nth simulation result for the factor in i-th level.

The results of 25 simulations are shown in Table 2. According to (11) and (12), the average value of each level and the range of each factor are calculated, Table 2 only shows the case that the target function is elongation for brevity, the full results are shown in Fig. 4. It can be concluded that the plasma elongation decreases monotonically with the increment of the currents in PF2 and PF4 and increases monotonically with the currents in PF6, PF8, PF10 and PF12. The plasma triangularity increases with the increment of the currents in PF2, PF4, PF6 and PF8 and decreases monotonically with the currents in PF10 and PF12. By comparing the range of each factor, R ₄ and R ₃ have the largest value when the target functions are elongation and triangularity, respectively, which indicates that among all the PF coil currents, the disturbances in PF6 and PF8 influence the triangularity and elongation dramatically. Generally speaking, the main role of a PF coil in plasma control is decided by its position and turn number, it is easy to qualitatively analyze the influence of a PF coil on the plasma, while the quantitative assessment can be carried out by linear fitting the experiment data. According to the simulation results, the response matrix that describes the linear relationships between elongation, triangularity and PF coil currents is calculated, the elements are values of PF coil currents (unit: A) per unit of elongation or triangularity, as shown in Table 3.

Fig. 4.

The dependence of elongation and triangularity on PF coil currents.

Table 2

Orthogonal Array with 6 Factors and 5 levels

Column	1	2	3	4	5	6	Elongation	Triangularity
Factor	I _pf2	I _pf4	I _pf6	I _pf8	I _pf10	I _pf12	𝜅	δ
1	−200	−200	−200	−200	−200	−200	1.758	0.347
2	−200	−100	−100	−100	−100	−100	1.789	0.369
3	−200	0	0	0	0	0	1.827	0.518
4	−200	+100	+100	+100	+100	+100	1.765	0.516
5	−200	+200	+200	+200	+200	+200	1.738	0.510
6	−100	−200	−100	0	+100	+200	1.815	0.509
7	−100	−100	0	+100	+200	−200	1.778	0.515
8	−100	0	+100	+200	−200	−100	1.771	0.516
9	−100	+100	+400	−200	−100	0	1.814	0.386
10	−100	+200	−200	−100	0	+100	1.843	0.407
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
25	+200	+200	+100	0	−100	−200	1.778	0.544
k ₁	1.801	1.801	1.786	1.760	1.788	1.795
k ₂	1.794	1.798	1.790	1.781	1.790	1.796
k ₃	1.790	1.796	1.796	1.792	1.793	1.798
k ₄	1.786	1.795	1.800	1.805	1.796	1.799
k ₅	1.776	1.793	1.811	1.820	1.800	1.802
R	0.025	0.008	0.025	0.06	0.012	0.007

Table 3

Response matrix of elongation and triangularity to PF coil currents

A	PF2	PF4	PF6	PF8	PF10	PF12
𝜅	16000	50000	16000	6667	33333	57142
δ	8000	22222	5714	8163	13333	15384

4. Conclusion

To study the impact of PF current disturbance on tokamak plasma shape, a free boundary equilibrium solver is presented in this paper, it is an efficient and effective simulation tool, and to reduce the number of simulations, the orthogonal method is used to design the tests. The response curves of elongation and triangularity in relation to the PF coil currents are obtained, as for other plasma parameters, the proposed methodology can also be applied to calculate their responses to PF coil currents. In the future work, the presented calculation will be used to design the feedback gain matrix in the EAST shape controller.

Footnotes

Acknowledgements

The author is grateful to Z.P. Luo for the support to the simulation. This work is supported in part by the National Natural Science Funds of China under Grant 51677051 and in part by Institute of Plasma Physics Chinese Academy of Sciences.

References

Albanese

Bellesia

Palmaccio

Portone

Mateos

F.R.

Rubinacci

, Electromagnetic transient studies in the ITER tokamak, International Journal of Applied Electromagnetics and Mechanics 39(1-4) (2012), 65–71.

Yuan

Q.P.

Xiao

B.J.

Luo

Z.P.

Walker

M.L.

Welander

A.S.

Hyatt

Qian

J.P.

Zhang

R.R.

Humphreys

D.A.

Leuer

J.A.

Johnson

R.D.

Penaflor

B.G.

and Mueller

, Plasma current, position and shape feedback control on EAST, Nuclear Fusion 53(4) (2013), 043009.

Elahi

A.S.

and Ghoranneviss

, Magnetic-based measurements of tokamak plasma equilibrium parameters, IEEE Transactions on Plasma Science 41(2) (2013), 334–340.

Ambrosino

and Albanese

, Magnetic control of plasma current, position, and shape in tokamaks: A survey of modeling and control approaches, Control Systems IEEE 25(5) (2005), 76–92.

Albanese

Artaserse

Mattei

Tudisco

and JET - EFDA Contributors , Plasma reconstruction in tokamaks with linearized approaches, International Journal of Applied Electromagnetics and Mechanics 26(3,4) (2007), 191–199.

J.J.

LWang

J.G.

Qian

and the EAST team , Effect of plasma disruption on superconducting magnet in EAST, Fusion Engineering and Design 88(9,10) (2013), 1513–1517.

Albanese

Ambrosino

Calabr’O

Castaldo

Crisanti

DeTommasi

Liu

Luo

Z.P.

Mele

Pironti

Xiao

B.J.

and Yuan

Q.P.

, A MIMO architecture for integrated control of plasma shape and flux expansion for the EAST tokamak, IEEE Conference on Control Applications (2016), 611–616.

Chen

S.L.

Villone

Sun

Y.W.

and Xiao

B.J.

, Comparison of different linearized plasma response models on the EAST tokamak, Nuclear. Fusion 58(7) (2018), 076009.

Qiu

Q.L.

Xiao

B.J.

Yong

Liu

Xing

and Humphreys

D.A.

, Simulation of EAST vertical displacement events by tokamak simulation code, Nuclear. Fusion 56(10) (2016), 106029.

10.

Qian

J.P.

Wan

B.N.

Lao

L.L.

Shen

Sabbagh

S.A.

Sun

Y.W.

D.M.

Xiao

B.J.

Ren

Q.L.

Gong

X.Z.

and Li

J.G.

, Equilibrium reconstruction in EAST Tokamak, Plasma Science and Technology 11(2) (2009), 142–145.

11.

Banerji

, Orthogonal array approach for test case optimization, International Journal of Advanced Research in Computer and Communication Engineering 1(9) (2012), 613–621.