Simulation of 3D fluid jets with application to the Muon Collider target
design

**J. Glimm**${}^{1}$**, H. Kirk**${}^{2}$**, X. L. Li**${}^{1}$**, J.
Pinezich**${}^{1}$**, R. Samulyak**${}^{2}$ **& N. Simos**${}^{2}$

${}^{1}$Department of Applied Mathematics and Statistics, State
University of New York at Stony Brook, USA.

${}^{2}$*Brookhaven National Laboratory, U.S. Department of Energy,
Upton, USA.*

# Abstract 

#  

The FronTier code was developed to model the propagation of three
dimensional fluid jets. FronTier is based on the method of front
tracking, a numerical technique for solving systems of conservation laws
in which the evolution of discontinuities is determined through the
solution of the associated Riemann problem. The code was used for the
numerical simulation of the Muon Collider target. The target is in the
form of a mercury jet interacting with high energy proton beams in the
presence of a strong magnetic field. The numerical simulation of the
target evolution driven by strong pressure waves is important for the
optimal target design. Theoretical estimations of the target behavior in
the magnetic field, based on the linear stability analysis, are also
presented. In addition, the response of the mercury target interacting
with the pulsed proton beam was studied using the ANSYS general purpose
code for the thermal and pressure wave propagation.

## 1 Introduction

In order to understand the most basic structure of matter and energy, an
advance in the energy frontier of particle accelerators is required.
Advances in high energy particle physics have been, and are expected to
continue to be, dependent on and paced by advances in accelerator
facilities. The majority of contemporary high-energy physics experiments
utilize colliders. A study group was organized at Brookhaven National
Laboratory to explore the feasibility of a high energy, high luminosity
muon-muon collider. Such a collider has numerous potential advantages
over traditional electron-positron machines (linear colliders), however,
several challenging technological problems remain to be solved. One of
the most important is to create an effective target able to generate the
high-flux muon beam. The need of operating high atomic number material
targets in particle accelerators that will be able to withstand the
intense thermal shock has led to the exploration of free liquid jets as
potential target candidates for the proposed Muon Collider \[12\]. The
target will be designed as a pulsed jet of mercury (high Z-liquid)
interacting in a strong magnetic field with a high energy proton beam. A
global numerical study and an analytical estimation of the behavior of
the mercury jet under such energy deposition is the main goal of the
present research.

The numerical simulation was performed using the computer code FronTier,
developed to model the propagation of three dimensional fluid jets.
FronTier is based on front tracking, a numerical method for solving
systems of conservation laws in which the evolution of discontinuities
is determined through solutions of the associated Riemann problems. For
the thermal deposition and pressure wave propagation analysis ANSYS, a
general purpose code, was also used.

The paper is organized as follows. In Section 2 we describe some details
of the Muon Collider target design and present analytical estimates of
the behavior of the mercury jet in the magnetic field based on linear
stability analysis of the magnetohydrodynamical system of equations.
Section 3 is devoted to the description of the main ideas of the method
of front tracking and their realization in the FronTier code. We present
results of the numerical simulation of the target behavior obtained
using the FronTier code in Section 4. In this section, we also present
some numerical results of the thermal analysis of the target and the
pressure wave calculations obtained by using the ANSYS code. Finally, we
conclude in Section 5 with a discussion of the consequences of the work
in this paper and possible directions for future research.

## 

## 2 Muon Collider target

##  

**2.1 Design parameters**

The Muon Collider target is shown schematically in Figure 1. It will
contain a series of mercury jet pulses of about 1 cm in radius and 30 cm
in length. Each pulse will be shot at a velocity of 10-15 m/s into a 20
Tesla magnetic field at a small angle $(6^{0})$ to the axis of the
magnetic field. When the jet reaches the center of the magnet it is hit
with two 2 ns proton pulses, 150 ns apart; each proton pulse will
deposit about 40 kJ of energy in the mercury. The value and spatial
distribution of the energy deposition was estimated based on the
interaction of a beam of $5 \cdot \text{10}^{\text{13}}$ protons with
mercury.

![](media/image1.jpeg){width="4.326388888888889in"
height="2.4277777777777776in"}

Figure 1. Muon Collider target

