New method to design stellarator coils without the winding surface

TL;DR

This paper introduces a novel coil design method for stellarators that does not rely on the traditional winding surface assumption, allowing for more flexible and potentially easier-to-build coil configurations.

Contribution

The authors develop a new coil design approach that represents coils as arbitrary curves and does not require a winding surface, expanding the design space for stellarator coils.

Findings

The method successfully designs coils for a simple configuration.

Applications to W7-X and LHD plasmas demonstrate its effectiveness.

The approach simplifies coil design constraints.

Abstract

Finding an easy-to-build coils set has been a critical issue for stellarator design for decades. Conventional approaches assume a toroidal "winding" surface. We'll investigate if the existence of winding surface unnecessarily constrains the optimization, and a new method to design coils for stellarators is presented. Each discrete coil is represented as an arbitrary, closed, one-dimensional curve embedded in three-dimensional space. A target function to be minimized that covers both physical requirements and engineering constraints is constructed. The derivatives of the target function are calculated analytically. A numerical code, named FOCUS, has been developed. Applications to a simple configuration, the W7-X, and LHD plasmas are presented.

Equations

$$R$$

$$Z$$

$${\left[\begin{array}{c}\mathbf{x}\\ \mathbf{t}\\ \mathbf{n}\\ \mathbf{b}\end{array}\right]}^{\prime }=\left[\begin{array}{ccc}1& 0& 0\\ 0& \kappa & 0\\ -\kappa & 0& \tau \\ 0& \tau & 0\end{array}\right]\left[\begin{array}{c}\mathbf{t}\\ \mathbf{n}\\ \mathbf{b}\end{array}\right],$$

$$x(t)$$

$${\chi }^2(\mathbf{X})=\sum _j{w}_j{\left(\frac{f_j(\mathbf{X})-f_{j,o}}{f_{j,o}}\right)}^2,$$

$$f_B(\mathbf{X})\equiv {\int }_S\frac{1}{2}{\left(\frac{\mathbf{B}\cdot \mathbf{n}}{\left|\mathbf{B}\right|}\right)}^{2}ds.$$

$${\mathbf{B}}_V(\bar{\mathbf{x}})=\frac{{\mu }_0}{4\pi }\sum _{i=1}^{N_C}{I}_i{\int }_{C_i}\frac{d{\mathbf{l}}_i\times \mathbf{r}}{r^3}.$$

$$\delta \mathbf{B}(\bar{\mathbf{x}})=\frac{{\mu }_0}{4\pi }{I}_i{\int }_0^{2\pi }\left[\frac{3\mathbf{r}\cdot {\mathbf{x}}_i^{\prime }}{r^5}\mathbf{r}\times \delta {\mathbf{x}}_i+\frac{2}{r^3}\delta {\mathbf{x}}_i\times {\mathbf{x}}_i^{\prime }+\frac{3\mathbf{r}\cdot \delta {\mathbf{x}}_i}{r^5}{\mathbf{x}}_i^{\prime }\times \mathbf{r}\right]dt.$$

$$\delta f_B(\mathbf{X})={\int }_S(\mathbf{B}\cdot \mathbf{n})\frac{{\mu }_0}{4\pi }{I}_i{\int }_0^{2\pi }\left[\frac{3\mathbf{r}\cdot {\mathbf{x}}_i^{\prime }}{r^5}\mathbf{n}\times \mathbf{r}+\frac{2}{r^3}{\mathbf{x}}_i^{\prime }\times \mathbf{n}+\frac{3{\mathbf{x}}_i^{\prime }\times \mathbf{r}\cdot \mathbf{n}}{r^5}\mathbf{r}\right]\cdot \delta {\mathbf{x}}_{i}dtds.$$

$$f_{\Psi }(\mathbf{X})\equiv \frac{1}{2\pi }{\int }_0^{2\pi }\frac{1}{2}{\left(\frac{{\Psi }_{\zeta }-{\Psi }_o}{{\Psi }_o}\right)}^{2}d\zeta ,$$

$${\Psi }_{\zeta }(\mathbf{X})$$

$$\mathbf{A}(\bar{\mathbf{x}})=\frac{{\mu }_0}{4\pi }\sum _{i=1}^{N_C}{I}_i{\int }_{C_i}\frac{d{\mathbf{l}}_i}{r}.$$

$$\delta f_{\Psi }(\mathbf{X})$$

$$\delta \mathbf{A}(\bar{\mathbf{x}})$$

$$f_L(\mathbf{X})=\frac{1}{N_C}\sum _{i=1}^{N_C}\frac{e^{L_i}}{e^{L_{i,o}}},$$

$$L_i(\mathbf{X})={\int }_0^{2\pi }|{\mathbf{x}}_i^{\prime }|dt,$$

$$\delta f_L(\mathbf{X})=\frac{1}{N_C}\frac{e^{L_i}}{e^{L_{i,0}}}{\int }_0^{2\pi }\frac{({\mathbf{x}}_i^{\prime }\cdot {\mathbf{x}}_i^{\prime }){\mathbf{x}}_i^{\prime \prime }-({\mathbf{x}}_i^{\prime }\cdot {\mathbf{x}}_i^{\prime \prime }){\mathbf{x}}_i^{\prime }}{({\mathbf{x}}_i^{\prime }\cdot {\mathbf{x}}_i^{\prime }{)}^{3/2}}\cdot \delta {\mathbf{x}}_{i}dt.$$

