%      
%\documentclass[11pt,dvips]{article}
%\usepackage{epsfig} 

%\begin{document}

%\newcommand\grad{\nabla}
%\newcommand\f{\frac}
%\newcommand\divg{\nabla\cdot}   %\newcommand\curl{\mbox{curl}}
%\newcommand\curl{\nabla\times}   %\newcommand\divg{\mbox{div}}
%\newcommand\p{\partial}
%\newcommand\B{\mathbf{B}}
%\newcommand\E{\mathbf{E}}
%\newcommand\J{\mathbf{J}}
%\newcommand\vel{\mathbf{u}}
%\newcommand\normal{\mathbf{n}}
%\renewcommand\H{\mathbf{H}}
%\newcommand\be{\begin{equation}}
%\newcommand\en{\end{equation}}


\section{Numerical Simulation of the Muon Collider Target}

The numerical simulation of hydro- and magnetohydrodynamics processes in the Muon Collider 
target is needed for better understanding of the physics involved and for answering 
 engineering design questions.
The Muon Collider target  will contain a series 
of mercury jet pulses of about 0.5 cm in radius and 60 cm in length. Each pulse will be
shot at a velocity of 30-35 m/sec into a 20 Tesla magnetic field at a small angle (0.1 rad)
to the axis of the field. When the jet reaches the center of the magnet, it is hit with 
3 ns proton pulses arriving with 20 ms time period; each proton pulse will deposit about
100 J/g of energy in the mercury. 

The main issues of the target design addressed in our numerical studies
are the distortion of the jet due to eddy currents
as  it propagates through the magnetic coil, the deformation of the jet surface due to 
strong pressure waves caused by the proton pulses and the probability of the jet breakup. 
Studying the state of the target during its interaction
with proton pulses will help to achieve  the maximal proton production rate
and therefore an optimal target performance. 


\subsection{Magnetohydrodynamics of Free Surface Liquid Flows}


The basic set of equations describing the interaction of a compressible conducting fluid 
flow and a magnetic field is contained in Maxwell's equations and in the equations of 
fluid dynamics suitably modified \cite{Chandra,Landau}. Namely, the systems contains the mass, 
momentum and energy conservation equations for the fluid which have hyperbolic 
nature and a parabolic equation for the evolution of the magnetic field. 
\begin{eqnarray}
\frac {\p \rho}{\p t} & = & - \nabla\cdot(\rho\vel),
\label{eq1}\\ 
\rho \left(\frac \p{\p t} +\vel \cdot\nabla\right)\vel
& = & - \nabla P + \rho \mathbf{X} + \frac1c(\J \times \B),
\label{eq2}\\
\rho \left(\frac \p{\p t} +\vel \cdot\nabla\right)U
& = & -P\nabla\cdot\vel + \frac1\sigma \J^2 - \frac1c\vel\cdot(\J \times \B),
\label{eq3}\\ 
\frac{\p \B}{\p t}
& = & \nabla\times(\vel \times \B) - 
\curl (\frac{c^2}{4\pi\sigma} \curl \B),
\label{eq4} \\ 
\divg \B & = & 0,
\label{eq5} 
\end{eqnarray}
Here $\vel,\, \rho$ and $U$ are the velocity, density and total energy of the fluid, 
respectively, $P$ is the total stress tensor, $\mathbf{X}$ includes external forces of 
non-magnetic origin, $\B$ is the magnetic field induction, $\J$ is the
current density distribution and $\sigma$ is the fluid conductivity. 
The magnetic field $\H$ and magnetic induction $\B$ are related by the 
magnetic permeability coefficient $\mu$: $\B=\mu\H$.

The system (\ref{eq1}-\ref{eq4}) must be closed with an equation of state. 
We are especially interested in fluid behavior under extreme thermodynamics conditions.
Some related equation of state models are discussed below.

The following boundary conditions must be satisfied at the jet surface:

i) the normal component of the velocity field is continuous across the 
material interface.

