Computational

An interface-tracking numerical algorithm for the simulation of magnetohydrodynamic multiphase / free surface flows in the low-magnetic-Reynolds-number approximation of (Samulyak R., Du J., Glimm J., Xu Z., J. Comp. Phys., 2007, 226, 1532) is described. The algorithm has been implemented in multi-physics code FronTier and used for the simulation of MHD processes in liquids and weakly ionized plasmas. In this paper, numerical simulations of a liquid mercury jet entering strong and nonuniform magnetic field and interacting with a powerful proton pulse have been performed and compared with experiments. Such a mercury jet is a prototype of the proposed Muon Collider / Neutrino Factory, a future particle accelerator. Simulations demonstrate the elliptic distortion of the mercury jet as it enters the magnetic solenoid at a small angle to the magnetic axis, jet-surface instabilities (filamentation) induced by the interaction with proton pulses, and the stabilizing effect of the magnetic field.


Introduction
The main driver for computational magnetohydrodynamics is the research in magnetically confined nuclear fusion which has recently been boosted by the International Thermonuclear Experimental Reactor project (ITER, [1]).Numerical algorithms for nuclear-fusion-simulation research are optimized for highly conductive, fully ionized plasmas.In simplified studies, even the infiniteconductivity approximation is widely used (ideal MHD approximation [2]).But even in thermonuclear fusion devices such as tokamaks, low-conductivity, weakly ionized plasma may still be present under special conditions despite very high temperature (the nuclear fusion ignition temperature is of the order of 10 keV or 10 8 K).In tokamaks, weakly ionized plasma is found in the ablated clouds of cryogenic fuel pellets.The injection of such frozen deuterium-tritium pellets is considered the most efficient technique for the fueling of tokamaks [3].Numerical algorithms designed for highly conductive plasma have very limited applicability to weakly ionized gases.Another large area of application of low-conductity MHD span liquid metals and liquid salts.Liquid metals, such as lithium, can be used in tokamaks as a plasma-facing component protecting tokamak walls.In this paper, we are exploring another application area of algorithms for liquid-metal MHD: the study of liquid-mercury-jet targets for future particle accelerators.
In many applications, liquid metals or weakly ionized plasmas have either free surfaces or moving, geometrically complex interfaces with either non-conductive or highly conductive media.The presence of interfaces imposes a major challenge on the development of high-quality numerical algorithms.The majority of numerical studies of free-surface MHD flows are based on semi-analytical treatment of simplified flow regimes [4,5].Simplified models have successfully been used for the description of self-organized filaments in dielectric barrier-glow discharges [6,7] and the numerical modeling of micro-plasma instabilities [8].But analytical models have a limited applicability for complex systems involving strongly coupled multi-physics phenomena.
The numerical resolution of interfaces can be performed using various methods developed for multiphase fluid dynamics.In our recent work [10], a free-surface MHD algorithm was developed based on explicit tracking of material interfaces.To the best of our knowledge, the only other fully numerical treatment of general, free-surface, incompressible liquid flows is implemented in the HIMAG code [11] using the level-set algorithm for fluid interfaces, the electric-potential formulation for electromagnetic forces, and the incompressible-fluid-flow approximation.The key feature of our algorithm is the use of the method of front tracking [12] for the propagation of vapor interfaces.The FronTier code is capable of tracking and resolving topological changes of large numbers of interfaces in two-and three-dimensional spaces [13].In the method of front tracking, the interface is a Lagrangian mesh moving through a volume-filling rectangular mesh according to the solution of the corresponding Riemann problem.High-resolution, shock-capturing, Godunov-type solvers are used to update hyperbolic states in the interior away from interfaces.The explicit treatment of interfaces typical of the method of front tracking greatly reduces the numerical diffusion and is especially advantageous for multi-physics problems involving phase transitions.It allows not only to solve accurately the Riemann problem for the phase boundary, but also to apply different mathematical approximations in the regions separated by interfaces to account for different material properties and, if necessary, eliminate fast time scales in numerical simulations.The three-dimensional, fronttracking algorithm has been recently upgraded with a so-called locally-grid-based tracking [14] that enabled robust resolution of very complex topological changes in the FronTier code.The comparison of front tracking and other numerical methods resolving interfaces, in particular the level-set method, can be found in [13].The main features of front tracking such as the conservative property and the ability to resolve complex interfaces are also present in particle methods such as the smoothed particle hydrodynamics [9].However the particle-based codes generally require longer computational time compared to grid-based codes.
Since most of our applications involve strong hydrodynamic waves, we solve the system of compressible Navier-Stokes equations coupled to electromagnetic forces.A numerical algorithm for the reduced MHD system by using the incompressible hydrodynamic approximations will be described in a forthcoming paper.
The front-tracking MHD algorithm has been used for the simulation of tokamak fueling by the ablation of cryogenic pellets [15,16] and the dynamics of laser-ablation plumes [17].In this paper, we apply the numerical MHD algorithms for the simulation of liquid mercury jet targets for the Neutrino Factory/Muon Collider.In particle accelerators, liquid metal targets convert intense proton beams into neutrons, pions or other particles used in fundamental and applied studies.The aim of the targetry group of the international, multi-institutional Neutrino Factory/Muon Collider Collaboration [18] is to explore the feasibility of high power targets for future particle accelerators and, in particular, for the proposed Neutrino Factory and Muon Collider.The target will contain a series of mercury jet pulses of about 1 cm in diameter and 30 cm in length.Each pulse will be shot at a velocity of 20 m/s into a 15-20-T magnetic field at a small angle to the axis of the magnetic field.When the jet reaches the center of the magnet, it interacts with a series of proton pulses depositing about 100 J/g of energy in the mercury.The aim of numerical simulations is to describe the hydrodynamic response of the target interacting with proton pulses in magnetic fields, the understanding of which is of major importance for reliable target design.
We have already performed mathematical modeling, software development and simulations of liquid mercury jet targets interacting with high-power proton beams in magnetic fields [19,20] and made predictions for the targetry experiment at CERN, so-called MERIT.The simulation predicted strong instabilities and cavitation of the mercury jet interacting with proton pulses at zero magnetic field and their reduction / stabilizing by a longitudinal magnetic field.Simulation predictions were qualitatively confirmed by the MERIT experiment conducted in CERN in the fall of 2007 [21].In this paper, we use exact input parameters of the MERIT experiment and quantitatively compare the simulations with experimental data.
The paper is organized as follows.The main governing equations are presented in section 2.
In section 3, we describe the front-tracking numerical algorithms for free-surface MHD and the algorithm for fluid cavitation.Numerical simulations of the mercury-jet target for the Neutrino Factory / Muon Collider are described in section 4. We conclude this paper with the summary of our results and perspectives for the future work.