$$M_i(\mathbf{X})=\sum _{n=1}^{N_F}{n}^p({X_{c,n}^i}^2+{X_{s,n}^i}^2+{Y_{c,n}^i}^2+{Y_{s,n}^i}^2+{Z_{c,n}^i}^2+{Z_{s,n}^i}^2),$$

$${\chi }^2(\mathbf{X})=w_B{f}_B(\mathbf{X})+w_{\Psi }{f}_{\Psi }(\mathbf{X})+w_L{f}_L(\mathbf{X}).$$

$$\frac{\partial f_B}{\partial X_{c,n}^i}$$

$$\approx \sum _{j=1}^{N_{\theta }}\sum _{k=1}^{N_{\zeta }}{\left[\left(B_x{n}_x+B_y{n}_y+B_z{n}_z\right)\left(\frac{\partial B_x}{\partial X_{c,n}^i}{n}_x+\frac{\partial B_y}{\partial X_{c,n}^i}{n}_y+\frac{\partial B_z}{\partial X_{c,n}^i}{n}_z\right)\sqrt{g}\right]}_{j+\frac{1}{2},k+\frac{1}{2}}\frac{2\pi }{N_{\theta }}\frac{2\pi }{N_{\zeta }}.$$

$$\frac{\partial B_x}{\partial X_{c,n}^i}$$

$$\frac{\partial \mathbf{X}}{\partial \tau }\equiv -\frac{\partial {\chi }^2}{\partial \mathbf{X}}-{\mathbf{X}}^{\prime }\left(\frac{\partial M}{\partial \mathbf{X}}\cdot {\mathbf{X}}^{\prime }\right).$$

$$\frac{\partial {\chi }^2}{\partial \tau }=\frac{\partial {\chi }^2}{\partial \mathbf{X}}\cdot \frac{\partial \mathbf{X}}{\partial \tau }=-\frac{\partial {\chi }^2}{\partial \mathbf{X}}\cdot \frac{\partial {\chi }^2}{\partial \mathbf{X}}-\left(\frac{\partial {\chi }^2}{\partial \mathbf{X}}\cdot {\mathbf{X}}^{\prime }\right)\left(\frac{\partial M}{\partial \mathbf{X}}\cdot {\mathbf{X}}^{\prime }\right)=-{\left(\frac{\partial {\chi }^2}{\partial \mathbf{X}}\right)}^2.$$

$$\frac{\partial M}{\partial \tau }=\frac{\partial M}{\partial \mathbf{X}}\cdot \frac{\partial \mathbf{X}}{\partial \tau }=-{\left(\frac{\partial M}{\partial \mathbf{X}}\cdot {\mathbf{X}}^{\prime }\right)}^2.$$

$$\mathrm{link}(C_i,Plasma)=\frac{1}{4\pi }{\oint }_{C_i}{\oint }_{axis}\frac{{\mathbf{r}}_1-{\mathbf{r}}_2}{|{\mathbf{r}}_1-{\mathbf{r}}_2{|}^3}\cdot d{\mathbf{r}}_1\times d{\mathbf{r}}_2.$$

$$R=$$

$$Z=$$

$$x$$

$$y$$

$$z$$

Full text

1 Introduction

The difficulties in designing and fabricating current-carrying coils to produce the magnetic fields required for confining plasmas for the purpose of creating fusion energy has been a critical problem since the beginning of research into magnetically confined plasmas in the 1950s. This problem remains crucial, as is evidenced by the recent NCSX [1] stellarator, which was partially built at Princeton Plasma Physics Laboratory but, because of budget constraints, was ultimately cancelled; and by the construction delays of the W7-X experiment [2] recently completed in Germany. The construction of the coils is only one component of modern experiments; but, realizing that it is the currents in the coils that produce the “magnetic bottle” that confines the plasma, it is easy to understand that designing and accurately constructing suitable coils is paramount.

Tremendous efforts have been made to find optimal coils that meet both the “physics” requirements of producing the desired magnetic field as precisely as possible; and the “engineering” requirement that it is realistically possible to construct the coils using existing technology. Coils that are well-separated, with ample space for the vacuum vessel and other components of the experiment, and with simple shapes are easier and cheaper to build. The task of identifying suitable coil configurations has a direct impact of the cost of magnetic confinement experiments, and thus also on fusion research, and indeed ultimately on the cost of the electricity [3].

Early stellarator designs, which pre-date the advent of modern computational architectures, usually began with the coils being described analytically, like the general “winding law” [4] by which the current filaments on a given toroidal surface are represented using trigonometric functions. By varying the coil parameters, the geometry of the coils can be varied to simultaneously obtain optimal plasma properties and to satisfy engineering constraints. However, the winding law represents a very restricted class of coil configurations.

