Simulation of Multiphase Magnetohydrodynamic Flows for Nuclear Fusion Applications Roman Samulyak

1. The Role of High Performance Computation in Economic Development Workshop at Rensselaer Polytechnic Institute, October 22 -24, 2008 Simulation of Multiphase Magnetohydrodynamic Flows forNuclear Fusion Applications Roman Samulyak AMS Department, Stony Brook University and Computational Science Center Brookhaven National Laboratory Collaborators: Tianshi Lu (BNL), Paul Parks (General Atomics) Wurigen Bo, Lingling Wu (Stony Brook) James Glimm, Xiaolin Li (Stony Brook) Kirk McDonald (Princeton), Harold Kirk (BNL)

2. Talk Outline • Brief description of wide range of applications • Numerical algorithms for multiphase hydro and MHD flows • Large scale computing (BlueGene) • Application specific mathematical models and numerical algorithms (for one selected applications - tokamak fueling) • Numerical simulations: main results and feedback to applications • Summary and conclusions

3. Fusion Energy. ITER project: fuel pellet ablation • ITER is a joint international research and development project that aims to demonstrate the scientific and technical feasibility of fusion power • ITER will be constructed in Europe, at Cadarache in the South of France in ~10 years Our contribution to ITER science: Models and simulations of tokamak fueling through the ablation of frozen D2 pellets Collaboration with General Atomics Laser driven pellet acceleration

4. New Ideas in Nuclear Fusion: Magnetized Target Fusion (MTF)

5. Fission Energy. DOE Nuclear Energy Research Initiative (NERI) • Stony Brook University, together with Rensselaer Polytechnical Institute, Columbia University, and Brookhaven National Laboratory, joined a DOE funded consortium within NERI initiative • Project title: Deployment of a Suite of High Performance Computational Tools for Multiscale Multiphysics Simulation of Generation-IV Reactors

6. Numerical Studies within NERI Consortium • Simulation of material relocation in a Core Disruptive Accident (CDA) in Gen. IV reactors • Numerical investigation of the surface instability and liquid entrainment in the expansion of an CDA bubble. • Studies of some fundamental problems of multiphase flows • Subgrid models for multiphase flows

7. Spallation Neutron Source: Cavitation induced erosion in mercury target Left: pressure distribution in the SNS target prototype. Right: Cavitation induced pitting of the target flange (Los Alamos experiments) • Injection of nondissolvable gas bubbles has been proposed as a pressure mitigation technique. • Numerical simulations aim to estimate the efficiency of this approach, explore different flow regimes, and optimize parameters of the system.

8. Mercury Jet Target for Neutrino Factory / Muon Collider • 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 disruptions Target schematic

9. Numerical Algorithms • Common features of all applications: • Multiphase / free surface hydro and MHD flows interacting with external energy sources • Phase transitions • They require: • Explicitly resolved multiphase (free surface) flows • Our choice is Front Tracking • Numerical algorithms for coupled hyperbolic - elliptic systems in geometrically complex domains • Large scale computing

10. MHD equations and approximations Full system of MHD equations Low magnetic Re approximation

11. 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

12. 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

13. The FronTierCode (SciDAC ITAPS Software) FronTier is a parallel 3D multiphysics code based on front tracking • Physics models include • Compressible fluid dynamics • MHD • Flow in porous media • Elasto-plastic deformations • Realistic EOS models, phase transition models • Exact and approximate Riemann solvers • Adaptive mesh refinement Turbulent fluid mixing. Left: 2D Right: 3D (fragment of the interface)

14. FronTier-MHD numerical scheme Elliptic step Hyperbolic step Point Shift (top) or Embedded Boundary (bottom) • Propagate interface • (solve Riemann problem for contact or phase transition equations) • 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

15. Normal propagation of interface points Contact discontinuity

16. Phase boundary problem Interface jump conditions Balance equations Temperature and pressure at the interface Subgrid model for temperature

17. 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 • Using finite difference for full cell and linear interpolation for cut cell flux calculation • Advantage: robust, readily parallelizable, compatible with FronTier grid-based interface tracking algorithm.

18. High Performance Computing • We are interested in the development of software for parallel distributed memory supercomputers • Efficient parallelization • Scalability to thousands of processors • State-of-art parallel visualization • Code portability (for use of local and external (NERSC) computational resources)

19. Application Specific Models and Simulations

20. Pellet Ablation for Tokamak Fueling: Main Models • Equation of state with atomic processes • Kinetic model for the interaction of hot electrons with the ablated gas • Surface ablation model • Cloud charging and rotation models • New conductivity model (ionization by electron impact) Schematic of processes in the ablation cloud

21. Relation to Other Projects • Macroscale Model: AMR simulation of ablation flow in plasma (Ravi Samtaney) • Focus on plasma flow • Ablation physics not resolved (simplified analytical models) • Local (“Microscale”) Model: Front Tracking simulation of pellet ablation • Focus on detailed ablation physics, ablation rates etc. • Far field plasma evolution not resolved

22. Equation of State with Atomic Processes. Saha equation for the dissociation (ionization) fraction

23. 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:

24. 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

25. 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.

26. Influence of Atomic Processes on Temperature and Conductivity

27. Electron Energy Deposition In the cloud: On the pellet surface:

28. Pellet charging model

29. 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

30. 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

31. 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. Clarified the role of directional electron beam heating on the ablation rate.

32. Performed first systematic studies of pellet ablation rates in magnetic fields 3 microseconds 1 microsecond 5 microseconds 9 microseconds 3 microseconds

33. Cloud Rotation • Revealed new propertied of the ablation flow: Supersonic rotation of the ablation channel • Supersonic rotation widens the ablation channel, re-distributes density, and changes the ablation rate • Resolution of this phenomenon greatly improves the agreement with experiments Isosurfaces of the rotational Mach number in the pellet ablation flow

34. 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 Schematic of the plasma pedestal. Plasma pedestal provides critical influence to the pellet ablation rate

35. Work in Progress • 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)

36. Neutrino Factory / Muon Collider target • Simulation tasks: • Studies of the jet formation in the nozzle • Entrance of the jet into the magnetic field • Mercury jet - proton pulse interaction. Jet cavitation and disruption, stabilizing effect of the magnetic field • Predictions for the CERN targetry experiment called MERIT (Fall of 2007)

37. Jet entering 15 T solenoid

38. Aspect ratio of the jet cross-section B = 15 TV0 = 25 m/s These studies has resulted in the change of MERIT design. The chamber was re-designed to allow much smaller angle with of the jet with the magnetic field.

39. Experimental data V = 15 m/s, B = 10T V = 20 m/s, B = 10T B = 15T B = 15T Simulations predicted slightly smaller value of the jet width in different view ports.

40. Liquid jet cavitation

41. Interaction of the mercury jet target with proton pulses 3D FronTier simulations of the disruption of Hg jet interacting with a proton pulse.

42. Summary of Target Studies • Accomplishments • Modeling and simulation of the mercury jet target • Change design parameters of the MERIT experiments • Successful predictions for MERIT: jet flattening in high gradient field, jet disruption by proton pulses, stabilizing effect of the magnetic field • Significance • Neutrino Factory Collaboration reached the conclusion that liquid mercury jet targets can successfully work with proton beams up to 8 MW