General Atomics July 14, 2009

1. General Atomics July 14, 2009 • Multiphase MHD at Low Magnetic • Reynolds Numbers • Tianshi Lu • Department of Mathematics • Wichita State University In collaboration with Roman Samulyak, Stony Brook University / Brookhaven National Laboratory Paul Parks, General Atomics

2. Motivation • Tokamak (ITER) Fueling • Fuel pellet ablation • Striation instabilities • Killer pellet / gas ball for plasma disruption mitigation Laser ablated plasma plume expansion Expansion of a mercury jet in magnetic fields

3. Talk Outline • Equations for MHD at low magnetic Reynolds numbers and models for pellet ablation in a tokamak • Numerical algorithms for multiphase low ReM MHD • Numerical simulations of pellet ablation

4. Equations for MHD at low magnetic Reynolds numbers Full system of MHD equations Low ReM approximation Elliptic Equation of state for plasma / liquid metal / partially ionized gas Ohm’s law Maxwell’s equations without wave propagation Parabolic

5. Models for pellet ablation in tokamak Global Model Local Model Courtesy of Ravi Samtaney, PPPL Tokamak plasma in the presence of an ablating pellet Pellet ablation in ambient plasma • Full MHD system • Implicit or semi-implicit discretization • EOS for fully ionized plasma • No interface • System size ~ m, grid size ~ cm • MHD system at low ReM • Explicit discretization • EOS for partially ionized gas • Free surface flow • System size ~ cm, grid size ~ 0.1 mm

6. Sheath Fluxes Sheath boundary Cloud Plasma f(z) Schematic of pellet ablation in a magnetic field Schematic of processes in the ablation cloud

7. Local model for pellet ablation in tokamak • Axisymmetric MHD with low ReM approximation • Transient radial current approximation • Interaction of hot electrons with ablated gas • Equation of state with atomic processes • Conductivity model including ionization by electron impact • Surface ablation model • Pellet penetration through plasma pedestal • Finite shielding length due to the curvature of B field

8. 1. Axisymmetric MHD with low ReM approximation Centripetal force Nonlinear mixed Dirichlet-Neumann boundary condition

9. 2. Transient radial current approximation f(r,z) depends explicitly on the line-by-line cloud opacity u. Simplified equations for non-transient radial current has been implemented.

10. 3. Interaction of hot electrons with ablated gas In the cloud On the pellet surface

11. 4. Equation of state with atomic processes (1) Dissociation and ionization fractions Saha equation for the dissociation and ionization Deuterium Ed=4.48eV, Nd=1.55×1024,ad=0.327 Ei=13.6eV, Ni=3.0×1021,ai=1.5

12. 4. Equation of state with atomic processes (2) High resolution solvers (based on the Riemann problem) require the sound speed and integrals of Riemann invariant type expressions along isentropes. Therefore the complete EOS is needed. • Conversions between thermodynamic variables are based on the solution of nonlinear Saha equations of (r,T). • To speedup solving Riemann problem, Riemann integrals pre-computed as functions of pressure along isentropes are stored in a 2D look-up table, and bi-linear interpolation is used. • Coupling with Redlich-Kwong EOS can improve accuracy at low temperatures.

13. 5. Conductivity model including ionization by impact Ionization by Impact

14. Influence of Atomic Processes on Temperature and Conductivity Temperature Conductivity

15. 6. Surface ablation model • Some facts: • The pellet is effectively shielded from incoming electrons by its ablation cloud • Processes in the ablation cloud define the ablation rate, not details of the phase transition on the pellet surface • No need to couple to acoustic waves in the solid/liquid pellet • The pellet surface is in the super-critical state • As a result, there is not even well defined phase boundary, vapor pressure etc. • This justifies the use of a simplified model: • Mass flux is given by the energy balance (incoming electron flux) at constant temperature • Pressure on the surface is defined through the connection to interior states by the Riemann wave curve • Density is found from the EOS.

16. 7. Pellet penetration through plasma pedestal