Governing equations
We are interested in the description of multiphase or multi-material systems involving conducting fluids or weakly ionized gases interacting with neutral fluids or gases in the presence of magnetic fields.Interfaces of the phase or material separation are assumed to be sharp (the thickness of the interface is negligible) and, in general, geometrically complex.The numerical simulation of liquid metals and liquid salts as well as weakly ionized plasmas, which are relatively weak electrical conductors, is difficult using the standard full systems of MHD equations [2].Fast diffusion of the magnetic field, caused by low value of electrical conductivity, introduces unwanted small time scales into the problem.If the time scale of the diffusion of the magnetic field is small compared to hydrodynamic time scale, the magnetic Reynolds number [4] Re where L is the typical length scale, u is the fluid velocity, and σ is the electric conductivity, is small.If in addition the eddy-current-induced magnetic field δB is small compared to the external field B, the full system of MHD equations can be simplified by neglecting the time evolution of the magnetic field.In this case, the generalized Ohm's law is used for the calculation of the current-density distribution instead of the Maxwell equation J = c 4π ∇ × H, where the magnetic field H and the magnetic induction B are related by the magnetic permeability coefficient µ: B = µH.The resulting system of MHD equations for compressible conductive fluids or gases in the low-magnetic-Reynolds-number approximation is where ρ, u, P , and e are the density, velocity, pressure, and specific internal energy of the fluid.Equation ( 5) is the consequence of Ohm's law, and the Poisson equation ( 6) follows from the local neutrality of fluids ∇ • J = 0 and the Ohm's law (5).Equation ( 4) is the equation of state (EOS) that describes material properties.Three material substances are present in our simulations: liquid mercury, ambient gas, and mercury vapor inside cavitation bubbles.The polytropic ideal gas equation of state is used for the ambient gas and mercury vapor, and the stiffened polytropic EOS [23] with the adiabatic exponent γ l = 3.2 and the stiffening constant P ∞ = 8 • 10 10 g/(cm s 2 ) is used for liquid mercury in this study.The stiffened polytropic EOS extends the range of pressure to negative values in order to account for the transient effect of tension in liquids.With the stiffened polytropic EOS, a thermodynamically consistent state can be achieved for P > −P ∞ as the sound speed becomes imaginary below this point: c 2 = γ(P − P ∞ )/ρ.The cavitation algorithm described in the next section prevents the pressure in liquids from falling below some specified critical pressure at which the liquid breaks in the form of cavitation bubbles that expand and relieve the tension.
The following boundary conditions must be satisfied at the interface Γ of a conducting fluid with a dielectric medium: i) the normal component of the velocity field is continuous across the interface; ii) the pressure jump at the interface is defined by the surface tension S and main radii of curvature: iii) the normal component of the current density vanishes at the interface giving rise to the Neumann boundary condition for the electric potential where n is a normal vector at the fluid free surface Γ.