ii) the normal and tangential components of the magnetic field at the 
material interface are related as
\begin{equation} 
\normal\cdot(\B_2 - \B_1) = 0, \label{eq6}
\end{equation}
\begin{equation}
\normal\times(\H_2-\H_1) = \frac{4\pi}c \mathbf{K}, \label{eq7}
\end{equation}
where ${\bf K}$ is the surface current density. The above jump conditions
define the refraction of magnetic lines on the material interface.
We can assume $\mu = 1$ for most fluids. Notice that the surface current density
$\mathbf{K}$ corresponds to a current localized in a thin fluid boundary layer 
($\delta$-functional current) which is non-zero only for superconducting materials.
The current density in fluids at normal conditions is distributed in the 3D volume 
and  $\mathbf{K}=0$. Therefore, the equations (\ref{eq6},\ref{eq7}) simply 
require the continuity of the normal and tangential components of the 
magnetic field.


The behavior of a fluid in the presence of electromagnetic fields is governed to a large 
extent by the magnitude of the conductivity. 
For fluid at rest (\ref{eq4}) reduces to the diffusion equation
\begin{equation}
\f{\p \B}{\p t}=\f{c^2}{4\pi\mu\sigma}\Delta\B \label{eq8}
\end{equation}
This means that an initial configuration of magnetic field will decay with typical
diffusion time 
\[
\tau=\f{4\pi\mu\sigma L^2}{c^2},
\]
where $L$ is a characteristic length of the spatial variation of $\B$. Despite being 
good enough conductors, most of liquid metals including mercury are characterized by
small diffusion times (33 microseconds for a mercury droplet of 1 cm radius) compared 
to some solid conductors (1 sec for a copper sphere of 1 cm radius). Therefore the 
magnetic field penetration in such liquid conductors can be considered as an 
instantaneous process.

Another crucial phenomena for MHD flows of compressible conducting fluids is the 
propagation of Alfven waves. For mercury at room temperature the Alfven velocity
\[
\mathbf{v}_A = \f{\B_0}{\sqrt{4\pi\rho_0}},
\]
where $\B_0$ and $\rho_0$ are unperturbed (mean) values of the magnetic induction
and density of the fluid, respectively, is $[B_0\mbox{(Gauss)}/13.1]$ cm/sec. 
This is a small number compared with the speed of sound of $1.45\times 10^5$ cm/sec 
even for the magnetic field of 20 T. In many cases, however, it is not desirable
to compute Alfven waves explicitly in the system. If, in addition, both the magnetic 
field diffusion time and the eddy current induced magnetic field are small, an 
assumption of the constant in time magnetic field can be made. The current density
distribution can be obtained in this case using Ohm's law
\begin{equation}
\J = \sigma\left(-\mbox{grad}\phi+\f1c\vel\times\B\right). \label{eq9}
\end{equation}
Here $\phi$ is the electric field potential.
The potential $\phi$ satisfies the following Poisson equation
\begin{equation}
\Delta\phi=\f1c\mbox{div}(\vel\times\B), \label{eq10}
\end{equation}
and the Neumann boundary conditions
\[
\left. \f{\p\phi}{\p \normal}\right|_\Gamma=
\f1c(\vel\times \B)\cdot \normal,
\]
where $\normal$ is a normal vector at the fluid boundary $\Gamma$. This approach 
is applicable for the study of a liquid metal jet moving in a magnetic field.

We shall use also the following simplification for the modeling of a thin jet
moving along the solenoid axis. Let us consider a ring of liquid metal of radius $r$
that is inside a thin jet moving with velocity $u_z$ along the axis of a solenoid 
magnet. The magnetic flux $\Phi=\pi r^2B_z$ through the ring varies with time because 
the ring is moving through the spatially varying magnetic field, and because the
radius of the ring is varying at rate $u_r=dr/dt$. Therefore, an azimuthal electric 
field is induced around the ring:
\begin{eqnarray*}
2\pi rE_\phi=-\f1c\f{d\Phi}{dt}
&=&-\f{\pi r^2}c \f{dB_z}{dt}-\f{2\pi ru_rB_z}c
\\
&=&-\f{\pi r^2u_z}c \f{\p B_z}{\p z}-\f{2\pi ru_rB_z}c.
\end{eqnarray*}