Another approach, introduced by Nührenberg & Zille [5], allowed the development of a “fully-optimized” stellarator by separating the problem of designing a high-performance plasma from the task of designing a suitable coils set. First, an attractive plasma equilibrium is identified by using a nonlinear optimization, where the equilibrium is typically calculated [6] using the macroscopic model of ideal magnetohydrodynamics, for which the geometry of the plasma boundary, ${\cal S}$, and various profiles such as the assumed pressure and current are required as boundary conditions. These boundary conditions are adjusted to discover an equilibrium that is stable with respect to small external perturbations and internal oscillations, reduces the transport and turbulent properties, and so forth.

The second stage of the design, which is the topic of this article, is to design a set of discrete current-carrying coils that creates the required external “vacuum” field for confining the “reference plasma” configuration. Hereafter, we will assume that a desired, reference plasma configuration is provided.

The vacuum field, ${\bf B}_{V}$, by which we mean the magnetic field produced by currents external to the plasma domain, must balance the magnetic field, ${\bf B}_{P}$, produced by the currents that may or may not be present in the plasma, so that the normal component of the total magnetic field, ${\bf B}={\bf B}_{V}+{\bf B}_{P}$, on ${\cal S}$ is zero.

For vacuum fields, the Ampere’s Law in ideal magnetohydrodynamics (MHD) equations reduces to $\nabla\times{\bf B}_{V}=0$. Together with the magnetic divergence constraint $\nabla\cdot{\bf B}_{V}=0$, we can easily derive $\nabla^{2}\phi=0$, where $\phi$ is the magnetic scalar potential. Thus, the coil determination problem is to solve Laplace’s equation with the boundary condition of $\mathbf{B}_{V}\cdot\mathbf{n}=-\mathbf{B}_{P}\cdot\mathbf{n}$ on ${\cal S}$ ($\mathbf{n}$ is the unit surface normal). The plasma boundary is usually specified by Fourier harmonics $R_{mn},Z_{mn}$ over the poloidal angle $\theta$ and toroidal angle $\zeta$ (restricting attention to stellarator symmetry [7] for convenience),

[TABLE]

Pioneering work in the field of coil design was performed by Merkel with the development of NESCOIL code [8], in which he assumed that the external magnetic field is produced by a surface current distribution on a closed toroidal surface surrounding the plasma. This toroidal surface constrains the location of the coils and is called the “current carrying surface” or the “winding surface”. The surface current density is expressed by a current potential $\Phi$ on the winding surface, $\mathbf{j}=\mathbf{n}\times\nabla\Phi$, where ${\bf n}$ is normal to the winding surface. A Green’s function method is then applied to solve the current potential distribution which minimizes the squared normal error $\displaystyle\epsilon^{2}=\oint_{\cal S}(\mathbf{B}\cdot\mathbf{n})^{2}ds$. Once the surface current potential is determined, a set of discretized coils can be obtained by selecting an appropriate number of contours of $\Phi$. This method was successfully applied to the design of the coils of W7-X [9], the world’s largest ever stellarator.

Merkel’s method leads to a Neumann condition problem that can be linearly solved. So it’s inherently fast and robust. However, due to the ill-conditioned nature of inverting the Biot-Savart integrals, the NESCOIL results for plasmas with complex shapes are usually not viable (e.g. unreasonably large amplitudes in high modes or impractically complex curvatures). A singular value decomposition (SVD) method [10] and, more recently, a Tikhonov regularization approach [11] were applied to provide improvements on NESCOIL. Nevertheless, these methods only indirectly control the geometry of the resultant coils set, and allow limited opportunities to optimize the engineering constraints.

A different approach that explicitly incorporates engineering constraints has been advanced by Drevlak with the extended NESCOIL code [12] and ONSET [13], by Strickler et al. with the code COILOPT [14], and later by Breslau et al. with COILOPT++ [15]. The coils are represented as “filaments”, one-dimensional curves, lying on a toroidal winding surface (pre-defined or optimized simultaneously). The magnetic field produced by $\delta$-function current-densities in the coils set is calculated using the Biot-Savart law. The geometry of the coils is varied using nonlinear optimization algorithms to minimize a “cost-function”, which represents a balance between the physics requirements (that the total normal magnetic field at the plasma boundary is as small as possible) and the engineering constraints (that the coils can be achieved by current engineering techniques). Coils of the compact stellarators NCSX [16] and ARIES-CS [17] are designed in such approaches.

For all the methods mentioned above, a toroidal winding surface is required to locate the coils (ONSET may need two constraining surfaces for interpolation). For NESCOIL and NESCOIL-like codes (NESVD, REGCOIL, etc.), the winding surface is strictly required, as it provides the surface where the current potential lies so that the currents can be explicitly represented and then be solved. For the discrete coils nonlinear optimization codes, the winding surface is not strictly required, but there are advantages. This decreases the number of freedoms and the plasma-coil separation, which is important for divertor/blanket access, can be easily enforced by the existence of winding surface.

However, on the other hand, the existence of winding surface also provides strong limitations. A “bad” winding surface directly results in the failure of finding acceptable coils set. Thus, before getting an acceptable coils set, a “good” winding surface should be obtained first. However, the standard of good or bad of a winding surface is actually evaluated by the final coil shapes. So nested optimizations targeting both the winding surfaces and the coils are required to be carried out back and forth.