Front-tracking algorithm for MHD equations
The existence of sharp, geometrically complex material interfaces presents a great challenge to numerical algorithms.We solve the coupled hyperbolic-elliptic system (1-8) in an operator-splitting fashion using the method of front tracking.The front-tracking technique, its implementation in the multiphysics code FronTier, and applications to computational fluid dynamics problems are described in detail in [12][13][14] and references therein.The front-tracking algorithm for multiphase MHD equations in the low magnetic Reynolds number approximation was developed in [10].In this paper, we summarize the main ideas of the algorithm and refer the reader to the above references for details.
Front tracking is an adaptive computational method in which a lower dimensional moving grid or interface (a triangulated surface in 3D space) represents material interfaces.The interface contains left and right physical states corresponding to materials on both sides and keeps the discontinuity sharp.The key feature of the method is greatly reduced numerical diffusion due to the absence of numerical finite differentiation across the interface.Front tracking is implemented in a multi-physics code FronTier [13].FronTier is capable of tracking and resolving the topological changes of a large number of interfaces in two-and three-dimensional spaces.
Using the operator splitting, we solve the hyperbolic equations first, followed by solving the elliptic boundary value problem (6) and (8).The time step starts with the advance of the interface.The propagation of the interface is accomplished by solving the MHD equations in the normal and tangential directions to each interface point.In the low-magnetic-Reynolds-number approximation, the MHD equations reduce to Euler equations with external forces that come from the elliptic step.The rectangular grid, interface, and states for the normal propagation of the interface are shown schematically in figure 1.After the propagation of the interface points, the new interface is checked for consistency of intersections.The untangling of the interface at this stage consists in removing the unphysical intersections, and rebuilding a topologically correct interface.The final phase of the hyperbolic time step update consists of computing new states on the rectangular spatial grid using one of the second order shock capturing solvers.
After the hyperbolic step, the elliptic problem for the electric potential (6, 8) is solved using the embedded-boundary method [24] (see [25] for the latest development of the embedded-boundary method for 3D moving elliptic-interface problems with front tracking).High-performance, parallelsoftware libraries of preconditioners and iterative solvers based on Krylov subspace methods such as PETSC [26] are used for solving the corresponding linear system of equations, and the electromagnetic terms are calculated.In the final stage of the time step, the states in the interior and on interfaces are modified using electromagnetic terms.

