THE OFFICE OF NONPROLIFERATION

1
CPPG Seminar,Princeton Plasma Physics
LaboratoryJuly 26, 2006, Princeton, NJ
Axisymmetric MHD Simulation of Pellet
Ablation Roman Samulyak Computational Science
Center Brookhaven National Laboratory
Tianshi Lu CSC/BNL, modeling, software
development, fusion applications Paul
Parks General Atomics, MHD theory, fusion
applications Hydrodynamic and MHD Multiphase
Flows James Glimm, Stony Brook University / BNL,
modeling, numerical algorithms Jian Du, Stony
Brook University, software development,
accelerator applications Zhiliang Xu, Xiaolin Li,
CSC/SBU, front tracking methods
2
Pellet Ablation in the Process of Tokamak Fueling
Detailed studies of the pellet ablation physics
(local models)
Studies of the tokamak plasma in the presence of
an ablating pellet (global models)
ITER schematic
3
Previous Studies Local Models

• Transonic Flow (TF) (or Neutral gas shielding)
  model, P. Parks R. Turnbull, 1978
• Provides the scaling of the ablation rate with
  the pellet radius and the plasma temperature and
  density
• 1D steady state hydrodynamics model,
  monoenergetic electron distribution
• Neglected effects MHD, geometric effects,
  atomic effects (dissociation, ionization)
• 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
• Non-spherical ablation flow (axial symmetry),
  proper treatment of scattering
• Kinetic calculation of the electron heat
  deposition, atomic physics processes
• MHD effects not considered

4
Previous Studies Global Models (examples)

• Simulations using MH3D code, H. Strauss W.
  Park, 1998
• Finite element version of the MH3D full MHD code
• Details of the ablation are not considered
• Pellet is given as a density perturbation of
  initial conditions
• Smaller values of density and larger pellet
  radius (numerical constraints)
• Simulations using MHD code based on CHOMBO AMR
  package, R. Samtaney, S. Jardin, P. Colella, D.
  Martin, 2004
• Analytical model for the pellet ablation moving
  density source
• 8-wave upwinding unsplit method for MHD
• AMR package significant improvement of
  numerical resolution

5
Improved model is needed

• Studies of the local pellet ablation physics
  were still missing
• MHD
• 3D effects
• Charging models for the ablation cloud
• Global plasma simulations in the presence of an
  ablating pellet need a better local model as
  input
• Studies of striation instabilities, observed in
  all experiments, are not possible without a 3D
  detailed physics model
• We are working on building and validations of
  such models

6
Pellet Ablation Model Based on Front Tracking

• Explicitly tracked interfaces resolution of
  material properties and multiple scales
• MHD equations in the low magnetic Reynolds
  number approximation and numerical methods for
  free surface flows
• Equation of state with atomic processes
• Kinetic model for the interaction of hot
  electrons with the ablated gas
• Surface ablation model

7
MHD equations and approximations
Full system of MHD equations
Low magnetic Re approximation
Assumptions that will be verified later. Near
the pellet
Downstream in the ablation channel
8
Simplification of the elliptic equation for the
pellet problem
Axial symmetry
Pellet cloud charging, polarization currents, and
the axial rotation of the ablation channel are
neglected in this study (implementation of the
charging and rotation of the channel is in
progress and will be discussed later)
Magnetic field is constant B (0,0,Bz)
Therefore,
and the only non-zero component of the current is
9
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

10
Stencil and equations for the interface point
propagate
11
Schematic of the interface point propagate
algorithm
12
Embedded Boundary Elliptic Solver

• Main Ideas
• Based on the finite volume discretization
• Domain boundary is embedded in the rectangular
  Cartesian grid, and the solution is treated as a
  cell-centered quantity
• The discretized operator is centered in centroids
  of partial cells
• Using finite difference for full cell and linear
  interpolation for cut cell flux calculation

