Analytical Beamlet Code 3D for neutral beam injectors: principles and applications

Neutral beam injection is the main heating system on a variety of fusion devices, and will be the main heating system of ITER. Especially during high-power operation in long pulse devices, it is important that transmission losses and power deposition in the beamline are well quantified. This paper presents the analytical beamlet code (ABC)3D. ABC3D is based on analytical Gaussian beamlets, and includes shadowing effects. The code is able to calculate emission spectra in the beamline, duct transmission losses, and neutral beam power loads on 3D components. Examples of its application are shown that highlight key results in beam emission spectroscopy, beamline design, and geometry verification in experiments. Limitations to the computational approach are discussed on the basis of a comparison with more advanced calculations.


Introduction
Neutral beam injection (NBI) has been successful in producing high temperature, high fusion performance plasmas in all of the large tokamaks operating worldwide [1]. It is also the main heating system of the LHD stellarator, and will be the heating system with the highest peak power on stellarator W7-X [2,3]. The physical process of beam transport, ionisation, and energy transfer is well understood. However, the technical implementation of a high-power, long pulse NBI system remains challenging.
In NBI sources, positive or negative ions are extracted from a plasma by a grid system with many small apertures, that each produce a beamlet. The beamlets, each with a certain divergence and direction, merge into one beam. A crucial component for any NBI system is the beam duct, which the beam transits from the sources to the fusion plasma. Design efforts on existing machines have led to increasingly complex 3D wall geometries, optimized to minimize local power loads while satisfying geometry requirements imposed by e.g. the magnetic field coils [4]. In a fusion power plant, space constraints on the beam duct will be even more stringent due to shielding geometry and breeding blanket considerations. Calculations of the transmission and local heat loads are important to determine the effect of design choices, and to optimize the beamlet directions.
This paper presents the analytical beamlet code (ABC) 3D, which is based on analytical Gaussian beamlets. Shadowing effects, i.e. components shielding power loads from other components, are included with library routines normally used in 3D computer graphics. ABC3D is able to calculate power loads on and transmission through a 3D duct geometry. The code directly reads exported CAD geometry files without the need for further data preparation, so that a batch of many different geometries can be calculated without added effort. Power loads on a 3D geometry can be visualized and exported for further modeling. Directly importing CAD geometries, calculating shadowing, power density projection on surfaces, subsequent 3D visualization and export for further modeling, are all features not present in the previously used in-house code DENSB [5]. Beam emission spectra can be calculated from the power density distribution and the viewing geometry given a few simplifying assumptions. This useful to determine beam parameters such as the divergence from measured beam emission spectra. The computational approach neglects background gas interaction and the effects of electromagnetic fields, which requires a different modelling approach based on particle tracking. One of the benefits of the analytical approach is that the drastically reduced computation time enables the calculation of detailed power density profiles, for example to test a range of geometries.
In section 2 the computational method is described. Since the code was initially developed for ASDEX-Upgrade and W7-X, section 3 briefly introduces these experiments, for which example applications are shown in section 4 to illustrate the versatility of the code. Beam spectra are calculated to determine beam parameters. Power loads on a duct geometry are calculated to illustrate the use in leading edge detection. Measured infrared beam footprints are compared with calculated heat loads to verify the beam geometry. Power loads on an existing duct geometry are calculated to show the limitations of the computational approach.

