Simulation of hydrodynamic tunneling induced by high-energy proton beam in copper by coupling computer codes

The design of machine protection systems for high-energy accelerators with high-intensity beams requires analyzing a large number of failures leading to beam loss. One of the most serious failures is an accidental impact of a large number of bunches at one location, for example, due to a deflection of the particle beams by the extraction kicker magnets with the wrong strength. The numerical assessment of such an event requires an iterative execution of an energy-deposition code and a hydrodynamic code, in case the hydrodynamic tunneling effect plays an important role in the beam-matter interactions. Such calculations have been performed for the CERNLargeHadron Collider (LHC), since the energy stored in the LHC beams exceeds previous accelerators by 2 orders of magnitude. This was done using the particle shower code FLUKA and the hydrodynamic code BIG2 [Tahir et al., Phys. Rev. STAccel. Beams 15, 051003 (2012)]. Later, simulations for a number of cases for other accelerators at CERN and the Future Circular Collider (FCC) were performed [Tahir et al., Phys. Rev. Accel. Beams 19, 081002 (2016)]. These simulations showed that the penetration depth of the beam in copper or graphite could be an order of magnitude deeper when considering the hydrodynamic tunneling, compared to a static approximation. BIG2 is very efficient and used at GSI. It is always advantageous to cross-check the simulation results with different codes, because dedicated experiments are very complex or even impossible to carry out. For this purpose, it was decided to study the use of a commercial tool (Autodyn) for such calculations and to compare the results with previous work. This paper reports a benchmarking study against beam experiments performed at the HiRadMat (High-Radiation to Materials) facility using beams at 440 GeV from the Super Proton Synchrotron. Good agreement has been found between the simulation results and the experiments as well as previous simulations with FLUKA and BIG2 [Tahir et al., Phys. Rev. E 90, 063112 (2014)], particularly in terms of the penetration depth of the beam in copper. This makes the coupling of FLUKA and Autodyn an alternative solution to the simulation of the hydrodynamic tunneling, at least in the parameter range of this case study. Studies with other parameters are planned for FCC and other high-beam-power accelerators.


