Two-dimensional simulation of pellet ablation with atomic processes

A two-dimensional hydrodynamic simulation code CAP has been developed in order to investigate the dynamics of hydrogenic pellet ablation in magnetized plasmas throughout their temporal evolution. One of the properties of the code is that it treats the solid-to-gas phase change at the pellet surface without imposing artificial boundary conditions there, as done in previous ablation models. The simulation includes multispecies atomic processes, mainly molecular dissociation and thermal ionization in the ablation flow beyond the pellet, with a kinetic heat flux model. It was found that ionization causes the formation of a quasistationary shock front in the supersonic region of the ablation flow, followed by a ‘‘second’’ sonic surface farther out. Anisotropic heating, due to the directionality of the magnetic field, contributes to a nonuniform ablation ~recoil! pressure distribution over the pellet surface. Since the shear stress can exceed the yield strength of the solid for a sufficiently high heat flux, the solid pellet can be fluidized and flattened into a ‘‘pancake’’ shape while the pellet is ablating and losing mass. The effect of pellet deformation can shorten the pellet lifetime almost 3 from that assuming the pellet remains rigid and stationary during ablation. © 2004 American Institute of Physics. @DOI: 10.1063/1.1769376 #


I. INTRODUCTION
Fueling tokamaks by pellet injection has achieved notable success in the last two decades. 1 Pellet penetration well inside the separatrix has enabled the plasma to surpass the empirical density limit, 2 which is clearly of importance for next step devices such as the International Tokamak Experimental Reactor ͑ITER͒. 3Additionally, peaking the density profile with strong central fueling by pellet injection significantly improves the energy confinement time 4 -6 and increases the fusion power density at fixed ͗␤͘ and B t . 7Conventional gas puffing has also been successful in building up the plasma density, however in stellarator devices such as the large helical device, 8 the particle source from gas puffing is strongly localized at the plasma surface.In a large fusion reactor using tritium rich fuel, gas puffing has to be replaced by pellet injection simply to avoid high edge recycling and excessive tritium retention in the plasma facing components. 3here are many important physics processes in connection with pellet fueling, not the least in the prediction of the pellet ablation rate and penetration distance in the plasma.Pellets ablating in hot plasma are surrounded by a cold, dense neutral ablation cloud which expands in all directions near the pellet, and then after the flow ionizes and loses pressure it becomes constrained to follow the magnetic lines of force.The inner ''inertially dominated'' region, where magnetic forces are negligible, largely controls the ablation rate: being a few pellet radii thick it provides most of the shielding of the pellet from the incident plasma electrons, which supply the heat flux driving ablation.From the quasi-steady transonic-flow ͑TF͒ ablation model 9 several extensions and modifications have been made over the years.In the TF model the actual Maxwellian distribution function of the incident plasma electrons was approximated by an equivalent monoenergetic distribution having the same heat flux and particle density.The model used an equivalent spherically symmetric, or isotropic, heat source to approximate the electronic energy deposition.However, the incident plasma electrons are constrained to follow the straight magnetic field lines in tight helices, so that spherical symmetry is broken.Pellet ablation therefore involves at least a twodimensional ͑2D͒ geometry, as sketched in Fig. 1͑a͒, while in Fig. 1͑b͒ we display the 1D geometry used in the construction of the isotropic heating approximation.Both the monoenergetic and spherically symmetric approximations were partially removed by Kuteev, 10 who used the actual Maxwellian distribution for the incident electrons, and attempted to account for the heating asymmetry imposed by the magnetic field.Although the ablation flow itself was forced to be spherically symmetric, and everywhere sonic ͑Mach number ϭ1͒, the pellet ablation rate was nevertheless in very good agreement with MacAulay's 2D numerical simulation of pellet ablation, where these geometrical effects were treated. 11oth investigators found that the scaling of the ablation rate with the plasma parameters and pellet size agreed fairly well with TF scaling, but the ablation rate was systematically ϳ2ϫ higher.Note that the TF model is well supported by experimentally measured pellet lifetimes, but this is partly the result of a fortuitous cancellation between the two approximations which tend to act in opposite directions. 10,11he question is then why do the improved models of Kuteev and MacAulay have a factor 2 discrepancy with experiments?
One possibility is due to the electrostatic shielding provided by the negatively charged pellet cloud, which would effectively reduce the incident electron density n eϱ and heat flux (ϳn eϱ ) by the Boltzmann factor e Ϫe/T eϱ , where T eϱ is the plasma electron temperature, and the normalized potential drop 12,13 across the cold cloud/hot plasma interface is e/T eϱ ϭ1.8-2.Hence the ''flux-limiter'' f appropriate for pellet ablation is e Ϫ1.8 to e Ϫ2 or 0.135-0.165.The Boltzmann correction to the density and heat flux is only relevant for a Maxwellian distribution for which the normalized energy distribution function is not altered by a potential drop.Since the ablation rate scales with n eϱ 1/3 , the ablation rate would simply be reduced by a factor e Ϫe/3T eϱ ϭ0.51-0.55,which could alone explain the factor of ϳ2 discrepancy.Similarly, the pellet surface pressure ϰn eϱ 2/3 would be reduced by e Ϫ2e/3T eϱ ϭ0.26-0.30.A direct extension of the TF model, including incident Maxwellian electrons, the geometrical effect of their heat deposition, and the above electrostatic correction, predicts for a pellet of radius r p the mass ablation rate 14