In this paper, we present a new method for designing stellarator coils that eliminates the need of a toroidal winding surface. The new coil representation, as well as the interpretations on the definition of objective functions and their analytical derivatives, are discussed in section II. The numerical implementation, a code named Flexible Optimized Coils Using Space curves (FOCUS), is introduced in section III. In section IV, we show some example calculations. Conclusions and future works are discussed in section V.

2 3D coil representations

2.1 The Fundamental Theorem of Curves

For convenience, we work within the filamentary approximation, by which we mean that each coil has zero cross-sectional area. (This approximation is not fundamental, and can easily be relaxed in future work.) We represent the geometry of the coils as closed, one-dimensional curves embedded in three-dimensional space.

The geometry of such curves is elegantly described by the Fundamental Theorem of Curves [18], which states that any regular curves in three-dimensional space with non-zero curvature can be determined, up to rigid displacements and rotations, by its curvature, $\kappa(s)$, and torsion, $\tau(s)$, where $s$ parameterizes arc-length. Given $\kappa$ and $\tau$, and boundary conditions (such as the starting point), the geometry, ${\bf x(s)}$, unit tangent vector, ${\bf t(s)}$, unit normal and bi-normal, ${\bf n(s)}$ and ${\bf b(s)}$, of the curve can be obtained by integration of the Frenet-Serret formulas [19],

[TABLE]

where the “prime” denotes derivative with respect to the parameter $s$, i.e., ${\bf x}^{\prime}\equiv d{\bf x}/ds$.

The benefit is that each curve is described by two independent functions of arc length, namely the curvature and torsion. Using the curvature function as an independent degree of freedom to implicitly define the coil geometry has its appeal, as the curvature of the coils is frequently an engineering constraint and any curvature constraints can potentially be enforced before the nonlinear optimization proceeds. However, solving 12 ODEs for each coil at each iteration of the nonlinear optimization turns out to be somewhat cumbersome. Furthermore, to exploit efficient nonlinear optimization algorithms, we shall later construct the derivatives of the coil geometry with respect to the independent degrees of freedom, and this requires integration of the derivatives of Eq.(1) with respect to the curvature and torsion, which becomes even more cumbersome. Besides, additional constraints, such as ensuring the curves closed and sufficiently smooth, must also be included.

2.2 Fourier representation

Instead, we employ an easier way of describing one-dimensional curves embedded in three-dimensional space. A curve is described directly, and completely generally, in Cartesian coordinates as $\mathbf{x}(t)=x(t)\,{\bf i}+y(t)\,{\bf j}+z(t)\,{\bf k}$. Three functions are required to specify the geometry, namely $x(t)$, $y(t)$ and $z(t)$, with the constraints that each function be periodic, e.g. ${\bf x}(t+T)={\bf x}(t)$ for some $T$. The curve parameter, $t$, at this point is arbitrary.

A variety of mathematical representations are possible. For purpose of illustration, and because our initial interest is in smooth coils, we use a Fourier representation:

[TABLE]

with $t$ varying between $[0,2\pi]$, and similarly for $y(t)$ and $z(t)$. The shape of a coil is then fully determined by the $3\times(2N_{F}+1)$ Fourier coefficients. No additional assumptions are made here, so this representation fits for all kinds of smooth coils, such as helical, modular, saddle, etc. (For coils with straight segments, for example the rectangular and “window pane” coils, the Fourier representation is not suitable. In future upgrades to FOCUS, additional representations for the curves, such as cubic splines and piecewise-linear, will be implemented.)

Hereafter, $\mathbf{X}$ will be used to represent the geometrical degrees of freedom of a set of $N_{C}$ coils, ${\bf x}_{i}$, and the coil currents, $I_{i}$, which are also considered as independent degree of freedom. The total number of degrees of freedom is $N_{D}=N_{C}\times[3\times(2N_{F}+1)+1]$.

3 Theory of the coil optimization problem

3.1 General objective functions

The coil parameters are to be varied to minimize a target function consisting of both “physics” and “engineering” objective functions,

[TABLE]

where $f_{j}({\mathbf{X}})$ is the value of the $j^{th}$ objective function, to be defined below, for a given set of coil parameters, $f_{j,o}$ denotes the desired value, and $w_{j}$ is a user-prescribed weight.

3.2 Normal field error on plasma boundary

The key criteria is to satisfy the condition of $B_{n}=\mathbf{B}\cdot\mathbf{n}=0$ on the pre-defined plasma boundary as much as possible with a finite number of coils. A smaller $B_{n}$ error usually implies a better approximation to the target equilibrium properties. This can be achieved by minimizing the squared relative error of the total normal field integrated over the boundary,

[TABLE]

The magnetic field at position $\mathbf{\bar{x}}$ produced by a coils set is calculated using the Biot-Savart law,

[TABLE]

