
Citation: | Ricardo FARENGO, Pablo L GARCÍA-MARTÍNEZ, Hugo E FERRARI. A simple and fast method to calculate the coil currents and the external poloidal flux and magnetic field for fixed boundary equilibria[J]. Plasma Science and Technology, 2024, 26(11): 115103. DOI: 10.1088/2058-6272/ad70e0 |
A simple and fast method to calculate the coil currents and the external poloidal flux and magnetic field consistent with a given internal magnetohydrodynamic (MHD) equilibrium and coil geometry is presented. The constrained optimization technique employed allows to fix the exact position of the X-points. The external poloidal flux and field are calculated using Green’s functions expressed in terms of elliptic integrals that are calculated using a fast and accurate routine. The calculation can be limited to the region between the plasma and the device wall. An example for a double null configuration is presented.
Alternating current (AC) operation of a tokamak reactor is an attractive way to generate a continuous output of electric energy without the need for a complicated non-inductive current driving system. In recent decades, AC operation experiments have been performed in a number of tokamak devices [1–11]. This approach has been demonstrated to be feasible and free from secular degradation in performance. AC tokamak operation is always associated with a magnetic equilibrium in which the current density reverses. It is well known that the first step in tokamak modeling and analysis is the reconstruction of MHD equilibrium and stability calculations. Some active theoretical research activities on the current reversal equilibrium configurations (CRECs) have been done analytically [12–14] and numerically [15–17]. However, previous research mainly focused on the existence of the CRECs [12, 13, 15] and the fixed-boundary CRECs [12, 14, 15, 17], which are important for the understanding of the physics of CRECs but invalid for practical application to experiments on the CRECs reconstruction calculations to determine the magnetic separatrix and X-point locations. In practical experimental equilibria, the position of the plasma boundary is determined by the interaction between the currents in plasmas and in the plasma equilibrium control coils. For the current reversal equilibrium configuration, the plasma magnetic surface is not nested, and the magnetic surface function is not suitable for normalization by the poloidal magnetic flux at the magnetic axis and plasma boundary. It is highly difficult to adjust the plasma current profile using the iterative method adopted by Johnson et al [18] when the macroscopic plasma parameters are known. Therefore, we must use the finite element method of unstructured mesh and transform the free boundary problem into a series of fixed boundary problems to deal with, that is, the Hagenow and Lackner method [19–21].
In this paper, we propose an improved Hagenow and Lackner method through which the current reversal equilibrium configuration and the external control coil current can be obtained from the prescribed plasma parameters in the tokamak configuration with plasma surface bounded by a given desired plasma shape or the X-point outside of the torus. In section 2 we describe the plasma equilibrium model and the method of solution. In section 3, the benchmark and some results related to the actual plasma shape are presented for equilibria with current reversal. The summary and conclusions are described in section 4.
The ideal MHD equilibrium in an unbounded domain can be described by a second-order partial differential equation, i.e. the so-called Grad-Shafranov equation. In axisymmetric geometry, it can be written in cylindrical coordinates (x,ϕ,z) as follows [18, 21]:
Δ∗ψ=xjϕ(x, ψ), | (1) |
with
jϕ(x,ψ)={xdβ(ψ)dψ+f(ψ)xdf(ψ)dψintheplasmaregionwithboundaryΓp∑Nci=1Iiδ(x−xci)δ(z−zci)invacuum, | (2) |
Δ∗=−x∂∂x1x∂∂x−∂2∂z2, | (3) |
where ψ=Ψ/B0r2a is the normalized poloidal magnetic flux; x=R/ra, z=Z/ra, β(ψ)=2μ0p(ψ)/B20, f(ψ)=F(ψ)/B0ra, with the magnetic field represented by {\boldsymbol{B}}= F\left(\Psi\right)\nabla\phi+ \nabla\Psi\times\nabla\phi ; \phi is an ignorable angle in the cylindrical coordinate system \left(R,\ \phi,\ Z\right); p\left(\psi\right) is the plasma pressure; B_{0} is the vacuum magnetic field evaluated at R=R_{0}; R_{0} and r_{a} are the major radius and the minor radius, respectively; j_{\phi}= \mu_{0}aJ_{\phi}/B_{0} , where J_{\phi} is the toroidal current density and \mu_{0} is the permeability of free space; and I_{i} and \left(x_{i}^{{\rm c}},z_{i}^{{\rm c}}\right) are the dimensionless currents and coordinates of the external coil respectively.
When the plasma is surrounded by a vacuum region, j_{\phi} is nonzero only in the interior of the plasma. We assume that the plasma is in a domain of the cross-section denoted by \Omega_{{\rm p}} , with boundary \Gamma_{{\rm p}} , and the outside of the plasma is a vacuum region denoted by \Omega_{{\rm v}} with a rectangular computational boundary \Gamma_{{\rm v}} . It is worth noting that for the free boundary current reversal equilibrium, the poloidal magnetic flux within the plasma region is non-nested and not monotonic, so the boundary of the plasma and vacuum cannot be determined by comparing the value of the poloidal magnetic flux with a given \psi_{\Gamma_{{\rm p}}} as was done in reference [18], which is widely used by EFIT and other codes etc. In addition, when we want to effectively optimize the plasma current profile to improve the convergence efficiency of the solver, the poloidal magnetic flux at this time cannot be normalized simply as done in reference [18], that is to say, their method does not work efficiently in this case. However, we extend the method of Hagenow and Lackner, and find that it may be an efficient method for the solution of the CREC of equation (1) in an unbounded domain. The solution process can be grouped in the following steps:
Step 1. Determine j_{\phi}\left(x,\ \psi\right) by solving equation (1) in the plasma region imposed the Dirichlet conditions along a given closed curve \Gamma_{{\rm p}} , with \psi on \Gamma_{{\rm p}} given as a constant \psi_{\Gamma_{{\rm p}}} .
Equation (1) is a nonlinear equation, and both \beta and f are general functions of \psi . Without the loss of generality, considering \psi_{\Gamma_{{\rm p}}}=0 , we take the profiles of \beta and f as the following polynomial of \psi :
\begin{align} \beta & =\beta_{{\rm min}}+a_{1}\psi+a_{2}\psi^{2}-\beta_{{\rm m}}, \end{align} | (4) |
\begin{align} f & =f_{{\rm e}}+b_{1}\psi+b_{2}\psi^{2}, \end{align} | (5) |
where f_{{\rm e}}=f\big|_{\Gamma_{{\rm p}}}\approx x_{0}=R_{0}/r_{a} . The value \beta_{{\rm min}} is the minimum of \beta , which can be found from the actual experimental measurements. However, for general theoretical numerical analysis, we can arbitrarily choose a suitable \beta_{{\rm min}} value. The value of \beta_{{\rm m}} is determined by the minimum value of \beta in the plasma domain in which \beta_{{\rm min}} is found. The coefficients a_i,\ b_i\left(i=1,\ 2\right) are not realistic experimental parameters, but they can be determined by a group of actual physical quantities. In this work, we consider that two of them are given arbitrarily and the other two are determined by the total toroidal current I_{{\rm p}} and the volume averaged \beta_{{\rm v}} . Here,
\begin{align} I_{{\rm p}} & =\frac{B_{0}a}{\mu_{0}}\int_{\Omega_{{\rm p}}}j_{\phi}{\rm d}S, \end{align} | (6) |
\begin{align} \beta_{{\rm v}} & =\frac{2\mu_{0}\langle p\rangle }{\langle B^{2}\rangle } \end{align} | (7) |
where \langle \cdots\rangle =\dfrac{\int_{\Omega_{{\rm p}}}x\left(\cdots\right){\rm d}S}{\int_{\Omega_{{\rm p}}}x{\rm d}S} , \Omega_{{\rm p}} is the plasma domain, and {\rm d}S is the plasma cross area element. In order to obtain the desired CRECs by adjusting a_{i} and b_{i} , we can estimate their initial values through the eigenvalues of the Grad-Shafranov operator. Using the eigenvalue problem [15]
{\text{Δ}}^{*}\psi =\lambda\psi, | (8) |
with the same boundary condition as in equation (1) with a toroidal current density j_{\phi}=\lambda\psi/x , and comparing with equation (2) we can obtain 2x^{2}a_{2}+b_{1}^{2}+2f_{{\rm e}}b_{2}\approx\lambda . An equilibrium solution with the nested flux surfaces \psi={\rm const} corresponds to the lowest eigenvalue and an eigenfunction without nulls inside the plasma domain. A lot of other eigenfunctions provide a series of equilibria with a different topology of the magnetic surfaces.
Step 2. Construct the Green’s function of the Grad-Shafranov operator to solve equation (1) in an unbounded domain.
The known Green’s function in a cylindrical geometry is given by the following [22]:
G\left(x, z;\; x^{'}, z^{'}\right) =-\frac{1}{2\pi}\frac{\sqrt{xx^{'}}}{k}\bigg[\left(\left(2-k^{2}\right)K\left(k\right)-2E\left(k\right)\right)\bigg], | (9) |
\begin{align} k^{2} & =\frac{4xx^{'}}{\left(x+x'\right)^{2}+\left(z-z'\right)^{2}}, \end{align} | (10) |
where G satisfies
\text{Δ}^*G\left(x,\ z;\ x^{'},\ z^{'}\right)=x{\delta}\left(x-x^{'}\right){\delta}\left(z-z^{'}\right), | (11) |
where K\left(k\right) and E\left(k\right) are complete elliptic integrals of the first and second kind, and \left(x^{'},z^{'}\right) is the observation point. Using Green's theorem gives
\nabla\cdot\left[\psi\frac{1}{x^{2}}\nabla G-G\frac{1}{x^{2}}\nabla\psi\right]=\frac{1}{x^{2}}\psi{\text{Δ}}^{*}G-\frac{1}{x^{2}}G{\text{Δ}}^{*}\psi. | (12) |
We integrate the above equation over all space and use Gauss’s theorem to convert the divergences to surface integrals which vanish at infinity, then the expression to obtain \psi outside or on the computational boundary is
\psi\left(x^{'},\; z^{'}\right)=\int_{\Omega_{\rm{p}}}^{ }G\left(x,\; z;\; x^{'},\; z^{'}\right)j_{\phi}{\rm{d}}S+\sum_{i=1}^{N_{\rm{c}}}G\left(x_i^{\mathrm{c}},\; z_i^{\rm{c}};\; x^{'},\; z^{'}\right)I_i. | (13) |
From the above equation (13), it is clearly found that \psi is composed of the magnetic fluxes induced by the plasma current (the first term on the right hand side (RHS)) and the external coil currents (the second term on the RHS).
Step 3. Compute the contribution from the plasma currents using the von Hagenow method.
Consider a function \psi_{\rm{g}}\left(x,\ z\right) which satisfies the same differential equation as equation (1) that \psi does in the interior of the rectangular computational boundary \Gamma_{{\rm v}} , but which vanishes on the boundary, i.e.:
{\text{Δ}} ^{*}\psi_{{\rm g}} ={\text{Δ}}^{*}\psi=xj_{\phi}\left(x,\psi\right), | (14) |
\psi_{{\rm g}}\big|_{\Gamma_{{\rm v}}}=0. | (15) |
The resulting flux-function distribution {{\psi_{{\rm g}}}} can be interpreted as the sum of the field induced by the plasma currents plus the field of the fictitious mirror currents in an unbound domain. In order to solve equation (14) efficiently, we can use interpolation to transform j\left(x,\ \psi\right) into j\left(x,\ z\right) , and equation (14) becomes a linear partial differential equation with very good convergence performance. In the same way as in step 2, we use the following form of Green's theorem:
\nabla\cdot\left[\psi_{{\rm g}}\frac{1}{x^{2}}\nabla G\right]-\nabla\cdot\left[G\frac{1}{x^{2}}\nabla\psi_{{\rm g}}\right] =\frac{1}{x^{2}}\psi_{{\rm g}}{\text{Δ}}^{*}G-\frac{1}{x^{2}}G{\text{Δ}}^{*}\psi_{{\rm g}}. | (16) |
Integrating the above equation over the computational domain \Omega_{{\rm c}}=\Omega_{{\rm v}}+\Omega_{{\rm p}} , using \psi_{{\rm g}}|_{\Gamma_{{\rm v}}}=0 , and evaluating the differential surface element related to the differential line element by {\rm d}S=2\pi x{\rm d}l , then we can obtain the identity as follows:
\begin{align} \int_{\Omega_{{\rm p}}}G\left(x,z;x^{'},z^{'}\right)j_{\phi}{\rm d}S & =\psi_{{\rm g}}\left(x^{'},z^{'}\right)+\int_{\rm b}G\left(x,z;x^{'},z^{'}\right)\frac{1}{x}\frac{\partial\psi_{{\rm g}}}{\partial n}{\rm d}l. \end{align} | (17) |
The second term on the RHS of the above equation represents the fictitious mirror currents in an unbound domain. Substituting equation (17) into equation (13) yields [20]
\begin{split} \psi\left(x^{'},\ z^{'}\right)= & \psi_{\rm{g}}\left(x^{'},\ z^{'}\right)+\int_{\rm{b}}^{ }\frac{{\rm{d}}l}{x}G\left(x,\; z;\; x^{'},\; z^{'}\right)\frac{\partial\psi_{\rm{g}}}{\partial n} \\ & +\sum_{i=1}^{N_{\rm{c}}}G\left(x_i^{\rm{c}},\; z_i^{\rm{c}};\; x^{'},\; z^{'}\right)I_i, \end{split} | (18) |
which is the solution of equation (1) at a point interior to the computational boundary.
Step 4. Determine the external coil currents to produce a given plasma shape.
To determine the necessary external currents, we solve the overdetermined problem of finding the least square error [21, 23-25] in the poloidal flux at the N_{k} plasma boundary points \left(x_k^{'},\ z_k^{'}\right) produced by N_{{\rm c}} coils at set locations \left(x_{i}^{{\rm c}},\;z_{i}^{{\rm c}}\right) with N_{k}>N_{{\rm c}} . Assuming that the plasma boundary points lie interior to the computational boundary, we use equation (18) to express the function to be minimized as
\begin{split} \epsilon =&\sum_{k=1}^{N_{k}}\left(\psi_{k}-\psi_{{\rm b}}\right)^{2}+\gamma_{1}\sum_{i=1}^{N_{{\rm c}}}I_{i}^{2}+\gamma_{2}\sum_{i=2}^{N_{{\rm c}}}\left(I_{i}-I_{i-1}\right)^{2} \\ =&\sum_{k=1}^{N_{k}}\Bigg(\psi_{{\rm g}k}+{{\int_{\rm b}}{\frac{{\mathrm{d}l}}{R}}}G\left(x,\;z;\;x_{k}^{'},\;z_{k}^{'}\right)\frac{\partial\psi_{{\rm g}}}{\partial n}+\sum_{i=1}^{N_{{\rm c}}}G\left(x_{i}^{{\rm c}},\;z_{i}^{{\rm c}};\;x_{k}^{'},\;z_{k}^{'}\right)I_{i}\\ &-\psi_{{\rm b}}\Bigg)^{2}+\gamma_{1}\sum_{i=1}^{N_{{\rm c}}}I_{i}^{2}+\gamma_{2}\sum_{i=2}^{N_{{\rm c}}}\left(I_{i}-I_{i-1}\right)^{2} \\ = &\sum_{k=1}^{N_{k}}\left(\sum_{i=1}^{N_{{\rm c}}}c_{ik}I_{i}-d_{k}\right)^{2}+\gamma_{1}\sum_{i=1}^{N_{{\rm c}}}I_{i}^{2}+\gamma_{2}\sum_{i=2}^{N_{{\rm c}}}\left(I_{i}-I_{i-1}\right)^{2}, \end{split} | (19) |
where
\begin{align} c_{ik} & =G\left(x_{i}^{{\rm c}},\;z_{i}^{{\rm c}},\;x_{k}^{'},\;z_{k}^{'}\right), \end{align} | (20) |
\begin{align} d_{k} & ={{\psi_{{\rm b}}-\int_{\rm b}}{\frac{{\rm d}l}{R}}}G\left(x,\;z;\;x_{k}^{'},\;z_{k}^{'}\right)\frac{\partial\psi_{{\rm g}}}{\partial n}-\psi_{{\rm g}}{}_{k}. \end{align} | (21) |
Here \gamma_{1} and \gamma_{2} are the regularization parameters which are introduced to stabilize the procedure against large unreal oscillations of the current I_{i} and the dramatic current change of adjacent coils {\text{Δ}} I_{i} ; \psi_{{\rm g}k}=\psi_{{\rm g}}\big|_{\Gamma_{{\rm p}}} ; and \psi_{{\rm b}} is the desired value of flux at the given boundary points \left(x_{k}^{'},\;z_{k}^{'}\right) . We minimize \epsilon by \partial\epsilon/\partial I_{i}=0 , and obtain
\begin{split} &\sum_{k=1}^{N_{k}}\left(\sum_{i=1}^{N_{{\rm c}}}c_{ik}I_{i}-d_{k}\right)c_{jk}+\gamma_{1}I_{j}\\ & + \begin{cases} -\gamma_{2}\left(I_{2}-I_{1}\right), & j=1\\ -\gamma_{2}\left(I_{N_{{\rm c}}}-I_{N_{{\rm c}}-1}\right), & j=N_{{\rm c}}\\ \gamma_{2}\left(2I_{j}-I_{j-1}-I_{j+1}\right) & N_{{\rm c}}>j>1 \end{cases}=0. \end{split} | (22) |
The above algebraic equation can be rewritten in a concise form as follows:
\begin{align} \sum_{i=1}^{N_{{\rm c}}}A_{ji}I_{i} & =B_{j}, \end{align} | (23) |
where
\begin{align} A_{ji} & =\begin{cases} \gamma_{1}+\gamma_{2}+\displaystyle\sum\nolimits_{k=1}^{N_{k}}c_{jk}^{2}, & i=j,\;j=1,\;N_{{\rm c}}\\ -\gamma_{2}+\displaystyle\sum\nolimits_{k=1}^{N_{k}}c_{ik}c_{jk}, & i=2,\;j =1\;{\rm and}\;i= N_{{\rm c}}-1,\;j=N_{{\rm c}}\\ \gamma_{1}+2\gamma_{2}+\displaystyle\sum\nolimits_{k=1}^{N_{k}}c_{jk}^{2}, & i=j,\;N_{{\rm c}}>j>1\\ -\gamma_{2}+\displaystyle\sum\nolimits_{k=1}^{N_{k}}c_{ik}c_{jk}, & i=j\pm1,\;N_{{\rm c}}>j>1\\ \displaystyle\sum\nolimits_{k= 1}^{N_{k}}c_{ik}c_{jk}, & i\neq j,\;j\pm1,N_{{\rm c}}\geqslant j\geqslant1 \end{cases} \end{align} | (24) |
\begin{align} B_{j} & =\sum_{k=1}^{N_{k}}d_{k}c_{jk}. \end{align} | (25) |
Finally, we summarize the equilibrium algorithm in figure 1. Initially, we guess a profile for the plasma current. We use this to give a first guess for the desired I_{\rm{p}0},\ \beta_{\rm{v}0} within the desired plasma boundary shape using equations (6) and (4). With these constraints held fixed, we solve equation (1) with fixed plasma boundary. When this converges, we use the current distribution so computed to obtain the poloidal magnetic flux \psi by solving equation (14) with the rectangular computational boundary and its corresponding adjoint Green's function. Secondly, we use the magnetic flux obtained above and the desired boundary flux value, by adjusting parameters \gamma_{1} and \gamma_{2} , hold the boundary flux fixed until convergence, and finally determine the coil current in the vacuum outside the plasma.
According to the above detailed solution process, we developed an Advanced Tokamak Equilibrium Code (ATEC), which is written in MATLAB language. The Grad-Shafranov (G-S) operator solver is from the PDE solver in the MATLAB toolkit. The computational domain is discretized by a triangular mesh finite element method. We directly call the function “initmesh” in the MATLAB mathematical function library to generate triangular mesh, and call the nonlinear PDE solver function “pdenonlin” to solve the nonlinear G-S equation. This code was benchmarked against the well-known equilibrium code EFIT [26] and TEQ [27] for the normal nested plasma equilibrium configuration. An example of the comparison between the equilibrium results from EFIT, TEQ and ATEC is shown in figure 2, where the plasma parameters, profiles and related equilibrium information all come from the EAST equilibrium data file g070754.003740 calculated by EFIT and the CFETR equilibrium data file g220905.7500_teq calculated by Tequation. The position, size, and turns of the poloidal coils come from table 1 in references [28] and [29]. The red lines in figure 2 are the results from ATEC using coil currents given by EFIT (a), ATEC (b) itself for EAST, and TEQ (c), ATEC (d) itself for CFETR. The blue dashed lines in figure 2 are the results from EFIT for EAST ((a) and (b)) and TEQ for CFETR ((c) and (d)). It can be found that the results of ATEC and the other two codes are in good agreement. Tables 1 and 2 show the comparison of the coil currents, the total coil currents, and the errors of the last closed flux surface (LCFS), X-point and magnetic axis calculated by the codes ATEC, EFIT, and TEQ for EAST and CFETR, respectively.
EAST | CFETR | |||||
Coils | I_{i} (EFIT) | I_{i} (ATEC) | Coils | I_{i} (TEQ) | I_{i} (ATEC) | |
PF1 | 5.4370e+05 | 5.0517e+05 | CS1U | −6.0194e+06 | −5.4287e+06 | |
PF2 | 5.4270e+05 | 5.0308e+05 | CS2U | 2.7232e+06 | 1.0234e+06 | |
PF3 | −1.4077e+05 | 7.3277e+04 | CS3U | 8.3238e+06 | 1.3514e+07 | |
PF4 | 1.7036e+05 | 1.9780e+05 | CS4U | 8.2293e+06 | −2.0077e+06 | |
PF5 | 3.8309e+05 | −1.2694e+05 | CS1L | −1.0959e+07 | −1.1195e+07 | |
PF6 | 5.7632e+04 | −2.0271e+04 | PF1U | 1.2451e+07 | 1.4823e+07 | |
PF7 | 3.5643e+05 | 1.0275e+06 | PF2U | −2.6124e+06 | −3.2489e+06 | |
PF8 | 3.6069e+05 | 8.9066e+05 | PF3U | −6.9199e+06 | −6.9212e+06 | |
PF9 | 1.6525e+06 | 1.1513e+06 | PF3L | −3.1210e+06 | −1.5712e+06 | |
PF10 | 1.6723e+06 | 1.0563e+06 | CS2L | 2.6956e+05 | 5.1617e+05 | |
PF11 | −4.9521e+05 | −4.9535e+05 | CS3L | 1.6669e+07 | 1.8536e+07 | |
PF12 | −4.8425e+05 | −4.3512e+05 | CS4L | 1.4653e+07 | 1.1879e+07 | |
PF13 | 4.3716e+04 | 2.9533e+04 | PF2L | −7.8468e+06 | −1.0188e+07 | |
PF14 | 4.4703e+04 | 5.4555e+04 | PF1L | 1.7482e+07 | 1.2930e+07 | |
DC1 | 4.9374e+06 | 8.9734e+06 |
EAST | CFETR | |||||
Coils | I_{i} (EFIT) | I_{i} (ATEC) | Coils | I_{i} (TEQ) | I_{i} (ATEC) | |
\sqrt{\sum\displaystyle I_{i}^{2}} | 2.6566e+06 | 2.3018e+06 | \sqrt{\sum\displaystyle I_{i}^{2}} | 3.7476e+07 | 3.7681e+07 | |
\left|\dfrac{\psi_{{\rm e}}-\psi_{{\rm e}}^{0}}{\psi_{{\rm e}}^{0}}\right| | 1.8720e−02 | 6.5512e−05 | \left|\dfrac{\psi_{{\rm e}}-\psi_{{\rm e}}^{0}}{\psi_{{\rm e}}^{0}}\right| | 8.6541e−03 | 2.7569e−04 | |
\left|\dfrac{\psi_{{\rm axis}}-\psi_{{\rm axis}}^{0}}{\psi_{{\rm axis}}^{0}}\right| | 1.2263e−02 | 8.3719e−06 | \left|\dfrac{\psi_{{\rm axis}}-\psi_{{\rm axis}}^{0}}{\psi_{{\rm axis}}^{0}}\right| | 3.2759e−03 | 2.0890e−04 | |
\dfrac{\left|X_{{\rm P}}^{0}-X_{{\rm P}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 6.6993e−03 | 6.6993e−03 | \dfrac{\left|X_{{\rm P}}^{0}-X_{{\rm P}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 6.2120e−03 | 6.2120e−03 | |
\dfrac{\left|X_{{\rm axis}}^{0}-X_{{\rm axis}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 1.9714e−03 | 1.9714e−03 | \dfrac{\left|X_{{\rm axis}}^{0}-X_{{\rm axis}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 5.7144e−03 | 5.7144e−03 |
In the following section, the second type of free boundary case is considered, and the plasma parameters for simulation are from the EAST data of shot # g070754.003740. The desired LCFS is also the shape reconstructed by magnetic probe diagnosis using EFIT. The main plasma parameters are as follows: R_{0}=1.86\;{\rm m}, r_{a}=0.43, B_{0}=1.80\;{\rm T} and \beta_{{\rm v}}=0.02 , where R_0,\ r_a,\ B_0 are the major radius, minor radius, and the vacuum magnetic field at the geometric center of the plasma, respectively.
Figures 3(a) and (b) are the equilibrium configurations when the plasma current is zero. Figure 3(a) is the high-low field current reversal equilibrium (HL-CREC), whose profile coefficients of plasma currents are a_{1}=2.467 e−01, a_{2}= 1.8000 e−02, b_{1}=−9.7532e−02 and b_{2}=3.1000{\rm e} −01 in equations (4) and (5). This equilibrium configuration has been observed experimentally [10]. As shown in figure 3(a), the plasma is separated into two magnetic islands balanced with opposite current directions in the left and right, which is somewhat similar to a single spherical mark with zero aspect ratio. The magnetic island triangle is severely distorted, resulting in the existence of three X-points, the top two X-points and the bottom one. In addition, the total plasma current is zero, namely the total loop voltage is zero, but in fact, the plasma current density is not locally zero. The high and low field plasma currents are opposite, i.e. there should be an electric potential difference inside the plasma, resulting in a current across the plasma cross section. Figure 3(b) shows the up-down current reversal equilibrium configuration (UD-CREC), and the corresponding profile coefficients of plasma current are a_{1}=9.1441{\rm e} −02, a_{2}=1.5000{\rm e} −02, b_{1}=1.3201{\rm e} −03 and b_{2}=2.2000{\rm e} −01. This UD-CREC and the above-mentioned HL-CREC are both the eigensolutions of the G-S equation. However, the UD-CREC has not been observed experimentally so far, so the existence and stability analysis of this equilibrium need to be further verified by theory and experiment. As shown in figure 3(b), just like that in HL-CREC, the magnetic island triangle is also severely distorted, resulting in the existence of three X-points. For the specific actual device, more deflector plates must be considered. Compared with the HL-CREC, the UD-CREC seems not to be a good choice for tokamak plasmas. Figure 3(c) is the center current reversal equilibrium configuration (C-CREC), whose profile coefficients of plasma current are a_{1}= −1.1207e−02, a_{2}=4.9015{\rm e} −02, b_{1}= −1.8829e−02 and b_{2}=1.0000{\rm e} −01. This configuration is not the eigensolution of the G-S equation. Here, we just adopt the discontinuous current density profile model, which is similar to reference [15] and artificially set the current density in the central region to be negative. It can be seen that there is only one X-point at the plasma boundary, just like the nested equilibrium configuration (NEC), but the difference is that there are multiple magnetic islands in the interior of the plasma. Although such a solution exists mathematically, no such configuration has been found in experiments so far. However, the zero-center current is found more frequently and is able to exist stably, and this is called the current-hole equilibrium configuration [30−32]. The possible reason is that this type of center reversal current equilibrium is very unstable and is very easily relaxed to the current-hole equilibrium [30]. Figure 3(d) is the normal or nested equilibrium configuration (NEC), whose profile coefficients of plasma currents are a_{1}=−4.9407e−02, a_{2}=−1.0345e−01, b_{1}= 2.0702e−02 and b_{2}=1.2000{\rm e} −01. Figure 4 shows that the profiles of the normalized poloidal magnetic flux in HL-CREC and UD-CREC are reversed with the reversal of the direction of current density, but the pressure profiles are still similar to that in the NEC, and are hollowed in the C-CREC. However, the toroidal magnetic fields in these current reversal configurations are all little changed. Table 3 shows that the total effective coil currents \sqrt{\sum I_{i}^{2}} of the HL-CREC and UD-CREC are a little bigger than those of the nested equilibrium under the same errors of the LCFS. This means that controlling the HL-, UD-, and C-CRECs is a little more costly than controlling the NEC. It is also shown in table 3 that the \beta_{{\rm p}} of the CRECs is a little smaller than that of the NEC. In addition, as for AC operation, one of the most critical problems is the maintenance of the plasma when the plasma current is zero. In such a case, people may mainly focus more on the existence and stability of the equilibrium configuration when the plasma current is zero, so if we do not care much about the exact shape of the LCFS and the performance of plasma parameters, we can adjust \gamma_{1} and \gamma_{2} to reduce by half or more the total effective control coil currents compared to those for the exact fits of the LCFS to maintain the CRECs with a rough desired shape (as shown in figure 5 and table 4). In order to better understand the equilibrium conditions and the efficiency of the coil current, we also show the absolute plasma current I_{{\rm pabs}}=\dfrac{B_{0}a}{\mu_{0}}\int_{\Omega}|j_{\phi}|{\rm d}S in table 4.
Coils | HL-CREC (A) | UD-CREC (A) | C-CREC (A) | NEC (A) |
PF1 | 1.4131e+06 | 8.7750e+05 | 1.4433e+06 | 1.3968e+06 |
PF2 | 2.5396e+06 | −1.6512e+05 | −2.6781e+05 | 1.0406e+06 |
PF3 | 5.3719e+06 | 3.1279e+06 | −1.7004e+06 | 1.0874e+06 |
PF4 | 1.7050e+06 | −2.3708e+05 | −1.2162e+06 | 1.2288e+06 |
PF5 | −1.5675e+07 | 4.8045e+06 | −1.1617e+06 | 1.8051e+06 |
PF6 | −5.3236e+06 | −8.2558e+05 | −6.3073e+04 | 2.0526e+06 |
PF7 | 3.0715e+07 | −3.8904e+05 | 4.6209e+06 | 5.9412e+05 |
PF8 | 1.5813e+07 | 1.1429e+06 | 2.0429e+06 | 2.2780e+05 |
PF9 | −2.8630e+07 | −4.5500e+06 | −1.1784e+06 | −1.6729e+06 |
PF10 | −1.5845e+07 | 1.9476e+06 | 9.8322e+05 | −1.4072e+06 |
PF11 | 2.2270e+06 | 2.9918e+06 | −1.2139e+06 | 3.7252e+05 |
PF12 | 2.1009e+06 | −2.9352e+05 | −8.7355e+05 | 5.4823e+05 |
PF13 | −1.3933e+06 | −9.3773e+05 | 3.8834e+05 | 2.2159e+05 |
PF14 | −1.3285e+06 | −7.8612e+05 | −1.8577e+05 | −1.6968e+04 |
\sqrt{\sum I_{i}^{2}} | 5.0907e+07 | 8.4196e+06 | 6.1792e+06 | 4.3440e+06 |
\beta_{{\rm p}} | 1.2630e+00 | 1.0259e+00 | 1.1347e+00 | 1.7049e+00 |
Coils | HL-CREC (A) | UD-CREC (A) |
PF1 | 2.2298e+06 | 1.9104e+06 |
PF2 | 2.2664e+06 | −6.6681e+05 |
PF3 | 5.8628e+05 | 2.9801e+06 |
PF4 | 7.3105e+05 | −2.3493e+05 |
PF5 | 3.7842e+04 | 9.5368e+05 |
PF6 | 1.0214e+05 | 3.0989e+05 |
PF7 | 1.2618e+05 | −7.8916e+05 |
PF8 | 9.6939e+04 | 9.7855e+05 |
PF9 | 1.9136e+05 | −8.3600e+05 |
PF10 | 1.6128e+05 | 8.8734e+05 |
PF11 | 1.1100e+06 | 1.6053e+06 |
PF12 | 1.2024e+06 | 4.7541e+05 |
PF13 | −1.0344e+06 | −6.1468e+04 |
PF14 | −9.9634e+05 | −1.3812e+05 |
\sqrt{\sum I_{i}^{2}} | 3.9784e+06 | 4.6709e+06 |
I_{{\rm pabs}} | 1.1572e+06 | 1.02595e+06 |
In summary, an advanced tokamak equilibrium calculation code ATEC is developed by using the improved von Hagenow method. This method can transform the free boundary problem into a series of fixed boundary problems, which greatly improves the efficiency of solving G-S equations under more general current profiles. This code can be used for numerical simulation of more general free boundary problems of tokamak plasmas including the current-hole and current-reversal equilibrium. The most critical plasma current zero-crossing equilibrium problem in the AC operation process is investigated using this code. The results show that the configuration control of zero plasma current equilibrium can be achieved under the condition in which the coil control current is of the same order as that of the normal plasma equilibrium. At the same time, the current reversal equilibrium configuration can maintain similar confinement efficiency and related plasma parameters as the conventional configuration under some appropriate current profile conditions. In AC operation mode, the code can provide useful basic numerical analysis tools for the optimization of the equilibrium control coil current, the design of divertor plates, stability analysis, transport research and so on. Therefore, it may lay an important foundation for the successful realization of advanced high-parameter and high-confinement-efficiency AC operation in the future.
This work was supported by the Comisión Nacional de Energía Atómica, CNEA (Controlled Nuclear Fusion Program) and CONICET (Grant PIP 11220200101929CO).
[1] |
Khayrutdinov R R and Lukash V E 1993 J. Comput. Phys. 109 193 doi: 10.1006/jcph.1993.1211
|
[2] |
Lukash V E, Khayrutdinov R R and Dokuka V 2012 https://www.slideserve.com/howie/common-description-of-dina-code.
|
[3] |
Jardin S C, Pomphrey N and DeLucia J 1986 J. Comput. Phys. 66 481 doi: 10.1016/0021-9991(86)90077-X
|
[4] |
Sadykov A D et al 2015 Nucl. Fusion 55 043017 doi: 10.1088/0029-5515/55/4/043017
|
[5] |
Hertout P et al 2011 Fusion Eng. Des. 86 1045 doi: 10.1016/j.fusengdes.2011.03.092
|
[6] |
Sychugov D Yu et al 2023 Plasma Sci. Technol. 25 035602 doi: 10.1088/2058-6272/ac9f05
|
[7] |
Liu J, Zhu P and Li H 2022 Phys. Plasmas 29 084502 doi: 10.1063/5.0096256
|
[8] |
Lackner K 1976 Comput. Phys. Comm. 12 33 doi: 10.1016/0010-4655(76)90008-4
|
[9] |
Johnson J L et al 1979 J. Comp. Phys. 32 212 doi: 10.1016/0021-9991(79)90129-3
|
[10] |
Jeon Y 2015 J. Korean Phys. Soc. 67 843 doi: 10.3938/jkps.67.843
|
[11] |
Hofmann F and Jardin S C 1990 Nucl. Fusion 30 2013 doi: 10.1088/0029-5515/30/10/003
|
[12] |
Fukushima T 2009 Celest. Mech. Dyn. Astron. 105 305 doi: 10.1007/s10569-009-9228-z
|
[13] |
Farengo R 2023 Plasma Phys. Control. Fusion 65 055003 doi: 10.1088/1361-6587/acc484
|
[14] |
Simpson J et al 2001 NASA/TM-2013-217919
|
[15] |
Aydemir A Y et al 2019 Phys. Plasmas 26 112503 doi: 10.1063/1.5124138
|
[16] |
Luxon J L and Davis L G 1985 Fusion Technol. 8 441 doi: 10.13182/FST85-A40083
|
[17] |
Leuer J A et al 2011 Nucl. Fusion 51 063038 doi: 10.1088/0029-5515/51/6/063038
|
[1] | H R MIRZAEI, R AMROLLAHI. Design, simulation and construction of the Taban tokamak[J]. Plasma Science and Technology, 2018, 20(4): 45103-045103. DOI: 10.1088/2058-6272/aaa669 |
[2] | Rui MA (马瑞), Fan XIA (夏凡), Fei LING (凌飞), Jiaxian LI (李佳鲜). Acceleration optimization of real-time equilibrium reconstruction for HL-2A tokamak discharge control[J]. Plasma Science and Technology, 2018, 20(2): 25601-025601. DOI: 10.1088/2058-6272/aa9432 |
[3] | Hailong GAO (高海龙), Tao XU (徐涛), Zhongyong CHEN (陈忠勇), Ge ZHUANG (庄革). Plasma equilibrium calculation in J-TEXT tokamak[J]. Plasma Science and Technology, 2017, 19(11): 115101. DOI: 10.1088/2058-6272/aa7f26 |
[4] | Shanwen ZHANG (张善文), Yuntao SONG (宋云涛), Kun LU (陆坤), Zhongwei WANG (王忠伟), Jianfeng ZHANG (张剑峰), Yongfa QIN (秦永法). Thermal analysis of the cryostat feed through for the ITER Tokamak TF feeder[J]. Plasma Science and Technology, 2017, 19(4): 45601-045601. DOI: 10.1088/2058-6272/aa57ec |
[5] | JIANG Chunyu (蒋春雨), CAO Jing (曹靖), JIANG Xiaofei (蒋小菲), ZHAO Yanfeng (赵艳凤), SONG Xianying (宋先瑛), YIN Zejie (阴泽杰). Real-Time Bonner Sphere Spectrometry on the HL-2A Tokamak[J]. Plasma Science and Technology, 2016, 18(6): 699-702. DOI: 10.1088/1009-0630/18/6/19 |
[6] | HONG Rongjie (洪荣杰), YANG Zhongshi (杨钟时), NIU Guojian (牛国鉴), LUO Guangnan (罗广南). A Molecular Dynamics Study on the Dust-Plasma/Wall Interactions in the EAST Tokamak[J]. Plasma Science and Technology, 2013, 15(4): 318-322. DOI: 10.1088/1009-0630/15/4/03 |
[7] | FU Jia (符佳), LI Yingying (李颖颖), SHI Yuejiang (石跃江), WANG Fudi (王福地), ZHANG Wei (张伟), LV Bo (吕波), HUANG Juang (黄娟), WAN Baonian (万宝年), ZHOU Qian (周倩). Spectroscopic Measurements of Impurity Spectra on the EAST Tokamak[J]. Plasma Science and Technology, 2012, 14(12): 1048-1053. DOI: 10.1088/1009-0630/14/12/03 |
[8] | LI Li(李莉), LIU Yue (刘悦), XU Xinyang(许欣洋), XIA Xinnian(夏新念). The Effect of Equilibrium Current Profiles on MHD Instabilities in Tokamaks[J]. Plasma Science and Technology, 2012, 14(1): 14-19. DOI: 10.1088/1009-0630/14/1/04 |
[9] | HE Zhixiong, DONG Jiaqi, HE Hongda, JIANG Haibin, GAO Zhe, ZHANG Jinhua. MHD Equilibrium Configuration Reconstructions for HL-2A Tokamak[J]. Plasma Science and Technology, 2011, 13(4): 424-430. |
[10] | GUO Wei, WANG Shaojie, LI Jiangang. Vacuum Poloidal Magnetic Field of Tokamak in Alternating-Current Operation[J]. Plasma Science and Technology, 2010, 12(6): 657-660. |
1. | Hao, G., Xu, J., Sun, Y. et al. Summary of the 11th Conference on Magnetic Confined Fusion Theory and Simulation. Plasma Science and Technology, 2024, 26(10): 101001. DOI:10.1088/2058-6272/ad5d8a |
1. | Hao, G., Xu, J., Sun, Y. et al. Summary of the 11th Conference on Magnetic Confined Fusion Theory and Simulation. Plasma Science and Technology, 2024, 26(10): 101001. DOI:10.1088/2058-6272/ad5d8a |
EAST | CFETR | |||||
Coils | I_{i} (EFIT) | I_{i} (ATEC) | Coils | I_{i} (TEQ) | I_{i} (ATEC) | |
PF1 | 5.4370e+05 | 5.0517e+05 | CS1U | −6.0194e+06 | −5.4287e+06 | |
PF2 | 5.4270e+05 | 5.0308e+05 | CS2U | 2.7232e+06 | 1.0234e+06 | |
PF3 | −1.4077e+05 | 7.3277e+04 | CS3U | 8.3238e+06 | 1.3514e+07 | |
PF4 | 1.7036e+05 | 1.9780e+05 | CS4U | 8.2293e+06 | −2.0077e+06 | |
PF5 | 3.8309e+05 | −1.2694e+05 | CS1L | −1.0959e+07 | −1.1195e+07 | |
PF6 | 5.7632e+04 | −2.0271e+04 | PF1U | 1.2451e+07 | 1.4823e+07 | |
PF7 | 3.5643e+05 | 1.0275e+06 | PF2U | −2.6124e+06 | −3.2489e+06 | |
PF8 | 3.6069e+05 | 8.9066e+05 | PF3U | −6.9199e+06 | −6.9212e+06 | |
PF9 | 1.6525e+06 | 1.1513e+06 | PF3L | −3.1210e+06 | −1.5712e+06 | |
PF10 | 1.6723e+06 | 1.0563e+06 | CS2L | 2.6956e+05 | 5.1617e+05 | |
PF11 | −4.9521e+05 | −4.9535e+05 | CS3L | 1.6669e+07 | 1.8536e+07 | |
PF12 | −4.8425e+05 | −4.3512e+05 | CS4L | 1.4653e+07 | 1.1879e+07 | |
PF13 | 4.3716e+04 | 2.9533e+04 | PF2L | −7.8468e+06 | −1.0188e+07 | |
PF14 | 4.4703e+04 | 5.4555e+04 | PF1L | 1.7482e+07 | 1.2930e+07 | |
DC1 | 4.9374e+06 | 8.9734e+06 |
EAST | CFETR | |||||
Coils | I_{i} (EFIT) | I_{i} (ATEC) | Coils | I_{i} (TEQ) | I_{i} (ATEC) | |
\sqrt{\sum\displaystyle I_{i}^{2}} | 2.6566e+06 | 2.3018e+06 | \sqrt{\sum\displaystyle I_{i}^{2}} | 3.7476e+07 | 3.7681e+07 | |
\left|\dfrac{\psi_{{\rm e}}-\psi_{{\rm e}}^{0}}{\psi_{{\rm e}}^{0}}\right| | 1.8720e−02 | 6.5512e−05 | \left|\dfrac{\psi_{{\rm e}}-\psi_{{\rm e}}^{0}}{\psi_{{\rm e}}^{0}}\right| | 8.6541e−03 | 2.7569e−04 | |
\left|\dfrac{\psi_{{\rm axis}}-\psi_{{\rm axis}}^{0}}{\psi_{{\rm axis}}^{0}}\right| | 1.2263e−02 | 8.3719e−06 | \left|\dfrac{\psi_{{\rm axis}}-\psi_{{\rm axis}}^{0}}{\psi_{{\rm axis}}^{0}}\right| | 3.2759e−03 | 2.0890e−04 | |
\dfrac{\left|X_{{\rm P}}^{0}-X_{{\rm P}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 6.6993e−03 | 6.6993e−03 | \dfrac{\left|X_{{\rm P}}^{0}-X_{{\rm P}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 6.2120e−03 | 6.2120e−03 | |
\dfrac{\left|X_{{\rm axis}}^{0}-X_{{\rm axis}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 1.9714e−03 | 1.9714e−03 | \dfrac{\left|X_{{\rm axis}}^{0}-X_{{\rm axis}}\right|}{\left|X_{{\rm P}}^{0}-X_{{\rm axis}}^{0}\right|} | 5.7144e−03 | 5.7144e−03 |
Coils | HL-CREC (A) | UD-CREC (A) | C-CREC (A) | NEC (A) |
PF1 | 1.4131e+06 | 8.7750e+05 | 1.4433e+06 | 1.3968e+06 |
PF2 | 2.5396e+06 | −1.6512e+05 | −2.6781e+05 | 1.0406e+06 |
PF3 | 5.3719e+06 | 3.1279e+06 | −1.7004e+06 | 1.0874e+06 |
PF4 | 1.7050e+06 | −2.3708e+05 | −1.2162e+06 | 1.2288e+06 |
PF5 | −1.5675e+07 | 4.8045e+06 | −1.1617e+06 | 1.8051e+06 |
PF6 | −5.3236e+06 | −8.2558e+05 | −6.3073e+04 | 2.0526e+06 |
PF7 | 3.0715e+07 | −3.8904e+05 | 4.6209e+06 | 5.9412e+05 |
PF8 | 1.5813e+07 | 1.1429e+06 | 2.0429e+06 | 2.2780e+05 |
PF9 | −2.8630e+07 | −4.5500e+06 | −1.1784e+06 | −1.6729e+06 |
PF10 | −1.5845e+07 | 1.9476e+06 | 9.8322e+05 | −1.4072e+06 |
PF11 | 2.2270e+06 | 2.9918e+06 | −1.2139e+06 | 3.7252e+05 |
PF12 | 2.1009e+06 | −2.9352e+05 | −8.7355e+05 | 5.4823e+05 |
PF13 | −1.3933e+06 | −9.3773e+05 | 3.8834e+05 | 2.2159e+05 |
PF14 | −1.3285e+06 | −7.8612e+05 | −1.8577e+05 | −1.6968e+04 |
\sqrt{\sum I_{i}^{2}} | 5.0907e+07 | 8.4196e+06 | 6.1792e+06 | 4.3440e+06 |
\beta_{{\rm p}} | 1.2630e+00 | 1.0259e+00 | 1.1347e+00 | 1.7049e+00 |
Coils | HL-CREC (A) | UD-CREC (A) |
PF1 | 2.2298e+06 | 1.9104e+06 |
PF2 | 2.2664e+06 | −6.6681e+05 |
PF3 | 5.8628e+05 | 2.9801e+06 |
PF4 | 7.3105e+05 | −2.3493e+05 |
PF5 | 3.7842e+04 | 9.5368e+05 |
PF6 | 1.0214e+05 | 3.0989e+05 |
PF7 | 1.2618e+05 | −7.8916e+05 |
PF8 | 9.6939e+04 | 9.7855e+05 |
PF9 | 1.9136e+05 | −8.3600e+05 |
PF10 | 1.6128e+05 | 8.8734e+05 |
PF11 | 1.1100e+06 | 1.6053e+06 |
PF12 | 1.2024e+06 | 4.7541e+05 |
PF13 | −1.0344e+06 | −6.1468e+04 |
PF14 | −9.9634e+05 | −1.3812e+05 |
\sqrt{\sum I_{i}^{2}} | 3.9784e+06 | 4.6709e+06 |
I_{{\rm pabs}} | 1.1572e+06 | 1.02595e+06 |