Models for phase transition and cavitation
We are interested in the simulation of liquids in the presence of strong rarefaction waves.If the liquid is not capable of sustaining tension, it breaks or cavitates [27].The capability of liquids to withstand tension depends on their purity.The homogeneous nucleation theory [27], applicable to liquids without impurities that serve as initial cavitation centers, gives the following probability of nucleation in a volume V during time interval τ : where is the critical energy necessary to create the surface of a cavitation bubble of the critical radius R C , and is related to the critical pressure (the critical strength of the tensile pressure in the liquid) as and S is the surface tension.The coefficient J 0 is where N is the number density of the liquid (molecules/volume) and m is the molecular mass.In real fluids, the existence of impurities significantly reduces the values of the critical pressure and energy compared to the values predicted by the homogeneous nucleation theory.Our DNS model of cavitation is implemented as follows: when the pressure falls below the critical pressure, cavitation bubbles are formed in some spatial distribution in the rarefaction wave by inserting numerically tracked bubble surfaces and replacing liquid states with vapor states inside.The values of the critical pressure and the initial number density of cavitation bubbles are selected individually for each numerical experiment by estimating liquid conditions from available experimental data.In most simulations, we assume that the amount of vapor in a bubble is constant and the expansion of cavitation bubbles is caused only by pressure gradients.In other words, we neglect the phase transition on the bubble surface.If the phase-transition-induced mass flux is essential for the interface dynamics, the phase-transition algorithm developed in [28] can be used.Fast vaporization or ablation is especially important in the pellet-fueling processes mentioned in the introduction.Initial cavitation bubble sizes in real liquids (for instance, R c = 1 micron for mercury at ∆P c = 10 bar) is close to numerically resolved limits.However in practice, we frequently use larger initial bubble sizes for coarser grid computations in large domains.This leads to some numerical errors including the loss of mass and momentum conservation.While we estimate that such numerical errors are insignificant in most of simulations, we are currently working on strict enforcement of the conservative properties using a multi-scale coupling with a stand-alone, resolved, one-dimensional simulation of a cavitation bubble in the ambient liquid.Similar cavitation algorithms have already been used in 2D [20] and 3D simulations [29].

Applications to accelerator targets
In this section, we apply the front-tracking algorithms for multiphase MHD to the simulation of the mercury target for the proposed Neutrino Factory / Muon Collider [18].The target will contain a series of 30-cm-long and 1-cm-diameter mercury jets entering a strong (∼ 15 Tesla) magnetic field at a small angle to the solenoid axis.When each jet reaches the center of the solenoid, it interacts with a powerful proton pulse penetrating the jet and depositing energy of the order of 100 J/g into mercury.The purpose of our numerical simulations is to evaluate the states of the target before and after the interaction with protons to optimize the target design.Preliminary simulations have been reported in [19,20,22].

Simulation of the mercury jet entering a solenoid magnet
The nonuniform transverse component of the magnetic field with respect to the jet trajectory, caused by a small angle between the jet and the magnetic solenoid axis, distorts the jet during the motion toward the solenoid center.To evaluate the state of the jet target before the interaction with the solenoid axis, we performed numerical simulations of the jet entering the solenoid using the real profile of the magnetic field along the jet trajectory in the MERIT experiment.Results are summarized in figure 3, which plots the transverse-distortion ratio of the jet or the maximum radius in the transverse cross section normalized by the unperturbed initial radius of the jet.We observe that the distortion strongly depends on the angle between the jet and the solenoid axis: the maximum distortion ranges from 1.2 at the angle of 0.05 rad to 2.75 at the angle of 0.15 rad.The latter distortion is unacceptable for the target: it transforms the jet into a thin sheet and significantly reduces the effective cross section with the proton pulse and, as a result, the pionproduction rate.To reduce the amount of distortion, the angle of 0.033 rad was used in the MERIT experiment.Simulations of mercury jets in transverse magnetic fields similar to those presented in figure 3 agree very well with theoretical calculations [10], other experiments [5], and simulations using the HYMAG code [30].However, the simulations underestimated the amount of distortion for the MERIT experiment.Numerical simulation of the mercury jet entering the 15-T solenoid at the angle of 0.033 rad predicted the distortion of 1.17 and the value of 1.8 was observed [21].The reason for this disagreement is related to deviations from ideal simulation conditions which are currently under investigation.

