Particle swarm optimization of divertor targets for heat load control

TL;DR

This paper presents a method using particle swarm optimization combined with FLARE-based heat load approximation to optimize divertor target geometry in fusion devices, considering plasma assumptions and heat load constraints.

Contribution

It introduces a novel integration of PSO with FLARE for divertor design optimization in fusion reactors.

Findings

Optimal divertor configurations depend on plasma assumptions.

The method effectively minimizes heat loads on divertor components.

Results demonstrate the importance of background plasma modeling.

Abstract

Divertor targets in magnetic confinement fusion devices must be designed to handle extreme heat loads. Fast approximation of heat loads with FLARE based on field line reconstruction from an unstructured flux tube mesh is utilized in particle swarm optimization (PSO) of the divertor target geometry. Optimization of the outer divertor target in ITER is evaluated with a constraint for the head loads onto baffles. The optimal configuration is found to depend on assumptions for the background plasma in the heat load proxy simulation.

Figures

Equations

$${\mathbf{v}}_i(t+1)=w{\mathbf{v}}_i(t)+c_1{\mathbf{r}}_1\odot ({\mathbf{p}}_i-{\mathbf{x}}_i(t))+c_2{\mathbf{r}}_2\odot (\mathbf{g}-{\mathbf{x}}_i(t)).$$

$${\mathbf{x}}_i(t+1)={\mathbf{x}}_i(t)+{\mathbf{v}}_i(t+1)$$

$$\mathbf{q}=-{\kappa }_{\parallel }{\nabla }_{\parallel }T-{\chi }_{\perp }n{\nabla }_{\perp }T$$

$$q_t=\gamma T_t{\Gamma }_t$$

$$f(\alpha ,s)=q_{\text{max}}^{\text{(OT)}}+f_{\text{penalty}}\left(q_{\text{max}}^{\text{(baffle)}}+q_{\text{max}}^{\text{(floor)}}\right),f_{\text{penalty}}=10$$

$${\Delta }_{\mathbf{g}}^{\ast (k)}=\Vert ({\mathbf{g}}_k-{\mathbf{g}}_{k-1})\oslash {\mathbf{v}}_{\text{max}}\Vert$$

$$r=\Vert \sigma (\mathbf{x})\oslash {\mathbf{v}}_{\text{max}}\Vert$$

Full text

Keywords: divertor targets, heat load approximation, particle swarm optimization

1 Introduction

Divertor targets are a key component in magnetic confinement fusion devices - especially in those aiming for burning plasmas [1, 2]. Power exhaust from the plasma can result in extreme heat loads, and divertor targets need to be designed such that they can tolerate those. Advanced materials can push the acceptable limit [3], while driving the plasma in front of the divertor targets into a “detached” state [4, 5, 6] can significantly reduce potential heat loads. Plasma detachment is typically achieved by gas puffing with seeded impurities (N, Ne, Ar, …) into the divertor volume for power dissipation.

An advantage of stellarator configurations over tokamaks is their potential for continuous operation and improved stability, but this comes at the price of a significantly more complex design. A particular challenge for stellarators is the toroidal localization of heat loads and the ability to accumulate neutral particles in the divertor volume for efficient removal [7]. High-fidelity numerical models for plasma detachment (such as EMC3-EIRENE [8]) can guide the design process. However, they are computationally expensive and are more suitable for design verification from scoping studies based on fast, low-fidelity models. Even though the latter can not (yet?) capture plasma detachment, they can be leveraged in an optimization loop for minimizing potential heat loads. While it may be sufficient to reduce the actual heat loads below the limit determined by the material properties (with an additional margin for model uncertainties), minimizing the potential heat loads (i.e. without taking into account power dissipation in the plasma) implies that fewer seeded impurities are required for detachment.

Reconstruction of magnetic field lines from a flux tube mesh provides considerable speedup [9]. It is exploited EMC3 for high-fidelity modeling of plasma detchment as well as for low-fidelity heat load approximation. Recently, an implementation based on an unstructured mesh layout has been integrated into FLARE which is suitable for fast heat load approximation in highly shaped island divertor and non-resonant divertor configurations [10]. Once the flux tube mesh is constructed, divertor targets can be adjusted as needed within a given magnetic casing (i.e. the first wall, or the boundary of available space for shape optimization).