Here, $I_{i}$ is the current in $i^{th}$ coil, $d\mathbf{l}_{i}\equiv{\bf x}_{i}^{\prime}\,dt$ is the differential line element, and $\mathbf{r}=\mathbf{\bar{x}}-\mathbf{x}_{i}$ is the displacement vector between the evaluate point on the surface and the differential element.

The plasma field, $\mathbf{B}_{P}$, produced by the reference plasma must be pre-computed using a suitable equilibrium code, and the following will assume that the plasma field does not change during the optimization of the coil geometry. The FOCUS code can accommodate an arbitrary $\mathbf{B}_{P}$, but for simplicity of illustration in this paper we restrict attention to vacuum fields, so that $\mathbf{B}=\mathbf{B}_{V}$, and we hereafter suppress the subscript $V$ for notational clarity.

If a small geometry deformation $\delta\mathbf{x}_{i}$ is applied to the $i^{th}$ coil, the corresponding change in $\mathbf{r}$ and $r$ are $\delta\mathbf{r}=-\delta\mathbf{x}_{i}$ and $\delta r=-\mathbf{r}\cdot{\delta\mathbf{x}_{i}}/r$, and the change in the magnetic field is

[TABLE]

The change in $f_{B}$ is, omitting the normalization term for notational clarity, given by

[TABLE]

What should be noted here is that the variation of $f_{B}$ (as well as other object functions) with respect to coil currents $\delta I_{i}$ can be derived in the same way. But we will not cover the details it in this article.

3.3 Constraint of the toroidal flux

To avoid trivial solutions, like when $I_{i}\rightarrow 0$ $\forall i$, $f_{B}\rightarrow 0$, it is sufficient to constrain the enclosed toroidal flux. If ${\bf B}\cdot{\bf n}=0$ on the boundary, then the toroidal flux through any poloidal cross-sectional surfaces is constant. We include an objective function defined as

[TABLE]

where the flux through a poloidal surface, ${\cal T}$, produced by cutting the boundary with plane $\zeta=const.$ is computed using Stokes’ theorem,

[TABLE]

Here $d\mathbf{l}$ is on the boundary curve of the poloidal surface and the total magnetic vector potential $\mathbf{A}$ is

[TABLE]

The variation of $f_{\Psi}$ resulting from $\delta\mathbf{x_{i}}$ is

[TABLE]

where $\displaystyle\delta\Psi_{\zeta}=\int_{\partial{\cal T}}\delta\mathbf{A}\cdot d\mathbf{l}$ and

[TABLE]

3.4 Engineering constraints

Besides the physical requirements, engineering constraints should be considered, such as the coil-coil separation, coil-plasma separation, minimum radius of curvature, electromagnetic forces, stored energy, etc.

Some of the engineering constraints are actually consistent with the magnetic field requirements. For example, to reduce the magnetic field ripple, the coils are usually placed away from the plasma, which is exactly of the same effect as maximizing the coil-plasma separation constraint for the capacity of divertors. And the coil-coil separation constraint is to prevent coils overlapping with each other. It also coincides with the need of lowering magnetic ripple. With finite number of coils, it’s better to equally place the coils in the toroidal angle; otherwise, the ripple would increase rapidly.

There are also engineering considerations that compete with the magnetic field requirement. Among them, the coil length constraint is the most straightforward. Without a constraint on length, the coils can become arbitrarily long to lower the ripple, and more “wiggles” can potentially be formed to better match the plasma shape. Besides, the total coil length is directly related to the usage of materials, i.e. cost, which is without any relation to the physical requirements. We include an objective function of the form

[TABLE]

where $L_{i}({\bf X})$ is the length of $i$-th coil,

[TABLE]

and $L_{i,0}$ is a user-specified normalization. The variations in $f_{L}$ resulting from a variation $\delta\mathbf{x}_{i}$ is

[TABLE]

3.5 Spectral condensation

The parameterization given in Eq.(2) is not unique. Consider the curve variation $\delta{\bf x}_{i}={\bf x}^{\prime}_{i}(t)\delta u$, where $\delta u$ is a small constant. Inserting this into Eq.(7), Eq.(11) and Eq.(15), we can easily get that $\delta f_{B}$, $\delta f_{\Psi}$ and $\delta f_{L}$ are all zero. This “variation” is tangential to the curve and leaves the geometry of curve unchanged, but it does effectively change the parameterization of the curve and thereby changes the Fourier harmonics.

This non-uniqueness of the curve parameters can be exploited for numerical efficiency. For each coil, a spectral width of the Fourier harmonics is defined as, [20]

[TABLE]

where $p\geq 1$. Curve parameterizations that minimize $M=\sum_{i}^{N_{C}}M_{i}$ are attractive for numerical purposes, but minimizing the spectral width $M$ should not compete with the physics or engineering properties. So when minimizing $M$, we only allow tangential variations in the geometry of coils. In the Fourier parameter space, the tangential variation is defined as $\delta{\bf X}\equiv{\bf X}^{\prime}\delta u$, where ${\bf X}^{\prime}$ is the tangential direction in parameter space and $\delta u$ is an arbitrary function. If we choose $\delta u\equiv-\partial_{\bf X}M\cdot{\bf X}^{\prime}$, the spectral width will be forced to decrease subject without changing the coil geometry. This shall be used below.

