Simulation of Pellet Ablation for Tokamak Fuelling

Abstract.

A magnetohydrodynamic numerical model and parallel software for the ablation of cryogenic deuterium pellets in the process of tokamak fuelling has been developed based on the method of front tracking of ITAPS Center. The main features of the model are the explicit tracking of material interfaces, a surface ablation model, a kinetic model for the electron heat flux, a cloud charging and rotation model, and an equation of state accounting for atomic processes in the ablation cloud. The software was used for the  first systematic studies of the pellet ablation rate and properties of the ablation channel in magnetic fields. Simulations revealed new features of the pellet ablation such as strong dependence of the radius of the ablation channel and  ablation rate on the ``warm-up'' time and supersonic spinning of the ablation channel.

Introduction.

The injection of small frozen deuterium-tritium pellets is an experimentally proven method of refuelling tokamaks [1].  Pellet injection is currently seen as the most likely refuelling technique for the International Thermonuclear Experimental Reactor (ITER). In order to evaluate the efficiency of this fuelling method, it is necessary to determine the pellet ablation rate.

The ablation of tokamak pellets and its influence on the tokamak plasma has been studied using several analytical and numerical approaches (see [2] and references therein). The inherent limitation of the previous ablation models has been the the absence of a rigorous inclusion of important details of physics processes in the vicinity of the pellet and in the ablation channel, and insufficient accuracy of numerical models due to extreme change of thermodynamics states on short length scales. The motivation of the present work is to improve both the accuracy of computational models by using the front tracking technology of the Interoperable Technologies for Advanced Petascale Simulations (ITAPS) Center [3], and physics modeling of the interaction of the pellet ablation channel with the tokamak magnetic field. The present work is a continuation of [5] which introduced sharp interface numerical MHD model for the pellet ablation, but omitted effects of charging and rotation of the pellet ablation channel. Both sharp interface numerical techniques and the resolution of complex physics processes in the pellet ablation channel are important not only for the calculation of pellet ablation rates and the fuelling efficiency, but also for the understanding of striation instabilities in tokamaks [6].

Main physics processes

In this section, we briefly summarize the main physics processes associated with the ablation of a cryogenic deuterium pellet in a tokamak magnetic field (Figure 1). Hot electrons travelling along the magnetic field lines hit the pellet surface causing a rapid ablation. A cold, dense, and neutral gas cloud forms around the pellet and shields it from the incoming hot electrons. After the initial stage of ablation, the most important processes determining the ablation rate  occur in the cloud. The cloud away from the pellet heats-up above the dissociation and then the ionization levels, and partially ionized plasma channels along the magnetic field lines. As it was shown in [5], this process can be described by  the system of MHD equations in the low magnetic Reynolds number approximation. The plasma cloud stops the incident plasma ions at the cloud / plasma interface, while the faster incident electrons penetrate the cloud where their flux is partially attenuated depending on the cloud opacity. We employ the kinetic model for the hot electron  - plasma interaction proposed in [9]. The tendency of the background plasma to remain neutral confines the main potential drop to a thin sheath adjacent to both end-faces of the cloud. Inside the cloud, the potential slowly changes along each field line. Since the cloud density and opacity vary radially, the potential inside the cloud varies from field line to field line causing E x B cloud rotation about the symmetry axis. The potential can be explicitly found using kinetic models for hot currents inside the cloud [6]. The fast cloud rotation widens the ablation channel, redistributes the ablated gas, and changes the ablation rate. We believe that it also causes  striation instabilities observed in pellet injection experiments [6].

Numerical method


In this section, we will describe numerical ideas implemented in the FronTier-MHD code. In general, the system of MHD equations in the low magnetic Reynolds number
approximation is a coupled hyperbolic - elliptic system in a geometrically complex moving domain. We have developed numerical algorithms and parallel software for 3D simulations of such a system
\cite{SamDu07} based on the ITAPS front tracking technology \cite{FT_lite}.