Power density calculation
Since the extraction grids have multiple apertures, the neutral beam consists of a sum of beamlets. The power density in each of the beamlets is assumed to have a single Gaussian distribution, which is a common assumption used in extensively validated codes [6,7]. Double Gaussian distributions, such as expected for ITER, and Lorentzian distributions, are implemented but not used in the shown examples [8]. The beamlets are parametrized by the divergence ò, which is half the opening angle of a cone with power densities that are 1/e of the on-axis power density. The normalized power density for a single beamlet as function of parallel (l P ) and perpendicular (l ⊥ ) distance with respect to the beamlet axis is In the experiment, reionisation due to background gas interaction leads to a decrease of the neutral beam power along l P . Since these losses are typically in the percent range, they are negligible with respect to the total beam power, and not included here. The total power density is the sum of the individual beamlets. This part of the calculation is implemented in a Python script which requires as input values: the beamlet origins, directions (which are not necessarily normal to the grids due to e.g. aperture offset steering), the divergence, total power, and the diagnostic locations. Optionally, a separate divergence and power can be specified for each individual beamlet. Additionally, pre-calculated shadowing information can be passed to this routine.
Beamline components are defined as surfaces, read from STL files. The name is an abbreviation of 'stereolithography', the format was invented in 1987 and has remained relatively unchanged since the first release. It describes triangular cells by specifying the vertices and surface normals, and can be generated from 3D models by all major CAD software packages.
To determine which beamlets are shadowed by the geometry for which diagnostic points, the lines of sight from all beamlets to all diagnostic points are checked for intersections. If there are triangle intersections present on the line from beamlet to the diagnostic point, that beamlet is considered shadowed and will not contribute power density to that diagnostic point. The highly optimized modified binary space partitioning tree algorithm is used, which is implemented in the Visualisation ToolKit (VTK) [9,10]. The shadowing calculation can be CPU intensive for geometries with many triangles and diagnostic points, but needs to be done only once per geometry. As a consequence, it is very efficient to optimize the transmission of a beamline by varying only the directions of the beamlets, since the shadowing does not change as long as the beamlet origins and the geometry elements stay the same. Power loads on surfaces are calculated by projecting the power density vector in the midpoint of the triangle onto the surface normal. The results are visualized with VTK [10]. Export of power loads for further modeling is implemented with VTK as well [11].
The analytical modeling approach outperforms particle tracking by orders of magnitude in terms of runtime when considering detailed heat loads in an NBI duct. Because a typical duct transmission is 90%, a naive particle tracking approach wastes 90% of the CPU time, since most of the tracked particles do not interact with the geometry. Resolving small details in the power density needs an extremely high number of tracked particles, since for a surface histogram, the optimum number of samples scales inversely with the bin width to the fourth power [12]. In the analytical approach, the number of operations scales linearly with the number of bins, i.e. inversely quadratic with the bin width for a surface. The analytical modeling approach enables detailed power load calculations on the order of seconds.

Beam emissions spectra calculation
Emission spectra are used to evaluate neutral beam characteristics by comparing them with calculated spectra. The wavelength of the observed light will be Doppler shifted due to the angle between the vector r ij  from beamlet i to observation volume j, and the line of sight of the optical system e LOS  according to where v is the speed of the emitter, c is the speed of light, and λ 0 is the wavelength of the emission line. The spread in particle velocity v is neglected here. If other broadening mechanisms are negligible, the measured spectral profile of the emission line is then only a function of the angular distribution of the beam particles and the instrument function of the spectrometer. Equation (3) describes the wavelength contribution of a single beamlet i to observation volume j. To calculate a spectrum, wavelengths and intensities at each observation volume for every beamlet are calculated. The intensity is given by the volume of the cell, multiplied by the power density at the center of the observation volume. A histogram is made from the wavelengths and intensities, which results in a synthetic spectrum. This spectrum needs to be convolved with the instrument function of the spectrometer in order to compare it with measured spectra. Note that variations in the background gas density and excited state fraction along the beamlet path are neglected here, but could be taken into account when known.
The calculated spectra are scaled and allowed a small shift to overlap with a measured spectrum. To estimate the relative goodness of fit, the coefficient of determination R 2 is calculated from the measured data d(λ) and the model contains the input parameters of the model, and N is the number of data points in the measured spectrum. An R 2 of 1 indicates that the model predictions perfectly fit the data.