The main goal of the present paper is to estimate the state of the
target after the interaction with the proton pulses. In particular, the
main issues are the distortion of the jet due to eddy currents as it
propagates into the magnetic field, the deformation of the jet surface
due to strong pressure waves caused by the proton pulses, the
possibility of the jet to break up into droplets and the pressure field
evolution.

2.  **Linear stability analysis**

Analytical estimation of the mercury jet behavior during its propagation
in the magnetic field provide important data about the state of the jet
before the interaction with the proton beam.

Classical linear stability analysis \[2,13,14,15\] of the inviscid
hydrodynamical system of equations leads to the following conclusions:

i\) The jet is stable for all non-axisymmetric deformations;

ii\) The jet is stable for axisymmetric deformations with wavelengths
exceeding the circumference of the jet;

iii\) The wave number of the most unstable wave is *k=*0.697/*R*, where
*R* is the radius of the jet, for which the coefficient of the
exponential wave growth is

$\sigma_{\text{max}} = 0\text{.}\text{3433}\sqrt{\frac{T}{R^{3}}\rho}$,
where *T* is the surface tension and $\rho$is the density of the liquid.
The characteristic time for the mercury jet breakup is about 0.5 sec
which is much greater than the propagation time of the target in real
experiments.

For the study of the magnetohydrodynamic stability of the jet we shall
follow the analysis presented in \[2\]. The equations of
megnetohydrodynamics applicable to the jet stability analysis are

$\begin{array}{r}
\rho\frac{du}{\text{dt}} = \text{div}\ P + \rho X + \mu J \times H, \\
\frac{\partial H}{\partial t} - \text{curl}(u \times H) = - \text{curl}\left( \eta\ \text{curl}H \right), \\
\text{curl}H = 4\pi J,
\end{array}$ (1)

where $\eta = \frac{1}{\sqrt{4\text{πμσ}}}$ is the resistivity, $\mu$is
magnetic permeability, $\sigma$is the conductivity, $P$is the total
stress-tensor and $X$includes the external forces of non-magnetic
origin. Then the linearized system of magnetohydroynamic equations for
an inviscid fluid of finite electrical conductivity is

$\begin{array}{r}
\frac{\partial u}{\partial t} = - \text{grad}\Pi + \frac{\mu H}{4\text{ηρ}}\frac{\partial h}{\partial z}, \\
\frac{\partial h}{\partial t} = \text{curl}(u \times H) - \eta\text{curl}^{2}h, \\
\text{div}h = \text{0,}\ \text{div}h = \text{0,}\ \Pi = \frac{\delta p}{\rho} + \frac{\mu\text{Hh}_{z}}{4\text{ηρ}},
\end{array}$ (2)

where $\delta p$and $h$are the pressure and magnetic field
perturbations, respectively. Conditions at the jet boundary are:

i\) The radial component of the velocity field on the boundary is
compatible with the form of the boundary, which we assume to be
harmonically perturbed:

$$w = R + \varepsilon_{0}\text{exp}\left( i(\text{kz} + m\varphi) + \sigma t \right);$$

ii\) The normal component of $h$is continuous across the boundary;

iii\) The normal component of the total stress tensor is continuous
across the boundary.

From these assumptions the following expression defining the growth rate
of the perturbation waves on the jet surface can be obtained

$\sigma^{2}\left\lbrack 1 + \frac{\Omega_{A}^{2}}{\sigma^{2}\Phi}\frac{{yI}_{0}(y)}{{xI}_{0}(x)} \right\rbrack = \frac{T}{R^{3}\rho}\left\lbrack 1 + \frac{\Omega_{A}^{2}}{\sigma^{2}\Phi}\frac{I_{1}(y)}{I_{1}(x)} \right\rbrack\frac{{xI}_{1}(x)}{I_{0}(x)}\left( 1 - x^{2} \right),$
(3)

where

$$\begin{array}{r}
\Phi = {yI}_{0}(y)K_{1}(x) + {xI}_{1}(y)K_{0}(x),\ x = \text{kR}, \\
y^{2} = x^{2} + \frac{R^{2}}{\text{ησ}}\left( \sigma^{2} + \Omega_{A}^{2} \right),\ \Omega_{A}^{2} = \frac{{\mu H}^{2}}{4\text{πρ}R^{2}}x^{2},
\end{array}$$

and $I_{n}(x),\ K_{n}(x),\ n = 0,1,$are the Bessel functions of order
*n* of the first and second kind, respectively.

Asymptotic analysis of equation (3) for different limiting values of the
resistivity makes it possible to conclude:

i\) In the limit of infinite resistivity the magnetic field has no
effect on the capillary instability of liquid jets;

ii\) In the limit of zero resistivity the magnetic field tends to
stabilize the jet, in particular, a field strength of 100 Gauss will
stabilize the jet for the varicose deformations of all wavelengths;

iii\) In the limit of high but finite resistivity (the case of mercury)
the stabilizing effect of the magnetic field is diminished. For a
mercury jet with radius *R*=1cm a magnetic field strength of order
$\text{10}^{4}\text{Gauss } = \text{ 1 Tesla}$ will be necessary to
stabilize the jet against the varicose deformations of all wavelengths.

## 3 Method of front tracking 

The FronTier code was developed to model the propagation of three
dimensional fluid jets. The FronTier code is based on front tracking
\[3,5,6\], a numerical method which allows the computation of
dynamically evolving interfaces. In general, interfaces can model many
types of discontinuities in a medium such as shock waves in gas
dynamics, boundaries between fluid-gas states, different fluids or their
different phases in fluid dynamics, component boundaries in solid
dynamics etc. The unifying feature of the models solvable by the method
of front tracking is the possibility to represent the governing system
of equations in the form of conservation laws and to determine the
evolution of the discontinuity through the solution of the Riemann
problems \[3,11\]. This method often doesn\'t require highly refined
grids and it has no numerical diffusion.

Front tracking represents interfaces as lower dimensional meshes moving
through a volume filling grid. The traditional volume filling finite
difference grid supports smooth solutions located in the region between
interfaces and the lower dimensional grid or interface defines the
location of the discontinuity and the jump in the solution variables
across it. The dynamics of the interface comes from the mathematical
theory of Riemann solutions, which is an idealized solution of a single
jump discontinuity for a conservation law. Where surfaces intersect in
lower dimensional objects (curves in three dimensions), the dynamics is
defined by a theory of higher dimensional Riemann problems such as the
theory of shock polars in gas dynamics. Nonlocal correlations to these
idealized Riemann solutions provide the coupling between the values on
these two grid systems.

The computation of a dynamically evolving interface requires the ability
to detect and resolve changes in the topology of the moving front. A
valid interface is one where each surface and curve is connected,
surfaces only intersect along curves and curves only intersect at
points. We say that such an interface is untangled. Two independent
numerical algorithms, grid-based tracking and grid-free tracking, were
developed \[5-7\] to resolve the untangling problem for the moving
interface. The advantages and deficiencies of the two methods are
complementary and an improved algorithm combining them into a single
hybrid method was implemented in the FronTier code and described in
\[7\]. Our numerical results presented in the next section were obtained
by using this hybrid method of front tracking.

2.  **Results of numerical simulation**

**3.1 Evolution of the target**

In this Section, we present the results of the Muon Collider target
simulation by FronTier. The initial jet, prior to interaction with the
first proton pulse, was obtained by the FronTier code in a natural way
using physical assumptions derived from the design of the jet gun. In
particular, it was assumed that mercury comes from the nozzle with a
prescribed velocity of 10 m/s, the outside pressure is 1 atm, the
pressure inside the jet differs by the value $\frac{T}{R}$, where *T* is
the surface tension and *R* is the radius of the jet and the temperature
of the ambient gas is equal to the temperature inside the jet. The state
of the ambient gas was modelled by the polytropic gas equation of state
with $\gamma = 1\text{.}4$while the liquid state was modeled by the
stiffened polytropic equation of state with the stiffening coefficient
$\gamma_{\text{st}} = 7\text{.}5$\[1,11\]. Since the viscous length
scale for the problem is negligibly small and the viscous force can play
a significant role only near the pinch-off point of the jet, the
viscosity effect was neglected. Thus the characteristic length scale for
the problem is the jet diameter. The propagation of the jet boundary was
determined through the Riemann solution to the compressible inviscid
fluid system of equations.

The state of the jet before the interaction with the first proton pulse
is shown in Figure 1a. We observe that the jet surface does not exhibit
any significant surface instabilities, in agreement with the linear
stability analysis predictions. The interaction of the jet with the
proton pulses was modeled by adding external energy to the interior
energy at points located near the axis of the jet at one specific moment
of time. This is a reasonable approximation since the energy deposition
time interval will be only 2 ns in real experiments. The energy
deposition was transformed within the code to a pressure perturbation by
the equation of state. Figures 2b and 2c depict the evolution of the jet
surface instabilities after the interaction with the first and second
proton pulses, respectively. Notice that the time scale of the numerical
experiments was increased to microseconds since the jet surface didn't
exhibit any essential changes within 150 ns. In order to study the
influence of the single proton pulse the time interval between pulses
was extended to $t = 4 \cdot \text{10}^{- 6}$sec.

![](media/image2.jpeg){width="1.99375in" height="2.8430555555555554in"}
![](media/image3.jpeg){width="2.411111111111111in"
height="2.8368055555555554in"}

a\) State of the jet before interaction b) State of the jet after
interaction with

with the first proton pulse, *t=*0 the first proton pulse,
$t = 4 \cdot \text{10}^{- 6}$sec

![](media/image4.jpeg){width="2.727777777777778in"
height="3.2715277777777776in"}

c\) State of the jet after interaction with

the second proton pulse, $t = 4 \cdot \text{10}^{- 6}$sec

Figure 2. Evolution of the mercury jet

# 3.2 Heat Generation and Diffusion

In order to address the question of thermal relaxation in the mercury
target between pulses, a transient heat diffusion analysis was performed
using the ANSYS code. Specifically, two beam pulses separated in time by
150 ns were implemented assuming that thermalization due to the energy
deposition occurs instantaneously. The ionization energy, the latent
heat and the radiation heat loss were neglected. Figure 3 depicts the
temperature distribution at the end of the first pulse. The temperatures
experienced by most of the jet exceed those of mercury vaporization.
There exists an outer shell layer, however, with temperatures still
below the vaporization temperature at the end of the first pulse. The
heat diffusion analysis that followed the end of the first pulse showed
that the thermal relaxation during the 150ns is insignificant. The
deposition of heat by the second beam pulse will result in temperatures
in the core of the target in excess of 2000 K. For a more realistic
thermal analysis the heat diffusion, ionization energy and latent heat
need to be implemented in the model along with better representation of
mercury properties as functions of temperature.

![](media/image5.jpeg){width="4.404861111111111in"
height="0.2791666666666667in"}

Figure 3. Temperature distribution in the target at the end of the

first proton pulse (only one half of the jet is shown)

# 3.3 Generation and propagation of pressure waves

The pressure distribution calculated according to the energy deposition
profile was implemented in the model as a non-equilibrium initial
condition and the solution to the corresponding acoustic wave equation
was obtained numerically by using the ANSYS code. The goal of this
simulation was to observe the propagation of the pressure waves due to
pressure gradients present in the initial conditions and to establish a
time frame in which pockets of negative pressure start to appear in the
wake of the wave front as well as a result of wave reflection at the jet
surfaces.

The first frame of Figure 4 depicts the initial pressure at a cross
section of the jet. The first simulation performed within the 150ns
interval representing the distance between two proton beam pulses
demonstrated that this time scale is negligibly small for any essential
advance of the pressure front. In order to observe the dynamics of the
process, the time scale of the simulation was increased to several
microseconds. The various frames of Figure 4 depict the pressure profile
within the jet as the front moves both axially and radially and reflects
at the boundaries. Some domains of negative pressure were observed at
later times (Figure 4d).

![](media/image6.jpeg){width="0.375in" height="3.61875in"}
![](media/image7.jpeg){width="0.375in" height="3.6222222222222222in"}
![](media/image8.jpeg){width="0.375in" height="3.613888888888889in"}
![](media/image9.jpeg){width="0.375in" height="3.613888888888889in"}

a\) $t = 0$ b) $t = 2\text{.}7 \cdot \text{10}^{- 6}$sec c)
$t = 3\text{.}7 \cdot \text{10}^{- 6}$sec d)
$t = 1 \cdot \text{10}^{- 5}$sec

Figure 4. Pressure waves caused by downward-projected proton beam.

3.  **Conclusions**

The FronTier code was developed to model the propagation of three
dimensional fluid jets and applied to the numerical simulation of the
Muon Collider target. The state of the mercury jet before its
interaction with the first proton pulse obtained by the FronTier code is
in a good agreement with the prediction of the classical linear
stability analysis. The linear stability analysis of the
magnetohydrodynamical set of equations demonstrate the stabilizing
effect of the magnetic field of the mercury jet.

The results of the numerical simulation demonstrate the evolution of the
surface instabilities of the mercury jet due to the proton pulse energy
deposition. The analysis however shows that the time interval between
the two proton pulses is not large enough for the development of any
significant surface instabilities. Therefore, as required for normal
muon production, the mercury jet will not immediately break into
droplets following the energy deposition process. Similar predictions
were obtained in \[8\]. However the observed surface instabilities will
probably lead to the jet break-up at later time when the jet is leaving
the magnetic field.

Numerical simulation of the thermal and pressure wave propagation
obtained using the ANSYS code confirmed the conclusion of the previous
analysis about the time scale of the problem. These results, however,
show that the temperature in the center of the jet reaches values above
the boiling point of mercury. Therefore, an improved analysis of the
target behavior is required, for which the FronTier code must be further
developed to accommodate a liquid-vapor mercury phase transition and a
modified equation of state for mercury near the critical point. This
will be the main effort of our future research.

**References**

1\. R. Courant and K. Friedrichs, *Supersonic Flows and Shock Waves*,

Interscience, New York, 1948.

2\. S. Chandrasekhar, *Hydrodynamic and Hydrodynamic stability*,
Clarendon

press, Oxford, 1961.

3\. I. R. Chern, J. Glimm, O. McBryan, B. Plohr and S. Yaniv, Front
tracking

for gas dynamics, *J. Comp. Phys*., **62** (1986), pp. 83-110.

4\. J. Eggers and T. Dupont, Drop formation in a one dimensional
approxima-

tion of the Navier-Stokes equation, *J. Fluid Mech*., **262** (1994),
pp. 205-221.

5\. J. Glimm, J. W. Grove, X. L. Li, K. M. Shyue, Q. Zhang and Y. Zeng,
Three

dimensional front tracking, *SIAM J. Sci. Comp*., **19** (1998), pp.
703-727.

6\. J. Glimm, M. J. Graham, J. W. Grove, X. L. Li, T. M. Smith, D. Tan,

F. Tangerman and Q. Zhang, Front tracking in two and three dimensions,

*Computers Math. Applic*., **35** (1998), pp. 1-11.

7\. J. Glimm, J. W. Grove, X. L. Li, and D. Tan, Robust computational

algorithms for dynamic interface tracking in three dimensions, Los
Alamos

National Laboratory Report LA-UR-99-1780.

8\. A. Hassanein, Simulation of intense heating and shock hydrodynamics

in free moving liquid targets, Nuclear Applications of Accelerator

Technology, ACCApp'99, Long Beach, CA, November 14-18, 1999.

9\. R. L. Holmes, J. W. Grove and D. H. Sharp, Numerical investigation
of

Richtmyer-Meshkov instability using front tracking, *J. Fluid Mech*.,

**301** (1995), pp. 51-64.

10\. J. Keller and M. Miksis, Surface tension driven flows, SIAM J.
Appl. Math.,

**43** (1983), pp. 268-277.

11\. R. Menikoff and B. Plohr, The Riemann problem for fluid flow of
real

materials, *Rev. Mod. Phys*., **61** (1989), pp. 75-130.

12\. R. B. Palmer, Muon Collider design, Brookhaven National Laboratory

report, BNL-65242, CAP-200-MUON-98C, 1998.

13\. Lord Rayleigh, On the instability of jets, *Proc. London Math.
Soc*., **10**

(1879), pp. 4 - 13.

14\. R. D. Reitz and F. V. Bracco, Mechanism of atomization of a liquid
jet,

*Phys. Fluids*, **25** (1982), pp. 1730 - 1742.

15\. G. I. Taylor, Generation of ripples by wind blowing over a viscous
fluid,

*Collected works of G.Taylor*, edited by G.Batchelor, Cambridge U.P.,
1958.