THE OFFICE OF NONPROLIFERATION

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
Ohms law
Maxwells 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
Schematic of pellet ablation in a magnetic field
Schematic of processes in the ablation cloud
7
Local model for pellet ablation in tokamak

1. Axisymmetric MHD with low ReM approximation
2. Transient radial current approximation
3. Interaction of hot electrons with ablated gas
4. Equation of state with atomic processes
5. Conductivity model including ionization by
   electron impact
6. Surface ablation model
7. Pellet penetration through plasma pedestal
8. 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)
Saha equation for the dissociation and ionization
Dissociation and ionization fractions
Deuterium Ed4.48eV, Nd1.551024,ad0.327 Ei13.6
eV, Ni3.01021,ai1.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
Supernova explosion
Richtmyer-Meshkov instability
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 Joules 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

Poly EOS Plasma EOS
Sonic radius 0.66 cm 0.45 cm
Temperature 5.51 eV 1.07 eV
Pressure 20.0 bar 26.9 bar
Ablation rate 112 g/s 106 g/s
32
2. Axially symmetric hydrodynamic simulation
Steady-state ablation flow
Temperature, eV
Pressure, bar
Mach number
33
3. Axially symmetric MHD simulation (1)
Plasma electron temperature Te 2 keV
Plasma electron density ne 1014 cm-3(standard) 1.6x1013 cm-3(el. shielding)
Warm-up time tw 5 20 microseconds
Magnetic field B 2 6 Tesla
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

• 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 ßltlt1 justifies the static B-field
  assumption

Channel radius Ablation rate ?B/B ? b/2
Non-rotating 2.3 cm 195 g/s 0.079
Rotating 2.8 cm 262 g/s 0.088
39

• 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