Since 18 of December 2019 uses Nucleus credentials. Visit our help pages for information on how to Register and Sign-in using Nucleus.
May 10 – 15, 2021
Virtual Event
Europe/Vienna timezone
The Conference will be held virtually from 10-15 May 2021

Coupling plasma and neutral kinetic models: Considerations and solutions

May 12, 2021, 2:00 PM
4h 45m
Virtual Event

Virtual Event

Regular Poster Magnetic Fusion Theory and Modelling P4 Posters 4


George Wilkie (Princeton Plasma Physics Laboratory)


Divertor detachment is a scenario characterized by the dominance of neutral interactions to mitigate the extreme plasma heat flux that would otherwise be incident upon solid walls of fusion reactors. Despite the critical role theory will play in predicting divertor performance, rigorous modelling of neutrals is plagued by the difficulty of directly solving the nonlinear Boltzmann equation. While several assumptions are appropriate in some contexts to alleviate this difficulty, these break down in the detached regime. Here we present new capabilities in the DEGAS2 neutral transport solver [1] that include a rigorous treatment of nonlinear collision operators. This is a unique capability among comprehensive kinetic simulation models and will be an important component in predicting the performance of advanced divertor scenarios in ITER and future magnetically-confined fusion reactors.

Neutral atoms and molecules, formed by recombination either at the solid wall or within the scrape-off layer plasma, play an important role in advanced divertor concepts. Their interaction with the plasma is purely through reactions, charge exchange, elastic and inelastic scattering. Because of their relatively long mean free path, neutrals are highly non-Maxwellian, as is the plasma they interact with in the scrape-off layer. Therefore, to rigorously simulate the dynamics of plasma and neutrals in the divertor region, kinetic models are needed. XGC is a suite of edge gyrokinetic solvers [2], with early versions fully coupled to DEGAS2. This identified important issues that need to be addressed for high-fidelity simulations of neutral-dominant regimes.

Cooling the plasma with a radiative divertor results in more neutrals generated from plasma recombination. In addition, the gas can become optically dense, where the finite mean free path of photons needs to be accounted for. A cooler plasma further leads to longer lifetimes for hydrogenic molecules. Elastic scattering between neutrals, as well as excitation of internal energy states also become important in these regimes. Beyond this additional complexity, an increased neutral density presents several fundamental computational challenges compared to the low-recycling regime. Firstly, since collisions between neutrals become more important, this necessitates the fully nonlinear Boltzmann collision operator. In DEGAS2, neutral-neutral interactions are currently modelled with an iterative BGK operator, but this is only an approximation and can be improved.

Another complication is the lack of conservation of energy and momentum in charge-exchange collisions when the non-Maxwellian nature of ions and neutrals is not respected. Although both XGC and DEGAS2 are fully kinetic solvers, the field particle species in collisions have historically been modelled as equivalent Maxwellians. Previous work identified this approximation as the cause of the lack of energy conservation when the physical charge exchange cross section is used [3]. Conservation can be recovered if a constant collision kernel is used instead [4]. If this approximation is to be relaxed for increased accuracy, more information about the neutral and ion distributions need to be accounted for in calculating the respective collision operators charge exchange, among other interactions.

One approach is to use DEGAS2 to directly calculate Monte Carlo estimates of the ion collision operator against neutrals. The charge-exchange collision operator is an integral operator, so for each spatial grid point, this requires $N_v^2$ Monte Carlo estimates, where $N_v$ is the number of velocity space grid points used for collisions in the XGC-1 (typically $N_v \sim 500$). This is prohibitively expensive and is indicative of the challenge in numerically solving even the linear Boltzmann collision operator.

In recent years, applied mathematicians have risen to this challenge and several methods have been developed to efficiently solve the Boltzmann equation. This work builds on one such method: the conservative spectral scheme of Gamba & Rjasanow [5]. The distribution function is expanded in an orthonormal Burnett basis:

$$ f\left(v, \theta, \phi\right) \approx \sum\limits_{i=1}^{N} f_{i} A_{k_i l_i} \left(\frac{v}{v_\mathrm{ref}} \right)^{l_i} e^{-v^2/2 v_\mathrm{ref}^2} L_{k_i}^{l_i+1/2}\left(\frac{v^2}{v_\mathrm{ref}^2} \right) Y_{l_im_i} \left(\theta, \phi \right), $$ where $i$ is a compound index encoding the triplet $(k_i, l_i, m_i)$, $L_k^{l+1/2}$ are generalized Laguerre polynomials, $Y_{lm}$ are spherical harmonics, $A_{kl}$ are normalization factors, and $v_\mathrm{ref}$ is a characteristic speed scale of the basis. The weak form of the Boltzmann equation is solved with test functions chosen to *manifestly* conserve collisional invariants. The discrete collision operator is a triply-indexed object of size $N^3$ where each element is an 8-dimensional integral over two velocities and two scattering angles. These integrals are precomputed with a combination of generalized Gaussian quadrature [6] in energy and Lebedev quadrature in solid angle, and these are stored in an online database for efficient application. A prototype framework for solving the Boltzmann equation with this method has been developed [7]. In the context of neutrals in divertors, the spectral scheme will be added to supplement the Monte Carlo scattering operators in DEGAS2. This is most usefully viewed as a *moment method*, wherein the collisional response of several moments (the coefficients $f_{i}$) are robustly calculated and evolved. For neutrals, these moments are calculated with Monte Carlo estimators in DEGAS2, the spectral collision operator finds their respective time derivatives, and a projection onto the velocity space grid provides a conservative collision operator for use in the total-f framework of XGC-1 [8]. It is found that Monte Carlo estimates of higher-order moments becomes difficult beyond $N \sim 27$ due to sampling noise. With such a few number of moments, the spectral scheme is remarkably efficient: requiring only about 0.1 ms per timestep of the nonlinear Boltzmann operator. As such, rigorous nonlinear neutral-neutral scattering can be included in the XGC1-DEGAS2 coupling with negligible computational overhead.

[1] D. P. Stotler and C. Karney. "Neutral Gas Transport Modeling with DEGAS 2". Contrib. Plasma. Phys., 34:392 (1994)
[2] C. S. Chang and S. Ku. "Spontaneous rotation sources in a quiescent tokamak edge plasma". Physics of Plasmas, 15:062510 (2008)
[3] D. P. Stotler et al. "Energy conservation tests of a coupled kinetic plasma–kinetic neutral transport code". Comput. Sci. Disc., 6:015006 (2013)
[4] R. M. Churchill et al. "Kinetic simulations of scrape-off layer physics in the DIII-D tokamak". Nuc. Mat. and Energy, 12:978 (2017)
[5] I. M. Gamba and S. Rjasanow. "Galerkin–Petrov approach for the Boltzmann equation". Journal of Computational Physics, 366:341 (2018)
[6] G. J. Wilkie. Microturbulent transport of non-Maxwellian alpha particles. PhD thesis, University of Maryland (2015)
[7] G. J. Wilkie. "Lightningboltz: A distributed framework for efficient solution of the boltzmann equation". In preparation.
[8] S. Ku et al. "A new hybrid-Lagrangian numerical scheme for gyrokinetic simulation of tokamak edge plasma". Journal of Computational Physics, 315:467 (2016)

Affiliation Princeton Plasma Physics Laboratory
Country or International Organization United States

Primary author

George Wilkie (Princeton Plasma Physics Laboratory)


Daren Stotler (Princeton Plasma Physics Laboratory) Randy Churchill (Princeton Plasma Physics Laboratory) Choong-Seock Chang (Princeton Plasma Physics Laboratory and KAIST)

Presentation materials