MHD of Multiphase Flows: Numerical Algorithms and Applications Roman Samulyak

CSC SeminarMarch 9, 2006, BNL MHD of Multiphase Flows: Numerical Algorithms and Applications Roman Samulyak Center for Data Intensive Computing, BNL Collaborators: Tianshi Lu, CSC/BNL, modeling, software development, fusion applications Jian Du, Stony Brook University, software development, accelerator applications Zhiliang Xu, Xiaolin Li, CSC/SBU, front tracking methods Paul Parks, General Atomics, MHD theory, fusion applications James Glimm, SBU/CSC, modeling, numerical algorithms

Talk outline • Motivation • MHD system for free surface multiphase flow at low magnetic Reynolds numbers; approximations • Numerical algorithms for MHD equations • Applications • Numerical simulations of the pellet ablation in a tokamak

Motivation: Studies of Tokamak Fueling • Problems • Pellet ablation • Striation instabilities • Laser driven pellet acceleration • Gyrotron driven pellet acceleration • Liquid jet for plasma disruption mitigation Laser driven pellet acceleration ITER schematic Fueling using a high speed gaseous jet

Motivation: Accelerator Design (Neutrino Factory/Muon Collider target) • The target has been proposed as a free mercury jet interacting with an intensive proton pulse in a 20Tesla magnetic field Exparinemtal data (BNL), no B field • Future experiment at CERN (2007) will operate with • 20 m/s mercury jet • 15 Tesla pulsed solenoid • proton beam line

MHD equations and approximations Full system of MHD equations Low magnetic Re approximation • MHD regimes, approximations, and discretization techniques • Astrophysics: full system, explicit time discretization, 7-wave Riemann problem • Magnetically Confined Fusion plasmas: full system, implicit or semi-implicit discretization • Liquid metal MHD: incomplessible fluid, low magnetic Re, implicit time discretization • Our approach: compressible fluid (shock waves), low magnetic Re, phase transitions, free surface flows, explicit time discretization

Numerical Approach • The low megnetic Re MHD is a copuled hyperbolic/elliptic system. Operator splitting. • The hyperbolic subsystem is solved on a finite difference grid in both domains separated by the free surface using flont tracking numerical techniques. • Implemented in FronTier code • Riemann problem for interface propagation • Complex interfaces with topological changes in 2D and 3D • High resolution hyperbolic solvers • Realistic EOS models • The elliptic subsystem is solved in geometrically complex domains • Dynamic grid generation, finite element diecretization of conforming grids • Embedded boundary finite volume discretization • Fast parallel linear solvers

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

Stencil and equations for the interface point propagate

Schematic of the interface point propagate algorithm

Comparison of techniques for elliptic equations • Point-shift grid generation • Compatible with mixed finite element formulation • The same order accuracy for the potential and gradients • Capable of generating grids for vector finite elements • Not robust (especially in 3D) • Colella’s embedded boundary method • Second order accurate for the potential • Robust • Simple implementation for parallel computing • Due to arbitrary shape elements, incompatible with mixed finite element formulation

Embedded Boundary Elliptic Solver • Main Ideas • Based on the finite volume discretization • Potential is treated as cell centered value, even if the center is outside the computational domain • 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

Embedded Boundary Elliptic Solver • Area is corresponding to the computational cell on which to make • integration • fi is the flux across the non-boundary cell edges and ff is the boundary edge flux given by Neumann Conditions. All fluxes are calculated at the middle point of the edge and dl is the related edge length

Algorithm • Reconstruct the interface, record the crossings, and set the component type for grid point • Simplify the control volume configurations, remove small volume cells (better condition number) • Count the local number of blocks (both full and partial cells), set the matrix and vector dimension for the linear system solver,set global index for the counted blocks • For partial cells, the necessary information such as cell edge type, position of cell edge center and partial cell edge length are recorded. • A nine-point stencil is set for each partial cell. With the above information, related stencil value from flux integration is set and addeded into the corresponding matrix row • Set up the right hand side for the linear system, which is also calculated from flux integration by divergence theorem. Note that for partial cell, the boundary flux from both left and right side are cancelled out by Neumann boundary conditions. • Solve the linear system Ax = b with some fast parallel linear solver (PETSc or HYPRE)