13
Validation of the Elliptical Problem
Let
be a solution of
Derive the R.H.S. and B.C., solve numerically and
compare with the exact solution.
2D problem
Mesh size Error Conv. Rate CPU time Iterations
64x64 9.09e-05 N/A 0.087 44
128x128 2.01e-05 2.175 0.389 98
256x256 4.80e-06 2.122 2.223 264
512x512 1.78e-06 1.893 15.445 500
3D problem
Mesh size Error Conv. Rate Iterations
32x32x32 1.32e-03 N/A 42
64x64x64 3.18e-04 2.050 76
128x128x128 8.05e-05 2.016 144
14
Validation of the MHD Code
A free mercury jet travels longitudinally along
the z axis in a magnetic field (0,By,0)
A satisfactory perturbation theory and
experimental data are available (S. Oshima et
al., JCME Int. J., 30 (1987), No. 261.
Solid line theory Dots simulations
15
Muon Collider target a brief summary of modeling
and simulation

• Simulation of the mercury jet target interacting
  with a proton pulse in a magnetic field
• Studies of surface instabilities, jet breakup,
  and cavitation
• MHD forces reduce both jet expansion,
  instabilities, and cavitation

Jet surface instabilities
Cavitation in the mercury jet and thimble
16
Equation of State with Atomic Processes.
Saha equation for the dissociation (ionization)
fraction
17
EOS with Atomic Processes
Incomplete EOS (known from literature)
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.
Using the second law of thermodynamics
we found the complete EOS and showed that the
compatibility with the second law of
thermodynamics requires
18
Complete EOS with Atomic Processes
Notations
We will define the sound speed in a form typical
for the polytropic gas
where the effective gamma is
the Gruneisen coefficient is
and the entropy is
19
Numerical Algorithms for EOS
For better numerical efficiency, FronTier
operates with three pairs of independent
thermodynamic variables

• For the first two pairs of variables, solve
  numerically nonlinear algebraic equation, and
  find T. Using , find the remaining
  state.
• Such an approach is prohibitively slow for the
  calculation of Riemann integrals (involves nested
  nonlinear equations).
• To speedup calculations, we precompute and store
  values on Riemann integrals as functions of the
  density and entropy. Two dimensional table lookup
  and bi-linear interpolation are used.

20
Influence of Atomic Processes on Temperature and
Conductivity
21
Real gas EOS (work in progress)
Redlich-Kwong EOS for the cold and dense gas

• We have derived an extension of the
  Redlich-Kwong EOS to include atomic processes
  (dissociation and ionization)
• The EOS contains three terms the partial
  pressure/energy of the molecular gas is written
  in the Redlich-Kwong form, and the partial
  pressure/energies of the dissociated and ionized
  components is written in the ideal EOS form.
• Complete EOS has been derived (expressions for
  entropy, sound speed etc.)
• The numerical implementation is in progress

22
Electron Energy Deposition
In the cloud
On the pellet surface
23
Physics Models for Pellet Studies Surface
Ablation

• Features of the pellet ablation
• 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

24
Main Simulation Parameters
Pellet radius rp 2 mm
Pellet density 0.2 g/cm3
Plasma electron temperature Te 2 keV
Plasma electron density ne 1014 cm-3(standard) 1.6x1013 cm-3(electrostatic shielding)
Length of the ablation channel Lc 15 cm
Warm-up time 5 20 microseconds
Magnetic field B 2 6 Tesla
25
Three types of simulation
1D hydrodynamic model spherically symmetric, no
JxB force 2D hydrodynamic model axially
symmetric, directional heating along magnetic
field lines, no JxB force 2D MHD model axially
symmetric, directional heating along magnetic
field lines, JxB force is applied
26
We will compare results with
TF model P. Parks and R. Turnbull, Phys. Fluids,
21 (1978), 1735. Kuteev B. V. Kuteev, Sov. J.
Plasma Phys, 11 (1985), 236. MacAulay A. K.
MacAulay, Nucl. Fusion, 34 (1994), 43. Parks
2000 P. Parks, W. Sessions, L. Baylor, Phys.
Plasmas., 5 (2000), 1968 Ishizaki R. Ishizaki,
P. Parks, N. Nakajiama, M. Okamoto, Phys.
Plasmas, 11 (2004), 4064
27
Spherically symmetric 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
28
Axially Symmetric Hydrodynamic Simulation
Temperature, eV
Pressure, bar
Mach number
Distributions of temperature, pressure, and Mach
number of the ablation flow near the pellet at 20
microseconds.
29
Axially Symmetric Hydrodynamic Simulation
Temperature, pressure, and Mach number of the
ablation flow near the pellet in the longitudinal
(solid line) and radial (dashed line) directions
at 20 microseconds.
30
Reduction of the ablation rate in 2D

• The ablation rate, 90 g/s, agrees with Kuteev
  and MacAulay, disagrees with Ishizaki
• The reduction of the ablation compared to the
  spherically symmetric case is 18.
• This disagrees with prevailing expectations of
  the factor of 2 reduction (Kuteev and Ishizaki)
• However Kuteev did not calculate the 1D ablation
  rate he compared with the FT model that used the
  monoenergetic electron distribution, and found a
  2.2 reduction
• Our explanation of the factor of 2.2 reduction
• Maxwellian heat flux increases the 1D ablation
  rate by 2.75 (Ishizaki)
• Directional heat flux in 2D reduces the ablation
  rate by 0.82 (this work)
• 2.75 x 0.82 2.25

31
(No Transcript)
32
Mach number distribution of the ablation flow
near the pellet in 6 Tesla magnetic field. Warm
up time is 20 microreconds.
3 microseconds
5 microseconds
9 microseconds
33
Temperature (eV) distribution of the ablation
flow near the pellet in 6 Tesla magnetic field.
Warm up time is 20 microseconds.
3 microseconds
5 microseconds
9 microseconds
34
Pressure (bar) distribution of the ablation flow
near the pellet in 6 Tesla magnetic field. Warm
up time is 20 microseconds.
3 microseconds
5 microseconds
5 microseconds
35
Pressure along the z-axis of steady state
ablation channel.
No shielding
Electrostatic shielding
Solid line 2 Tesla, dashed line 4 Tesla,
dotted line 6 Tesla. Warm up time is 10
microseconds.
36
Temperature along the z-axis of steady state
ablation channel.
No shielding
Electrostatic shielding
Solid line 2 Tesla, dashed line 4 Tesla,
dotted line 6 Tesla. Warm up time is 10
microseconds.
37
Mach number along the z-axis of steady state
ablation channel
No shielding
Electrostatic shielding
Solid line 2 Tesla, dashed line 4 Tesla,
dotted line 6 Tesla. Warm up time is 10
microseconds.
38
Normalized temperature and pressure in the
ablation channel
Solid lines Simulation in 2 Tesla field with
electrostatic shielding. Dashed lines
Theoretical parallel flow model for the ablation
channel (Parks 2000)
39
Radius of the ablation channel
Solid line 10 microseconds warm up time, ne
1.0e14 cm-1 Dashed line 10 microseconds warm
up time, ne 1.6e13 cm-1 Dotted line 5
microseconds warm up time, ne 10e14 cm-1
40
Density along the axis of symmetry and the
ablation rate
Solid line MHD model, B 6 Tesla,
ne 1.0e14 cm-1 Dashed
line MHD model, B 2 Tesla,
ne 1.6e13 cm-1 Dotted line 1D
spherically symmetric model
Solid line tw 10, ne 1.0e14 cm-1 Dashed
line tw 10, ne 1.6e13 cm-1 Dotted line
tw 5, ne 10e14 cm-1
41
Justification of the Electrostatic Approximation
B ne 2 Tesla 2 Tesla 4 Tesla 4 Tesla 6 Tesla 6 Tesla
B ne
1014cm-3 0.530 0.822 0.128 0.223 0.051 0.100
1.6x1013cm-3 0.110 0.180 0.029 0.055 0.015 0.026
The ratio of induced magnetic field to the
toroidal magnetic field.
42
Current work implementation of the pellet
charging model
43
Other Low Magnetic Reynolds Number Flows in
Tokamaks

• Laser driven pellet acceleration
• Gyrotron driven pellet acceleration
• Liquid jet for plasma disruption mitigation
• High density, low temperature gas jets for
  tokamak fueling (inefficient method?)

Laser driven pellet acceleration
Fueling using a high speed gaseous jet
44
Conclusions and future work

• Developed MHD code for free surface low magnetic
  Re number flows
• Front tracking method multiphase/multimaterial
  flows
• Elliptic problems in geometrically complex
  domains
• Phase transition models
• Performed numerical simulation of tokamak
  fueling through the injection of frozen deuterium
  - tritium pellets
• Computed ablation rates in hydro and MHD case
• Explained the factor of 2 reduction of the
  ablation rate
• Performed first systematic studies of the
  ablation in magnetic fields
• Future work
• 3D simulations of the pellet ablation
• Studies of striation instabilities
• Coupling our pellet ablation model as a subgrid
  model with a tokamak plasma simulation code
• Laser -- plasma interaction model with sharp
  absorption front, laser acceleration of pellets
  (possible)