Experimental layout
As ABC3D was originally developed to model IPP's neutral beam injectors, those systems are briefly described here. At the ASDEX Upgrade tokamak, there are two neutral beam boxes. Both beam boxes are equipped with 4 NBI sources, with 2.5 MW beam power arriving in the plasma each [13]. At the Wendelstein 7-X stellarator, there are two neutral beam boxes. In the experimental phase OP1.2b, one box is equipped with two sources [14].
Neutral beam injection at the ASDEX Upgrade tokamak and the W7-X stellarator make use of essentially the same beam sources. In each NBI source, protium or deuterium atomic and molecular ions (H + /D + , + + H D 2 2 , and + + H D 3 3 ), are extracted from a plasma by means of 3 grids with 774 apertures, where the second grid prevents back-acceleration of electrons [15]. In each source, the two grid halves are tilted vertically by 0.87°to focus the beam halves at 8.5 m from the source. In addition, the holes in the plasma facing grid halves have an offset that increases linearly with distance from the grid half center, to induce an inward steering of the outer beamlets. To correctly model the geometry, the single beamlets need to be taken into account.
Directly after the source the extracted ions interact with hydrogen background gas in a tapered pipe, the neutralizer, leading to dissociation of the molecular ions and neutralization. In addition, background gas interaction leads to excitation and light emission. The Balmer-α line, with an average air wavelength of 656.093 nm for deuterium, is used to diagnose the beam [16]. The n=3 excited state fraction, the upper level for the Balmer-α transition, is constant after approximately the first 10 cm into the neutralizer [17]. The emission spectrum is collected along a vertical line of sight with an angle of 50°to the source normal, at 1.5 m distance from the source. The light is coupled into a fiber by a 135 mm lens. For ASDEX-Upgrade and W7-X, it is a good approximation to treat the viewing system as perfectly collimated, i.e. e LOS  is identical for all the observation volumes. As previously shown, fine structure of the Balmer-α line due to the Lorentz, Stark, or Zeeman splitting is negligible in the ASDEX-Upgrade NBI beamlines [18]. A schematic view of the emission spectroscopy with the middle and outermost beamlets for each grid half is shown in figure 1.
The neutrals produced by the different extracted atomic/ molecular ions have a different speed, and thus Doppler shift. The intensity ratio between the different peaks can be used to deduce the species distribution over these energy components [18]. The spectral shape of the emission peak is used to determine the angular distribution of the beam.

Beam emission spectroscopy
Beam emission spectra, collected along the line of sight in the neutralizer, contain information about the angular distribution of the beam particles. By comparing measured and calculated spectra, the angular distribution of the individual beamlets can be deduced. In figure 2, the main energy component of the spectrum is shown, which is measured at ASDEX Upgrade at an extraction potential of 93 kV.
Several effects contribute to the shape of the measured spectrum: the angles between the individual beamlets and the collection geometry, the divergence of the individual beamlets, and the instrument function of the spectrometer. The two tilted grid halves combined with the beamlet steering and measurement location lead to the double peak structure.
Synthetic spectra are calculated as function of two input parameters, the beamlet steering factor δ which quantifies the amount of beamlet steering in degrees per mm offset of the apertures, and the divergence ò. The non-relativistic version of Formula 3 is used, since D at 93 keV moves at less than 1 percent of the speed of light. The synthetic spectra are convolved with the instrument function, and fit to the measured spectra by scaling the amplitude and allowing a wavelength shift on the order of 10 pm. The scaling is necessary because the optical system is not absolutely calibrated. The wavelength shift is needed because of submillimeter differences in the fiber alignment between the lines of sight.
To estimate the relative goodness of the different fits shown in figure 2, the R 2 for several combinations of beamlet steering and divergence is shown in figure 3. In general the fit quality is quite good, which shows that for this positive ion extraction system, a single Gaussian divergence is well justified. The overall quality of the fit is equally good for several combinations of the beamlet steering factor and the divergence, because both parameters factor into the broadening of the spectrum. As shown in figure 2, it is not possible to determine which of these synthetic spectra is best, an additional constraint on either the beamlet steering factor or divergence is needed. Since the divergence results from a combination of beamlet formation and space charge compensation which is hard to quantify, the beamlet steering factor was taken from IBSimu ion optics calculations [19,20]. With the calculated beamlet steering of 1.44°per mm, the divergence of the full energy component is 0.74 degrees, with error bars on the order of 0.1°. The half and third energy components have divergences within the error bar of the main energy components, but with larger uncertainties due to worse signal to noise ratios. The measured values compare well with the previously determined value of 0.8° [18]. With R 2 as goodness of fit indicator, the interdependencies between the different fit factors are clear.