17. 8. Finite shielding length due to the curvature of B field • The grad-B drift curves the ablation channel away from the central pellet shadow. To mimic this 3D effect, we limit the extent of the ablation flow to a certain axial distance. • Without MHD effect, the cloud expansion is three-dimensional. The ablation rate reaches a finite value in the steady state. • With MHD effect, the cloud expansion is one-dimensional. The ablation rate would goes to zero by the ever increasing shielding if a finite shielding length were not in introduced.

18. Talk Outline • Equations for MHD at low magnetic Reynolds numbers and models for pellet ablation in a tokamak • Numerical algorithms for multiphase low ReM MHD • Numerical simulations of the pellet ablation in a tokamak

19. Multiphase MHD Solving MHD equations (a coupled hyperbolic – elliptic system) in geometrically complex, evolving domains subject to interface boundary conditions (which may include phase transition equations) • Material interfaces: • Discontinuity of density and physics properties (electrical conductivity) • Governed by the Riemann problem for MHD equations or phase transition equations

20. Main ideas of front tracking Front Tracking: A hybrid of Eulerian and Lagrangian methods • Two separate grids to describe the solution: • A volume filling rectangular mesh • An unstructured codimension-1 Lagrangian mesh to represent interface • Major components: • Front propagation and redistribution • Wave (smooth region) solution • Advantages of explicit interface tracking: • No numerical interfacial diffusion • Real physics models for interface propagation • Different physics / numerical approximations in domains separated by interfaces

21. Level-set vs. front tracking method Explicit tracking of interfaces preserves geometry and topology more accurately. 5th order level set (WENO) 4th order front tracking (Runge-Kutta)

22. The FronTier code FronTier is a parallel 3D multi-physics code based on front tracking • Physics models include • Compressible fluid dynamics • MHD • Flow in porous media • Elasto-plastic deformations • Realistic EOS models • Exact and approximate Riemann solvers • Phase transition models • Adaptive mesh refinement Interface untangling by the grid based method

23. Main FronTier applications Rayleigh-Taylor instability Richtmyer-Meshkov instability Supernova explosion Targets for future accelerators Tokamak refuelling through the ablation of frozen D2 pellets Liquid jet break-up and atomization

24. FronTier – MHD numerical scheme Elliptic step Hyperbolic step Point Shift (top) or Embedded Boundary (bottom) • Propagate interface • Untangle interface • Update interface states • Apply hyperbolic solvers • Update interior hydro states • Calculate electromagnetic fields • Update front and interior states • Generate finite element grid • Perform mixed finite element discretization • or • Perform finite volume discretization • Solve linear system using fast Poisson solvers

25. Hyperbolic step Interior and interface states for front tracking • Complex interfaces with topological changes in 2D and 3D • High resolution hyperbolic solvers • Riemann problem with Lorentz force • Ablation surface propagation • EOS for partially ionized gas and conductivity model • Hot electron heat deposition and Joule’s heating • Lorentz force and saturation numerical scheme • Centripetal force and evolution of rotational velocity

26. Elliptic step Embedded boundary elliptic solver • Based on the finite volume discretization • Domain boundary is embedded in the rectangular Cartesian grid. • The solution is always treated as a cell-centered quantity. • Using finite difference for full cell and linear interpolation for cut cell flux calculation • 2nd order accuracy For axisymmetric pellet ablation with transient radial current, the elliptic step can be skipped.

27. High Performance Computing • Software developed for parallel distributed memory supercomputers and clusters • Efficient parallelization • Scalability to thousands of processors • Code portability (used on Bluegene Supercomputers and various clusters) Bluegene/L Supercomputer (IBM) at Brookhaven National Laboratory

28. Talk Outline • Equations for MHD at low magnetic Reynolds numbers and models for pellet ablation in a tokamak • Numerical algorithms for multiphase low ReM MHD • Numerical simulations of the pellet ablation in a tokamak