The present analysis explores the optimization of the divertor target geometry for heat load control. A brief introduction of the particle swarm optimization (PSO) method is given in section 2. The PSO method is chosen here because it is easy to implement, does not require information about the gradient of the problem being optimized, and has been empirically shown to perform well on many optimization problems. Then, in section 3, the specific optimization problem for this analysis is introduced. Here, rather than jumping into the deep end of stellarator divertor design, the focus is on a toy problem in order to highlight the implications of different parameters. Specifically, we look at the outer divertor target in the ITER and leverage the axisymmetry of the tokamak configuration in order to reduce the number of free parameters (which is not to say that ITER is a toy or that tokamak divertors are trivial). The purpose of this simplification is to analyze the potential of PSO for divertor heat load control, which is covered in section 4.

2 Particle swarm optimization

Particle swarm optimization (PSO) is a stochastic optimization method inspired by social behavior of bird flocking or fish schooling [11, 12]. In PSO, a population of candidate solutions (referred to as particles) are iteratively improved by moving around through search space. The movement is influenced by each particles best known position but also guided towards the global best known position (a variation of the algorithm is instead based on the local neighborhood best). Each particle $i$ has the following characteristics:

[TABLE]

The personal best known position ${\bf p}_{i}$ is either the current position ${\bf x}_{i}$ or one of its previous values - whichever yields the highest fitness value for that particle. For minimization of a multi-variate function $f({\bf x})$, higher fitness is determined by a lower function value. Furthermore, let ${\bf g}$ be the global best known position of the entire swarm of particles. Then, in the next iteration ($t+1$) of the swarm, the velocity of particle $i$ is updated as follows:

[TABLE]

The parameter $w$ in the first term is referred to as the inertia weight [13, 14]. The parameters $c_{1}$ and $c_{2}$ are often called cognitive coefficient and social coefficient. One particular combination that will be used later is $w=0.729$ and $c_{1}=c_{2}=1.494$ [14], but other combinations are possible (within certain limits for stability [15, 16]). The second and third terms guide the particle towards its personal best and global best known positions with uniform random weights ${\bf r}_{1}$ and ${\bf r}_{2}$ which are applied element wise ($\odot$). The position update is then:

[TABLE]

Particles are initialized with random locations and velocity. The best known positions ${\bf p}_{i}$ and ${\bf g}$ are updated when applicable, but not every iteration yields an improvement. Velocities are clampled to $\pm{\bf v}_{\textnormal{max}}$ but with a very generous ${\bf v}_{\textnormal{max}}\,=\,{\bf b}_{u}\,-\,{\bf b}_{l}$ where ${\bf b}_{u}$ and ${\bf b}_{l}$ are the upper and lower bounds of the problem space, respectively. Particles are allowed to fly beyond those boundaries, but $f({\bf x})$ is not evaluated there.

3 The head load optimization problem

Heat loads along the divertor targets and first wall are approximated from the solution of the linearized heat conduction equation $\nabla\,\cdot\,{\bf q}\,=\,0$ with

[TABLE]

where $\kappa_{\parallel}$, $\chi_{\perp}$ and $n$ are considered to be fixed model parameters. More precisely, the thermal conductivy along field lines $\kappa_{\parallel}\,=\,\kappa_{0}\,T^{5/2}$ is based on the physical constant $\kappa_{0}$ and an input background temperature $T$. On the target (here: on any surface), the Bohm boundary condition

[TABLE]

is applied with $\Gamma_{t}\,=\,n_{t}\,c_{st}$ and sound speed $c_{st}\,=\,c_{s}(T_{t})$. The sheath heat transmission coefficient $\gamma\,=\,\gamma_{e}\,+\,\gamma_{i}\,\approx 7$ is applied. We introduce $\chi_{\parallel}\,=\,\kappa_{\parallel}/n$ for convenience in order to relate the impact of cross-field transport to transport along field lines. The total power into the SOL $P_{\mathrm{SOL}}$ is a scaling parameter to which resulting heat loads can be normalized. This is the same model as in EMC3-Lite [17], but adapted in FLARE for an unstructured flux tube mesh. Both implementations are based on a Monte Carlo method which relies on fast reconstruction of field lines. Monte Carlo particles are generated along the inner mesh boundary (typically located just inside the last closed flux surface) and are tracked until they are deposited on the divertor targets or first wall.