Leading edge detection for duct optimisation
The amount of space available for an NBI duct on a fusion reactor is limited, since it has to fit between the field coils. In advanced stellarators such as W7-X this is especially challenging due to the complicated coil structure. In narrow sections, the duct components should be at an oblique angle to the beams to spread the scraped power over a large area while avoiding leading edges. In the stellarator Wendelstein 7-X, each beam duct contains two constrictions due to the field coils [4,5]. To avoid excessive power loads at these locations, the beams are gradually scraped by an inertially cooled copper bellows scraper, which has this name since it protects a bellows, and duct inner liner made of inertially cooled CFC tiles clamped onto CuCrZr heat sinks. It is crucial that the beam is not scraped too much, and that the power is spread evenly so that the local power load is not too high. Figure 4 shows the original duct geometry that was under consideration. The view is in beam direction, and shows the box exit scraper (beige), bellows exit scraper (orange/red), and duct inner liner (gray).
A calculation of the power load serves as a final check to see if the maximum power loads are not exceeding the material limits. A divergence of 1°is assumed in these calculations instead of the expected 0.8°to have a safety margin, as previously done in W7-X duct power load and transmission calculations [5]. All 4 sources produce beams in the calculation; a total beam power of 10 MW is assumed. Figure 5 shows the calculated power load for the whole geometry on the right, and a cutout of the duct inner liner tiles on the left. An excessive power load due to a leading edge at a carbon tile at the left side of the duct is clearly visible. Locally the power load at this tile exceeds 20 MW m −2 , which will likely lead to cracking. A chamfered tile was then used instead, which spreads the power more evenly. The power loads on the copper shield, and the duct constriction at the right are within material limits for 10 s pulses. Because of the interplay between source geometry and beam divergence, it is difficult to rigorously identify leading edges based on a visual inspection of the model. With a power load calculation, determination of leading edges is straightforward.

Use in geometry verification
Due to e.g. mechanical tolerances, there could be small deviations between the design NBI geometry, and the actual geometry in the experiment. With a combination of infrared measurements and power load calculations, the differences between design and experimental geometry were determined in ASDEX Upgrade. Infrared photon fluxes were measured with the ASDEX Upgrade infrared system during short NBI shots into the empty torus [21]. Only the footprints of the box 2 sources are in view of the camera. The diagnostic view was mapped onto the 3D construction design with the visualization software AUGpy, so that the measured IR footprints could be directly compared with calculated power density profiles [22,23]. Figure 6 shows an IR measurement overlaid with 3D models of the vessel components, together with a calculated power density profile seen from the same perspective. Many power density profiles were calculated, with varying horizontal and vertical angle of the source, and constant 0.8°divergence, 1.44°/mm steering factor, and 2.5 MW beam power. The closest match was determined by visually comparing the measured IR footprints to the calculated power deposition profiles. Although the beam size is on the order of 25 cm at the wall location, the beam location could be determined with an accuracy of approximately 1 cm. Corrections to the design geometry of at most 0.2°were found, which correspond to 4cm at the wall location. Althought the corrections are relatively small, they have a detectable influence on e.g. the radial fast ion distribution. The measured geometry was implemented in codes that simulate transport and interpret diagnostic output, leading to more consistent results [24,25].

Heat loads in shadowed areas and reionisation
A comparison of calculated and measured power loads at ASDEX Upgrade beam box 2 serves to illustrate the limitations of the computational approach. At the connection of beam box and duct, a box exit scraper is installed to shield components later in the beamline such as the DN800 gate valve. From 1996 to 2017 a D-shaped box exit scraper was used. Melt damage of the gate valve behind the scraper occured several times during this period. With the shadowing algorithm, the fraction of beamlets with a direct line of sight to the surface was calculated. The locations which were damaged in the experiment are rigorously shadowed, i.e. the damage must have come from reionized and subsequently deflected particles. A more in-depth study based on a particle tracking approach, which required additional modeling of the 3D neutral gas density profile and the magnetic field,   confirmed this hypothesis [26]. A new box exit scraper based on the W7-X design has been installed in each ASDEX Upgrade beam box.

Conclusion and outlook
The Analytical Beamlet Code 3D was developed to calculate power loads for 3D components including shadowing effects. The analytical approach combined with highly optimized raycasting routines results in the efficient calculation of detailed power density profiles, which are useful to detect leading edges or interpret measured infrared camera data. With a few simplifying assumptions, beam emission spectra can be calculated as well. 3D data visualisation and data export are implemented with VTK. The ABC3D code supersedes the in-house DENSB code, which is much more limited in terms of geometry importing and data visualization and exporting [5].
To illustrate the versatility of ABC3D, several example applications were shown which all require a substantial number of relatively detailed power density profiles. Limitations of the computational approach were illustrated with an example in which damage in the experiment was caused by reionization and subsequent deflection, which are not included in the calculation.
The computational approach is excellently suited for geometry optimization. This is especially crucial in future fusion reactors such as DEMO, where the NBI port size needs to be minimized in view of the neutron damage and tritium breeding ratio, but needs to be large enough to transmit a large fraction of the beam. Calculations are ongoing to optimize the beamline geometry for possible DEMO designs [27,28].