The numerical method uses the operator splitting method. We decouple the hyperbolic and elliptic parts of the MHD system for every time step. The mass, momentum, and energy conservation equations are solved first without the electromagnetic terms
(Lorentz force). We use  the front tracking hydro code FronTier with free interface support
for solving the hyperbolic subsystem. The electromagnetic terms are then found, in the general case,
from the solution of the Poisson equation for the electric potential in the conducting medium using
the embedded boundary method, as described in \cite{SamDu07}.
%In the current 2.5D axisymmetric model
%for the pellet ablation, the elliptic step is eliminated as the current
%density in each point of  the ablation cloud is a function of the hydrodynamic state and some integrated
%quantities along the corresponding magnetic field line. In our current work in progress on 3D pellet
%ablation, solving of a modified 3D Poisson problem is required.
At the end of the time step, the fluid states are integrated along every grid line in the longitudinal
direction in order to obtain the electron heat deposition and internal hot currents.
The heat deposition changes the internal energy and temperature of fluid states, and therefore
the electrical conductivity. The Lorentz force and the centrifugal force are then added to the momentum equation.

FronTier represents interfaces as lower dimensional meshes moving through a volume filling grid (\Fig{FT_structures}).
The traditional volume filling finite difference grid supports smooth solutions located
in the region between interfaces. The location of the discontinuity and the jump in the solution
variables are defined on the interface. A computational stencil is constructed at every interface
point in the normal and tangential direction, and stencil states are obtained through interpolation.
Then Euler equations, projected on the normal and tangential directions, are solved. The normal
propagation of an interface point employs a predictor - corrector technique. We solve the Riemann problem for
left and right interface states to predict the location and states of the interface at the next time step.
Then a corrector technique is employed which accounts for fluid gradients on both sides  of the interface.
Namely, we trace back characteristics from the predicted new interface location and solve Euler
equations along characteristics. After the propagation of the
interface points, the new interface is checked for consistency of intersections. The untangling of the
interface at this stage consists in removing unphysical intersections, and rebuilding a topologically
correct interface. The update of interior states using a second order conservative scheme is performed
in the next step. The tracked interface allows us to avoid the integration across large discontinuities
of fluid states, and thus eliminate the numerical diffusion. The FronTier code has been used for
large scale simulations on various platforms including the
IBM BlueGene supercomputer New York Blue located at Brookhaven National Laboratory.

References

\bibitem{Baylor00}
[1] L. Baylor et al., Improved core fueling with high field pellet injection in the DIII-D tokamak, Phys. Plasmas, {\bf 7} (2000), 1878-1885.

\bibitem{Pegourie}
[2] Pegourie B Review: Pellet injection experiments and modelling, Plasma Phys. Control. Fusion 49 (2007) R87- R160. 

\bibitem{ITAPS}
[3] http://www.tstt-scidac.org/

\bibitem{SamDu07}
[4] R. Samulyak, J. Du, J. Glimm, Z. Xu, A numerical algorithm for MHD of free surface flows at low magnetic Reynolds numbers, J. Comp. Phys., 226 (2007), 1532 - 1549.

\bibitem{SamLuParks07}
[5] R. Samulyak, T. Lu, P. Parks, A magnetohydrodynamic simulation of pellet ablation in the electrostatic approximation, Nuclear Fusion, 47 (2007), 103-118.

\bibitem{Parks96}
[6] P. Parks, Theory of pellet cloud oscillation striations, Plasma. Phys. Control. Fusion, 38 (1996) 571 - 591.

\bibitem{FT_lite}
[7] J. Du, B. Fix, J. Glimm, X. Li, Y. Li, L. Wu, A Simple Package for Front Tracking, J. Comp. Phys., 213, 613–628, 2006.

\bibitem{NGS}
[8] P. Parks, R. Turnbull,  Effect of transonic flow in the ablation cloud on the lifetime of a solid hydrogen pellet in plasma,  Phys. Fluids, 21 (1978), 1735 - 1741.

\bibitem{Ishizaki}
[9] R. Ishizaki, P. Parks, N. Nakajima, M. Okamoto, Two-dimensional simulation of pellet ablation with atomic processes, Phys. Plasmas, 11 (2004), 4064 - 4080.