This electric field leads to an azimuthal current density
\begin{equation}
J_\phi=\sigma E_\phi=-\f{\sigma r u_z}{2c}\f{\p B_z}{\p z}-\f{\sigma u_rB_z}c, \label{eq11}
\end{equation}
which defines the Lorentz force in the momentum equation (\ref{eq2}) and leads to the
distortion of the jet moving in a non-uniform magnetic field.

The linear stability analysis of thin conducting liquid jets moving along the axis of a uniform
magnetic field \cite{Chandra} and the corresponding analysis for the Muon Collider target
\cite{Gl4} show that an axial uniform field tends to stabilize the jet surface. The influence of 
a strong nonuniform field must be studied by means of the numerical simulation.


\subsection{Numerical Implementation}


In this section, we shall describe numerical ideas implemented in the FronTier MHD code.
FronTier represents interfaces as lower dimensional meshes moving through a 
volume filling grid \cite{Gl1}. The traditional volume filling finite difference grid supports smooth 
solutions located in the region between interfaces. The location of the discontinuity and 
the jump in the solution variables are defined on the lower dimensional grid or 
interface. 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
\cite{Gl1,Gl3} 
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 \cite{Gl3}. 

We solve the hyperbolic subsystem of the MHD equations, namely the equations (\ref{eq1}-\ref{eq3}),
on a finite difference grid in both domains separated by the free surface  using FronTier's 
interface tracking numerical techniques.
Some features of the FronTier hyperbolic solvers include the use of high resolution methods 
such as MUSCL, Godunov and Lax-Wendroff with a large selection of  Riemann solvers such as 
the  exact Riemann solver, the Colella-Glaz approximate Riemann solver,  
the linear US/UP fit (Dukowich) Riemann solver, and the Gamma law fit. We use 
realistic models for the equation of state such as the polytropic and stiffened 
polytropic equation of state, the Gruneisen equation of state, and the SESAME tabular 
equation of state.

The evolution of the free fluid surface is obtained through the solution of the 
Riemann problem for compressible fluids \cite{Gl1,Menikoff}. 
Notice that since we are primarily interested 
in the contact discontinuity propagation, we do not consider the Riemann problem
for the MHD system and therefore neglect elementary waves typical for 
MHD Riemann solutions.

We have developed the following numerical approaches for solving 
the equation (\ref{eq4}) for the magnetic field evolution and the Poisson equation (\ref{eq10}): 

{\bf Finite difference solvers based on regular grids.}
We have developed a parallel elliptic solver for the current density distribution 
equation (ref{eq10})
based on finite difference technique on
rectangular grid in regular domains. We have performed smoothening of discontinuities of 
physics fields
and applied the cubic spline interpolation of the fields within several grid blocks.
The corresponding parallel iterative solvers for linear systems of equations explore different 
Krylov subspace and preconditioning methods and are based on PETSC libraries \cite{PETSC}.
The solver is very effective for modeling multifluid MHD phenomena for fluids with smooth 
interfaces without small length scale perturbations. 
For MHD flows with unstable interfaces, we use FE solvers on
meshes dynamically conforming to evolving interfaces.


{\bf Finite element solvers.}
Finite element solvers are based on grids dynamically conforming to the moving interface.
The grid generation method is based on the point-shifted grid technique of McBryan 
\cite{McBryan}.
For the elliptic and parabolic systems considered here we favor the mixed-hybrid finite 
element formulation of Chavent and Jaffre \cite{Chavent}. We utilize Glowinski-Wheeler domain 
decomposition \cite{Glowinski}.  The linear systems obtained for the subdomain "interior problems" 
are solved using direct methods which preconditioned conjugate gradients for the global 
"wire basket" problem.  Preconditioners such as vector probing and balancing 
are used.



\subsection{Modeling  Thermodynamics Properties of Fluids}

Material properties strongly influence the structure and dynamics of waves
in any continuous mechanical system. In particular, elementary waves in compressible fluids
defined by Riemann problems depend both quantitatively and qualitatively on the property 
of an EOS. 

