The MH3D++ code has been used to simulate the 3D dynamics of pellet injection. MH3D++  is a modification of the PPPL MH3D resistive full MHD code , which replaces the original finite difference / spectral discretization, with an unstructured mesh and finite element , / spectral discretization. The unstructured mesh is made with triangular and quadrilateral cells. A sample mesh similar to that used in the computations, but with only 1/4 as many mesh points for clarity, is shown in Fig.1.
The MHD equations are discretized with a zero residual Galerkin approach, using piecewise linear basis functions. The differential operators in the equations become sparse matrices involving integrals of the basis functions and their derivatives.
The finite element unstructured mesh discretization has been incorporated into MH3D++ with an object oriented approach. The unstructured mesh objects generate an unstructured mesh and produce the sparse matrices which implement differential operators including gradient, curl, and divergence, as well as various Poisson solvers based on an Incomplete Cholesky Conjugate Gradient solver. Also included are mesh operations needed for line integration and contour plots.
An important feature of this approach is that most of the MH3D code is retained. This allows direct benchmarking of the two versions against each other. Equilibrium and stability calculations using the two versions have been compared, and agree well.
To represent the effect of electron free streaming along magnetic field lines, or large parallel temperature transport, the MH3D code uses the ``artificial sound" method . In this approach, an artificial parallel velocity is introduced, satisfying
The constants are chosen so that the temperature relaxes rapidly to a near equilibrium (5). The temperature diffusivity is small and improves numerical stability. In the absence of damping and dissipation, the system (12), (13) conserves the flux tube averaged pressure as well as an "artificial energy", so that in the presence of damping the system approaches an equilibrium (5). The system (12), (13) models electron free streamimg along field lines, which is the cause of temperature equilibration.
The following computations were done on a poloidal mesh of slightly more than 2000 grid points, and at least 9 toroidal Fourier harmonics in addition to the n=0 mode. The plasma is assumed bounded by a rigid conducting wall.
An initial equilibrium was prepared, starting with a prescribed initial, non equilibrium state, and evolving in 2D, removing kinetic energy, until an equilibrium is approached. In the following three cases, the equilibrium has a rotational transform at the magnetic axis initial peak and the aspect ratio R/a = 5. In addition, the pressure is boosted an additional constant amount with This additional constant pressure enhances the sound speed, while having no effect on the background magnetic equilibrium. This makes the total peak The initial equilibrium density is constant. The equilibrium flux surfaces are shown in Fig.2, and the profile in Fig.3.
The initial equilibrium was perturbed with a density blob representing the ionized pellet ablation cloud. In the following, the blob's peak density is up to 15 times the background density. The blob is cigar shaped, with a circular cross section in the poloidal plane, and a dependence on toroidal angle proportional to The density perturbation extends roughly 1/4 of the way around the torus. The model ablation cloud is considerably less localized than in experiments. This was done because of constraints on numerical resolution, that should be relaxed in future simulations. Similarly the density contrast is less than in experiments. Nonetheless, the effects demonstrated here are expected to be sufficiently robust to qualitatively account for experimental observations.
The temperature was initialized adiabatically, as described in Section II. Before the pellet is injected, the temperature is given by where is the magnetic flux function. After the pellet is inserted, the density is flux surface averaged, as in (8), and used to calculate a new temperature, (9). The new pressure is then obtained from (3). This gives an initial pellet state with a pressure highly nonuniform on flux surfaces, although the flux surface averaged pressure is unchanged. This causes a 3D motion towards a new equilibrium state. Since the profile is the same with and without the pellet, there is a tendency to return to the initial equilibrium state after the pellet cloud spreads out on magnetic field lines. However, when reconnection occurs, the profiles change and the pellet does not spread out to the initial equilibrium.
Simulations were also done with a simpler initialization in which the temperature was not modified by the adiabatic correction. These simulations were quite similar in behavior, as far as the maximum displacement of the pellet cloud is concerned. However, since in this case the initial profile was changed, the pellet had less tendency to return to the initial equilibrium.