Gϭ
1.36ϫ10 Ϫ8 A 2/3 r p where A is the atom mass in amu (Aϭ2.014 for deuterium pellets͒, I * ϭ7.5 eV is the mean excitation energy of hydrogenic atoms, T eϱ and n eϱ are the plasma parameters, and cgs eV units are used.The numerical value of the ablation constant represents a synthesis of refinements and modifica-tions contributed by other published papers, and it includes an order unity multiplier to eliminate the small discrepancy with the experimental ablation rate inferred from the multimachine International Pellet Ablation Database ͑IPADBASE͒ of measured pellet lifetimes. 15Pellet injectors usually produce cubic or cylindrical pellets.Such pellets are likely to become nearly spherical early on because passage through the guide tube and/or ablation in the plasma quickly smoothes out irregularities such as edges and corners.In modeling ablation, the usual practice is to use the radius of an equivalent equal-mass spherical pellet, i.e., r p ϭ(3/16) 1/3 D is used in the case of a right circular cylinder pellet with length L equal to diameter D. This adjusted definition of r p is understood for the pellet data archived in IPADBASE.To convert from G to dN/dt (atoms/s) one can use dN/dtϭN A G/A, where N A is Avogadro's number.
A definitive identification of other important effects on the ablation rate is still incomplete.For example, all of the above models have assumed that the pellet remains rigid and nearly spherical during ablation.Deviations from this assumption can arise because the heat flux and the areal ablation rate are larger on pellet surface elements normal to the magnetic field lines than those that are tangential.As pointed out by Kuteev, 16 the pellet would become progressively more oblate or ''lentil'' shaped.The heat flux asymmetry also imparts a nonuniform ablation recoil pressure at the pellet surface ͑rocket effect͒.Under high heat flux conditions, this pressure variation can exert a shear stress on the pellet which is larger than the yield strength.The solid pellet will then be fluidized and flattened on the surfaces facing the heat flux. 14o gain a fuller understanding of these two-dimensional effects, we developed a time-dependent code CAP in order to rigorously treat the dynamical coupling between interior pellet motions and the exterior ablation flow.The 2D axisymmetric geometry of the problem is illustrated in Fig. 1͑a͒.The CAP code uses an Eulerian approach based on the cubicinterpolated pseudoparticle ͑CIP͒ method. 17,18The CIP numerical scheme has two important features.One is that it is designed to accurately resolve fine-scale temperature and density structure near the interface between two disparate regions, for example, a shock surface.The other is that it can treat the nebulous boundary layer where the solid undergoes transformation to a gas and the pellet surface recedes.In that sense, the CIP method is particularly suitable for the study of pellet ablation involving the interaction between a compressible phase ͑abated material͒ and an incompressible phase ͑the cold pellet͒.This numerical scheme can simultaneously treat both regions in one program, and therefore avoids having to join solutions in two distinct regions by the use of boundary conditions at the interface of separation.
The code has been first benchmarked with the results of the previous 2D code 11 and good agreement was obtained for the ablation rate G. Dissociation and ionization in the ablation flow can also be treated by using the assumption of local thermodynamic equilibrium ͑LTE͒ as in Felber et al., 19 who included atomic processes in the original TF model.In keeping with Ref. 19, the CAP code finds that the energy sinks made available by the latent heats of dissociation and ionization lead to only fractionally lower ͑ϳ10%͒ ablation rates.FIG. 1. ͑Color online͒.͑a͒ Cylindrical axisymmetric coordinate system (r,z) used to study anisotropic heating by incident plasma electrons tied to the magnetic field B ¯ϭBz ˆ. ͑b͒ Spherically symmetric coordinate system r used for the isotropic heating approximation, in which the heat deposition at any r is the same as the heat deposition along the polar axis, ϭ0.
The present work confirms the Mach number depression near the point in the supersonic region of the flow where ionization is nearly completed. 19Reference 19 also found that when the pellet is subjected to high heat fluxes, the Mach number can dip below unity, i.e., the flow becomes singular preventing a continuation of the exterior flow solution beyond that point.We have found that the Mach number never dips below unity.Instead, the ablation flow evolves towards a quasistationary ''double-transonic'' flow solution with a shock discontinuity, and thus the flow singularity encountered in Ref. 19 is resolved.Then, the CAP code reinforces the validity of the quasisteady approach pertaining to the inertially dominated part of the flow domain.
Because the pellet behaves as almost a hydraulic fluid under high shear stress, it can deform into a pancake shape, flattened in the direction of the magnetic field and elongated in the perpendicular direction.As we shall see, such 2D deformations can cut the pellet lifetime by ϳ3-4ϫ compared with the lifetime of a spherical pellet that remains rigid, and ϳ1.5-2ϫ compared with a pellet ablated by an equivalent spherically symmetric ͑isotropic͒ heat source.This large increase in the ablation rate is much greater than was previously predicted by an analytical model. 14ur final goal is to rigorously simulate the smooth 2D transition from the inner nearly isotropic flow region to the outer field-aligned, or channel flow, region by including the magnetic JϫB force in the CAP code.The buildup of an extended ''plasma shield,'' 20 can provide additional shielding against the plasma heat flux directed along the lines of force.The effectiveness of the plasma shield on reducing the ablation rate needs to be fully quantified and compared with other shielding mechanisms.The physical variables of the plasma shield are also needed to accurately model the fast EϫB advection of the ionized pellet fuel across the plasma column, 21,22 and its homogenization within a magnetic flux shell. 23We have not yet included such hydromagnetic effects.One of the main purposes here is to demonstrate the reliability of the CAP code, and compare ablation rates and fluid profiles against results from other models which have also neglected hydromagnetic effects.
In Sec.II, basic fluid equations used in the code are presented, and a semianalytic kinetic heat flux model is also described.In Sec.III, details of the CIP numerical scheme used in the code are presented.In Sec.IV, we make the 1D spherically symmetric approximation for the heat flux deposition in order to isolate 1D and 2D effects.The full twodimensional effects are treated in Sec.V, with special emphasis on deformation of soft pellets by the nonuniform surface pressure distribution, including Kuteev's lentil effect.Conclusions and summary are given in Sec.VI.

A. One-fluid multispecies conservation equations
Due to the high collisionality in the cold, dense ablation cloud, all species: molecules, atoms, ions, and electrons can be treated as one fluid with common temperature T, fluid velocity v ¯, and mass density .Because of the dominance of the collisions over the radiative rates, and the fact that the time scales of collisional excitation and deexcitation are much shorter than the characteristic ablation flow time scale, the LTE condition is assumed to hold at all points in the ablatant.In contrast, the hot collisionless plasma electrons that penetrate the ablatant must be considered as part of a separate ''suprathermal distribution'' with temperature T eϱ .These electrons have very long-mean-free paths proportional to energy squared so that the energy transported by these electrons must be calculated by solving the Fokker-Planck equation.Other than providing the heat source that drives ablation, the external plasma otherwise has minimal effect on ablation.Electromagnetic effects are ignored in the present study: the only effect of the magnetic field is to provide directionality to the incident electron heat flux ͑see Fig. 1͒.In general, the dynamics in the ablation cloud as well as the solid ͑pellet͒ region are given by the following fluid equations of mass, momentum, and energy conservations: where d/dtϭ‫ץ/ץ‬tϩv ¯"" is the advective time derivative, e is the total specific internal energy, and q ¯is the heat flux due to the incident plasma electrons.The pressure tensor P ¯is decomposed into a scalar pressure p and a stress tensor S ¯, namely, P ¯ϭ pI ¯ϩS ¯where I ¯is a unit tensor.Notice that these equations are written in nonconservative form, which makes it possible to treat incompressible and compressible regions simultaneously by use of the CIP method outlined in Sec.III.When a solid is heated, it is transformed to gas ͑mol-ecules͒ by absorbing the sublimation energy (⑀ s ϳ0.01 eV).Subsequently, the molecules are transformed to atoms by absorbing the dissociation energy (⑀ d ϭ4.48 eV), and the atoms are transformed to ions and electrons by absorbing the ionization energy ⑀ i ϭ(1/2)m e (c/137) 2 ϭ13.6 eV.In the onefluid model, p and e are given by where ␥ m and ␥ are specific heat ratios for molecules and atoms, respectively, k is the Boltzmann constant, m is the mass of the atom ͑ion͒, p 0 is the pressure of the solid pellet phase, and the specific internal energy of the solid phase e 0 is negligibly small.Fractions for sublimation f s (T,), dissociation f d (T,), and ionization f i (T,) are defined by in which n t ϵ2n 0 ϩ2n g ϩn a ϩn i ϭ/m stands for the total number density of nuclei, and n 0 , n g , n a , and n i denote, respectively, the number densities of solid D 2 molecules, gas D 2 molecules, D atoms, and D ϩ ions.These three equations constitute a complete set of algebraic equations for the number densities of the four species, given T and .Note that the ideal gas equation of state ͑EOS͒ is assumed for all phases except the solid phase, which we justify shortly.
The fractional ionization can be found from the Saha equilibrium formula 24

͑4͒
in eV ͑cgs͒ units.The Saha relation assumes a mixture of atoms and ions coexisting in equilibrium without the presence of molecules.This is actually the case here since sublimation and dissociation processes are usually completed well before ionization begins ( f s ϭ1,f d Х1) so that n t Хn a ϩn i whenever f i becomes appreciable.Moreover, the Saha equation terminates the electronic partition function summation at the ground state level and thus overestimates the fractional ionization slightly.However, the resulting overestimate in e coming from terms involving f i is compensated by the lack of terms in e contributed by the electronic excitation energy.The excitation level cutoff corresponds to the idea that the electron ''prefers'' to be removed from an atom rather than occupy an excited level; particularly so in the case of a highly dense ablation plasma which effectively limits the actual number of excited levels. 24The dissociation fraction is expressed in analogy to the ionization fraction, namely, 19,24

B. The cryogenic solid-gas transition region near the pellet surface
We propose a somewhat ad hoc function for the sublimation fraction f s , but before we do so, some discussion involving the properties of the hydrogens at low temperatures and moderately high pressures near the pellet surface is needed to justify our approach.The f s function is used simply for computational convenience: it dictates at what temperature the incompressible solid is supposed to be converted to the compressible ideal gas phase, and it does this by mocking up a continuously smooth transition function which joins the two regions, in similitude with the f d and f i functions.In a thin transition layer between the two regions, f s varies from zero ͑all solid phase, with no gas͒ to unity ͑all gas phase, with no solid͒.Of course f d and f i are both zero in the layer.So the relative contribution of the pressure in the layer coming from the gas phase and the solid phase is quantified by the local value of f s .In actuality, however, the layer often consists of a real gas or supercritical fluid, neither ideal gas or solid.Nevertheless, we will now show that this does not matter because the pressure difference across the layer is negligible.Hence, the pressure on the pellet side of the layer p 0sur will be the same as the pressure on the ideal gas side of the layer p ig .
Pellet surface pressures cannot be easily measured, or indirectly inferred from experimental data.However, from theoretical predictions 14 surface pressure magnitudes of interest lie within the range p 0sur ϳ0.2-20 MPa.The upper value here corresponds to plasma temperatures ϳ6 keV for pellet sizes and plasma densities expected in ITER: pellets usually cannot penetrate to distances in the plasma where the temperature is much higher than 6 keV.Now the real gas EOS for hydrogen at cryogenic temperatures can be written as pϭZkT/(2m), where the compressibility factor Z is a measure of the departure from the ideal gas EOS, where Z ϭ1.By restricting ourselves to this pressure range, approximate ideal gas behavior will be valid for TϾϳ60 K.This can be illustrated in Fig. 2͑a͒, where we used the hydrogen EOS data tables 25 to plot Z for a family of isotherms.The initial dip below unity is due to the van der Waals intermolecular attraction forces.At still higher pressures, the ''hardcore'' repulsive forces begin to dominate and Z rises sharply above unity for pϾ10 MPa.The minimum temperature indicated on the plot is just slightly above the critical point temperature for hydrogen.Above the critical point temperature, unique liquid and gas phases do not exist for any pressure, there is only one phase-a supercritical ''fluid''-with an inflection point (‫ץ‬p/‫ץ‬V) T ϭ(‫ץ‬ 2 p/‫ץ‬V 2 ) T ϭ0 at the critical point on a Ϫp phase diagram as sketched in Fig. 2͑b͒.The critical point parameters for deuterium ͑hydrogen͒ are 26 T c ϭ38.3 (33.2) K, and p c ϭ1.66 (1.315) MPa.Often pellet ablation lies in the supercritical regime, p 0sur Ͼ p c ; once again, no liquid-to-gas phase transition exists for temperature increases along a supercritical isobar, as illustrated by the complementary ϪT phase diagram in Fig. 2͑c͒.In that case, one can ask how the sublimation energy ⑀ s is defined, and how is the solid-to-fluid phase transition captured computationally if we confine ourselves to the ideal gas approximation outside the solid phase?
To address these questions, we first note that the transition layer is very thin compared with the radius of the pellet.This is because the temperature on the fluid side of the layer where ideal gas behavior is valid, T ig ϳ60 Kϭ0.0052 eV, is much less than the temperature T * ϳ1 eV at the sonic ͑Mach number M ϭ1) surface of the ablation flow, which lies more than one-half pellet radii away from the pellet surface. 9,14ssuming that temperature in the transition layer changes while the pressure remains constant, as depicted in Fig. 2͑c͒, the (T) dependence within the layer can be determined from the real gas EOS.It then follows that the specific internal energy of the supercritical fluid phase, expressed as e f ͓T,͔, becomes a function of temperature only, e f ͓T,(T)͔.Hence, we can define the heat of transition or effective energy of sublimation per molecule ⑀ s as where 0 is the heat capacity of the solid phase.The precise distinction between the solid and the fluid phase still remains at our typical ablation pressures, i.e., a melt line exists with T m being the melting temperature.The heat of fusion associated with crossing the melt line is h f us ϳ200 J/mol for deuterium. 26Employing the SESAME EOS tables, 25 we obtain ⑀ s Ϸ1260 J/mol, or 0.013 eV/molecule for deuterium.Interestingly enough, these supercritical transition energies are essentially the same as the energy of sublimation at the triple point where the corresponding deuterium pressure is much lower, p tp ϭ0.017 MPa.Thus ⑀ s ϳ0.013 eV may be treated as a constant for all parameters of interest.
Motivated by our simple heuristic discussion on the solid-to-gas phase transition, we suggest the mock-up model