Simulation of the mercury jet interacting with proton pulses.
When the mercury jet reaches the solenoid center, it interacts with a proton pulse, depositing energy of the order of 100 J/g into mercury.Because of the short time scale of the interaction,  we assume that the increase of the internal energy and pressure is an isochoric process.In computations, we modify the internal energy and pressure states of the jet during the initial time step according to predictions of atomistic Monte-Carlo simulations of the mercury-proton pulse interaction using the MARS code [31].In most of simulations, we assumed that the initial shape of the jet is distorted by the transverse component of the magnetic field, and the cross section of the jet is an ellipse with the long and short radii of 0.8 and 0.4 cm, correspondingly.But in some simulations we also used the unperturbed cylindrical jet.The profiles of pressure after the interaction with 24-GeV, 10-teraproton and 14-GeV, 10-teraproton pulses are shown in figure 4. For the mercury jet with the elliptic cross section, we monitored the jet surface velocity in four radial directions, or points A, B, C and D shown in figure 4(b).After the energy deposition, the high-pressure wave propagates outward and reflects from the mercury-air interface as a strong rarefaction wave.Since the mercury is incapable of sustaining such a large tension, it cavitates and the growth of cavitation bubbles causes a reduction of tension, rapid jet expansion and surface instabilities.A snapshot of cavitation bubbles in the jet is shown in figure 5.If the jet cavitation and expansion occurs in a longitudinal magnetic field, the radial motion of the fluid induces vortices of azimuthal current.Then, the Lorentz force opposes the fluid motion.Hence the magnetic field reduces the amount of cavitation and tends to stabilize the mercury jet.Snapshots of the jet surfaces at 100 µs after the interaction with the proton pulse in magnetic fields ranging from 0 to 15 Tesla are shown in figure 6.
Figures 7-10 show the evolution of velocities of elliptic-jet surface filaments in four radial directions A, B, C, and D, as explained above.The maximum velocity of filaments ejected in the direction of the short axis reaches 35 m/s.In all directions, the velocity decreases with the increase of the magnetic field strength.We would like to emphasize that the formation and evolution of filaments cannot be attributed solely to classical fluid interface instabilities such as the Rayleigh-Taylor and Richtmyer-Meshkov instabilities.We have shown that the formation of mercury jet filaments critically depends on the presence of cavitation.In figure 11, we compare the velocity of jet-surface filaments with experimental data, which is sparse; it is available for different values of the beam intensity but not always for the 10-teraproton beam which was used for numerical simulations.The simulated filament velocity at 5-T magnetic field agrees well with the experiment but the experimental value must be assigned a large error bar because of inconsistency of measurements for 10-, 15-and 20-teraproton beams.The simulated velocity at 10-T field agrees well with the interpolated value using the experimental data points at 5 and 15 teraproton.There is only one measurement available for the 15-T magnetic field which shows no velocity reduction compared to the 10-T field.The corresponding simulation indicates a further reduction of the filament velocity compared to 5-and 10-T fields.Because of the sparseness of experimental data and limited number of spikes used in the data analysis, it is difficult to make a conclusion on the last case.