{\bf SESAME equation of state for mercury.} To model accurately the interaction of the mercury target with
proton pulses, we created a tabulated equation of state for mercury
in a wide temperature - pressure domain which includes the
liquid-vapor phase transition and the critical point. The FronTier code uses
SESAME Equation-of-Setate Library \cite{SESAME}, developed at Los Alamos National Laboratory, 
to work with such tabulated EOS format. Necessary
data describing thermodynamic properties of mercury were obtained courtesy of T. Trucano 
of Sandia National Laboratory. The behavior of isotherms of the specific internal energy, pressure 
and entropy of mercury obtained using this equation of state are shown in Figures 
\ref{figure_energy} - \ref{figure_entropy}. The computational domain included the two phase
region for mercury and the critical point. 


%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/energy.eps,width=0.6\textwidth}
\caption[Isotherms of mercury specific internal energy \textsl{vs.} density]  ]{Isotherms of mercury specific internal energy as a function of density in the mixed liquid - vapor domain.}
\label{figure_energy}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/pressure.eps,width=0.6\textwidth}
\caption[Isotherms of mercury pressure \textsl{vs.} density]{Isotherms of mercury pressure as a function of density in the mixed liquid - vapor domain.}
\label{figure_pressure}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/entropy.eps,width=0.6\textwidth}
\caption[Isotherms of mercury entropy \textsl{vs.} density]{Isotherms of mercury entropy as a function of density in the mixed liquid - vapor domain.}
\label{figure_entropy}
\end{figure}

{\bf Isentropic EOS for cavitation flows.} The tabulated EOS model allows an accurate description 
of equilibrium thermodynamic processes in mercury. However processes which occur in the 
mercury jet interacting with proton pulses go beyond the assumptions of equilibrium 
thermodynamics. Namely a regions of negative pressure (liquid tension) appear during the 
interaction which lead to the formation of cavities. To resolve such problem of 
negative pressures, we have developed an analytic two phase EOS model for an isentropic 
regime of the flow. The model is based on the polytropic EOS for the gas phase, the 
stiffened polytropic EOS for the liquid phase \cite{Menikoff} and the following
EOS type model (pressure - density relation):
\be
P = P^{sat}_l + P_{gl}\mbox{log}\left[\f{\rho_ga_g^2(\rho_l+\alpha(\rho_g-\rho_l))}
{\rho_l(\rho_ga^2_g-\alpha(\rho_g a^2_g-\rho_la^2_l))}\right],
\en
where
\[
P_{gl} = \f{\rho_ga^2_g\rho_la^2_l(\rho_g-\rho_l)}{\rho_g^2a_g^2 - \rho_l^2a_l^2},
\]
Here $\rho_g, \, p_g, \, a_g,\, \rho_l, \, p_l, \, a_l$ are the density, pressure, and sound
speed of the saturated gas and liquid respectively,   
and the void fraction $\alpha$ is
\[
\alpha = \f{\rho - \rho_l}{\rho_g - \rho_l}.
\]
We have developed the corresponding software library for FronTier. 


\subsection{Numerical Simulation Results}

{\bf MHD studies of the Muon Collider target}
We shall present here numerical simulation results of thin jets of 
conducting fluid moving in highly nonuniform magnetic fields. In this numerical experiment, 
a 1 cm radius liquid jet is sent into a 20 T solenoid with the velocity 90 m/sec along the
solenoid axis. The density of the liquid is $1 \,g/cm^3$, the electric conductivity
is $10^{16}$ in Gaussian units, and the initial 
pressure in the liquid is 1 atm. 
The electrically and magnetically neutral gas outside the jet has density 0.01 g/cm and 
the same initial pressure. The thermodynamic properties of the ambient gas were modeled 
using the polytropic equation of state \cite{Courant} with  the ratio of specific heat
$\gamma = 1.4$ and the ideal gas constant $R = 1$.  The properties of the liquid jet 
were modeled using the stiffened polytropic equation of state 
with the Gruneisen exponent $ = 5$ and the stiffened gas constant 
$P_\infty = 3\cdot 10^9g/(cm\cdot sec^2)$. The field of a magnetic coil of rectangular 
profile and $8\times8\times20$ cm size was calculated using exact analytical expressions. 

A set of images describing the evolution of the liquid jet as it enters 
and leaves the solenoid
is depicted in Figure \ref{fast_jet}. The strong nonuniform magnetic field near the 
solenoid entrance squeezes and distorts the jet. The  magnetic field 
outside the solenoid stretches the jet which results in the jet breakup. 
Notice that 
these simulations use the expression for the current density distribution
(\ref{eq11}) which looses the quantitative accuracy when the jet is close
to the breakup. Numerical results obtained using accurate numerical solutions
of equation (\ref{eq10}) are given below.