͑6͒
To ensure that f s does not reach unity before T reaches T ig , we shall arbitrarily set T tran ϳ58.7 K and ⌬Tϭ1.163K.The solid-gas interface ͑pellet surface͒ movement can be traced by obtaining the recession velocity, or velocity of the ablation front uϭϪ(‫ץ‬T tran /‫ץ‬t)/ٌ͉T tran ͉.
To justify that a nearly constant pressure exist within the transition region 0Ͻ f s Ͻ1, we consider solutions of the conservation equations ͑1͒ in plane geometry with local outward normal coordinate x.All quantities are assumed time independent in the frame moving with the ablation front.Transforming to the coordinate ϭxϪut, such that ϭ0 coincides with the TϭT tran surface, we have As noted earlier, the solid pellet phase is assumed to be incompressible.Because of the high sound speed in solid deuterium, c s0 ϭ1700 m/s, pressures above 40 MPa are needed to cause significant ͑Ͼ10%͒ density compressions. 26Additionally, the ablation front advances inside the pellet at subsonic velocities (u/c s0 ϳ10 Ϫ3 -10 Ϫ2 ), hence the solid phase remains at rest, except for slow internal motions associated with deformation, which have no consequence in what follows.Accordingly, we may integrate Eqs.͑7a͒ and ͑7b͒ from ϭϪϱ with boundary conditions, v→0, → 0 , and p →p 0sur to ϭ, and obtain a pressure-density relation within the layer

͑8͒
In the language of shock physics, if u were the shock velocity, this relation would correspond to the Rayleigh line.Referring to Fig. 2͑c͒ again, we see that on thermodynamic grounds, ‫͉‪T‬ץ/ץ‬ p Ͻ0 for a real gas fluid.Hence, the transition is accompanied by a fluid expansion, since the temperature has to be continuously increasing away from the pellet surface, dT/dxϾ0.Furthermore, the fluid never gets compressed to a density higher than the solid density even at high pressure; the high compressibility factor Z ensures this, the ideal gas law does not, and for this reason it can easily fail in the transition region.It then follows that the last term in Eq. ͑8͒ is positive, so that dp/dxϽ0 at any moment within the layer.We note parenthetically that pellet ablation in tokamaks is quite different from the radiation driven ablation problem connected with inertial confinement fusion.In this situation, the piston action of ablation drives a shock front ahead of itself. 27An analogous model put forward by Lengyel,28 in which the ablation layer is accompanied by a shock front (dp/dxϾ0) is therefore not relevant for pellet ablation in plasmas.Continuing, we substitute the ideal gas parameters on the opposite side of the transition layer, f s ϵ1, ϭ ig , and pϭp ig ϭ ig kT ig /2m, into Eq.͑8͒ and find two solutions for the density ratio across the layer ig / 0 .Introducing the notation c 0 2 ϭp 0sur / 0 , and c ig 2 ϭkT ig /2m ϭ1.24ϫ10 5 m 2 /s 2 , the physical solution can be written as ig / 0 ϵc 0 2 /c ig 2 Ϫu 2 /c 0 2 ϩu 2 /c ig 2 , and thus the pressure drop across the layer ⌬pϭ p 0sur Ϫp ig is given by ⌬p in which we have assumed that u 2 /c 0 2 Ӷc 0 2 /c ig 2 .For example, taking r p ϭ2 mm, 0 ϭ200 kg/m 3 ͑deuterium͒, n e ϭ10 20 m Ϫ3 , and T e ϭ3 keV, Eqs.͑2͒ and ͑3͒ of Ref. 14 give uϭϪ7.83m/s, p 0sur ϭ7.47 MPa, c 0 2 ϭ3.73ϫ10 4 m 2 /s 2 , and therefore ⌬ p/p 0sur ϭ0.0038.This example verifies that the pressure within the layer hardly changes, as stated above.Accordingly, we can integrate the energy equation, Eq. ͑7c͒, across the layer while holding pressure fixed, and use the boundary conditions e→e 0 , → 0 , and q→0 at →Ϫϱ, and e→e ig , → ig and q→q sur on the other side, with q sur signifying the attenuated electronic heat flux incident on the layer.Then if we use our earlier definition for the sublimation energy per unit mass e ig Ϫe 0 ϵk⑀ s /2mϭ3.11ϫ10 5 m 2 /s 2 , the energy equation can be expressed in the familiar ''ablation-rate'' form Here, the second term in the denominator represents the specific work done by expanding the pellet gas/fluid material against the pressure force.Using our above example, we see that this term is smaller than the sublimation energy per unit mass: c ig 2 Ϫc 0 2 ϭ8.65ϫ10 4 m 2 /s 2 .Since both specific energy reservoirs are small, as well as q sur , the details of the physics within the transition layer are of little consequence.Hence, any nonideal gas behavior in the thin layer, in conjunction with the use of our artificial sublimation function f s (T) cannot have any significant effect on the surface recoil pressure p 0sur that develops by heating and expanding the ablated material, or the surface recession rate u, which is determined by the ablated mass flow through some surface enclosing the pellet, and not directly by Eq. ͑10͒.
The fact that the surface transition details are somewhat immaterial for pellet ablation under relevant plasma conditions is reinforced by the fact that arbitrarily changing the sublimation energy by factors of order 2 or more has virtually no effect on the ablation rate computed by the CAP code.Hence, there is no need for clustering the grid points in the high-gradient region near the pellet surface in order to accurately resolve the structure of the transition layer.More will be said about the computational aspects of the problem in Sec.III, and the pellet deformation problem will be discussed in detail in Sec.V.