29. Previous studies • Transonic Flow (TF) (or Neutral Gas Shielding) model, P. Parks & R. Turnbull, 1978 • Scaling of the ablation rate with the pellet radius and the plasma temperature and density • 1D steady state spherical hydrodynamics model • Neglected effects: Maxwellian hot electron distribution, geometric effects, atomic effects (dissociation, ionization), MHD, cloud charging and rotation • Claimed to be in good agreement with experiments • Theoretical model by B. Kuteev et al., 1985 • Maxwellian electron distribution • An attempt to account for the magnetic field induced heating asymmetry • Theoretical studies of MHD effects, P. Parks et al. • P2D code, A. K. MacAulay, 1994; CAP code R. Ishizaki, P. Parks, 2004 • Maxwellian hot electron distribution, axisymmetric ablation flow, atomic processes • MHD effects not considered

30. Our simulation results • Spherical model • Excellent agreement with TF model and Ishizaki • Axisymmetric pure hydro model • Double transonic structure • Geometric effect found to be minor • Plasma shielding • Subsonic ablation flow everywhere in the channel • Extended plasma shield reduces the ablation rate • Plasma shielding with cloud charging and rotation • Supersonic rotation widens ablation channel and increases ablation rate Spherical model Axis. hydro model Plasma shielding

31. 1. Spherically symmetric hydrodynamic simulation Polytropic EOS Plasma EOS Normalized ablation gas profiles at 10 microseconds • Excellent agreement with TF model and Ishizaki. • Verified scaling laws of the TF model

32. 2. Axially symmetric hydrodynamic simulation Steady-state ablation flow Temperature, eV Pressure, bar Mach number

33. 3. Axially symmetric MHD simulation (1) Main simulation parameters: Velocity distribution Channeling along magnetic field lines occurs at ~1.5 μs

34. 3. Axially symmetric MHD simulation (2) Mach number distribution Double transonic flow evolves to subsonic flow

35. Dependence on pedestal properties • Critical observation • Formation of the ablation channel and ablation rate strongly depends on plasma pedestal properties and pellet velocity. Simulations suggest that novel pellet acceleration technique (laser or gyrotron driven) are necessary for ITER. -.-.- tw = 5 ms, ne = 1.6  1013 cm-3 ___ tw = 10 ms, ne = 1014 cm-3 ----- tw = 10 ms, ne = 1.6  1013 cm-3

36. 4. MHD simulation with cloud charging and rotation (1) Supersonic rotation of the ablation channel Density redistribution in the ablation channel Steady-state pressure distribution in the widened ablation channel Isosurfaces of the rotational Mach number in the pellet ablation flow

37. 4. MHD simulation with cloud charging and rotation (2) Pellet ablation rate for ITER-type parameters G, g/s

38. 4. MHD simulation with cloud charging and rotation (3) Normalized potential along field lines Potential in the negative layer Channel radius and ablation rate • Grot is closer to the prediction of the quasisteady ablation model Gqs = 327 g/s • Magnetic β<<1 justifies the static B-field assumption

39. Striation instabilities • Current work focuses on the study of striation instabilities • Striation instabilities, observed in all experiments, are not well understood • We believe that the key process causing striation instabilities is the supersonic channel rotation, observed in our simulations Striation instabilities: Experimental observation (Courtesy MIT Fusion Group)

40. Plasma disruption mitigation Pressure distribution without rotation Gas ball R = 9 mm Killer pellet R = 9 mm

41. Plasma disruption mitigation Mach number distributions in the gas shell

42. Conclusions and future work • Developed MHD code for free surface low magnetic Re number flows • Front tracking method for multiphase flows • Elliptic problems in geometrically complex domains • Phase transition and surface ablation models • Axisymmetric simulations of pellet ablation • Effects of geometry, atomic processes, and conductivity model • Warm-up process and finite shielding length • Charging and rotation, transient radial current • Ablation rate, channel radius, and flow properties • Tracking of a shrinking pellet • Future work • 3D simulations of pellet ablation and striation instabilities • Asymptotic ablation properties in long warm up time • Natural cutoff shielding length • Magnetic induction • Systematic simulation of plasma disruption mitigation using killer pellet / gas ball • Coupling with global MHD models