Conclusions and future work
We have described a numerical algorithm for the simulation of free-surface magnetohydrodynamic flows at low magnetic Reynolds numbers.The corresponding governing equations constitute a coupled hyperbolic-elliptic system in a geometrically complex and evolving domain.The numerical algorithm includes the interface-tracking technique for the hyperbolic problem, a Riemann problem for the material interface, discretization of elliptic equations in irregular domains with interface constraints using the embedded-boundary method, and high-performance, parallel solvers such as MUSCL-type schemes for hyperbolic problems and iterative solvers implemented in the PETSc package.An extensive theoretical analysis of the method of front tracking for hyperbolic systems of conservation laws has already been performed in earlier works, and the method has been validated and tested on problems of Rayleigh-Taylor and Richtmyer-Meshkov surface instabilities.The embedded-boundary technique for irregular geometric domains within the method of front tracking was developed in [10] and generalized to include elliptic problems with discontinuities along interfaces in [25], and shown to be second order accurate for both the electric potential and its gradient.
The front-tracking MHD algorithm has been used for the simulation of tokamak fueling by the ablation of cryogenic pellets [15,16] and the dynamics of laser ablation plumes [17].In this paper, we applied the numerical MHD algorithms for the simulation of liquid mercury jet targets for the Neutrino Factory/Muon Collider.Simulations of the mercury jet entrance into a non-uniform magnetic field at small angles to the solenoid axis showed that the transverse component of the magnetic field distorts the jet.The jet cross section becomes an elongated ellipse with the long axis parallel to the direction of the transverse field.Simulations of the mercury jet interacting with proton pulses in magnetic fields have also been performed.After the proton-pulse energy deposition, the pressure in the center ranges from 10 to 30 kbar, depending on the proton-pulse properties.Then the high-pressure wave propagates towards the jet boundary and causes tension that leads to severe cavitation in mercury and the formation and growth of surface filaments.The fastest filaments reach velocities of 35 m/s.The cavitation and surface filamentation is reduced by longitudinal magnetic fields.Simulations are in reasonably good agreement with the results of the mercury target experiment MERIT [21] conducted at CERN in the fall of 2007.
In our future work, we will develop a front-tracking MHD algorithm for multiphase or freesurface incompressible fluids in the low-magnetic-Reynolds-number approximation.The algorithm will enable the FronTier code to perform long-time-scale simulations of liquid metals in magnetic fields for applications which do not involve strong waves.Such are liquid-metal blankets on the interface with plasma that can be used for the protection of tokamak walls, mercury-handling components in accelerator targets and other applications.We are also working on the free-surface algorithm for the opposite MHD regime: the ideal MHD approximation for highly conducting plasmas.Such an algorithm will allow us to simulate the behavior of magnetized plasma targets compressed by liners with the purpose of achieving nuclear fusion ignition.Hybrid magneto-inertial fusion methods [32] have recently gained large attention due to potential capability of overcoming difficulties of traditional approaches.

Figure 1 .
Figure 1.Rectangular grid, interface, and states for the method of front tracking.States contain density, momentum, and energy density of the fluid, and references to the EOS model and other parameters.

Figure 2 .
Figure 2. Schematic of the mercury jet target for Neutrino Factory / Muon Collider.

Figure 3 .
Figure 3. Normalized elliptic deformation of the mercury jet entering a 15-T solenoid at the angle of 0.15 rad (top line), 0.1 rad (middle line), and 0.05 rad (bottom line) to the magnetic axis.

Figure 5 .
Figure 5. Cavitation bubbles in the mercury jet at 20, 130, 200, and 250 µs after the interaction of the cylindrical jet with the proton pulse.

Figure 6 .
Figure 6.Filaments on the surface of initially cylindrical jet at 150 µs after the interaction with the proton pulse in magnetic fields ranging from 0 to 15 Tesla.

Figure 7 .
Figure 7. Expansion velocity of jet surface filaments at zero magnetic field.

Figure 8 .
Figure 8. Expansion velocity of jet surface filaments in 5 Tesla magnetic field.

Figure 9 .
Figure 9. Expansion velocity of jet surface filaments in 10 Tesla magnetic field.

Figure 10 .
Figure 10.Expansion velocity of jet surface filaments in 15 Tesla magnetic field.

Figure 11 .
Figure 11.Comparison of experimental and simulated values of the jet surface filament velocity.