4 Numerical implementations

4.1 Analytically calculated derivatives

Our target function is now composed of the normal field error, toroidal flux error and the length constraint, and is represented as,

[TABLE]

The task of designing optimal coils is now reduced to a single-object minimization problem. A lot of minimization algorithms can be applied here, like the famous Levenberg-Marquardt (LM) method [21] which was successfully used in COILOPT and COILOPT++. However, unlike other nonlinear coil-designing codes where the derivatives are computed numerically, we can offer the analytically calculated derivatives from the Eq.(7), Eq.(11) and Eq.(15).

Calculating the derivatives involves line and surface integrals. Each coil is treated as a single filament with $N_{s}$ segments. And the predefined plasma surface is discretized into $N_{\theta}\times N_{\zeta}$ surface elements. We shall illustrate how we implement the numerical calculations in the code by one example. The derivative of the normal field error $f_{B}$ with respect to an arbitrary cosine harmonic of the $i^{th}$ coil denoted as $X_{c,n}^{i}$ is computed as,

[TABLE]

Here, the normalization term is omitted and the integral is evaluated by the mid-point rule. And the calculation of $\partial{B_{x}}/\partial{X_{c,n}^{i}}$ , which involves a line integral, can be approximated like,

[TABLE]

Similarly, the derivatives for toroidal flux error and coil length constraint can also be calculated (fast and accurately).

4.2 Steepest descent method

Analytically calculated derivatives have great advantages in accelerating the computation speed and accessing more applicable algorithms. For the proof of concept, we will apply the steepest descent algorithm here. Defining an artificial “time” $\tau$, the descent direction for the DoF vector is given by

[TABLE]

Here, in the RHS, the first term is the usual gradient descent and the second is defined in tangential direction specialized for the spectral condensation. So the relationship between the target function and time is:

[TABLE]

The first derivative $\partial\chi^{2}/\partial\tau$ is always non-positive. After reaching the minima ($\partial\chi^{2}/\partial{\bf X}=0$), the spectral term has the same monotonicity,

[TABLE]

By solving Eq.(20), the steepest descent method will keep minimizing the target function, which covers the physical requirements (the squared normal field error and the toroidal flux error) and engineering constraint (the length constraint), as well as the numerical non-uniqueness eliminator (the spectral condensation term) . The ODEs in Eq.(20) are integrated using a fixed order Runge-Kutta method with the subroutine D02BJF in the NAG Library [22].

The steepest descent method is continuous. This means that the coils cannot pass through the plasma, as this would result in a singularity in the normal field. We can define the linking number of each coil and the plasma using the Gauss linking integral [23] with respect to the coil and the magnetic axis:

[TABLE]

With the steepest descent algorithm, this linking number is implicitly conserved during the optimization.

5 Applications

To demonstrate the feasibility of this method, three cases are investigated. A simple rotating elliptical plasma boundary is used for investigating the convergence properties of the code. Illustrations for the existing stellarators W7-X and LHD are presented to explore the code’s capacity of designing various coils configuration for realistic stellarators.

5.1 Rotating ellipse

A two-field period, rotating elliptical plasma boundary is specified as,

[TABLE]

where $N_{P}=2$. As one of the simplest configurations for stellarators, the rotating ellipse can be achieved with conventional helical coils (like the LHD machine [24]). Helical coils, with certain pitch angles, can efficiently produce desired rotational transforms, particularly for classical stellarators with helix structures. However, the access to the plasma and ease of maintenance are restricted badly. Modular coils, on the other hand, can easily overcome these difficulties and they can closely match macroscopic magnetic parameters, but differences in the local structures of the fields do exist [25]. Here, we’ll try to use the FOCUS to design modular coils for it, just for demonstrations.

To initialize the coil-design calculation in FOCUS, we must supply a suitable initial guess. We anticipate that, generally, codes such as NESCOIL and REGCOIL will be used. But here, a very simple initial guess for the coils is made by toroidally placing 16 circular coils, as shown in Fig.1. The radius of the coils is $0.75\ m$, and the chosen Fourier resolution is $N_{F}=4$.

To test the different objective functions separately, we first include only the physical constraints, i.e. we set the length constraint weight to zero, $w_{L}=0$. The limit on the time integration is set to an arbitrary value, e.g. $\tau_{end}=10^{3}$.

As illustrated in Fig.3, $f_{B}$ significantly decreases from $3.63\times 10^{-2}$ to $1.26\times 10^{-5}$ at $\tau=10^{3}$. The final shape of the optimized coils are shown in Fig.1. We can see that all the coils are getting longer and further away from the plasma to better match the plasma and reduce the ripple. A Poincar$\acute{\rm e}$ plot of the flux surfaces produced by the optimized coils, Fig.3, indicates that a satisfactory magnetic field has been obtained.

Note that in this optimization, the only targets were to minimize $w_{B}\;f_{B}+w_{\Psi}\;f_{\Psi}$ ($w_{B}/w_{\Psi}=100/1$). Engineering constraints such as coil-coil separation, coil-plasma separation, even the coil length constraint were not included. The coils are not constrained to lie on a winding surface and are completely free to move in space. The fact that we get reasonable coils is remarkable.