%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/jet_0200.ps,width=0.8\textwidth}
\centering\epsfig{file=roman/jet_1000.ps,width=0.8\textwidth}
\centering\epsfig{file=roman/jet_2000.ps,width=0.8\textwidth}
\caption{Liquid metal jet in a 20~T solenoid.}
\label{fast_jet}
\end{figure}


The numerical simulation  demonstrates the influence
of strong magnetic field gradients on the jet distortion. To avoid such instabilities
of the mercury target during its propagation in a 20 T magnetic field, the 
nozzle of the mercury target was placed inside the main 20 T resistive magnetic coil
(see  Figure \ref{magnet-e951}). Therefore the jet will not meet 
strong magnetic field gradients before the interaction with the proton pulses. 

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/magnet-e951.ps,width=0.9\textwidth}
\caption{Schematic of the Muon Collider magnet.}
\label{magnet-e951}
\end{figure}

We have performed numerical simulations of an off-axial  mercury jet moving in 20T, 15T, and 10T 
solenoids using accurate numerical solution for the current density distribution 
(see Figures \ref{target_mhd_20T} -- \ref{target_mhd_10T}).
Our results show that the magnetic field influence will not lead to significant jet deformations. 
Currently, the jet is not modeled as coming form the nozzle but is initialized as a long cylinder
(see Figure \ref{magnet-jet}).
The perturbation of the tail part of the jet (Figure \ref{target_mhd_20T}) is caused
by the nonuniform field behind the solenoid. 

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/magnet-jet.ps,width=0.9\textwidth}
\caption{Initialization of the mercury jet.}
\label{magnet-jet}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/jet1-000.eps,width=0.8\textwidth}
\centering\epsfig{file=roman/jet1-1000.eps,width=0.8\textwidth}
\caption{Evolution of the mercury target in 20T solenoid}
\label{target_mhd_20T}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/jet1-000.eps,width=0.8\textwidth}
\centering\epsfig{file=roman/jet2-1250.eps,width=0.8\textwidth}
\caption{Evolution of the mercury target in 15T solenoid}
\label{target_mhd_15T}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/jet1-000.eps,width=0.8\textwidth}
\centering\epsfig{file=roman/jet3-1250.eps,width=0.8\textwidth}
\caption{Evolution of the mercury target in 10T solenoid}
\label{target_mhd_10T}
\end{figure}


{\bf Interaction of mercury jet with proton pulses}

Numerical simulations presented in this section 
aid in understanding the behavior of the target under the influence of the proton pulse, and in
estimating the evolution of the pressure waves and surface instabilities. 
We have neglected the influence of the 
external magnetic field on the hydrodynamic processes driven by the proton energy deposition.

The influence of the proton pulse was modeled by adding the proton energy density distribution
to the internal energy density of mercury at a single time step. The value and spatial
distribution of the proton energy was calculated using the MARS code \cite{Ozaki}.  
The stiffened polytropic EOS model was used for these calculations with coefficients calculated 
using the tabulated EOS for mercury.

The evolution of the mercury target during 120 mks due to the interaction with a proton
pulse is shown in Figure \ref{jet_beam}. Simulations show that
velocities of surface instabilities are in the range 20 - 60 m/sec. The 
pressure field inside the jet developed small regions of negative pressure which indicates 
a possibility for the formation of cavities. Detailed studies of the cavitation phenomena
in such a target will be done using the developed two-phase equation of state for mercury. 


%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\centering\epsfig{file=roman/sp3_0000.eps,width=0.6\textwidth}
\centering\epsfig{file=roman/sp3_0200.eps,width=0.6\textwidth}
\centering\epsfig{file=roman/sp3_0600.eps,width=0.6\textwidth}
\centering\epsfig{file=roman/sp3_1000.eps,width=0.6\textwidth}
\caption[Evolution of the target driven by the proton energy deposition]{Evolution of the target driven by the proton energy deposition. Snapshots 
are shown at the initial time and at 27 mksec, 78 mksec, and 120 mksec.}
\label{jet_beam}
\end{figure}