C. Electronic heat flux model
The heat source Ϫ"•q ¯comes from energy deposition by hot, long-mean-free path, plasma electrons streaming into the ablation cloud from both sides along the magnetic field lines.For this reason the CAP simulation uses a cylindrical axisymmetric coordinate system (r,z), with origin at the center of the pellet and z axis parallel to the magnetic field ͑see Fig. 2͒.The parallel heat flux moment can be written as q ¯ϭq ϩ z ˆϪq Ϫ z ˆwhere the subscript ''ϩ'' ͑Ϫ͒ refers to the right ͑left͒-going plasma electrons.The heat source can thus be calculated field line by field line.Assuming incident Maxwellian electrons, a kinetic treatment using a Bethe-like collisional stopping power formula leads to the ''Bessel function'' heat flux moment 22 ͑see the Appendix for derivation͒: where K 2 is a standard modified Bessel function, and is the heat flux incident on the jet from ''infinity.''The argument u Ϯ ϵ Ϯ / ϱ is a dimensionless opacity, where are the respective line-integrated densities penetrated by right-going, left-going electrons arriving at the point (r,z) from infinity.Note that ϩ (r,z)ϭ Ϫ (r,Ϫz) on account of reflection symmetry.The projected range in the z direction for Maxwellian electrons with temperature T eϱ is given by The projected range is exactly 1/2 the total ''path length'' traveled by a gyrating electron with initial energy T eϱ as it slows down to rest ͑neglecting the weak energy dependence in the Coulomb logarithm ln ⌳͒.
A composite Coulomb logarithm is used to take into account slowing down interactions involving both free ͑f͒ and bound ͑b͒ electrons in the ablation cloud in which ⌳ eb ϭE/I * , I * ϭ7.5 eV, ⌳ e f ϭ0.2E/ប pe , បϭPlanck's constant, and pe ϭ(4n e e 2 /m e ) 1/2 is the plasma frequency.The density-effect correction due to polarization of the plasma medium by the incident electron is included in ⌳ e f .We suppress the energy dependence in the Coulomb logarithms by substituting Eϭ2.52T eϱ , which represents an averaged energy-flux-weighted value, so in eV ͑cgs͒units, we have Finally, taking Eqs.͑11͒-͑16͒ together, we arrive at 22 Ϫ"•q ¯ϭ q ϱ n t ͑ r,z

III. NUMERICAL SCHEME
The ablation cloud is a compressible material, because the expansion velocity may exceed the sound speed.On the other hand, the pellet behaves like an incompressible material.If the ablation pressure exerted on the pellet surface varies by more than the yield strength for solid deuterium, a few tenths of a MPa, 26 then the behavior of the pellet is like that of a hydraulic fluid.Otherwise, the pellet remains rigid while its surface recedes.We do not consider the material behavior in the transition through the yield point because of the uncertain nature of the viscoelastic-plastic behavior of the solid.Now in general the numerical scheme for a compressible fluid must be treated differently than one for an incompressible fluid, although their equations of motion are the same.In the present work, the CIP method is used that can treat compressible and incompressible fluids simultaneously without employing any special boundary conditions at the pellet surface as done in the previous steady-state ablation models.
Originally, the CIP method is one of the numerical schemes used to solve hyperbolic equations, such as wave equations. 17The set of hyperbolic Eqs.͑1͒ can be separated into an advection phase and a nonadvection phase by using time splitting.In the advection phase, the spatial profile of some hydrodynamic quantity f is interpolated within each grid with a cubic polynomial shifted in space according to the purely advection term in the hydrodynamic equation.In that way, f is updated from time step label n to a provisional time step label *.In the nonadvection phase, a conventional explicit method is used to update f from * to nϩ1, with the ''source'' terms being the terms on the right-hand side of Eqs.͑1͒.In determining the interpolated profile, both f and its spatial derivative grad f are used.The key feature of the CIP method is that the hydrodynamic equations for f are solved simultaneously with the corresponding equations for grad f.Thus, not only function values but also their spatial derivative are known at each grid point.In that way, overshooting during interpolation can be minimized for any hydrodynamic variable possessing a steep gradient.The CIP method has already been successfully employed in problems dealing with steep gradients, such as in situations where shock wave phenomena is expected.
Recently, the CIP method was extended by Yabe in order to solve incompressible and compressible regions simultaneously. 18The key point is that the pressure is determined by using the following equation instead of the EOS ͓that is, pϭkT/(2m) if ideal͔ in the nonadvection phase where c 2 ϭ(‫ץ‬p/‫)ץ‬ T is the isothermal sound speed squared, P TH ϭT(‫ץ‬p/‫ץ‬T) is the effective pressure that is reduced to the conventional pressure in an ideal fluid, C V ϭ(‫ץ‬e/‫ץ‬T) is the specific heat at constant volume, HϭϪ"•q ¯is the energy source, and ⌬t is a time increment.This equation is derived by combining Eqs.͑1b͒ and ͑1c͒ with the total differential dpϭ(‫ץ‬p/‫)ץ‬ T dϩ(‫ץ‬p/‫ץ‬T) dT, and making use of the above thermodynamic definitions.The essence of Eq. ͑18͒ is that it can automatically distinguish between incompressible and compressible regions by introducing a sound speed.In order to verify this feature, Eq. ͑18͒ is considered when Hϭ0 for simplicity, namely, the second term on the right-hand side is zero.In an incompressible region, each term in Eq. ͑18͒ has the following ordering: where wϵ⌬x/⌬t, with ⌬x being the grid width in the code, and we used the ordering p/ϳv 2 , and "ϳ1/⌬x ͑steep gradient near pellet surface͒.Simultaneously, the speeds satisfy the inequality vϽwӶc because the sound speed is regarded as infinity, and the Courant-Friedrichs-Lewy ͑CFL͒ stability condition should be satisfied.Then, the first term on the right-hand side of Eq. ͑19͒ becomes negligible compared with the other terms.Thus, Eq. ͑18͒ becomes " •("p nϩ1 /*)ϭ"•v ¯*/⌬t.Substituting the equation of motion in the nonadvection phase into that equation, we obtain In other words, Eq. ͑18͒ determines p nϩ1 such that "•v ¯nϩ1 becomes zero, namely, the incompressible condition is satisfied.
On the other hand, in the compressible ideal region, the ordering of Eq. ͑18͒ tends to ͑ c/w ͒ 2 ϳ1ϩv/w, ͑20͒ since p/ϳc 2 in the compressible ideal gas region.Now taking into account the CFL condition cϳvϽw, the lefthand side of Eq. ͑20͒ becomes negligible compared with the other terms.Especially, since the first term on the right-hand side of Eq. ͑18͒ becomes (p nϩ1 Ϫ p*)/͓␥p*(⌬t) 2 ͔ in an ideal fluid, Eq. ͑18͒ is reduced to temporal evolution for pressure in the nonadvection phase: p nϩ1 Ϫ p*ϭϪ␥p*" •v ¯*⌬t.This corresponds to an equation derived by substituting an ideal EOS into an energy equation with Hϭ0.In other words, p nϩ1 is determined such that the ideal gas EOS is satisfied in the compressible region.As a result, Eq. ͑18͒ can treat the entire region without calculating individually incompressible and compressible regions.