An example heat load distribution along the divertor targets and first wall is shown in figure 1 (a,b) where two strike lines can be seen: one on the inner target and one on the outer target. Model parameters are $n\,=\,10^{20}\,\mathrm{m}^{-3}$, $T\,=\,160\,\mathrm{eV}$ and $\chi_{\perp}\,=\,2.0\,\mathrm{m}^{2}\,\mathrm{s}^{-1}$. The simulation has been conducted with $M\,=\,400000$ Monte Carlo particles and took about 1 minute on a desktop computer with 12 CPUs (using MPI for parallelization). It should be noted that there are a number of factors determining the run time: the number of Monte Carlo particles $M$, the level of cross-field transport $\chi_{\perp}$ relative to transport along field lines $\chi_{\parallel}$, the position of the inner mesh boundary (here: at $\psi_{N}=0.99$ normalized poloidal flux), the connection length of field lines in the SOL, and (to some degree) the number of surface patches used to model the divertor targets and first wall. Figure 1 (c) shows the heat load profiles on the outer divertor target for a few different choices of model parameters and $P_{\mathrm{SOL}}=100\,\mathrm{M}\mathrm{W}$ in comparison to an EMC3-EIRENE simulation. The latter includes Ne seeding at a low level of D2 gas puffing and is consistent with the corresponding SOLPS-ITER simulation presented in [18]. It can be seen that the heat load proxy based on upstream values (red) is rather peaked while choosing model parameters based on typical downstream values (green) substantially overestimates the broadening of heat loads. The intermediate case (orange) suggests that a good choice of model parameters should be biased towards upstream conditions.