\section{Proposed Research}

In the future, we propose the following research:

- Improve the robustness of the dynamic grid generator and implement new
elliptic/hyperbolic solvers.

- Develop a general analytic EOS for cavitation flows.

- Numerical simulation of nozzle effects.
 
- Further studies of liquid metal jets in magnetic fields.

- Numerical simulation of the mercury target interacting with a high
energy proton beam in the presence of a strong magnetic field. 

- Numerical simulation of turbulent jets.

%
% ---- Bibliography ----
%
\begin{thebibliography}{30}
%

\bibitem{Chandra} Chandrasekhar, S.: Hydrodynamics and Hydrodynamic stability. Clarendon Press, 
Oxford (1961)

\bibitem{Chavent} G. Chavent and J. Jaffre, Mathematical Models and Finite Elements for Reservoir 
Simulation.  North Holland, 1986.

\bibitem{Chern} Chern, I.R.,  Glimm, J., McBryan, O., Plohr, B., Yaniv, S.: Front tracking for Gas 
Dynamics. J. Comp. Phys. {\bf 62} (1986) 83-110

\bibitem{Courant} Courant, R., Friedrichs, K.: Supersonic Flows and Shock Waves. Interscience, New 
York (1948)


\bibitem{Gl1} Glimm, J., Grove, J.W., Li, X.L., Shyue, K.M., Zhang, Q., Zeng, Y.: Three  
dimensional front tracking. SIAM J. Sci. Comp. {\bf 19} (1998) 703-727  

\bibitem{Gl3} Glimm, J., Grove, J.W., Li, X.L., Tan, D.: Robust computational algorithms for 
dynamic interface tracking in three dimensions. Los Alamos National Laboratory 
Report LA-UR-99-1780

\bibitem{Gl4}  Glimm, J.,  Kirk, H., Li, X.L., Pinezich, J., Samulyak, R., Simos, N.: Simulation of
3D fluid jets with application to the Muon Collider target design. Advances in
Fluid Mechanics III (Editors: Rahman, M., Brebbia, C.), WIT Press, Southampton Boston (2000) 191 - 200

\bibitem{Glowinski} Glowinski and M. Wheeler, Domain Decomposition and Mixed Finite Element 
Methods for Elliptic Problems.  In R. Glowinski, G. Golub, G. Meurant, and J. 
Periaux, Eds., Proceedings of the First International Symposium on Domain 
Decomposition Methods for Partial Differential Equations, pp. 144-172.  SIAM, 
Philadelphia, 1988.

\bibitem{SESAME} K. S. Holian, T-4 handbook of material properties data bases, 
Vol. 1c: Equations of state, Los Alamos Report LA-10160-MS, 1984.

\bibitem{Kirk} Kirk, H., et al.: Target studies with BNL E951 at the AGS. 
Particles and Accelerators 2001, June 18-22 (2001) Chicago IL

\bibitem{Landau}  Landau, L.D.,  Lifshitz, E.M.: Electrodynamics of Continuous Media. Addison -    
Wesley Publishing Co., Reading Massachusetts (1960)

\bibitem{McBryan} O. McBryan, Fluids, Discontinuities and Renormalization Group Methods.  In Brittin, 
Gustafson, and Wyss, Eds., Mathematical Physics VII, pp. 481-494.  North-Holland 
Publishing Company, Amsterdam, 1984.

\bibitem{Menikoff} Menikoff, R., Plohr, B.: The Riemann problem for fluid flow of real materials. 
Rev. Mod. Phys. {\bf 61} (1989) 75-130

\bibitem{PETSC} S. Balay, W. Gropp, L. McInnes and B. Smith", PETSc Users Manual,
ANL-95/11 - Revision 2.1.1, Argonne National Laboratory, 2001.

\bibitem {Ozaki}
Ozaki, S., Palmer, R., Zisman, M., Gallardo, J. (editors)
Feasibility Study- II of a Muon-Based Neutrino Source.
BNL-52623 June 2001; available at http://www.cap.bnl.gov/mumu/studyii/FS2-report.html

\end{thebibliography}
%\end{document}




