IV. PELLET ABLATION IN THE ONE-DIMENSIONAL SPHERICALLY SYMMETRIC HEATING APPROXIMATION
The utility of the CIP method is demonstrated by comparing numerical predictions with known solutions of pellet ablation models in dimensional ͑1D͒ spherically symmetric geometry, as shown in Fig. 1͑b͒.Thus, Eqs.͑1͒ are written in spherically symmetric form with radial coordinate r, and the origin coincides with the center of the pellet.We also replace Ϫ"•q ¯by the equivalent 1D approximation, namely, Ϫdq/dr.
In the first test of the CAP code, we compared the ablation rate and fluid profiles with the original TF model. 9For that comparison we used the monoenergetic heat source instead of the Maxwellian heat source.We also neglect dissociation and ionization by setting f d ϭ f i ϭ0, and take ␥ϭ7/5 appropriate for a molecular gas.After a few microseconds of exposure to a constant heat flux, the ablation flow becomes fully developed, and the ablation rate settles down to a steady state just as predicted by the TF model.We choose a deuterium pellet with initial radius r p ϭ2 mm, plasma electron temperature T eϱ ϭ2 keV, and number density n eϱ ϭ10 14 cm Ϫ3 .In all the cases in this paper, T eϱ and n eϱ are assumed to be constant throughout the temporal evolution.For this example, the steady-state ablation rate was found to be Gϭ41 g/s by simulation.Remarkably, the simulation reproduces the TF result G TF ϭ41.27 g/s given by Eq. ͑31͒ of Ref. 9 with Qϭ1.The CAP code also reproduced identical dimensionless fluid profiles.
Using the same pellet/plasma parameters, the ablation rate increased from 41 g/s to 113 g/s when the Maxwellian heat source was employed in the simulation.This significant increase was expected, since most of the heat flux transported through the dense ablation cloud comes from the long-range suprathermal electrons in the Maxwellian distribution; 10 the ablation cloud density and thickness is much larger for the Maxwellian heat source than the monoenergetic one.As a result, the energy deposited per particle ablated is less, and the temperature and flow velocity become uniformly lowered.For example, at the radius rϭ0.5 cm, the temperature of the cloud with the monoenergetic heat source is 8.2 eV, but it is only 3.5 eV for the Maxwellian heat source.
Next we included the effect of atomic processes in the ablation cloud and used the Maxwellian heat source.Only a small ϳ6.2% reduction in the ablation rate, from 113 g/s to 106 g/s, was found.The most significant effect of atomic processes was the lowering of the ablation temperature.At the first sonic surface, T * dropped from 3.44 eV to 1.11 eV.Such large ϳ3ϫ temperature drops were also noted in Ref. 18 using the monoenergetic heat flux approximation.Figures 3͑a͒ and 3͑b͒ show radial profiles of various quantities at 10 s without and with atomic processes, respectively.Figure 3͑a͒  are the pressure, temperature, and radius, respectively, at the first sonic (M ϭ1) point.The pressure decreases, and the temperature and Mach number increase in the radial direction outside the pellet because the ablation cloud expands while heated.On the other hand, the pressure is constant inside the pellet because the pellet is rigid ͑for a spherically symmetric pressure distribution͒ and incompressible.Figure 3͑b͒ shows p/p * , T/T * , M, fractional dissociation f d , and ionization f i as functions of r/r * .It is clear that the ablation phenomena is changed by including the atomic processes.The fact that radial profile in Fig. 3͑b͒ has a jump structure is the most essential point.There are two features in that fact.The first feature is that the flow is supersonic ahead of the jump structure and subsonic downstream.This means that the jump structure obviously represents a shock.The stationary shock front is found in the region where ionization is nearly completed f i ϳ1, and it is due to the removal of the ionization energy sink which tends to cause a rollover in the Mach number profile, possibly driving the Mach number below unity in some cases.However, the Mach number never dips below unity which would lead to a flow singularity as in Felber et al.; 19 instead, the ablation flow develops a stationary shock front, across which the fluid undergoes the transition from supersonic to subsonic state.Farther downstream, a ''second'' sonic surface appears.The second feature lies in the difference between the effects of dissociation and ionization on the structure of the Mach number profile.Note that a stationary shock front does not appear in the region where dissociation is dominant (0Ͻ f d Ͻ1) because dissociation is always completed in the subsonic region of the flow.Figure 4 displays the temporal evolution of the Mach number profile during the initial transient period (0.2рtр3 s).At tϭ0.2 s, we see a local depression in the Mach number near the pellet surface, but no shock.This is because during that moment the flow is fairly supersonic when the energy sink associated with the latent heat of dissociation has ''dried up.''However, the Mach number profile is becoming less supersonic in time, which contributes to the formation of a definite shock front at tϭ0.4 and 0.6 s in the region near the pellet where dissociation is completed.Farther out, in the region where ionization is completed, the Mach number experiences only small depressions because the Mach number is still well above unity during those times.However, the shock structure associated with dissociation becomes gradually flattened by the sound wave since the region where dissociation dominates is becoming more and more subsonic in time.On the other hand, the shock formed by ionization has been maintained because in the approach to steady state the flow always remains supersonic where f i ϳ1.Therefore, the ablation cloud has two sonic surfaces and one shock front in the case with atomic processes, though it has one sonic surface in the case without atomic processes.It is not clear to us why the simulations in Ref. 11 with atomic processes revealed no shocklike structure in the flow parameters.Also Eq. ͑40͒ in Ref. 11, which is a fit to the numerical results for T * without atomic processes, seems to predict unrealistically high values, casting doubt on the results when these effects are included.
Using the same pellet/plasma parameters as in Fig. 3, we display in Fig. 5 the (r,t) behavior of the pellet, sonic, and shock surface radii over the long time scale of the pellet lifetime.Figure 5͑a͒ shows the case without atomic processes, and Fig. 5͑b͒ shows the case with atomic processes.Note that the pellet lifetime in the case with atomic processes is slightly longer in proportion to the 6.2% reduction in the ablation rate.The ablation cloud without atomic processes is fully developed after ϳ2 s, and subsequently, the pellet and sonic surface radii are reduced by the surface recession.It is clear that the flow structure shown in Fig. 3͑a͒ changes slowly, i.e., the ablation flow is quasi-steady-state.However, with atomic processes, the time to fully develop the flow and establish the first sonic, shock, and second sonic surfaces is longer, ϳ10 s.Since the shock front associated with dissociation is erased early on, it does not appear on the long time scale of Fig. 5.The shock front induced by ionization remains almost stationary.It is not carried inwards while the pellet shrinks, like the first sonic surface is, because its inward motion is compensated by an increase in the velocity of the upstream flow disappearing into the shock.Interestingly, a new dissociation shock appears near the end of the pellet's life at ϳ80 s.This is because the dissociation dominant region is crossed by the first sonic surface as it follows the inward progression of the pellet surface.The crossing results from the fact that the ablation temperature at the first sonic surface is becoming cooler as the pellet shrinks consistent with the TF scaling laws. 9Hence, the region where dissociation occurs has to move up the temperature gradient, i.e., it has to expand into the supersonic region.Consequently, when the dissociation energy sink is removed, f d ϳ1, a shock front may appear.A new sonic surface induced by dissociation approaches zero as the pellet surface radius approaches zero similar to one in the case without atomic processes shown in Fig. 5͑a͒.
As anticipated, the double-transonic flow solution with a shock discontinuity avoids the flow singularity often encountered in Ref. 19, and generalizes the validity of the quasisteady assumption under a wider set of plasma-pellet parameters.We carried out simulations over a nominal parameter range in order to compare the ablation rates between the cases without (G iso w/o ) and with (G iso with ) atomic processes.Our results, displayed in Table I, show that the difference is consistently less than 10%.Felber et al. 19 found that atomic processes can extend pellet lifetimes by 10%-20%.This dis-  crepancy is presumably due to the use of the monoenergetic heating approximation in Ref. 19, which, as noted earlier, uniformly heightens the ablation temperature.Thus, the ionization and dissociation heat sinks are moved closer to the pellet surface where the influence on the ablation rate becomes stronger. 19The most important effect of the atomic processes combined with the Maxwellian heat source is that the ablation temperature at the first sonic radius is consistently reduced by about a factor of ϳ3.Using the parameters in Fig. 3 again, we find that without atomic processes the temperature at the sonic surface is T * ϭ3.44 eV, and with atomic processes T * ϭ1.11 eV.This leads to a significant reduction in the electrical conductivity of the ablation cloud, and would clearly influence how the ablation flow is influenced by the induction currents generated by the component of the flow across the magnetic field.This is a problem left for a future study.