Fig.3 also indicates that the target function is not converged at $\tau=10^{3}$. Actually, it turns out that even at $\tau=10^{6}$, the target function is still not converged. We’ll explore the convergence properties by varying different parameters.

We now vary the Fourier harmonics number, $N_{F}$, to explore its relations with the achieved normal field error. As shown in Fig.4, when the number of Fourier harmonics is restricted to $N_{F}=1$, which means that the coils can only be planar, the normal field error cannot decrease as much amplitude as other cases. The optimized coils, as illustrated in Fig.5, are elliptical and rotated to follow the pattern of the plasma shape. And for $N_{F}\geq 4$ cases, there are no distinct differences. In other words, for this elliptical configuration, the coils can be adequately represented by $N_{F}=4$.

The number of initial coils, $N_{C}$, can also be changed. In principle, the more coils, the more accurate magnetic field we can get. However, more coils imply higher cost and less separations between each coils which are vital for diagnostic ports. Usually, we have to determine the most balanced number of coils for a practical machine. Fig.4 shows the convergence curves with respect to different number of coils. As expected, with more coils, the lower $f_{B}$ has been achieved. Also, we observe that there exists a big gap between the curves $N_{C}=32$ and $N_{C}=36$, which implies that this could be the most effective number. Without enough coils, the shape of the finally optimized coils would be unrealistically complicated, like the coils in Fig.5. The four coils have large twists to produce the desired magnetic field.

For the calculations shown in Fig.4 and Fig.4, the initial coils are circular coils with a arbitrarily chosen radius of $R_{C}=0.75m$. For different initial radii (of course differing in a small range), as shown in Fig.4, we found that the final $B_{n}$ error are trending to be identical. And all the optimized coils are having close geometries (details not shown here), which indicates that they are at the same minima. So, the optimizing process is not very sensitive to the initial radius, at least for this equilibrium.

Now, we explore the effect of length constraint weight. We assign $w_{L}$ with different values ranging from $10^{-6}$ to $1.0$. The results in Fig.4 show that the constraint on coil length are suppressing the decrease of the normal field error. And the evolutions now converge faster as the length weight increases. This can be used to control the convergence of the code and make a balance between the physical accuracy and engineering feasibility.

With the experience of the previous explorations, we can come up with a suitable coils set for the rotating ellipse configuration by enabling the length constraint. All the three constraints are activated with normalized weights so that at the initial state $w_{B}f_{B}=1.0$, $w_{\Psi}f_{\Psi}=0.5$ and $w_{L}f_{L}=0.2$. Starting with 16 circular coils with a radius of $0.75m$, we get a satisfactory coils set at $\tau=10^{3}$ as shown in Fig.5, which makes a reasonable tradeoff between all the constraints. The coils shape now looks much more simple with a comparably smaller size.

5.2 Designing modular coils for the W7-X

The rotating ellipse is just a simple fictitious classical stellarator configuration. It would be more challenging to design coils for the W7-X, due to its complicated geometries. The coils system of W7-X mainly consists 50 non-planar modular coils and 20 planar coils in five-field periods [26]. Modular coils produce the required magnetic field, while the planar coils are mainly used for the flexibility of different rotational transform ($\iota$) scenarios. The last closed flux surface (LCFS) in the “standard case” configuration [27] is used here as the target plasma boundary. To produce this configuration, the planar coils are not used (zero currents) and the currents in modular coils are $1.62$ MA (108 turns $\times$ 15.0 kA/turn).

First of all, some preliminary calculations are performed to verify FOCUS’ availability for the W7-X. The actual W7-X coils (single-filamentary CAD models of the as-built W7-X coils set) are fitted with $N_{F}=6$ Fourier series and taken as the initial coils. One of the simplest verifications is perturbing only one parameter and remaining all the others unchanged to see if FOCUS can recover the perturbed parameter (Here we are assuming that the actual coils are already at/close the minimum.) The current in an arbitrary modular coil (marked as modular A for notation) is changed to $5.0$ MA, and as shown in Fig.7, FOCUS can recover the current to the correct value of $1.62$ MA. Another verification is taken out by optimizing all the $2000$ parameters of the actual 50 modular coils to lower $B_{n}$ error. The original $f_{B}$ is $1.86\times 10^{-4}$. And after a quick FOCUS optimization, the $f_{B}$ is reduced to $3.56\times 10^{-5}$. Although, the differences between the optimized coils and the original ones are visually subtle, as shown in Fig.7.

Now, we’ll start the optimization processes with an arbitrary initial guess to better illustrate the capability of FOCUS. The initial guess comes from a set of circular coils ($R_{coil}=1.25$m) that are positioned surrounding the plasma surface with equal toroidally-interval angles, i.e., the $i^{th}$ coil is placed in the poloidal plane of $\zeta=2\pi(i-1)/{N_{C}}$. The total number of coils ${N_{C}}$ varies from 50 down to 20 to obtain various solutions. All the three constraints introduced above are taken into account with artificially chosen weights.