The ITER divertor includes a dome below the X-point of the magnetic separatrix (as indicated in figure 2. This is important for efficient pumping of netural particles that are released from the divertor targets after recombination of the incident ion and electron fluxes. For the purpose of this optimization study, however, we neglect the dome structure and the transport of neutral particles. Furthermore, we focus on the outer divertor target which we approximate as a plate of fixed length as shown in blue in figure 2. The base point of the plate (black dot) is anchored along the outer divertor leg of the magnetic separatrix and defined by the variable $s\,\in\,[0,1]$ going from the lowest point ($s=0$) up to the X-point ($s=1$). The orientation of the plate is determined by the angle $\alpha$ between the surface normal and the magnetic separatrix: i.e. $\alpha=0\,\deg$ corresponds to a normal incident of the poloidal projection of the separatrix onto the divertor target (note that this is not the field line incident angle - field line incident is at a grazing angle in toroidal direction). Positive values $\alpha\in(0,90)\,\deg$ imply a reflection of neutral particles into the private flux region, as is the case for the actual ITER divertor. Negative values $\alpha\in(-90,0)\,\deg$ imply that neutral particles are reflected outwards into the far SOL. Finally, the divertor volume is closed by a floor connecting the inner and outer targets (indicated in green) and an outer baffle (indicated in orange) connecting to the main chamber first wall.

This is a crude approximation of the actual divertor geometry, but the advantage is that it can be described by two parameters such that the entire configuration space can be sampled. Figure 3 (a) shows the dependence of the peak heat load $q_{\textnormal{max}}^{\textnormal{(OT)}}$ on the outer divertor target on $\alpha$ and $s$, and the black dot indicates the actual divertor setup in ITER. It can be seen that a more grazing incident of the poloidal projection of the magnetic separatrix is more favorable, regardless of the direction in which neutral particles would be reflected (and which are neglected here). This is to be expected since $\alpha$ is linked to the field line incident angle. Furthermore, it can be seen that larger $s$ are more favorable for heat loads than smaller $s$. This is a result of the flux expansion around the X-point. Obviously, configurations very close to $s=1$ should be considered with caution (and be excluded) as the target plate may stick into the plasma or shadow the inner target depending on the sign of $\alpha$.

The peak heat load on the outer divertor target should be viewed in context of the other exposed surfaces: reducing $q_{\textnormal{max}}^{\textnormal{(OT)}}$ below the material limit is necessary, but minimizing $q_{\textnormal{max}}^{\textnormal{(OT)}}$ itself can have unintended side effects. In the following we will consider the objective function

[TABLE]

which penalizes heat loads to the outer baffle and to the floor and is shown in figure 3 (b). The penalty factor is chosen to reflect different material limits, e.g. $10\,\mathrm{M}\mathrm{W}\,\mathrm{m}^{-2}$ for divertor targets and $1\,\mathrm{M}\mathrm{W}\,\mathrm{m}^{-2}$ elsewhere. It can be seen that configurations with $|\alpha|\approx 90\,\deg$ are heavily penalized, and so are configurations with $s\approx 1$. As a result, it is found that there is a main valley of good configurations along the left side of the plot and a secondary valley on the right side towards smaller $s$. The latter includes the actual divertor configuratoin in ITER.

4 Performance analysis

In order to explore the potential of particle swarm optimization (PSO) for divertor heat load control, we take the results from figure 3 and construct a continuous objective function via interpolation (using RegularGridInterpolator from Python’s scipy.interpolate package). Snapshots of a particular swarm’s evolution are highlighted in A where it is shown that a good solution can be found after only a few iterations. As PSO is a stochastic algorithm, we will run simulations for a number of swarms for statistical analysis. Figure 4 shows the results from 10 swarms with 12 particles each. It can be seen in figure 4 (a) that all swarms find a good solution, but not all swarms find the global minimum. Instead, one swarm finds a nearby local minimum as highlighted in figure 4 (b). Nevertheless, most swarms get within $10\,\%$ of the global minimum after only a few iterations, and converge to within $1\,\%$ of their final value (taken after 400 iterations) after 10-40 iterations.

Progress is not always steady: it may take several iterations before a swarm finds a better best known solution. Let us define ${\bf g}_{k}$ as the k-th best known position and $i_{k}$ the iteration number in which it was found. Then the improvement over the last best known position ${\bf g}_{k-1}$ can be given as the $L^{2}$ norm

[TABLE]

where the division ($\oslash$) is taken element wise so that the differences are normalized to the dimensions of the problem space. The improvement $\Delta_{{\bf g}}^{\ast\,(k)}$ at iteration $i_{k}$ is shown in figure 4 (c). The swarms achieve $\Delta_{{\bf g}}^{\ast}\approx 1\,\%$ within 12-40 iterations and $\Delta_{{\bf g}}^{\ast}\approx 0.1\,\%$ within 36-65 iterations. Improvements tend to get smaller with advancing iterations, however, a small $\Delta_{{\bf g}}^{\ast}$ is not a guarantee for convergence. It is entirely possible that a swarm finds a slightly better solution while it is not yet close to a minimum. This happens e.g. to the gray swarm in figure 4 (c) in the 13-th iteration where $\Delta_{{\bf g}}^{\ast}=3.4\cdot 10^{-5}$, but then a better solution is found a few iterations later further away with $\Delta_{{\bf g}}^{\ast}=9\cdot 10^{-3}$. This is perhaps good enough for divertor optimization given the uncertainty of low-fidelity models, but it is worth noting that stopping the swarm once $\Delta_{{\bf g}}^{\ast}$ undercuts a given tolerance may lead to premature convergence.

Another quantity that is related to convergence is the swarm’s radius. Particles will buzz around their personal best known positions and the swarm’s best known position within a volume that tends to shrink around the latter. Thus, the swarm’s normalized radius

[TABLE]

can be an indicator for the swarm’s potential to further explore problem space. Here, $\sigma({\bf x})$ is the standard deviation of the particle positions ${\bf x}_{i}$. One can even go one step further and define the swarm’s energy $E\,=\,\sum_{i}\,\|{\bf v}_{i}\|$ in order to avoid the situation where the swarm is momentarily collapsed to a small volume. However, a weaker - but perhaps sufficient - condition to determine if no further improvement can be expected is the requirement of $r<r_{\textnormal{tol}}$ for two consecutive iterations. Figure 4 (d) shows $r$ remains about an order of magnitude larger than $\Delta_{{\bf g}}^{\ast}$. In particular, it can be seen that one swarm remains at a radius of a few $\%$ for a few 100 iterations before it finally coalesces. At the same time, $\Delta_{{\bf g}}^{\ast}$ is already below $10^{-6}$.

Robustness is evaluated in figure 5 for different swarm sizes. This time, 1000 swarms are simulated for each swarm size, and the fraction of swarms is evaluated that find the global minimum value within $4\,\%$ accuracy for a given number of iterations. It can be seen that a large fraction of swarms reach a good solution within about 20 iterations regardless of the swarm size. Larger swarms tend to find a good solution within fewer iteration and with higher probability, but this is offset by a larger amount of required evaluations of the objective function. Figure 5 (b) shows that larger swarms are not necessarily more efficient. A swarm size of 12 appears to be working well.

As the swarm is guided by its previously known best positions ${\bf p}_{i}$ and ${\bf g}$, one aspect to consider is how noisy values of the objective function can impact the optimization results. Here, the heat load proxy is based on a stochastic method. Under unlucky circumstances, the peak heat load in a simulation may become the new best known value even if the true value is larger. This may guide the swarm to converge to a position that is not even a local minimum. Increasing the number of Monte Carlo particles will mitigate the noise of peak heat load values, but this implies increased computational cost. Thus, it is important to evaluate the performance of PSO in the presence of noise.

The data in figure 3 that went into the interpolated objective function already had some noise in it (this is likely the reason for the other, nearby local minimum found in figure 4 (a,b)). However, evaluation of the objective function was deterministic. In the following we will introduce artificial noise to the problem in order to analyze the performance of PSO. First, though, data values are smoothed with the gaussian_filter from scipy_ndimage in order to remove residual noise. Then, a new objective function $\tilde{f}(\alpha,s)$ is constructed by returning a random sample from a normal distribution with mean $f(\alpha,s)$ and standard deviation $\sigma$. This noisy objective function is then used in another batch of PSO simulations. Results from 1000 simulations are shown in figure 6 for $\sigma/f=4\,\%$ and $16\,\%$. It can be seen that the swarms no longer converge to the global (or local) minimum, but rather converge to points which lie within the noise range of the minimum value $f_{\textnormal{min}}$. Target angles of -60 to -35 deg are found to be optimal at a noise level of $4\,\%$, and solutions are spread over quite a chunk of problem at a noise level of $16\,\%$. It is unclear how this translates to more complex geometries, but it suggests that the potential impact of noise should not be ignored.

Finally, we circle back to the matter of model parameters for the heat load proxy. Figure 1 (c) showed that quite different peak values and profile shapes are found for different approximations of the background plasma. The implications of those choices are revealed in figure 7. It can be seen that a more tangential incident of the magnetic separatrix is preferred for more peaked profiles (red), and that a target location further away is preferred for broader heat load profiles (green). The resulting shape is quite different from the actual setup, but this is likely caused by the focus on heat loads while neglecting requirements for particle exhaust. In any case, the candidate solution(s) found by the present procedure would need to be validated by high-fidelity modeling which includes power dissipation by seeded impurities and recycling of neutral particles.

5 Conclusions

A procedure for the design of divertor targets for heat load control based on particle swarm optimization has been introduced. This relies on fast approximation of heat loads with FLARE enabled by reconstruction of field lines from an unstructured magnetic flux tube mesh. The potential of PSO has been explored based on a simplified divertor geometry for ITER. Even though that it is not guaranteed that PSO finds the global minimum, the procedure appears to be robust enough to find a solution with good enough accuracy. However, if this translates to more complex geometries still needs to be verified and compared to other global optimization methods.

The candidate design is found to depend on the choice of model parameters which are necessary for the approximation of heat loads. It is therefore recommended to generate a few candidate designs for different sets of parameters witin an adequate range, which can then be benchmarked against each other with high-fidelity modeling. An extension of the procedure to multiobjective optimization which accounts for pumping of neutral particles is in preparation.

Acknowledgments

This work was supported by the U.S. Department of Energy under award No. P-240001537.

Appendix A Swarm evolution

The evolution of one instance of a swarm of 12 particles is shown in figure 8. It only takes a few iterations to find a fair approximation of the global minimum, but this is perhaps related to the low dimension of the problem space. It can be seen that a good approximation of the global minimum is found after 20 iterations with further refinement after 50 iterations.