V. ABLATION WITH ANISOTROPIC HEATING
In this section, we examine 2D effects resulting from anisotropic heating imposed by the magnetic field.We use the cylindrical axisymmetric coordinate system (r,z), and again we use the Maxwellian heat source model in this geometry as described in Sec.II C. Again, the readers are reminded that no effect of the magnetic field, except for setting the direction of the electron heat flux, is included in our present simulations.The first set of simulations is designed to isolate just the geometrical effect of anisotropic heating on the ablation rate without including the effect of pellet deformation.This is accomplished by arbitrarily increasing the density of the solid pellet by a factor of 100 so that the ablation cloud can evolve to a steady-state ablation rate after a few microseconds, while the pellet remains effectively rigid or ''hard.''Obviously, this artificial density increase has no effect on the ablation rate or the ablation flow quantities since they have no density dependence. 9,11Figure 6 compares the ablation rates between the CAP code simulation and the 2D simulation by MacAulay. 11Figure 6 also shows the modified TF formula 14 given in the Introduction, which is consistent with the experiments.In Figs.6͑a͒-6͑c͒, the dependencies of the ablation rate on plasma temperature, density, and pellet radius, respectively, are separated.All results have essentially the same TF scaling behavior.However, the multiplicative constant in the ablation rate differs among those results.As noted in the Introduction, the ablation rate from Macaulay's simulation is systematically a factor of ϳ2 higher than that given by the modified TF model.The ablation rate from the CAP simulation lies in between those.By comparing the ablation rates for cases 1-7 displayed in Fig. 6 with the ablation rates displayed in Table I for the isotropic heating model, we find that the ratios of the ablation rates are fairly consistent, namely, G aniso/hard G iso ϳ0.51-0.57,͑21͒ where G iso and G aniso/hard stand for the ablation rates in isotropic heating and anisotropic heating with a hard pellet, respectively.The nearly factor of 2 reduction in the ablation rate caused by anisotropic heating arises because shielding of the pellet is primarily due to the high-density layers of the cloud adjacent to the pellet.For a detailed explanation see Ref. 14.The main idea is that the effective area intercepting the electron heat flux, i.e., the projected area of the pellet normal to the B field is the main factor that determines the ablation rate. 10,11,14For a spherical pellet, the projected area is half the surface area of the pellet, hence the factor of 2 reduction.A different result, factor of ϳ0.68 reduction, was found in MacAulays's simulation, 11 although it should be mentioned that the simulation in 1D ͑isotropic heating͒ was not possible.To study these geometrical effects on the ablation rate systematically, MacAulay simulated 2D ablation with the monoenergetic heat flux model, and compared the simulation result with the original TF ͑isotropic heating͒ ablation rate expression. 9The reduction factor stated in Ref. 11 was systematically 0.52.However, a more careful comparison with Eq. ͑31͒ of Ref. 9 indicates that the reduction factor is closer to ϳ0.68 ͑see also Ref. 14͒.
In the remaining simulations, we use the actual mass density of solid deuterium 0 ϭ0.2 g/cm 3 in order to include pellet interior flows and deformation during the ablation process, and we also include atomic processes.Figure 7 shows the isodensity contours in the (r,z) plane at 2 s, using the pellet/plasma parameters of Fig. 3.The central solid ͑red͒ region represents the pellet.The contour lines surrounding the pellet identify the oval-shaped first and second sonic surfaces, and the shock front nested in between those.Figures 8͑a͒ and 8͑b͒ shows, respectively, the fluid profiles along the z direction ͑parallel to polar axis and through rϭr p /2) and along the r direction ͑equatorial midplane at zϭ0) at 2 s.The temperature at the sonic point along z direction is T * ϭ1.14 eV, very close to the previous value of 1.11 eV for the isotropic heating case, as expected.However, T * ϭ1.51 eV on the equatorial midplane; it is higher because there is less shielding ͑line-integrated density͒, and also because the incident plasma electrons are approaching from both directions outside the ''shadow'' of the pellet.We see in Fig. 9 that the pellet surface pressure p 0sur at the two poles labeled by z ϩ and z Ϫ is a factor of 2 higher than the pressure at the two equatorial points r ϩ and r Ϫ .The absolute pole pressure is 12 MPa.Taking into account the correction due to electrostatic shielding mentioned in the Introduction, the actual pressure becomes 0.3ϫ12ϭ3.6MPa.This result is in close agreement with Eq. ͑3͒ of Ref. 14.The surface pressure differential is of course due to the fact that the heat flux following magnetic field lines intercepting the pellet near the polar axis is entirely blocked by the pellet and thus fully absorbed, while the heat flux following field lines just grazing the pellet surface is partially transmitted across the equatorial midplane, so there is less absorption at the equatorial points.Since the areal ablation rate is more enhanced at the poles than at the equatorial points, the ablation cloud becomes thicker in the polar direction than it is along the radial coordinate in the equatorial midplane.This leads to the ovalshaped density contours elongated in the polar direction as shown in Fig. 10͑a͒ at different times: 10, 20, 30, 40, 50, and 60 s.The contour of the pellet density ͑which is uniform inside the pellet͒ has the opposite shape.The anisotropic pressure distribution deforms the pellet into an oval shape flattened in the polar direction and elongated in the radial direction.After a longer time, the pellet becomes flattened into a thin pancake while losing mass by ablation.The isopressure contours corresponding to these time frames are shown in Fig. 10͑b͒.
Figure 10͑a͒ also illustrates the temporal evolution of the sonic and shock surfaces, whose shape can be understood as follows.As we shall see below, when the pellet is allowed to deform, the ablation rate actually increases at first.This is contrary to what the conventional ablation models indicate, namely, that the ablation rate should decrease while its ra-dius is shrinking, since Gϰr p 4/3 .Since the density is increasing with the increasing ablation rate, the fluid velocity and the temperature must decrease.However, the sound speed ϳT 1/2 decays more slowly than the velocity.Eventually, the subsonic region has expanded to such an extent that the first sonic surface with the adjacent shock front no longer exist in the flow domain.However, the second sonic surface still survives, but it lies beyond the image box as shown by the image at 40 s.Later, however, the ablation rate is coming down along with the cloud density.This in turn uniformly increases the temperature causing the reemergence of the first sonic surface and shock front.All three surfaces start to collapse towards the pellet, and we see their reappearance in the image box at 60 s.
Figure 11 shows the temporal behavior of pressure at the pellet surface, where p A is defined on the equator A and p B is defined at the point B where projection on the plane with z ϭ0 corresponds to point C which is located at half radius of the equator as shown in Fig. 12.Here we have included atomic processes.Initially, the difference ⌬ pϭp B Ϫp A reaches 11 MPa; adjusting for electrostatic shielding, ⌬ p ϭ3.3 MPa.This difference is somewhat greater than the yield stress of solid deuterium ϳ0.5 MPa. 26 Hence, pellets with r p ϳ0.2 cm should become fluidized, and will begin to flatten out when T eϱ տ2 keV and n eϱ տ10 14 cm Ϫ3 .Figure 13 shows the time evolution of the major equatorial radius r pellet , and half thickness z pellet of the disklike pellet shape of the pellet.Here, r pellet corresponds to the equatorial radius OA in Fig. 12, while z pellet corresponds to the line segment CB defined in Fig. 12 instead of OF since pellet erosion along OF is completed sooner.As we see, pellet deformation promotes the increases in r pellet at first, but later on mass erosion intensifies at the edges of the disk, causing r pellet to rapidly fall off.Apparently, the motion of r pellet is determined by a balance between the pellet deformation and the surface recession due to ablation.The rollover in r pellet begins when the changing pellet shape starts to bring p A and p B closer together, which in turn lessens the shear stress causing deformation.On the other hand, the thickness z pellet continually decreases due to ablation ͑lentil effect 16 ͒.Figures 14͑a͒ and 14͑b͒ show the (z,t) and (r,t) trajectories, respec- tively, of the pellet, sonic, and shock surfaces.The solid lines are plots of the pellet surface chords, z pellet (CB) and r pellet (OA).The dotted lines mark the first sonic surface chords, z 1st (CE) and r 1st (OD), and the second sonic surface chords, z 2nd (CI) and r 2nd (OJ).The dot-dashed lines mark the corresponding chords for the shock surface, z shock (CG) and r shock (OH).It is found that the dimensions of all three surfaces increase at first, and decrease later on.It is clear that the shock surface in the r direction disappears sooner than one in the z direction.That is due to the fact that the temperature in the r direction is higher than it is in the z direction, as we mentioned above.Then, because the subsonic length in r direction is longer than the length in the z direction, the shock surface in the r direction disappears sooner.Regarding the cloud thicknesses z 2nd Ϫz pellet (BI) and r 2nd Ϫr pellet (AJ) as the distance between the pellet surface and second sonic surface, it is clearly understood that those thicknesses are comparable.This is due to the fact that expanding motion of the cloud is unaffected by the magnetic field even if the pellet is deformed.Of course, these results, in particular, are likely to be much different when the JϫB force is included in our simulation.
Figure 15 plots the time dependence of the volume and area of the pellet projected normal to the magnetic field lines r pellet 2 , i.e., the effective area which absorbs the incoming heat flux.Although the volume decreases monotonically with time due to surface recession, the area increases at first, reaches a peak, and then falls to zero later on.
One of the key results of our study is that pellet deformation can shorten the life span of a pellet.Figure 16  ϳ10 s as shown in Fig. 15, because the total heating power available for driving ablation is simply proportional to this area.The dashed line in Fig. 16 plots the ablation rate for the 1D spherically symmetric isotropic heating model used in the simulation.After the sharp rise to the peak value of 113 g/s mentioned in Sec.IV, whereupon the ablation flow becomes fully developed, the ablation rate steadily decreases in accordance with the monotonically decreasing pellet size.We can compare this ablation rate curve with the corresponding analytical result ͑thin solid line͒ obtained by a simple integration of the TF ablation rate formula, where M 0 is the pellet mass, G 0 is the ablation rate at t ϭ0, and the pellet lifetime is t p ϭ(9/5)M 0 /G 0 .The two curves match remarkably well, adding again to the validity of the code.Notice that the pellet lifetime for the 1D isotropic heating simulation is 103 s.The corresponding pellet lifetime predicted by the modified TF model 14 is 309 s (G 0 ϭ39 g/s).As noted in the Introduction, the modified TF model is consistent with IPADBASE. 15However, the 2D model predicts a pellet lifetime of 70 s; a factor of 4 smaller than the modified TF model.Now according to Eq. ͑21͒, the 2D simulation of the ''hard pellet'' that is not allowed to make a deformation, would predict a pellet lifetime of ϳ103/͑0.54͒ϭ191s.Hence, there is a large difference-a factor of 2.7-in the lifetimes between hard and ''soft'' pellets!An obvious consideration is that we did not take into account electrostatic shielding.As noted in the Introduction, electrostatic shielding cuts the ablation rate by a factor of 0.51-0.55.If we choose 0.53, the readjusted pellet lifetimes for the three simulation cases become t p iso ϭ191 s, t p aniso/hard ϭ360 s, t p aniso/so f t ϭ132 s, t p TF ϭ309 s.
Clearly the modified TF lifetime is intermediate between the soft and hard pellets' lifetimes.This might suggest that the pellet undergoes deformation only in the beginning of ablation, when the surface pressure variation ⌬p is at a maximum ͑see Fig. 11͒, and ⌬pϾ.Later on ⌬p diminishes due to the deformation, until at some point ⌬ pϽ.Thereafter, the pellet would ablate as a rigid body, and consequently the lifetime would be extended.We might also add that experimental determination of G(t) and pellet lifetime ͑or penetration distance͒ is done in plasmas with realistic plasma profiles such that T eϱ and n eϱ are varying with time in the reference frame of the pellet.The experimental curves for G(t) always fall off much more sharply at the end of the pellet's life than ablation curves calculated by all sphericalpellet ablation models.The abrupt decline in the 2D ablation curve in Fig. 16 as compared with the 1D spherically symmetric curve strongly suggests an explanation that lies in the changing pellet shape during ablation.Finally, we notice that the 2D ablation rate curve in Fig. 16 forms a pedestal at about 2 s.The ablation rate at the pedestal is found to be about one half the value of that for 1D isotropic heating since the pellet still keeps a spherical shape in the beginning as the ablation flow becomes fully developed ͑see Fig. 7͒.However, when a pellet starts to deform, the ablation rate is enhanced due to increase of effective area absorbing the heat flux, which results in the shorter lifetime.
In order to investigate the pellet lifetimes more systematically, several runs are carried out for plasma/pellet parameters shown in Table I.The results presented in Fig. 17 are not readjusted for electrostatic shielding.The abscissa of the plot corresponds to the lifetimes predicted by the 1D simulation without atomic processes for the various case numbers.The ordinate indicates the corresponding lifetimes of the other models: 1D simulation with atomic processes ͑circles͒, 2D soft pellets without atomic processes ͑crosses͒.MacAulay's simulation ͑squares͒, and the modified TF model ͑triangles͒.It is clear that the lifetime for the 2D soft pellet is always lower than the MacAulay and TF models, however, the initial ablation rate for the hard pellet always lies in between the MacAulay and TF models ͑see Fig. 6͒.
Our present simulations do not include viscosity and  I, respectively.1D simulation with atomic processes ͑circles͒, 2D soft pellets without atomic processes ͑crosses͒, MacAulay's simulation ͑squares͒, and the modified TF model ͑triangles͒.
elastic-plastic effect in the pellet.Therefore, calculation by an advanced code including those effects should be required in order to obtain more detailed information about the deformation.