After the optimizations, Poincar$\acute{\rm e}$ plots of the flux surfaces produced by the magnetic field of the optimized coils, as well as the actual coils, are produced by field-line tracing, as shown in Fig.8. It’s clear that FOCUS can produce reasonably close agreements at the target plasma boundary, even with only 20 coils. The main visible difference is that the $n/m=5/5$ island chain outside the LCFS in the standard configuration is not visual in the Poincar$\acute{\rm e}$ plots of FOCUS coils. The magnetic island content is determined by the rotational-transform profile and the magnetic perturbation spectrum. There are clearly small differences in both between the FOCUS coils and the actual W7-X coils; however, this is to be expected. As the number of coils used in FOCUS is reduced, the ability to exactly reconstruct a given boundary is also reduced.

The “magnetic perturbation” as defined in this paper is ${\bf B}\cdot{\bf n}$ on the given reference boundary, and it is this quantity that is minimized during the optimization. However, realizing that it is only the resonant part of this perturbation that is associated with the formation of magnetic islands, in future work we intend to (i) identify the resonant harmonics of the perturbation spectrum that are more important to control; and (ii) include additional weights in our target function that explicitly penalize the resonant harmonics of the perturbation spectrum. We expect that this will allow greater control of the island content.

This paper is primarily concerned with the introduction of a new coil optimization algorithm, and the W7-X calculation is used for an initial illustration. We have not performed a comprehensive experimental design study of W7-X and we do not wish to elaborate upon this calculation; however, there is one calculation we present that, in a practical sense, quantifies the difference between the actual W7-X coils and the coils set presented above. Ultimately, the magnetic field provided by the coils is required to confine a plasma equilibrium. So, we compute the free-boundary, MHD equilibrium using the widely used VMEC [28] code with zero pressure and zero plasma currents. Fig.9 shows the differences in the nested flux surfaces of the $N_{C}=30$ coils set and the actual coils, together with a comparison of the coil shapes. The results indicate that the two coils sets have comparatively close plasma boundaries.

5.3 Designing helical coils for the LHD

Besides the W7-X device, another existing large stellarator/heliotron machine is the Large Helical Device (LHD) [29], in Japan. Unlike the modular coils used in W7-X, the main coils of LHD are a pair of helical windings. The approximated LHD coils, including two helical windings and three pairs of vertical field coils, are shown in Fig.10. Helical curves on a torus can be represented by the Fourier series precisely:

[TABLE]

where $R$, $r$ are the major, minor radius of the torus and $N=10/2$ for the LHD case.

When fitting the actual LHD coils with $N_{F}=6$ Fourier series, truncating errors are introduced. Consequently, the fitted coils have slightly different geometries from the actual coils, which ultimately results in the relatively large errors for the produced magnetic field ($f_{B}=1.062\times 10^{-1}$). Fig.11 (a) shows that the initial fitted coils have deteriorate magnetic flux surfaces. There are no closed flux surfaces at the target boundary and the magnetic axis is shifted outwards badly. When we use FOCUS to optimize the fitted coils, the ${\bf B}\cdot{\bf n}$ error $f_{B}$ is decreased to $7.112\times 10^{-3}$ at at $\tau=10^{2}$. The Poincar$\acute{\rm e}$ plots in Fig.11(b) indicate that the optimized coils fix the errors and can produce remarkably good flux surfaces that match the target plasma boundary.

6 Conclusions

In this paper, we have described a new method for designing stellarator coils without the winding surface and introduced the numerical tool, the FOCUS code.

Theoretical interpretations show that the winding surface is not essential. Fourier series is chosen to represent the three-dimensional coils. And for the first time, we have analytically derive the relations between coil parameters and the target functions. This offers a better physical guidance for how to minimize the target function and is then used to analytically calculate the first derivatives. The steepest descent minimization algorithm is successfully employed to optimize a well-constructed target function that includes both physical requirements and engineering constraints. Numerical applications have been performed. The convergence studies show that the FOCUS code is quite effective and flexible in finding appropriate coils. Results of the W7-X and LHD cases indicate that the code is also applicable to practical complicated configurations, with both modular coils and helical coils. And by getting rid of the winding surface, the FOCUS code has the potential to find out more feasible coil solutions.

In the future, we will keep constructing more constraints regarding to different needs, like optimizing the resonant perturbation harmonics and other engineering constraints. And also, we intend to implement more advanced algorithms. The steepest descent is not usually the fastest optimization method, and it’s a local minimizer, but it is used here for simplicity of illustration. The capacity of handling configurations that have nonzero plasma currents will also be carried out, such that the code can be used for more general cases, even for designing 3D coils in tokamaks.

7 Acknowledgements

This work was supported by the Chinese Scholarship Council (CSC) with No. 201506340040 and U.S. Department of Energy with under DE-AC02-09CH11466. The authors gracefully acknowledge the fruitful discussions with S. Lazerson, J. Breslau, D. Gates and N. Pomphrey. Special appreciations are expressed to J. Geiger and J. Loizu from IPP for providing the W7-X data, and Y. Suzuki from NIFS for providing the LHD data. The author C.Z. would also like to thank the PPPL Theory Department for hosting his visiting with this work.