Algorithm Stencil Setting Figure 3. Stencil for partial cells

Validation Test Function: Computational Domain Domain Boundary: Sin wave perturbed circle Convergence Rate:

Validation Illustration of flux error (X direction)

3D implementation Parallel 3D implementation has been completed and fully tested • Same principle as 2D • Bilinear interpolation of flux

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

Mercury jet entering magnetic field.Schematic of the problem. Magnetic field of the 15 T solenoid is given in the tabular format

Incompressible steady state formulation of the problem

Results: Aspect ratio of the jet cross-section B = 15 TV0 = 25 m/s

Consequences of the jet distortion • The cross-section of the mercury jet interaction with the proton pulse is significantly reduced • This reduces the pion production rate • In order to avoid these undesirable consequences, the angle between the magnetic field and the solenoid axes was reduced to 0.3 rad. This implies some new constraints on the hardware design • Another possible solution was the use of an elliptic nozzle to compensate the MHD distortion. This option has not been explored due to time constraints on the final design.

Processes in the ablation cloud

Pellet ablation in tokamaks: previous works. Local models • Neutral gas shielding (NGS) 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 • Neglected effects: MHD, geometric effects, atomic effects (dissociation, ionization) • Theoretical studies of MHD effects, P. Parks • 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

Pellet ablation in tokamaks: previous works. Global models • 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

Improved model is needed: • Studies of the local pellet ablation physics are still missing • MHD • 3D effects • Model for currents in the ablation cluod • 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

Physics Models for Pellet Studies : Electron Energy Deposition In the cloud: On the pellet surface:

Physics Models for Pellet Studies : Surface Ablation • 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.

Physics Models for Pellet Studies : Equation of State with Atomic Processes. 1 Saha equation for the dissociation (ionization) fraction

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. Second law of thermodynamics: We found the complete EOS and showed that compatibility with the second law of thermodynamics requires:

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.

Conductivity of weakly ionized plasmas

EOS with Atomic Processes. Numerical Results

Real gas EOS (future work) • Gas (deuterium vapor) near the pellet surface is cold and dense • Ideal EOS model is not accurate • A very good approximation is given by the Redlich-Kwong EOS: • We propose an extension of the Redlich-Kwong EOS to include atomic processes (dissociation and ionization) • The EOS will contain 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. • Challenges to derive the complete EOS.

Simulation results B = 0 = 1 = 2.5 = 5 Tesla Formation of pellet ablation channels in magnetic fields ranging from 0 to 5 Tesla

Simulation results B = 2.5 Tesla B = 5.0 Tesla P T P T Developed ablation flow channels in magnetic fields of 2.5 and 5 Tesla.

Numerical results. Spherical model G = 118 g/sec

Numerical results. 2D axisymmetric model Pressure distribution at t = 2 microseconds Polytropic EOS Plasma EOS

Numerical results. 2D axisymmetric model

Striation instabilities • Striation instabilities have been observed in all pellet experiments; Important problem for the plasma confinement • Fluctuations of emitted light due to the cloud pellet separation; fingering of the ablated material • Several theories have been proposed, but they don’t satisfactory explain the phenomenon • Idea: Striation instabilities are caused by the Rayleigh-Taylor instability due to the cloud rotation. Linear model of the E x B rotation by Parks (1996). • Nonlinear studies (numerical simulations) are necessary

Let’s add missing physics 3D detailed physics model is critical for the study of cloud rotation, pellet – cloud separation, and striation instabilities. Work on the 3D model is in progress.

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 • Numerical simulations of the Muon Collider target in magnetic fields • 2D numerical simulations of the pellet ablation • Axially symmetric ablation flow • Kinetic calculation of the electron heat deposition, atomic physics processes • MHD effects; formations of ablation channels • Calculated pellet ablation rates, channel properties • Future work • 3D simulations of the pellet ablation and striation instabilities (3D detailed physics model is critical for the study of cloud rotation, pellet – cloud separation, and striation instabilities) • Laser -- plasma interaction model with sharp absorption front • Laser acceleration of pellets

MHD of Multiphase Flows: Numerical Algorithms and Applications Roman Samulyak