VI. CONCLUSIONS
A high-resolution hydrodynamic code CAP has been developed to investigate the pellet ablation problem.The code uses the CIP numerical scheme, which can accurately capture the solid to gas phase change without having to impose special boundary conditions at the pellet surface.For nominal surface ablation pressures, the use of the ideal gas EOS was found to be a valid approximation near the pellet surface where the solid is transformed to a dense supercritical fluid.It is shown that a stationary shock wave is formed in the ablation flow where thermal ionization is nearly completed.The shock wave resolves the flow singularity encountered in previous steady-state ablation models, and thus allows the flow to reach a quasi-steady-state early on.No shock front is associated with molecular dissociation process because dissociation is completed in the subsonic flow region, except near the end of the pellet's life when a small shock appears.The ablation rate is reduced only by a few percent by the energy sinks made available by the latent heats of ionization and dissociation.However, these atomic processes can reduce the Spitzer electrical conductivity of the cloud by a factor of ϳ5 which would have significant impact on the magnetohydrodynamic forces acting on the flow, a problem left for future study.
We have shown that the ablation processes quickly reach quasi-steady-state for isotropic 1D heating and for anisotropic 2D heating, if the pellet does not deform.An initially spherical pellet can deform into a pancake shape, elongated in the direction perpendicular to the B field, due to the uneven surface ablation pressure, when the pressure variations exceed the yield strength of solid deuterium by a few times.In light of the surface pressure scaling p 0sur ϳn eϱ 2/3 T eϱ 1.54 , this study concludes that the threshold condition for such deformations are obtained for millimeter sized pellets only if (n eϱ /10 14 ) 2/3 (T eϱ /10 3 ) 1.54 Ͼ3; for example, when T eϱ տ2 keV and n eϱ տ10 14 cm Ϫ3 .Thus, it is likely that ablation pressure driven deformation is not significantly affecting the pellet lifetime and penetration distance in present-day devices.We also found that for constant plasma conditions, a quasi-steady-state for such soft pellets does not exist.Instead, the ablation rate has to be determined at each moment in order to predict the pellet lifetime.Furthermore, the deformation is self-healing: the pressure variation ⌬ p weakens in time.At some point, it is likely that ⌬p will drop below the yield strength, stopping further deformation.In the case of pellet injection into a fusion device with a high edge ''pedestal'' temperature Ͼ4 keV, such deformations can become worrisome.We find under constant plasma conditions that the lifetime of a soft pellet is shortened by a factor of 2.7 from that of a hard pellet, and its lifetime is shortened by a factor of 2.3 compared with that predicted by the modified TF model, 14 which is in agreement with experiments.
Though not explored in this study, the penetration of pellets into plasmas with realistic profiles could be done with the CAP code.One must, however, take into consideration that the pellet has both soft and hard features that depend on the relation between the surface pressure asymmetry at any given moment along the trajectory and the yield strength of the solid.Future penetration predictions for next-step fusion experiments will include the influence of plasma profile effects on pressure driven pellet deformations.We also intend to use the CAP code to study the formation and structure of the extended field-aligned plasma shield caused by the nozzling effect of the magnetic field in order to better quantify its ability to shield the pellet and reduce the ablation rate.

ACKNOWLEDGMENTS
One of us ͑R.I.͒ acknowledges fruitful discussions with Professor K. Itoh.This is a report of work supported in part by the U.S. Department of Energy under Grant No. DE-FG03-95ER54309, and in part by grant-in-aid from the Ministry of Education, Culture, Sports, Science and Technology, Japan.

