# Presentation: 178689

## MHD of Multiphase Flows: Numerical Algorithms and Applications

**Roman Samulyak**
*Center for Data Intensive Computing, BNL*

**CSC Seminar****March 9, 2006, 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

- Fueling using a high speed gaseous jet

- ITER schematic

## Motivation: Accelerator Design

**Motivation: Accelerator Design**
**(Neutrino Factory/Muon Collider target)**

- Exparinemtal data (BNL), no B field

-  The target has been proposed as a free mercury jet interacting with an intensive proton pulse in a 20Tesla magnetic 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

**MHD  equations and approximations**

**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

**Full system of MHD equations**

**Low magnetic Re approximation**

## 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

**Hyperbolic step**

**FronTier-MHD numerical scheme**

**Elliptic step**

-  Propagate interface
-  Untangle interface
-  Update interface states

-  Apply hyperbolic solvers
-  Update interior hydro states

-  Generate finite element grid
-  Perform mixed finite element discretization
- or
-  Perform finite volume discretization
-  Solve linear system using fast Poisson solvers

-  Calculate electromagnetic fields
-  Update front and interior states

- Point Shift (top) or Embedded Boundary (bottom)

## Stencil and equations for the interface point propagate 

**Stencil and equations for the interface point propagate **

## Schematic of the interface point propagate algorithm

**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
* **f* i   is the flux  across the non-boundary cell edges and  *f* f  is the boundary edge flux given by Neumann Conditions.  All fluxes are calculated at the middle point of the edge and d*l *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:**

**Domain Boundary:       Sin wave perturbed circle**

**Convergence  Rate:**

| Mesh Size |    error (fx) |       R |
| --- | --- | --- |
| 32x32 | 2.496e-2 | N/A |
| 64x64 | 8.841e-3 |      1.497 |
| 128x128 | 3.093e-3 |      1.506 |
| 256x256 | 1.096e-3 |      1.503 |

**Computational Domain**

## Validation

- Illustration of flux error (X direction)

## 3D implementation

-  Same principle as 2D
-  Bilinear interpolation of flux

- Parallel 3D implementation has been completed and fully tested

## Muon Collider target: a brief summary of modeling and simulation 

**Muon Collider target: a brief summary of modeling and simulation **

- Jet surface instabilities

- Cavitation in the mercury jet and thimble

-  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

## Mercury jet entering magnetic field.Schematic of the problem.

**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

**Incompressible steady state formulation of the problem**

## Results: Aspect ratio of the jet cross-section

**Results: Aspect ratio of the jet cross-section**

**B = 15 T****V0 = 25 m/s**

## Results: Aspect ratio of the jet cross-section

**Results: Aspect ratio of the jet cross-section**

**B = 15 T****V0 = 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

**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: 

**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 :

- Saha equation for the dissociation (ionization) fraction

**Physics Models for Pellet Studies :**
**Equation of State with Atomic Processes. 1**

## EOS with Atomic Processes 

**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 

**Numerical Algorithms for EOS **

-  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

**Conductivity of weakly ionized plasmas**

## EOS with Atomic Processes. Numerical Results 

**EOS with Atomic Processes. Numerical Results **

## EOS with Atomic Processes. Numerical Results 

**EOS with Atomic Processes. Numerical Results **

## Real gas EOS (future work) 

-  Ideal EOS model is not accurate
-  A very good approximation is given by the Redlich-Kwong EOS:

**Real gas EOS (future work) **

-  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 

**Simulation results **

- B = 0                   = 1           = 2.5        = 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.

**Simulation results **

## Numerical results. Spherical model 

**Numerical results. Spherical model **

- G = 118 g/sec

## Numerical results. Spherical model 

**Numerical results. Spherical model **

- G = 110 g/sec

## Numerical results. Spherical model 

**Numerical results. Spherical model **

- G = 105 g/sec

## Numerical results. 2D axisymmetric model 

- Polytropic EOS          Plasma EOS

**Numerical results. 2D axisymmetric model **

- Pressure distribution
- at t = 2 microseconds

## Numerical results. 2D axisymmetric model 

**Numerical results. 2D axisymmetric model **

## Striation instabilities 

-  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

**Striation instabilities **

## Let’s add missing physics 

**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