I. INTRODUCTION
In the Large Hadron Collider (LHC), the energy stored in one of the two counterrotating proton beams reaches 362 MJ, under nominal beam parameters, i.e., 2808 bunches at 7 TeV with a bunch intensity of 1.15 × 10 11 [1].This energy is sufficient to melt about 500 kg of copper.To increase its annual integrated luminosity by a factor of 10, the LHC will be upgraded to the High Luminosity LHC (HL-LHC) by around 2025 [2].The bunch intensity will increase to 2.2 × 10 11 , benefiting from the LHC injectors upgrade (LIU) project [3].In addition, a higher-energy LHC (HE-LHC) is under study to make the colliding proton energy as high as 13.5 TeV (center-of-mass energy 27 TeV) [4,5].Compared with the LHC, the energy stored in the beams will be doubled (for HL-LHC) or more (for HE-LHC).Moreover, the conceptual design of the Future Circular Collider (FCC) is being carried out in a global collaboration hosted by CERN [6].For the proton-proton collider (FCC-hh), the beams will be accelerated up to 50 TeV in a 100 km tunnel.The nominal number of bunches in one beam is 10 400 and the bunch intensity is 1.0 × 10 11 , leading to a beam energy of 8.3 GJ, about 20 times higher than that of the LHC.In these high-energy colliders, any uncontrolled release of the beam energy could potentially result in severe damage to accelerator components.Thanks to the machine protection system [7,8], the LHC has been running safely with high availability and reliability and hence impressive luminosity performance [9].At present, the operating beam energy in the LHC has reached 6.5 TeV, very close to the designed top energy.
For the design and operation of the machine protection system, it is important to classify different beam loss scenarios according to the beam lifetime and then apply proper strategies to avoid equipment damage [10][11][12].For the LHC, under normal operating conditions, the beam lifetime is typically longer than 10 hr.Steady beam losses are unavoidable due to machine imperfections and protonproton collisions.To protect sensitive components from beam impact and magnets from quench, a multistage collimation system is installed to continuously clean beam halos via betatron and momentum collimations defining aperture limitations.Under some particular operating conditions, e.g., during the change of optics, the beam lifetime may be reduced considerably, and the minimum allowed beam lifetime is 12 min by design.The corresponding beam loss power will be 500 kW.In such situations, safe operation of the LHC is possible only with the help of efficient collimation systems.Fast equipment failures, e.g., magnet powering failures, may reduce the beam lifetime to seconds or even milliseconds.Once a failure or an abnormal beam parameter is detected by the monitoring system, a beam dump request will be triggered and transmitted to the dumping system via the interlock system.The beam will be dumped in time before the collimators and other equipment can be damaged.Before reaching the beam dump block, the beam is diluted by transverse dilution kicker magnets.The dump block is designed to absorb the entire beam energy without degradation.Singlepassage beam loss can possibly happen during injection or extraction.In this case, protection relies on dedicated beam absorbers.
It is of fundamental importance to study beam impacts on the beam intercepting devices (BIDs) such as collimators, absorbers, dump blocks, and other targets that interact with the beam frequently [13][14][15].The main purposes are to evaluate their damage threshold and to develop advanced materials suitable for the protection task.The first step is to simulate the energy deposition distribution in the target using a Monte Carlo code such as FLUKA [16][17][18], MARS [19], and Geant4 [20].The energy deposition is then used as the input of internal energy for mechanical analysis.As illustrated in Ref. [21], the dynamic response induced by a particle beam in material can be classified into the elastic regime, the plastic regime, and the shock wave regime, depending on the deposited power density and the duration of the interaction.In the elastic regime and the plastic regime, the problem can be solved by standard finite element (FEM) tools like ANSYS [22] using implicit time-integration schemes.In the shock wave regime, where the deposited energy typically exceeds 10 kJ=cm 3 , wave propagation codes or hydrocodes with explicit timeintegration schemes must be employed, such as ANSYS- Autodyn [23], LS-Dyna [24], and BIG2 [25].Complex material constitutive models including equations of state (EOS), strength models, and failure models are required, since the target will usually undergo phase transitions or even be physically displaced in the most impacted region.In case a large number of high-energy bunches are successively lost at the same point, energy deposited by the heading bunches produces an outgoing radial shock wave which reduces the density along and around the beam axis.The subsequent bunches and their secondary hadronic shower will penetrate deeper and deeper into the target, while the energy deposition per proton in the upstream part is reduced, because the inelastic interaction is density dependent for given atomic and mass numbers.This effect is called hydrodynamic tunneling [26][27][28].
For most cases, the number of particles lost in the aperture is relatively low.Therefore, a one-way coupling between the energy deposition code and the hydrodynamic code is sufficient.It means that the energy distribution calculated under nominal density is used for the overall beam impact analysis, without taking into account the density change and hence the dose redistribution in the target, as illustrated in Refs.[14,[29][30][31].This is valid as long as the density change due to the beam impact is less than a few percent.In this paper, we focus on worst-case failures where the entire beam or a bunch train with many bunches is lost at the same point.Such a beam accident could happen during injection and extraction, if, for example, the kickers or septum magnets deflect the beam by a wrong angle due to failures in the energy-tracking system or in the magnets themselves [32].In 2004, the full Super Proton Synchrotron (SPS) beam (288 bunches, 3.4 × 10 13 protons, 450 GeV) was extracted with a wrong angle due to a switch-off of the extraction septum [33].The stainless steel vacuum chamber and one magnet in the transfer line were severely damaged and had to be replaced.Interlocking systems and operational procedures were modified afterwards to avoid similar incidents.Another possible failure mode is when the beam is extracted towards one spot (without energy dilution) on the beam dump block due to a dilution kicker malfunction [34].In general, an accidental loss of the whole beam on a single spot is very unlikely, but it is still worth performing detailed studies concerning the consequences.The aforementioned hydrodynamic tunneling effect will most likely play an important role in the beam-matter interactions for such failures.Therefore, an iterative coupling of the two kinds of numerical codes is necessary.
Before we present the calculation of the hydrodynamic tunneling with novel tools, a brief summary of existing methods is made.Two-and three-dimensional (2D and 3D) hydrodynamic calculations were performed to understand the behavior of the Superconducting Super Collider (SSC) graphite beam dump under abnormal beam aborts at 20 TeV [26].Since the designed proton-beam pulse length of 300 μs was comparable to the characteristic time of expected hydrodynamic motions, the energy deposition code MARS and the hydrodynamic codes MESA [35] and SPHINX [36] were combined.Using Autodyn, the structure response of the LHC beam dump absorber block in the event of a total beam dilution failure was simulated [37].However, 2D loading of the beam was simply modeled by depositing the specific internal energy as a function of the radial and axial position, and local density, using a discretized numerical model.Two case studies coupling FLUKA and LS-Dyna iteratively have been reported.In one case, it was assumed that a cylindrical tungsten target was facially irradiated by 30 LHC bunches at 7 TeV [14].In another case, the impact of 60 LHC bunches at 7 TeV on a tungsten collimator was simulated [38].A voxel-based model (3D pixel) was used for the parallelepiped target.By coupling FLUKA and BIG2, some of the worst possible beam loss scenarios have been studied for the SPS, the LHC, and the FCC.The targets were homogeneous cylinders made of copper, graphite, and tungsten.2D hydrodynamic calculations were performed using BIG2.At early stages, the energy deposition was normalized with axial line density (analytic approximation) without running new FLUKA simulations [28,[38][39][40][41][42].Since 2012, the two codes have been run iteratively [43].BIG2 uses the energy deposition data from FLUKA as the energy input, while FLUKA uses the modified target density distribution from BIG2 to calculate a new energy-deposition map which again serves as the input of BIG2.Within each iteration, the maximum allowed density drop is 10%-15% to reach a compromise between computation time and accuracy.Case studies showed that the penetration depth of a full nominal LHC beam with an rms beam size of σ x;y ¼ 0.2 mm was about 35 m in copper due to the hydrodynamic tunneling effect [28,42], not 2 m as indicated by the static simulation using FLUKA.In graphite, the penetration depth reached 25 m with σ x;y ¼ 0.5 mm compared to the static value of 3 m [43].Recent simulations illustrated that the 50 TeV FCC beam would penetrate 350 m in copper when σ x;y was 0.2 mm [44] and 1.3 km in water if σ x;y was 0.4 mm [45].The water target was studied for the FCC beam to examine the possibility of a water beam dump without the need for dilution kickers.The study suggested that the beam size must be increased from 0.4 mm to centimeters to make the water tank shorter and allow the survival of a beam window separating the beam transfer line and water.Otherwise, a beam with small beam sizes could easily melt the beam window.Dedicated experiments at the High-Radiation to Materials (HiRadMat) facility with the 440 GeV SPS beams have validated the numerical simulation method coupling FLUKA and BIG2 [46][47][48][49].
In these studies, no suitable scaling law has been suggested to predict the hydrodynamic tunneling range for varying beam parameters such as beam energy, beam size, bunch intensity, and time structure of the beam, since this is a complex nonlinear multiphysics problem.Specific numerical or experimental studies are needed on a case-bycase basis.Dedicated experiments are very complex or even impossible to perform, e.g., for the FCC beams.In many cases, we rely on sophisticated numerical techniques.
As the demand still grows, including the feasibility study of the water beam dump for FCC under various beam parameters, and many other cases for the LHC, the HL-LHC, the HE-LHC, etc., it was decided to develop a new way for the simulation of the hydrodynamic tunneling using different codes.The availability of different numerical solutions allows us to carry out a direct comparison among them.This provides more reliable data for the design and operation of advanced machine protection systems.This will also serve as an alternative solution when other codes are unavailable or not suitable; e.g., BIG2 is a 2D hydrodynamic code with limitations to calculate the case of nonround beams.Therefore, we have developed a new way to simulate hydrodynamic tunneling using FLUKA and ANSYS-Autodyn.Autodyn is an explicit analysis tool for modeling nonlinear dynamics of solids, fluids, and gases as well as their interactions in two or three dimensions especially in short-duration, highly nonlinear transient events.
In this paper, we report our benchmarking case study where beam and target parameters in the HiRadMat experiment mentioned earlier are used.Our results from the coupling of FLUKA and ANSYS-Autodyn are in good agreement with the test results and with the results from the simulation using FLUKA and BIG2.
In Sec.II, the beam and target parameters used in the experiment and simulation will be described.The simulation strategy and algorithm improvement of the benchmarking study will be presented in Sec.III together with corresponding results.For an easy comparison with previous studies, we decided to show the results in a similar way as in previous papers [47,49].All the results in Sec.III shown in Figs.2-11 were obtained using FLUKA and Autodyn.In Sec.IV, a comparison of the results between different numerical codes and the beam test will be shown, mostly in terms of the penetration depth of the beam into the target.Conclusions drawn from this work will be summarized in Sec.V.

II. BEAM AND TARGET PARAMETERS
As the second-largest machine in CERN's accelerator complex, the SPS receives protons or ions from the Proton Synchrotron and then accelerates them to provide beams not only for the LHC, but also for a number of present and past fixed-target research programs including the COMPASS experiment, the CNGS experiment, and the AWAKE experiment.Among others, the unique beam test facility HiRadMat is designed to study beam impacts on materials and accelerator components using the SPS beams.It is able to deliver up to 288 proton bunches at 440 GeV with a maximum bunch intensity of 1.7 × 10 11 , i.e., a total of 4.9 × 10 13 protons storing an energy of 3.4 MJ.The typical bunch length is about 0.5 ns, while the transverse beam size can be tuned from 0.1 to 2 mm (rms value) and the bunch spacing can be adjusted to 25, 50, 75, or 150 ns.Full details about the HiRadMat facility can be found in Ref. [50].
In order to confirm the existence of the hydrodynamic tunneling effect and validate the numerical approach coupling FLUKA and BIG2, dedicated experiments were performed at the HiRadMat facility in 2012.Three copper targets were facially irradiated by beams from the SPS respectively.Each target comprised 15 copper cylinders separated by a 1 cm gap between each other to allow for visual inspection after the beam test.Each cylinder had a radius of 4 cm and a length of 10 cm.Information on the experiment assemblies can be found in Ref. [46].In the experiments, the proton energy was 440 GeV, the bunch intensity 1.5 × 10 11 , the bunch length 0.5 ns, and the bunch spacing 50 ns.Target 1 was irradiated by a beam of 144 bunches with σ x;y ¼ 2 mm.Two beams with the same rms beam size of σ x;y ¼ 0.2 mm irradiated targets 2 and 3, while the number of bunches were 108 and 144, respectively.The protons were delivered in three or four packets (for 108 or 144 bunches).The packets were spaced by 250 ns, and each packet consisted of 36 bunches.
After eight months of cool down, the targets were opened for visual inspection in February 2013.Droplets and splashes of molten and evaporated copper were found on the copper cylinders, the aluminum housing, and the front aluminum caps.Later, the targets were dissected into finer pieces for visual and microscopic inspections to establish the precise penetration depths of protons and their hadronic shower.Moderate and significant tunneling was found in target 2 and target 3, respectively.The measured melting fronts reached 58, 79.5, and 85 cm, respectively, for targets 1-3.To have a direct observation, Fig. 1 reproduces a few pictures presented in Ref. [48] showing the front and back faces of the first, fifth, and ninth cylinders of target 3 after irradiation.
For the simulation with FLUKA and Autodyn, we chose the case of target 3 as a benchmarking study, since the hydrodynamic tunneling is more visible than in target 1.The case of target 2 was covered during the simulation when 108 bunches were delivered.In the simulations, we used exactly the same beam parameters as described including the time structure.As a simplification, an extended copper target with a radius of 4 cm and a total length of 150 cm was applied.The 1 cm gap was not taken into account, since the temperature and pressure gradients are much stronger in the radial direction than in the axial direction, as suggested in the previous study [47,49].

III. SIMULATION PROCEDURE AND RESULTS
First, the detailed simulation procedure coupling FLUKA and ANSYS-Autodyn iteratively is described.In this study, target reconstruction in FLUKA is different compared to previous work.Numerical simulation results are then reported.

A. Simulation methodology
FLUKA is a multipurpose Monte Carlo simulation tool for the calculation of particle transports and interactions with matter.Primary particles irradiate the target and then produce a particle shower.These secondary particles deposit part of their energy inside the target, while some of them escape from the target with residual energy.During the simulation, FLUKA takes a primary particle from the predefined beam, starts the transport, and repeats until the predetermined number of primary histories is reached.Finally, the results are normalized to one primary proton for easy scaling.In our case, the distribution of dose (energy deposition per unit mass per proton) in units of GeV=ðg pÞ has been scored under a circular symmetry condition, i.e., in the r-z geometry, since a round beam is irradiating perpendicularly the front face of the copper target along the target and the beam axis.
ANSYS-Autodyn is a commercially available 2D and 3D software for explicit analysis of nonlinear dynamics.The dynamics of continuous media are described by a set of differential equations derived from the principles of conservation of energy, mass, and momentum.Analytical or tabular EOS correlate the pressure (stress), density (specific volume), and temperature (internal energy) for a given material over a wide range in different physical phases.In Autodyn, a tabular SESAME EOS is implemented directly [51,52].The deviatoric behavior of a material is usually expressed by the strength model, while a dynamic failure model is typically used to determine the structural limits.For the solid copper target, the empirical Johnson-Cook strength and failure model is adopted, which is particularly suitable for metals and ductile materials [53,54].From the known state of the material, Autodyn calculates the state of the material after a certain time step using an explicit time integration scheme.Lagrange and Euler formulations are two fundamental approaches to the spatial discretization of continuous media.The Eulerian description uses a grid, which remains spatially fixed in time.The energy, mass, and momentum may flow across cell boundaries to avoid mesh tangling when large material deformations occur.Also, note that the Eulerian mesh does not foresee material failure.In our case, the Lagrangian description that uses a set of points attached to the material is employed for the copper target taking advantage of its computational efficiency.
With FLUKA and Autodyn, we modeled the target as a 150-cm-long copper cylinder with a radius of 4 cm.The calculations were carried out under the axisymmetric condition in a cylindrical coordinate (r-φ-z) system, since the beam was round in transverse.Note that 3D simulation can also be executed coupling the two codes for a nonround beam.The workflow is as follows.
(i) A first Monte Carlo simulation was performed using FLUKA with a nominal density (8.94 g=cm 3 ) of solid copper.The locally deposited heat was estimated by multiplying the dose (per proton) by the total number of protons.The most impacted zone was thus predicted, where the material was expected to undergo density FIG. 2. Two-dimensional distributions of the dose in units of GeV=ðg pÞ in the r-z plane calculated by FLUKA for bunches delivered at different times: (a) using a nominal solid-copper density of 8.94 g=cm 3 at t ¼ 0 μs, when bunch 1 started impacting the target, (b) using the density distribution provided by Autodyn at t ¼ 2 μs, just before bunch 37, (c) using the density distribution provided by Autodyn at t ¼ 4 μs, just before bunch 73, and (d) using the density distribution provided by Autodyn at t ¼ 6 μs, just before bunch 109.Target axis (cm) FIG. 3. Dose along the target axis for bunches delivered at different times t ¼ 0, 2, 4, and 6 μs.
changes later.For the present study, the zone was defined from r ¼ 0 to 0.5 cm and from z ¼ 0 to 100 cm, where the target was divided into 2500 fine regions.Each fine region had a radial thickness of Δr ¼ 0.1 mm (σ x;y =2) and a length of 2 cm.Each region was assigned later a density corresponding to the modified density distribution from Autodyn simulation.(ii) A map of the dose distribution from FLUKA was imported into Autodyn as an internal energy load via a user subroutine to calculate the thermodynamic and hydrodynamic processes in the target.Using the same dose map, a certain number of bunches was successively simulated until the maximum density drop along the target axis reached 10%-15%.An updated density distribution just before a subsequent bunch arrived was then exported for the next FLUKA modeling.
(iii) Before executing a new FLUKA calculation, the target geometry was rebuilt by assigning densities to the fine regions based on the updated density map.This was the most crucial procedure.To do so, two steps of data processing were needed.First, density data within one region were merged by means of averaging.When a mesh distortion due to material deformation was evident, a linear interpolation was applied to estimate the density value at the desired location.Second, the merged density values were categorized into 100 predefined density levels ranging linearly from 0.1 to 10.0 g=cm 3 .Note that an individual material must be defined for each density level, and in maximum 700 materials can be handled in FLUKA by default.The upper limit (700) was not used, since it is not  meaningful to distinguish two densities with a difference of less than 1% according to our FLUKA experience.It is worth mentioning that an extension of the zone with fine regions will be necessary if the density on the boundaries deviates from the nominal density by more than 10%, even though this did not happen in the present study.Note that some of the fine regions might need to be split into smaller ones to limit the density gradient between neighboring regions for a better simulation accuracy.Also, note that the number of regions should not exceed 20 000 in FLUKA by default.Therefore, regions with the same density were merged into one region (regardless of whether or not they were adjacent to each other) via Boolean operations when necessary.
(iv) After the construction of the new target geometry, a new FLUKA simulation was carried out to provide a modified energy deposition distribution for the following mechanical simulation of Autodyn.
(v) FLUKA and Autodyn were hence run iteratively until all bunches had been delivered onto the target.From the Autodyn results, the physical state of the target during and after the beam impact can be derived.
Compared with the previous method coupling FLUKA and BIG2 [55], the general principle is similar, but this way FIG.6. Two-dimensional distributions of temperature in the r-z plane calculated by Autodyn (a) at t ¼ 1.8 μs, after 36 bunches have been delivered, (b) at t ¼ 3.8 μs, after 72 bunches, (c) at t ¼ 5.8 μs, after 108 bunches, (d) at t ¼ 7.8 μs, after 144 bunches, (e) at t ¼ 10 μs, 2.2 μs after the delivery of the beam, and (f) at t ¼ 20 μs, 12.2 μs after the delivery of the beam.
of performing the calculations with FLUKA has two main advantages.First, predefined density levels are employed in all iterations instead of defining discrete density levels in each iteration.In this way, the simulation accuracy of FLUKA is not influenced, whereas data analysis and target modeling are simplified significantly, since, for any new density level, a new material together with a new compound has to be defined in FLUKA.The definition process of new materials and compounds is rather time consuming.Second, we are able to assign one material to different regions, so that the number of regions could be up to 20 000.In contrast, as reported previously [55], one material could only be assigned to one region, which limited the total number of regions to be less than 700 due to the fact that FLUKA does not support more than 700 materials at once by default.As a result, regions with similar (or the same) densities must be merged accordingly, which is another time-consuming process.
In FLUKA simulations, the equivalent number of proton primaries was around 150 000 to control the relative statistical error below 5% at the end of FLUKA run.The spatial dose distribution was scored in a regular binning.The angular coordinate from 0 to 2π was spanned in one bin, taking advantage of the circular symmetry condition.A radial bin size of 0.05 mm (σ x;y =4) and an axial bin size of 2 mm were adopted.The FLUKA simulations showed that 60% of the particle energy was absorbed in the copper target, while 40% escaped.The total beam energy of 144 bunches is 1.5 MJ, which means that the total energy deposited in target 3 is about 0.9 MJ.
In Autodyn, the time structure of the beam was taken into account by using a calculation time step of 0.5 ns.As a reasonable approximation, the energy deposition for each proton was instantaneously loaded onto the target as internal energy, considering that the proton and its hadronic shower would deposit their energy throughout the target only in a few nanoseconds.During the bunch length of 0.5 ns, the energy load was ramped from 0 to the energy deposited by the entire bunch.The Lagrangian mesh might distort with the material it models due to forces from neighboring elements.For the present study, the mesh nodes moved 0.1 mm in maximum in the most impacted area in the last iteration step.The density gradient was small between neighboring meshes.Therefore, a linear interpolation was performed to estimate the density of concerned spatial positions where the mesh distortion was evident.This was considered a reasonable approximation.For higher-energy and higher-intensity beam impact on materials, where a significant element distortion will likely occur leading to a large energetic error when performing data interpolation, the Eulerian mesh will be preferred, which can be employed by Autodyn alternatively.
The beam pulse lengths are 5.8 and 7.8 μs for 108 and 144 bunches, respectively.It was found that the maximum density drop reached 13% after the first 12 bunches (600 ns) had been delivered.Therefore, 12 bunches were simulated within each iteration, and in total 12 iterations were executed coupling FLUKA and Autodyn.In this way, the maximum dose change from one iteration to the next was also limited to the order of 10%.The overall simulation results are reported in the following subsection, in terms of dose, specific energy, target temperature, pressure, and density.

Energy deposition per proton
Spatial distributions of dose in units of GeV=ðg pÞ for bunches delivered at different times are shown in Figs. 2 (in the r-z plane) and 3 (along the target axis).At time t ¼ 0 μs, when bunch 1 starts impacting on the target, the nominal density of solid copper (8.94 g=cm 3 ) has been used in the FLUKA simulation.As can be seen in Fig. 2(a) and curve 1 in Fig. 3, the maximum dose is 3.8 GeV=ðg pÞ around z ¼ 12 cm.The energy deposition distributions presented in Fig. 2(b) and curve 2 of Fig. 3 were calculated by FLUKA, using the density distribution provided by Autodyn at t ¼ 2 μs, which equals the time just before bunch 37 arrived at the target.The broadened range of energy deposition indicates that the hadronic shower penetrates deeper into the target.The maximum dose has been reduced to 2.6 GeV=ðg pÞ due to the density reduction in the most heated area where the outgoing radial stress wave has been generated.At t ¼ 4 μs (before bunch 73), the range of energy deposition is further broadened as shown in Fig. 2(c), while a double peak behavior is observed with a maximum dose of about 1.8 GeV=ðg pÞ as shown in curve 3 in Fig. 3.The double peak behavior becomes more visible in Fig. 2(d) and curve 4 in Fig. 3, where the energy deposition has been simulated by FLUKA using the density data obtained from Autodyn at t ¼ 6 μs (before bunch 109).The second peak with a dose of 1.7 GeV=ðg pÞ at z ¼ 45 cm is even higher than the first one.The reason is that the density in the previously heated  region has decreased so much that the shower passes through this region with very little energy loss and deposits the majority of its energy deeper into the target.

Accumulated specific energy
The doses obtained from FLUKA were converted into specific energy as internal energy in Autodyn to simulate the mechanical behaviors of the target.The protons irradiated the target successively following the mentioned time structure of the beam.The specific energy in units of kJ=g was obtained by multiplying the dose by a certain number of protons.Figures 4 and 5 present the accumulated specific energy distributions at different times in the r-z plane and along the target axis (r ¼ 0 cm), respectively.Four points in time have been selected including t ¼ 1.8 μs (after 36 bunches have been delivered), t ¼ 3.8 μs (after 72 bunches), t ¼ 5.8 μs (after 108 bunches), and t ¼ 7.8 μs (after 144 bunches).As shown in Figs.4(a)-4(d) and curves 1-4 in Fig. 5, the beam-heated region gradually expands with time due to the successive delivery of the proton bunches.The maximum specific energy increases from about 3 to 7 kJ=g as time grows from 1.8 to 7.8 μs.However, the increase rate of the maxima slows down with  time; i.e., the accumulated specific energy is not a linear function of the bunch number.This is because of the redistribution of the energy deposition of the subsequent bunches due to the hydrodynamic tunneling effect; namely, the hadronic shower penetrates deeper into the target and heats a larger volume of material as time grows.In a static approximation, the total specific energy induced by the beam is obtained by multiplying the dose per proton at t ¼ 0 μs by the total number of protons, neglecting the hydrodynamic tunneling.For instance, curve 5 in Fig. 5 plots the total specific energy of the 144 bunches along the target axis.Compared to curve 4 in Fig. 5, in the static case as shown in curve 5, the specific energy becomes much more concentrated with an almost doubled maxima of 13 kJ=g and a much narrower shape.Taking into account the latent heat, a melting energy of 0.67 kJ=g is needed to melt copper from 300 K and 6.3 kJ=g to evaporate it.As an estimation, both curves 4 and 5 indicate that part of the copper target can be evaporated, while curve 4 shows a much deeper melting depth of about 90 cm than curve 5 (67.5 cm in the static approximation).This implies that the hydrodynamic tunneling plays an important role in the studied beam-matter interactions.

Target temperature
The high level of internal energy load results in strong heating of the target material.In Figs.6(a)-6(f) and 7, we present temperature distributions at different times in the r-z plane and along the target axis obtained from Autodyn.In addition to the four points in time used in Figs. 4 and 5, two more times are shown in Figs. 6 and 7, namely, t ¼ 10 μs (2.2 μs after the 144 bunches have been delivered) and t ¼ 20 μs (12.2 μs after the 144 bunches).The purpose is to display evolutions of temperature both during and after the beam impact.From Figs. 6(a)-6(d), a continuous extension of the beam-heated region can be clearly seen, both in the radial direction due to the accumulating specific energy and in the longitudinal direction due to the energy accumulation as well as the hydrodynamic tunneling effect.Correspondingly, curve 1 in Fig. 7 shows that a maximum temperature of 5200 K has been generated at t ¼ 1.8 μs, while curve 2 shows a maximum value of 7000 K for t ¼ 3.8 μs.Curves 3 and 4 in Fig. 7 imply that there is no significant increase in the maximum temperature as more bunches have deposited their energy, but the profiles become wider with time.In Fig. 7, there is a flat region on the right part of each curve around 1360 K, which corresponds to the melting temperature of copper.During the delivery of the beam, this plateau moves deeper in the beam direction due to the hydrodynamic process.As shown in curve 3 in Fig. 7, the melting region lies between z ¼ 77 and 83 cm after 108 bunches have been delivered (corresponding to target 2).The melting region ranges from z ¼ 89 to 95 cm after 144 bunches (target 3) as can be seen in curve 4. Figures 6(e) and 6(f) show the evolution of the temperature after the beam impact.It is seen that the high-temperature area spreads only in the radial direction due to the pressure wave propagating radially, without further longitudinal extension on the simulated time scale, because the thermal transfer is much slower than the pressure wave propagation.Consider the time constant over an element of 1 cm long, it is of the order of 10 μs for the pressure wave propagation that depends on the speed of sound in copper, while 1 s for the thermal propagation that depends on the thermal diffusivity.The radial melting range reaches a radial position of about r ¼ 0.45 cm at 20 μs.The corresponding curves 5 and 6 in Fig. 7 show that the temperature decreases in the range of z ¼ 0-70 cm, while it does not change after z ¼ 70 cm on the simulation timescale.Since the high temperature is diluted quickly in the radial direction after the beam impact, the melting front stops to penetrate deeper after the entire beam has deposited its energy in the target and after the simulated 20 μs. Figure 6, especially Fig. 6(f), shows that the target is damaged via melting and evaporation not only at the axis but also in a certain volume around the axis; i.e., a hole is generated in the target by the high-energy and high-intensity beam.

Target pressure
The high temperature induced by the energy deposition produces very high pressure in the beam-heated region.The high pressure generates a strong radial stress wave as well as a weaker axial stress wave, which is initially compressive (positive pressure).The outgoing stress wave, especially in the radial direction, leads to continuous density reduction at and around the target axis.The pressure distributions from Autodyn are shown in Figs. 8 (in the r-z plane) and 9 (along the target axis) for the same six points in time as used in Figs. 6 and 7.As shown in Fig. 8(a), a maximum pressure of 2.4 GPa is generated at t ¼ 1.8 μs, while the radial high-pressure stress wave has Target axis (cm)  already spread from the axis to a radial position of r ¼ 0.7 cm.During the extension of the high-pressure area, the maximum pressure decreases; e.g., a maximum value of 3.2 GPa has been reached at 600 ns (after the first 12 bunches).In other words, the cylindrical spreading involves a decrease in the wave amplitude along its propagation.The reason is that the same energy is spread over a volume, which is bigger and bigger in time.In the cylindrical case, this decrease in amplitude is logarithmic close to the impact point and proportional to the square root of the radius when the distance from the impact is large (i.e., the wave has propagated at a distance r → ∞).From Figs. 8(a)-8(d), it can be seen that the high pressure generated by previous bunches travels outwards with time, while the new high-pressure area generated by subsequent bunches is located around the axis but at a deeper axial position due to the hydrodynamic tunneling.and turned into tensile waves, because the shock impedance of the surrounding medium (air) is lower than that of the target.Corresponding to Figs. 8(a)-8(d), curves 1-4 in Fig. 9 show the pressure along the target axis.The maximum pressure decreases from 2.3 to 1.5 GPa as the time grows from 1.8 to 3.8 μs, while the high-pressure front moves deeper into the target.Curve 3 in Fig. 9 shows a double peak behavior in accordance with the energy distribution presented in curve 4 in Fig. 3.In Fig. 9, curve 4 has a second peak appearing deeper with a similar maximum pressure of 1.4 GPa compared to curve 3.This is due to the development of the hydrodynamic tunneling in the target.The pressure distributions after the beam impact are shown in Figs.8(e) and 8(f), as well as in curves 5 and 6 in Fig. 9, for the times at t ¼ 10 and 20 μs.
It can be seen in Fig. 8(e) that negative pressures have been created near the axis due to the radial stretching of the material and close to the target surface due to the tensile stress, which keeps the target boundary at the original position.A positive pressure region can still be seen in the middle of the target with a reduced maximum value of about 0.4 GPa.The radial compressive wave generated by the beam impact is immediately followed by a tensile wave generated at the axis, which also propagates radially.Once the beam impact stops, the tensile wave becomes important and unloads the sample, which at the beginning was held under compression due to mass inertia.Figure 8(f) shows that at t ¼ 20 μs the majority of the target is filled with negative pressure, implying that the inner and outer negative pressure areas have combined.The lowest negative pressure remains about −0.5 GPa during the simulated period.

Target density
Modified density maps from Autodyn were used to carry out new FLUKA simulations during the iterations.10(a), there is a compressed region with a density of about 9.0 g=cm 3 at a radial location between r ¼ 0.6 and 0.7 cm.This agrees with Fig. 8(a), which shows that the high-pressure stress wave has just arrived at this area at t ¼ 1.8 μs.The compressed area disappears in Figs.10(b)-10(d), because the high pressure has spread out of the radius of r ¼ 1.0 cm (maximum plotted radius in Fig. 10) at and after t ¼ 3.8 μs.In Fig. 10 (f), we clearly see a region at r ¼ 0-0.25 cm and z ¼ 0-75 cm where the density is close to zero, implying that a slim cavity has been generated in the irradiated target at t ¼ 20 μs.The phenomenon is in accordance with the experimental result.In Fig. 11, curves 1-4 show the density distributions along the target axis from t ¼ 1.8 to 7.8 μs.The minimum density in the most impacted region decreases from 6.0 to 0.6 g=cm 3 with time, while the depletion front moves deeper into the target because of the hydrodynamic tunneling.Corresponding to Figs. 10(e) and 10(f), curves 5 and 6 indicate that the density remains unaltered at axial positions of z > 75 cm after the irradiation, while the density before z ¼ 75 cm continues to reduce until a steady state is reached.

IV. RESULT COMPARISON
A comparison between the above simulations and those performed with FLUKA and BIG2 [47] has been made, in terms of the accumulated specific energy, the evolutions of the target temperature, pressure, and density.Results just after the delivery of the beam (at t ¼ 7.8 μs) are shown in Fig. 12. Differences of the order of 10% can be seen, which are acceptable, considering that different FLUKA scoring (or binning), material constitutive model, and mesh have been used for the present study in comparison with Ref. [47].As listed in Table I, different from BIG2, the tabular SESAME EOS and the empirical Johnson-Cook model have been used in Autodyn.Meanwhile, the Lagrangian mesh has been adopted in Autodyn, taking advantage of its computational efficiency.Note that, for Autodyn, either analytical or tabular EOS, and other material constitutive models and meshes can be applied easily on a case-by-case basis, according to beam parameters and materials.
Target damage caused by the beam heating can be evaluated by the range of target melting.Table II shows a comparison of the molten length in targets 2 (108 bunches) and 3 (144 bunches) between simulations and measurements [46,48].For simulations, the melting ranges which appeared as the flat regions in Fig. 7 are listed.In addition, the results of the static approximation (linear scaling of the initial FLUKA data without considering the ) Target axis (cm) density change) are listed to show the importance of the hydrodynamic tunneling effect in the studied beam-impact case.It can be concluded that our results agree well with the measurements and the simulations performed using FLUKA and BIG2, with a difference of the order of 10%.A significant lengthening of the damage distance can be seen when comparing these results with the results from a static simulation.This shows that it is necessary to make an iterative coupling simulation when the hydrodynamic tunneling is expected to be significant, e.g., in case the maximum reduction of density exceeds 10%.The parameters (temperature, pressure, and density) reached in the simulations show that different phase states of copper have been generated including liquid and gas, and the most heated part has been converted into high-energy-density (HED) matter.
For machine protection, a conservative design is required.Therefore, numerical errors in the simulations are analyzed.We consider that the errors mainly arise from the FLUKA statistical errors (<5%, simulating 150 000 primaries for each run of the 12 iterations), a finite iteration step (defined by a density drop of 10%-15%), errors of the material constitutive models (of the order of 10%), and the simulation accuracy of the hydrocodes (error of the order of 5%).On the other hand, there are also deviations between   the simulated beam parameters and the actual values, due to the errors from the beam diagnostics.The main concern is the beam size with an error of the order of 10%.

V. CONCLUSIONS AND FUTURE WORKS
Hydrodynamic tunneling is an important process to be considered when a number of high-energy bunches are lost at the same place.It has been investigated already for SSC, and in detail during the past 15 years by coupling FLUKA and BIG2.Here we present a new approach to the simulation of hydrodynamic tunneling, where the energy deposition code FLUKA and the commercial code ANSYS-Autodyn are run iteratively.
One critical step is the modelling of the target in FLUKA that has also been improved.The modified density distribution from Autodyn must be used by FLUKA to update the energy distribution.This has been simplified by predefining a certain number of density levels for all iterations instead of defining density levels in each iteration as before.
The new approach was benchmarked against the beam experiment performed at the HiRadMat facility where the SPS beams were used to irradiate solid copper cylinders.The comparison between our simulations and measurements as well as with previous calculations showed good agreement in terms of the melting depth of the target after the irradiation.In order to compare the two hydrodynamic codes, Autodyn and BIG2, in our simulations, the small gaps between neighboring cylinders were neglected, using the same simplification as in the previous study coupling FLUKA and BIG2.The segmentation in the experimental target would surely affect the pressure wave propagation.The simulated pressure contours presented in Fig. 8 would have reflections at these interfaces.Since the temperature and pressure gradients are much stronger in the radial direction than in the axial direction, the influence of the target simplification on the density change and tunneling depth is considered very limited.It is worth mentioning that the simulations also indicated that a significant part of the target was converted into HED matter, which means that the SPS beam can be used to study HED physics via the HiRadMat facility.
For the long term, the option of performing hydrodynamic-tunneling simulations using codes available at CERN is of prime importance.The numerical method presented in this paper is an alternative solution to such simulations.It can be applied to many other practical case studies relating to the machine protection of high-energy colliders, e.g., beam impact on a graphite dump block during the failure of the dilution kickers at the LHC and the HL-LHC.Potential applications in the FCC (design of a possible water dump), the HE-LHC, and other accelerators with high-energy and high-intensity beams are foreseen.

FIG. 1 .
FIG. 1. Target 3 after irradiation, front face (left) and back face (right) of (a) the first cylinder, (b) the fifth cylinder, and (c) the ninth cylinder.Each cylinder is 10 cm long with a radius of 4 cm.

14 FIG. 5 .
FIG.5.Accumulated specific energy along the target axis at t ¼ 1.8, 3.8, 5.8, and 7.8 μs.Note that the total specific energy at t ¼ 7.8 μs in the static approximation is also shown in curve 5.The energy lines implying evaporation and melting are shown in the figure.

FIG. 10 .
FIG.10.Two-dimensional distributions of density in the r-z plane calculated by Autodyn (a) at t ¼ 1.8 μs, after 36 bunches have been delivered, (b) at t ¼ 3.8 μs, after 72 bunches, (c) at t ¼ 5.8 μs, after 108 bunches, (d) at t ¼ 7.8 μs, after 144 bunches, (e) at t ¼ 10 μs, 2.2 μs after the delivery of the beam, and (f) at t ¼ 20 μs, 12.2 μs after the delivery of the beam.Note that a different color scale is used in (a) to show the compressed region.
Figures 10 and 11  plot the density distributions at six different points in time, namely, 1.8, 3.8, 5.8, 7.8, 10, and 20 μs.Figures 10(a)-10(d) show that during irradiation the low-density region extends continuously in the radial as well as the longitudinal directions with a reducing minimum value of density, while Figs.10(e) and 10(f) indicate that the extension occurs mainly radially after the beam heating.As shown in Fig.

FIG. 12 .
FIG. 12.Comparison between present simulations (FLUKA and Autodyn) and previous simulations (FLUKA and BIG2), in terms of (a) accumulated specific energy, (b) temperature, (c) pressure, and (d) density along the target axis, just after the delivery of the entire beam of 144 bunches (at t ¼ 7.8 μs).

TABLE II .
Comparison of the molten depth between simulations and measurements.

TABLE I .
Material constitutive models and meshes employed by the simulations.