APPENDIX: DERIVATION OF SEMIANALYTIC KINETIC HEAT FLUX EXPRESSION
In this appendix, we give an explicit derivation of the semianalytical kinetic formula for the transport of electron energy flux through the ablation cloud, namely, the normalized energy flux used in Eqs.͑7͒ and ͑13͒ in Ref. 22. Attenuation of the heat flux by friction with the cloud electrons is mainly due to electron slowing down on free and bound electrons.In either case, the slowing down force for an energetic electron with energy E has the Bethe-like form, dE/dsϭ(2n t e 4 /E)ln ⌳(E), where ds is the differential path length and n t is the number density of nuclei, and the Coulomb logarithm is defined in Eq. ͑15͒ for a partially ionized hydrogenic gas.The derivation begins by neglecting pitch-angle scattering; each electron is assumed to retain its original cosine pitch angle while it slows down.Since the one way electron distribution function is a half-space Maxwellian, and thus isotropic over the interval ͑0,1͒ to begin with, the diffusion of pitch angles within and out of this interval will not affect the energy flux to any great extent because the energy flux moment involves a convolution of over this interval.A further approximation, especially for energetic electrons with EϾ300-500 eV, is to neglect the weak energy dependence in the Coulomb logarithm.The solution of the kinetic equation for the distribution function of electrons traveling in the right ϩz ͑left Ϫz) directions that have penetrated a line-integrated density ϩ ( Ϫ ) is given by 10,29,30 where where q ϱ ϭ(2/m e ) 1/2 n eϱ T eϱ 3/2 is the incident heat flux.Making a variable change yϭ(x 2 ϩ 2 ) 1/2 gives q ¯Ϯϭ ͵ 0 1 d ͵ ϱ ͑ y 2 Ϫ 2 ͒ 1/2 y exp͑Ϫy ͒dy, ͑A5͒ where ϭ(u Ϯ /2) 1/2 .To carry out the integral over y, we consider the auxiliary integral I(a,), which is a standard integral, 31 I͑a, ͒ϵ ͵ ϱ ͑ y 2 Ϫ 2 ͒ 1/2 exp͑Ϫay ͒dyϭ a K 1 ͑ a ͒, ͑A6͒ with K n being a modified Bessel function of order n.Then by differentiating both forms of I(a,) with respect to a, and setting aϭ1, it follows that q ¯͑u ͒ϭ ͵ 0 1 ͑ u/2 ͒K 2 ͓͑u/2 ͒ 1/2 ͔d, ͑A7͒ where q ¯ϭq ¯Ϯ and uϭu Ϯ .The angle integration in Eq. ͑A7͒ cannot be done analytically, nevertheless, we can supply a justification for a simple approximate form having great accuracy.Consider first the heat deposition rate per unit volume in the case of no pitch-angle scattering ͑NS͒.This is simply Ϫ"•q ¯NS ϭϪ n t q ϱ ϱ ͩ dq ¯ϩ du ϩ dq ¯Ϫ du ͪ .

͑A8͒
Making use of the relation we obtain When the ablation cloud is fully ionized and has no bulk motion, the heating rate given by Eq. ͑1c͒ and the definition of e in Eq. ͑2b͒ is given by ‫ץ‬T ‫ץ‬t ϭ Ϫ"•q NS 3n t .͑A11͒ In the thin cloud limit (uӶ1), we see that g(u)→Ϫ1/4.This is the limit that corresponds to the situation where a cold Maxwellian plasma is being heated up by an intermingling hot Maxwellian plasma with temperature T eϱ ӷT, and uniform density n eϱ .Then, Eq. ͑A11͒ reduces to a similar result obtained by Spitzer 32 ‫ץ‬T ‫ץ‬t ϭ T eϱ 2t eq ͑ in limit uӶ1 ͒, ͑A12͒ where t eq is Spitzer's temperature equilibration time for T eϱ ӷT, t eq ϭ 3m e 2 8͑2 ͒ 1/2 n eϱ e 4 ln ⌳ ͩ T eϱ m e ͪ 3/2 .͑A13͒ Obviously, the reason that the factor of 1/2 appears in Eq. ͑A12͒ and not in Ref. 32 is because the energy deposited on the cold electrons is shared with an equal number of cold ions having the same temperature T. It should be mentioned that the reduction to the Spitzer-like heating rate in the small u limit was obtained for a fully ionized plasma.However, our expression for "•q ¯is still valid for the partially ionized case, provided the ''composite'' Coulomb logarithm given by Eq. ͑15͒ is used.The last step is to replace the assumed isotropic distribution in pitch angle with an equivalent ␦-function distribution, i.e., we replace 1d in Eq. ͑A10͒ by ␦( Requiring this to match the exact g(u) in the small u limit yields ¯ϭ1/2.Thus for all u an appropriate expression for g becomes Integrating Eq. ͑A15͒ gives the corresponding normalized energy flux q ¯ϭ 1 2 uK 2 ͑ u 1/2 ͒ ͑A16͒ given by Eq. ͑11͒.
The Bessel function form can also be found from Eq. ͑12͒ of Ref. 29 after setting pϭ1 and Z a ϭ0 ͑zero pitchangle scattering͒.Note also that Ref. 29 used a ␦-function distribution for to evaluate the heat flux, but no justification for choosing ¯ϭ1/2 was provided, as done in this appendix.The Bessel function heat flux formula is in remarkable agreement with the heat flux moment of the electron distribution function obtained by a numerical solution of the Fokker-Planck calculation which included pitch-angle scattering. 11The reason for the agreement is that the cumulative effect of pitch-angle scattering on the energy flux moment is very weak for low-Z hydrogen, as supported by the rather small value of reflectivity coefficient 0Ͻ␦ R Ͻ␦ R max ϭ8%.11 Here, the reflectivity ␦ R is defined as the ratio of the backscattered energy flux to the forward energy flux at the surface of the cloud, and it can be taken into account in our formulas by simply replacing the incident heat flux q ϱ with the net flux q ϱ (1Ϫ␦ R ).The maximum reflectivity ␦ R max occurs for infinitely thick targets, e.g., on field lines intersecting the pellet.Such a small correction has a negligibly small ͑ϳ3%͒ effect on the ablation rate.

FIG. 2 .
FIG. 2. ͑Color online͒.͑a͒ The compressibility factor Zϭ͓2mp/(pkT)͔ for hydrogen as a function of pressure for a family of isotherms above, but near to the critical temperature T c ϭ33.2 K. ͑b͒ A sketch of the -p phase diagram for a family of isotherms near the critical point, which is marked by a circle.͑c͒ A sketch of the -T phase diagram for a family of isobars near the critical point.The dashed curve is the boundary of the two-phase liquidvapor coexistence region, and the notations Ϫ, c, ϩ, indicate a value below, at, and above the critical.
shows normalized pressure p/p * , temperature T/T * , and Mach number M as functions of normalized radius r/r * , where p * , T * , and r *

FIG. 3 .
FIG. 3. ͑Color online͒.Normalized fluid profiles at 10 s in 1D spherical symmetry of ͑a͒ ablation without atomic processes and ͑b͒ with atomic processes.Here, r p ϭ2 mm, T eϱ ϭ2 keV, and n eϱ ϭ10 14 cm Ϫ3 .The solid, dashed, and dot-dashed lines are M, p/p * , and T/T * as functions of r/r * , respectively.Dotted lines in ͑b͒ show f d and f i profiles.

FIG. 5 .
FIG. 5. Pellet, sonic, and shock surface radii as a function of time for ͑a͒ ablation without atomic processes and ͑b͒ with atomic processes.Pellet, sonic, and shock surface radii are shown by solid, dashed, and dot-dashed lines, respectively.

FIG. 7 .
FIG. 7. ͑Color͒.Isodensity contours on the r-z plane in the 2D cylindrical axisymmetric coordinate system at 2 s.Black lines indicate first sonic, shock, and second sonic surfaces.

2 afterFIG. 12 .
FIG.12.Sketch of the deformed pellet and surrounding ablation cloud.Solid, dashed, and dot-dashed ovals are the pellet, sonic, and shock surfaces, respectively.O is the center point.A, D, H, and J are points on the equators of the pellet, first sonic, shock, and second sonic surfaces, respectively.Point C marks the half radius of the pellet equator (COϭCA).B is a point on the pellet surface whose projection on the zϭ0 plane corresponds to point C. Point E, G, and I mark the intersections of the line extending from CB to the first sonic, shock, and second sonic surfaces.Point F marks the intersection between the pellet surface and the z axis.

FIG. 16 . 3 .FIG. 17 .
FIG. 16.Time dependence of the pellet ablation rate for a constant heat flux over the lifetime of the pellet.The solid curve pertains to 2D cylindrical axisymmetric geometry and the dashed curve is for the 1D spherically symmetric geometry.The thin solid line is the ablation rate from the transonicflow model with the well-known scaling law Gϰr p 4/3 .

TABLE I .
Parameters used in the simulation.rp , T eϱ , and n eϱ are pellet radius, electron temperature, and density in the bulk plasma.G iso w/o and G iso with are ablation rates without and with atomic processes.Cases are corresponding to numbers, respectively, in Figs.3 and 11.
FIG. 6. ͑a͒ T eϱ , ͑b͒ n eϱ , and ͑c͒ r p dependences of ablation rates in 2D cylindrical axisymmetric geometry.The circles with numbered labels correspond to simulation results listed in Table I.The solid lines are results from the modified TF model of Ref. 14.The dashed lines are results from Ref. 11.