The evolution of the temperature field during cavity collapse in liquid nitromethane. Part I: Inert case

We study the effect of cavity collapse in non-ideal explosives as a means of controlling their sensitivity. The aim is to understand the origin of localised temperature peaks (hot spots) which play a key role at the early stages of ignition. Thus we perform 2D and 3D numerical simulations of shock induced gas-cavity collapse in nitromethane. Ignition is the result of a complex interplay between fluid dynamics and exothermic chemical reaction. To understand the relative contribution between these two processes we consider in this first part of the work the evolution of the physical system in the absence of chemical reactions. We employ a multi-phase mathematical formulation which accounts for the large density difference across the gas-liquid interface without generating spurious temperature peaks. The mathematical and physical models are validated against experimental, analytic and numerical data. Previous studies identified the impact of the upwind side of the cavity wall to the downwind one as the main reason for the generation of a hot-spot outside of the cavity; this is also observed in this work. However, it is apparent that the topology of the temperature field is more complex than previously thought and additional hot spots locations exist, arising from the generation of Mach stems rather than jet impact. To explain the generation mechanisms and topology of the hot spots we follow the complex wave patterns generated and identify the temperature elevation or reduction generated by each wave. This allows to track each hot spot back to its origins. We show that the highest hot spot temperatures can be more than twice the post-incident shock temperature of the neat material and can thus lead to ignition. By comparing the maximum temperature observed in the domain in 2D and 3D simulations we show that 3D calculations are necessary to avoid belated ignition times in reactive scenarios.


Introduction
Many explosives used in mining are insensitive and do not detonate under typical shock loading conditions, to ensure their safe transportation to mining sites. Before being detonated in blastholes these explosives have to be sensitised by including artificial impurities such as gas cavities or glass micro-balloons.
The mechanism by which the sensitivity of the explosives is increased is a combination of hydrodynamic, thermodynamic and chemical effects when these cavities collapse under the influence of an incident shock wave. Upon collapse, regions of locally high pressure Email addresses: lm355@cam.ac.uk (L.Michael), nn10005@cam.ac.uk (N. Nikiforakis) and high temperature are generated, which are often referred to as hot spots. If these regions are large enough and maintain high temperatures for a sufficient amount of time they can lead to the local ignition of the material and subsequently thermal runaway. In such cases, the time to ignition is observed to be considerably shorter than the time needed for the ignition of the neat explosive, thus the cavities are considered to sensitise the material.
To this end, experimental and numerical studies have been performed, aiming to better understand the physical processes behind the shock-induced cavity collapse and the sensitisation mechanism, which are imperative to the optimisation of the performance of non-ideal, inhomogeneous explosives. Studying inhomogeneous explosives at the macro-scale, for example in the form of a rate stick simulation or experiment, provides insight to the overall effect the inhomogeneities have on the behaviour of the bulk of the explosive. However, the length-scales and time-scales governing the collapse of the cavities and the initiation and transition to detonation of the explosive bulk are massively disparate. As a result, it is very difficult for a macro-scale simulation or an experiment to simultaneously and adequately capture both processes.
Hence, a lot of focus has turned to understanding the meso-scale physics by means of resolved simulations or experiments of cavity collapse to thereafter upscale to bulk. To reduce the complexity of the process induced by the chemical reactions, several researchers concentrated on the shock-induced cavity collapse phenomenon in isolation, as a purely hydrodynamical process, emphasising on the mechanical and the induced thermodynamic effects.
For example, extensive experimental studies have been performed to identify the process governing the cavity collapse and determine the mechanical effects behind hot spot generation. A series of experiments of air cavities (diameter 5 − 8 mm) collapsing under the influence of shocks of various strengths (0.3 − 5 GPa) in gelatinous materials have been undertaken by Field and co-workers [1,2]. Bourne and Field [2] studied the collapse process of isolated cavities of various shapes (circular, triangular, semi-circular, elliptical) in an inert medium. Dear and Field [1] considered multi-cavity systems collapsing in an inert material and observed phenomena like jet deviation, asymmetric collapse and double jetting in columns of cavities and shielding effects and pressure amplification in rows of cavities. More recently, Swantek and Austin [3,4] studied collapsing cavity arrays under stress loading and analysed the velocity fields around them. Similar studies were done by numerical means, looking at the shock-cavity interaction process in various configurations. A lot of numerical studies include collapse in gaseous media or focus on pressure fields alone. Here, we discuss selectively only a small number of papers that involve non-gaseous media and concentrate on the ones demonstrating results on temperature fields, as these are directly relevant to this work. Mader [5] was one of the first to simulate the collapse of a cavity in nitromethane and calculate the nitromethane temperature field around the cavity, to compare against the postshock, neat material temperature. However, this work was limited by resolution, as the computers at the time were not yet capable of handling the heavy computation needed for adequately resolving such processes. Bourne and Milne [6,7] modelled the temperature within the collapsing cavity, for circular, triangular and hemispherical isolated cavities and cavities in arrays. Ball et al. [8] and Hawnker and Ventikos [9] simulated a single cavity collapse in water and computed the water temperature field, while Lauer et al. [10], Michael et al. [11] and Betney et al. [12] simulated the collapse of arrays in water and focused on the effect of the collapse on the pressure field and subsequent pressure amplification. Similarly, the collapse of circular and ellipsoidal cavities collapsing in an inert, hydrodynamically-modelled solid material was studied by Ozlem et al. [13] who also focused on the pressure fields recovery and Kapila et al. [14] who presented temperature fields in an exemplary explosive. Limited results on the shock-cavity collapse in nitromethane were presented by Michael and Nikiforakis [15], in HNS by Kittel and Yarrington [16] and in HMX by Menikoff [17]. In an elastoplastic framework, Tran and Udaykumar [18,19] and Kapahi and Udaykumar [20] studied the response of inert and reactive HMX with micron-sized cavities, while Rai et al. [21,22] considered the resolution required for reactive cavity collapse simulations and the sensitivity behaviour of elongated cavities in HMX.
The studies mentioned above have contributed to our understanding of the cavity-collapse mechanisms. However, there are still challenges in performing a complete, physical simulation of the initiation of a condensed phase explosive by shock-induced cavity collapse. An important step forward would be to address and rectify the current simulation challenges and work towards a complete simulation of the phenomenon. Such challenges include the non-trivial use of complex equations of state for describing the materials involved in the simulation, retaining at least 1000:1 density difference across the cavity boundary, maintaining oscillation free interfaces (in terms of pressures, velocities and temperatures), obtaining realistic temperature fields in the explosive matrix and heavy computations required for well-resolved, three-dimensional simulations.
A major challenge worthy of special attention is the accurate and oscillation-free recovery of temperature fields in the explosive matrix. This is of critical importance as the ignition process of the energetic material is a temperature-driven effect, thus the accurate prediction of ignition relies on physically meaningful temperatures. An essential prerequisite for recovering accurate temperature fields is the characterisation of the explosive material by means of a realistic equation of state, for example in Mie-Grüneisen rather than in polytropic form, which increases the complexity of the underlying governing equations and hence the difficulty of their numerical integration. Moreover, the large density differ-ence across the material interface (at least 1000:1 for a gas bubble in a liquid) generates additional problems, compared to a uniform material simulation. Assuming that a non-oscillatory, interface capturing numerical scheme is employed, thin but purely numerical mixture zones are generated at the material boundaries, which can give rise to velocity, pressure and temperature oscillations [23,24,25]. Even if the velocity and pressure oscillations are eliminated (e.g. by an energy correction step [23,25] or by augmenting the Euler system of equations with equations evolving parameters of the equation of state [26,27,28,29]), the temperature oscillations may persist, as demonstrated in previous work by the authors (Fig. 4 in [25]). This could be attributed to the fact that the specific heat of the artificial mixture cannot be physically derived and hence must be constructed mathematically. Independent of their origin, the temperature oscillations would be destructive in any ignition studies, as they would spuriously activate the temperature-dependent reaction rate, thus generating false ignition. It is, thus, important to ensure that spurious numerical artefacts are either removed or do not arise in the first place.
Most work found in the literature that present temperature fields is either focusing on the temperature inside the cavity only, done in two-dimensions or use idealised equations of state. In this work, we overcome the difficulties in numerically simulating the cavity collapse mentioned above (including the spurious temperature oscillations) by exercising a mathematical model proposed by the authors in previous work [25], allowing for a complete shock-cavity interaction simulation. We simulate the three-dimensional collapse of isolated air cavities in nitromethane, using the Cochran-Chan equation of state for the explosive and demonstrating the validity of this equation of state in the regime of pressures found in the simulations by comparing to experiment. The temperature fields recovered are also free from spurious oscillations and the 1000:1 density difference between the cavity and ambient material is captured directly. In the first instance we consider the induced temperature field in the explosive in the absence of chemical reactions so that hydrodynamical effects are elucidated. However, the oscillation-free temperature fields will allow for the use of a temperature-dependent reaction rate law in future work. In this work, we move beyond solely reporting potential hot spot loci and identifying high temperature peaks; we look in depth how the mechanical effects (e.g. the generation and propagation of waves) lead to the thermodynamical effects (e.g. elevation of temperatures) that generate the localised hot spots. Each wave generated in the complex cavity col-lapse scenario is labelled and its effect on the temperature field is thereafter determined. In experiments such detail of study is limited by instrumentation, the locus of luminescence downstream the cavity is noted as the main ignition site and the jet impact as the principal ignition mechanism (e.g. [2,30,31]). In numerical studies, where more detailed wave patterns can be obtained, the shock waves produced upon jet impact are considered as a single wave moving outwards of the cavity remains (e.g. [6]). We demonstrate, however, that two distinct shock wave components can be identified upon collapse; one travelling downstream and one travelling upstream. It is also demonstrated that it is the upstream travelling wave that leads to the highest temperature locus in collapse scenarios of the incident shock pressure considered here. We identify the regions with temperatures that are higher than the post-shock temperature and are candidates for ignition loci in reactive simulations. In this process, we report a new hot spot location, in the Mach stem generated at late stages of the collapse. This hot spot is omitted in previous studies (e.g. [6]) but is considered here in detail as its effect could be similar to that of the principal hot spot along the cavity centreline. By comparing to the equivalent two-dimensional cavity collapse we demonstrate also the importance of performing this study in a three-dimensional framework as the two-dimensional simulations are shown to underpredict the temperature fields in some stages and overpredict them in others, something that could lead to no or false ignition sites.
In Part II of this work, that will follow in the future, the reactions in nitromethane will be activated and the ignition of the material can then be traced back to the elevation of the temperature fields and accountable waves described here.
The rest of the paper is outlined as follows; the next section presents the underlying formulation used in this work in terms of governing PDEs, the equations of state that close the system and the method for temperature recovery. A section on validation follows, where the mathematical system and numerical methods are tested against known solutions. The mathematical and physical models are also validated against theoretical and experimental temperatures of shocked nitromethane. In the results section that follows we present the threedimensional collapse of an air cavity in inert, liquid nitromethane and follow the events leading to the generation of locally high temperatures, study the temperature field topology and generated hot spots, the evolution of pressure and temperature fields on constant latitude lines and the maximum temperatures in the 2D and 3D simulations. Finally a section with the conclusions of this work is presented.

Mathematical and physical model
A formulation which rectifies the issues commonly presented in shock-bubble simulations discussed in the previous section was proposed by Michael and Nikiforakis [25] (MiNi16). This formulation considers the cavity as an inert phase (phase 1) and the surrounding material as a reacting phase (phase 2). This work focuses on the collapse of cavities in an inert medium, hence we neglect the reaction source terms and the products of reaction and as a result we use a reduced form of the MiNi16 formulation. The reactive terms are reinstated in the follow-up, Part II of this paper. Considering the gas inside the cavity as phase 1 and the liquid nitromethane around the cavity as phase 2, as shown in Fig. 1, the governing equations for this system take the form ∂z 1 ρ 1 ∂t + ∇ · (z 1 ρ 1 u) = 0, where for i = 1, 2, ρ i are the densities for the air and nitromethane, z i are their corresponding volume fractions (z 1 + z 2 = 1), ρ is the total density given by ρ = z 1 ρ 1 + z 2 ρ 2 , u is the velocity vector and p is the total pressure. The total specific energy is given by E = 1 2 u 2 + e, where e is the total specific internal energy. The mixture rule for the total internal energy is given by ρe = ρ 1 z 1 e 1 +ρ 2 z 2 e 2 , where e i for i = 1, 2 are the specific internal energies for the air and nitromethane, given by their corresponding equations of state. A mixture rule for ξ = 1 γ−1 , where γ is the total adiabatic index, is also required and in this case is given by ξ = z 1 ξ 1 + z 2 ξ 2 . The sound speed for the total mixture is given by where y i is the mass fraction of phase i, given by y i = ρ i z i ρ , for i = 1, 2.

Equations of state
To close the system, the Cochran-Chan equation of state [32] is employed to describe the liquid nitromethane. This is an equation of state of Mie-Grüneisen form and is given by with reference pressure given by reference energy given by and Grüneisen coefficient Γ(ρ) = Γ 0 . gas inside the cavity is modelled by the ideal gas equation of state, which is of Mie-Grüneisen form as well, with p ref = 0 and e ref = 0 . The parameters for the equations of state of the two materials are given in Table 1.

Recovery of temperature
The multi-phase nature of the model allows for separate temperature fields to be computed for each material as As a result, the nitromethane temperature (T NM = T 2 ) is computed explicitly from the equation of state and can be used directly in the reaction rate law.
Computing the temperature of a general condensed phase explosive (T CF = T 2 ) can involve completing the equation of state starting from the basic thermodynamic ΓT v and by integrating to obtain a reference temperature (T ref CF ) such that:  When the reference curve is an isentrope, dS /dv = 0 and hence we can simply compute T ref CF = T 0 ( ρ ρ 0 ) Γ . When the reference curve is a Hugoniot curve, the basic thermodynamic law cannot be integrated directly and often the Walsh Christian technique and numerical odeintegration techniques are used to compute the reference Hugoniot temperature.
For this work, substitution of the parameters of the equation of state for nitromethane and imposing an initial temperature of 298K for ρ = ρ 0 gives T 0 = 0. This is in line with other work using the Cochran-Chan equation of state, where T 0 takes zero or very small values [33,34,35,36]. The form (6) gives temperatures that match experiments as demonstrated in the validation section, but for other materials or other equations of state (e.g. shock Mie-Grüneisen) care should be taken as a different reference curve (as per Eq. 7) would be necessary.
The ideal gas equation of state results in the overheating of the gas inside the cavity since high pressures are reached within the cavity during the collapse process. However, at these timescales heat transfer would not physically take place and thus the cavity temperature would not affect the ambient nitromethane temperature. Since temperatures inside the cavity are not of interest for this work, they are not presented henceforth.
Note that the two-phase nature of the model allows for large (>1000:1) density gradients to be sustained across material boundaries and both the density and temperature fields are maintained oscillation-free.

Validation
System (1) is integrated numerically using a high resolution shock-capturing numerical scheme, namely the MUSCL-Hancock finite volume method with an underlying HLLC Riemann solver. Hierarchical, structured, adaptive mesh refinement (AMR) is used to dynamically increase the resolution locally [37].
The aim of this section is to validate the resulting code using test problems with known solutions and to assess the suitability of the physical models to predict realistic temperature fields.

Air pocket collapsing in water
The shock-induced collapse of an air bubble in water is used as a validation example for the mathematical and numerical methodologies in this work, due to its close relation to the application we are considering. We follow the set up considered, among others, by Terashima and Tryggrason [38] and Xu and Liu [39]. As the solution to this problem is self-similar and to compare directly with results from both [38] and [39], we consider this test in non-dimensional variables. Denoting the dimensional variables by over-bars, the non-dimensional variables in consideration are given as: where ρ 0 = 1000 kg m −3 , u 0 = 10 m s −1 and l 0 = 1 mm for [38] and l 0 = 1 m for [39].
A domain of size (x, y) ∈ [−6, 6] × [0, 6] is considered with base resolution 250×125 and three levels of AMR, doubling the resolution with each refinement level, leading to an effective resolution of 2000 × 1000 cells and an air bubble of initial diameter D = 6. The initial conditions for this test are given in non-dimensional form in Table 2.
The air is modelled as an ideal gas and the water as a stiffened gas. Both equations of state can be written in the Mie-Grüneisen form (3) In Fig. 3a the evolution of the bubble boundary under the effect of the incident shock wave is shown at time instances t = 1, 1.5, 2, 2.5, 3, 3.5 and 4. The top half of the image (y > 0) shows results from this work and the  (a) The evolution of the bubble boundary due to its interaction with the incident shock wave, at instances t = 1, 1.5, 2, 2.5, 3, 3.5, 4. The top half (y > 0) shows results from this work and the bottom half (y < 0) shows results from [39].
(b) Comparison of the evolved, non-dimensional height (H/D) and width (W/D) between this work and results from [38] and [39]. bottom half (y < 0) shows results from [39]. The two sets of results match reasonably well and small differences are attributed to the fact that a diffused interface method (shock-capturing) is used in this work whereas a shock-tracking method is used in [39]. In Fig. 3b we present the evolved, non-dimensional, normalised with respect to the bubble diameter width and height of the bubble. We compare the width and height obtained from this work to the results of [38] and [39].

Shocked nitromethane temperature
To validate the physical models used in this work we compare the temperature of the shocked, neat, nonreacting nitromethane computed from our simulations with experimentally measured values, as well as with the temperature calculated from analytic theory and molecular simulation, for various post-shock pressures; this compilation is presented in Fig. 4. The set of ex- perimental data most commonly used as a set of reference values is given by Lysne and Hardesty [40] and have been deduced numerically from experimentally determined parameters. Other experiments measuring post-shock temperatures include the work by Winey et al. [41] and Hervouët et al. [42], as well as work by Delpuech and Menil [43] and Dufort [44]. Analytic work includes Winey et al. [41] and Cowperthwaite and Shaw [45], while results from microscopic simulations are given by Soulard [46], Jones [47], Liu et al. [48], Hervouët et al. [49] and Desbiens et al. [50].
The solid circles ( ) in Fig. 4 present the temperature of neat nitromethane predicted from our simulations, calculated by expression (6). It is important to note the wide spread of the results from different sources and the fact that there is no indication of experimental error. Nevertheless, the post-shock temperatures predicted from our simulations are well within the bounds of this compilation and they are closest to the ones by Lysne and Hardesty [40]. This ensures the validity of the temperature field for the explosive material, before the inclusion of cavities.

The collapse of a single cavity in liquid nonreacting nitromethane
In this section, an isolated air-filled cavity of radius 0.08 mm 1 collapsing in non-reactive liquid nitromethane due to a 10.98 GPa incident shock wave (ISW) is considered. The air inside the cavity is modelled by the ideal gas equation of state and the nitromethane is modelled by the Cochran-Chan equation of state as described in Sec. 2.1, with parameters as listed in Table 1. As the ignition and thermal runaway in an explosive are attributed to the complex interaction between non-linear gas-dynamics and chemistry, it is intuitive to consider in the first instance the induced temperature field in the explosive in the absence of chemical reactions. This will allow the purely gasdynamical effects to be elucidated. The simulation is performed in three-dimensions, with effective grid size dx = dy = dz = 0.3125 µm, selected after the convergence study presented in Appendix Appendix A. The initial conditions for the simulation are given in Table  3, corresponding to the sketch of Fig. 1.
The aim of this section is to follow in detail the generation of complex wave pattern during the collapse process and identify the effect of each wave on the nitromethane temperature field. This will lead to identifying which waves are generating the high-temperature regions in the explosive and the magnitude of the temperature field that could account for the local ignition in a reactive simulation.
4.1. Early stages of the three-dimensional cavity collapse In order to study the complex wave pattern emanating from the cavity collapse process, xz-slices through the centre of the cavity are taken. The xy-plane plots show the same results and are thus omitted. As the cavity collapse is symmetric about the x-centreline, the lower part of the cavity is omitted. However, to give the threedimensional perspective of the setup, we first present in Fig. 5 a couple of indicative three-dimensional plots of the shock-cavity collapse in terms of the nitromethane temperature field. Fig. 6 illustrates the mock-schlieren plots on this plane, for the early stages of the collapse 1 typical size of large cavities used in mining applications [Pa] x < 100µm shocked nitromethane 2. In Fig. 6a, the transmitted air shock (S 2 ) and reflected rarefaction wave (R 1 ) generated by the interaction of the ISW (or S 1 ) with the cavity interface are depicted. In Figs. 6b-c, the non-uniform pressure and velocity fields around the cavity are seen to lead to the formation of the intruding jet.
Before the jet is fully formed, the air shock S 2 traverses a small part of the cavity and its upper part moves downwards due to the high pressure around the upper part of the cavity. This results in self-focusing of the shock, leading to a separation of the upper and lower parts of the wave, labelled S 2A and S 2B in Fig. 6c, and a kink where the two parts are connected. Two shocks, S 3 and S 4 , emanate from the kink (Fig. 6d). The air shock travels faster than the ISW and the upper part of S 2 forms an inward moving shock in the air, S 5 , and an outward moving shock in the nitromethane, S 6 ( Fig.  6d). Waves S 4 and S 5 interact inside the cavity, while a secondary kink and two more shock waves are generated along S 5 , as seen in Figs. 6e-f. In Fig. 6f, the lower part of the air shock (S 2B ) is undergoing a transmission/reflection process 2 upon reaching the front cavity face, while the upper part (S 2A ) has already completed this process. The most interesting waves in this case are the ones emanating from the interaction of S 2B with the downstream cavity wall. These are the shock labelled S 7 , travelling downstream in the ambient nitromethane and the shock wave labelled S 8 travelling upstream in the cavity (Fig. 6g).
The last row of Fig. 6 presents the stages close to the collapse time of the cavity. Fig. 6g in particular, illustrates the waves in the system just before the jet reaches the downstream cavity wall. Upon reaching this wall 3 2 When a shock wave travelling in a low impedance material reaches a material boundary ahead of which a material of higher impedance is encountered, two new waves are generated; a downstream-travelling shock wave in the high impedance material and an upstream-travelling shock wave in the lower impedance material. In short, we call this a transmission/reflection process. 3 we call this the time of collapse two new shock waves are generated in the nitromethane: a downstream-travelling shock wave labelled as S 9 and an upstream-travelling shock wave labelled as S 10 , both shown in Fig. 6h. Other numerical studies do not distinguish between these two parts of the waves and present it as a single 'hammer' shock (e.g. [6]). However distinguishing between the two will be proven important in the following sections when the temperature field is studied. In the 3D-space, the waves S 7 , S 9 and S 10 can be considered as spherical caps with different centres, radii and heights. Upon collapse, the cavity forms a toroid of trapped air. On each two-dimensional slice this phenomenon is seen as splitting of the cavity in two lobes. Any gas that was in the region between the jet and the downstream cavity wall before the collapse is pushed into the toroid (or lobes), where it is compressed further as time passes. This violent compression results in another shock wave labelled as S 11 in Fig. 6i travelling upwards in the lobe [51] and as this passes over the front cavity wall it generates the shock wave S 11B travelling in nitromethane. In the same figure, a liquid infusion zone Z 1 is illustrated. This is a zone created by the liquid penetrating into the gas region due to the high pressure acting on the liquid in this area and obtains a toroidal shape in 3D-space. As this zone is treated as a mixture zone by the mathematical model, its thermodynamic details are open to question (similarly in [51]).

Late stages of the three-dimensional cavity collapse
The wave patterns emerging at the late stages of the collapse process are complex, hence planar mockschlieren plots remain the most convenient means to visualise them (Figs. 7-8). In Fig. 7a, the air shock S 8 is seen to reach the downstream side of the lobe. Its interaction with the material interface results in a shock wave traversing the lobe, S 13 , and a shock travelling in the liquid, S 12 , (see Figs. 7b-c) obtaining a spherical cap shape in the 3D space. The shock wave S 11 travels upwards in the lobe and interacts obliquely with wave S 13 , resulting in the appearance of a new pair of weak shock waves, labelled S 13A in Fig. 7b. As S 11 contin- ues its upward travel generating S 11B it also interacts with the ISW (S 1 ) and this complex process leads to the formation of S 14 . While moving upwards S 11 also interacts with S 13 and other weak waves that are present in the lobe forming a front labelled S 15 . Upon reaching the end of the lobe, it undergoes a transition/reflection process again resulting in a spherical-cap wave (S 16 ) in the liquid and a different spherical-cap wave, S 17 , (Fig. 7e) in the lobe. The transmitted wave S 16 combines quickly with S 14 and the two move together outwards from the lobe, thus labelled as S 14,16 in Fig. 7e. This wave then combines with S 12 , as seen in Fig. 7f and is thus labelled as S 12,14,16 . Meanwhile, several reflection/transmission processes take place at the lobenitromethane interface, which lead to the lobe closing up due to its rapid compression. As a result, such processes are difficult to track from this point onward (e.g. wave S 18 in Fig. 7e). In general, the waves transmitted from the lobe continue to expand outwards and the ones traversing the lobe compress the cavity gas further.
The spherical-cap nature of the waves S 12,14,16 and S 9 , together with the curvature of the ISW, lead to the superposition of these three waves with S 7 and S 11B in zone Z 2 , as can be seen in Fig. 7f.   Fig. 8a illustrates waves S 19 and S 20 , which have been generated during late-time transmission/reflection processes at the lobe boundary. The downstreamtravelling wave generated upon the collapse (S 9 ) interacts with the ISW and changes shape, thus labelled now as S 9A and S 9B (Fig. 8b). S 19 interacts with S 9A and it is similarly split in two parts, S 19A and S 19B , as illustrated in Fig. 8b. New waves, S 21 and S 22 , are generated by more transmission/reflection processes at the lobe boundary and S 22 is also split into two parts, S 22A and S 22B , due to its interaction with the upstreamtravelling wave generated upon cavity collapse, namely S 10 (Fig. 8c). Waves S 23A and S 23B are generated by the interaction of the lower part of the ISW (S 1B ) with the lobe. At this point, the ISW has already been split in two parts (S 1A and S 1B ) by its interaction with S 9 (Figs. 8c,d). S 23C emerges from the lower side of the lobe and merges with S 23A and S 23B . In Fig. 8f two new waves are seen to emanate from the lobe, labelled S 24 and S 25 . As the ISW continues to interacts with S 9A and S 9B , a Mach stem region Z 3 is formed (Fig. 8e). This new region, often omitted in the literature, will be seen in the following sections to play an important role in the potential ignition of nitromethane.
In Fig. 8f wave S 26 is the equivalent of the S 12,14,16 wave emanating from the lower lobe of the cavity (which is not shown here). Similarly, wave S 27 is the equivalent of waves S 22A and S 22B generated at the lower lobe.

Temperature distribution during the early stages of the collapse
Having studied in detail the generation of the several waves in the collapse process, we will look now into the effect they have on the nitromethane temperature field. When nitromethane is shocked to 10.98 GPa it attains a  temperature of 1263 K, as illustrated in Fig. 9a. Note that since we are only interested in the nitromethane temperature and to ease the interpretation of the results, the cavity region is masked out in the temperature plots of this work. In the early stages of the collapse, the rarefaction wave R 1 generated by the interaction of the ISW with the cavity boundary (see Fig. 6a), lowers the temperature of nitromethane upstream of the cavity. Specifically, the temperature is distributed according to the distribution of density in the rarefaction wave. The air shock transmitted into the liquid, S 7 (Fig. 6g), heats the traversed nitromethane region slightly, raising its temperature from 298 K to about 310 K seen in the inlet of Fig. 9b.   Fig. 10 shows the temperature field in nitromethane and the isotherms of over-temperature, during stages close to the collapse time. Over-temperature is defined as the local nitromethane temperature minus the nitromethane post-shock temperature: T NM −1263 K in this case. This definition is in fact a play on the term 'over-pressure' often used in the literature (e.g. [14]). Fig. 10a illustrates the temperature distribution just before collapse. The highest temperature at this stage (1263K) is found behind the ISW. The over-temperature contours are not plotted at this stage as they have negative values. As soon as the cavity collapses, the new shock waves generated (waves S 9 and S 10 in Fig.  6h) heat the liquid nitromethane above the post-shock temperature, as seen in Fig. 10b. Waves S 9 and S 10 initially generate over-temperatures of the range 0 K-400 K downstream and of the range 0 K-1700 K upstream of the cavity, in regions we will call front hot- spot (FHS) and back hot spot (BHS). In a similar manner we label the waves generating these hot spots; we call the wave S 9 front collapse shock wave (FCSW) as it is the wave generating the FHS and the wave S 10 back collapse shock wave (BCSW) as it is the wave generating the BHS. As the BCSW travels outwards, it traverses a region of nitromethane that has already been shocked by the ISW and as a result, it leads to a generally higher temperature in the BHS than the temperature in the FHS, which is a region that has not already been compressed. As shown in Fig. 10c, over-temperature of the range 0−1897 K is observed upstream of the cavity, in the BHS, while over-temperature of the range 0−723 K is observed downstream of the cavity, in the FHS. The inner contours are the ones with the higher over-temperature value but enclose smaller areas.
To visualise the generation of the high-temperature regions we turn to Fig. 11. On the left side of each image, a colour plot of the density field is seen on a plane through the centre of the cavity. The white contour is the T NM = 1264 K contour, enclosing the regions with positive over-temperatures. On the right side of each image, a three-dimensional plot of the z = 0.5 contour is seen in blue, denoting the cavity boundary. A three-dimensional equivalent of the T NM = 1264 K contour is seen in a the white-orange-black scale. The three-dimensional surface is cut such that it only covers one quadrant to allow better visualisation. At the early stages that are close to the collapse time, the threedimensional format of the FHS and BHS are seen in Fig. 11(top). At this stage, both high temperature regions have a spherical-cap form. The flow going into the jet creates another region with positive but low overtemperature seen by the lower grey part of Fig. 11(top).

Temperature distribution during the late stages of the collapse
The temperature distribution in nitromethane during the late stages of the collapse process is presented using two-dimensional slices in Fig. 12 and Fig. 13 and three-dimensional contours in Fig. 11. With the outward propagation of the FCSW and the BCSW(S 9 and S 10 ), the regions of positive over-temperature grow in volume. The maximum over-temperature, however, decreases. This is seen to happen eventually in all the over-temperature regions during the late stages of the collapse. The maximum over-temperature in Fig. 12a is 1826 K and it is observed in the BHS. The maximum over-temperature in the FHS is 823 K. This is contrary to the usual literature reporting of the highest temperature to be ahead of the collapsed cavity (e.g. [2,6,7]). The transmitted waves S 7 (Fig. 6g) and S 11B (in Fig. 7b) do not significantly increase the nitromethane temperature. They result in an over-temperature of only 50 K-80 K.
In Fig. 12a, the temperature around the lobe is seen to have been affected by the RW. At this point, sample temperatures around the lobe reach values of about 650 K-780 K. In Fig. 12b the shock wave S 12 exits the lobe, gets linked with S 14,16 and heats the nitromethane as it travels outwards of the lobe. At this point, sample temperatures in the region traversed by the wave S 12,14,16 are in the range of 679 K−745 K.
Meanwhile, region Z 2 (Fig. 7f) is formed. This zone does not play a significant role in increasing the temperature (Fig. 12c), yet, as nitromethane is heated only up to about 760 K in this region. The maximum overtemperature is 1570 K and is found in the BHS. The maximum over-temperature in the FHS reaches 750 K. Figure 11: Demonstration of over-temperature regions in 3D space. On the left side of each image is a colour plot of the total density field (ρ) on a plane through the centre of the cavity and projected backwards for illustration purposes. On each 2D plot, the black contour represents the cavity boundary (z = 0.5) and the white contours (T N M = 1264) enclose the regions of (positive) over-temperature. On the right side of each image the three-dimensional cavity boundary is seen as a blue surface (contour z = 0.5). The over-temperature regions are demonstrated only in one quadrant to aid visualisation, using an 3D isosurface of T N M = 1264. Note that the collapse process is shown here to move from bottom to top rather than left to right. The images here can be considered as rotated by 90-degrees counter-clockwise compared to Figs. 6-10 and 12-13. This is done solely for illustration purposes.
In Fig. 12c, the growth of wave S 12,14,16 generates a new over-temperature region upstream of the cavity, part of which is shown to be bounded by the green overtemperature contour. An isotherm of over-temperature 0 K enclosing a region of maximum temperature 1345 K is seen in the same figure.
As more shock waves emanate from the lobe, new regions of over-temperature are formed around it (see Fig.  12d). As the wave S 12,14,16 is propagated downwards, it leads to the growth of the over-temperature regions, like the one illustrated by the fuchsia contour in Fig. 12d.
Even though only one contour of this region is included in this figure, the nitromethane temperature is increased throughout the whole region traversed by this wave, as inferred by the values of over-temperature in this figure.
Meanwhile, the ISW and the FCSW (S 1 and S 9 ) generate zone Z 3 , which we call the Mach stem hot spot (MSHS). This is a high temperature region which, at this instance, encloses over-temperature contours of values 0 K−230 K, as seen in Fig. 12d. The highest temperature (2773 K), however, is still observed within the over-temperature contours generated by the BCSW, up- stream of the cavity.
However, as zone Z 3 grows, the local nitromethane temperature increases, as seen in Fig. 13a and encloses over-temperatures ranging between 0 K and 1380 K. This is a region that is omitted in previous studies found in the literature (e.g. [6,7]). Its importance is clear as so high temperatures are found in this region that could lead to ignition in reactive simulations.
The maximum over-temperature (1421 K) is still located within the BHS. However, under other initial conditions (e.g. weaker ISW) not presented in this work, the authors have observed the maximum temperature to move in zone Z 3 at this stage. The highest overtemperature in the FHS is 640 K. Fig. 11(bottom left) shows in 3D the high temperature regions of FHS, BHS and the very early stages of the generation of the MSHS. Behind the cavity, a large region of positive over-temperatures is seen with relatively low over-temperature values, though, compared to the other three regions.
Waves S 12,14,16 and S 22B now also travel downwards, enlarging the over-temperature regions, as seen in Fig.  13a. As wave S 12,14,16 is superimposed with the BCSW (S 10 ) it generates the over-temperature regions enclosed in dark blue and turquoise, with over-temperature values 361 K and 446 K respectively. The spherical nature of waves S 12,14,16 and S 22B leads to the growth of the overtemperature regions generated in a three-dimensional sense, all around the lobe.
Wave S 26 is generated from the reflection of the wave S 12,14,16 at the bottom boundary (Fig. 13b). As a result, the liquid in the region the wave S 26 traverses is shocked for the third time, leading to a temperature of about 1622 K. Part of S 26 is superimposed with the BCSW, leading to a temperature of about 1748 K. The highest temperature in the BHS is now 2600 K, in the FHS 586 K and in the MSHS 2500 K.
In Fig. 13d, the wave S 22B is also reflected at the bottom boundary, generating wave S 27 , which is then superimposed with the BCSW and waves S 22A and S 22B . This leads to temperatures of the range 1580 K-1650 K. Fig. 11(bottom left) shows in 3D the growth of the four high temperature regions and also the additional increase of temperature in the back cloud of small overtemperature due to the traversing and superposition of the waves emanating from the lobes along the cavity centreline. The highest temperatures are noted in the BHS and the MSHS, as seen in Fig. 14 and are identified as the contours enclosing temperatures above 2000K. These are noted as the most possible regions for observing ignition in reactive simulations. In this figure the three-quarters of the three-dimensional temperature contour is seen.

Evolution on constant latitude lines
It is informative to consider the evolution of the flow field on lines of constant latitude on slices like the ones presented in Figs. 9,10,12,13. We consider the pressure and temperature of nitromethane along the lines y = 201 µm (Figs. 15a,15b) giving insight on the FHS and BHS and y = 290 µm (Figs. 15c,15d) giving insight on the MSHS.
In Fig. 15a (first red) the ISW is seen before interacting with the cavity. The corresponding line in Fig. 15b shows the increase of the temperature in nitromethane due to the shock wave to 1263 K. Upon the interaction of the ISW with the cavity, the pressure behind the cavity is lowered due to the RW generated ( Fig. 15a (blue)), translated into a decrease of nitromethane temperature ( Fig. 15b (blue)). As the jet is formed and vortices are generated at the back of the cavity, the RW and the jet flow are working against each other, the former towards lowering the pressure and temperature and the latter towards raising them. The raising of the pressure and temperature prevails as the RW travels upstream while the jet becomes more profound. This is seen in Figs 15a and 15b (pink-black).
The FCSW and BCSW are seen in Fig. 15a (orange) and correspond to a downstream-travelling and an upstream-travelling thermal front respectively seen in Fig. 15b. These fronts are the ones expected to lead to ignition in reactive simulations. As the temperature rise from the BCSW is higher than that from the FCSW, we expect the BHS to lead to faster reaction than the FHS. In between the two thermal fronts of 15b (orange, grey) the nitromethane temperature is seen to reach 0. This is because this region is still occupied by the air of the collapsed cavity hence nitromethane does not reside there. The superposition of waves emanating from the two lobes along the centreline of the cavity is seen as two new, upstream pressure and temperature peaks in the second green graph of Figs. 15a and 15b.
Figs. 15c and 15d show the evolution of the pressure and nitromethane temperature along the y = 290 µm line, at the same time instances as presented in 15a and 15b. As the RW takes longer to reach this line, the pressure and temperature seem to maintain their post-shock values for longer (Figs. 15c, 15d (first red-blue)). The effect of the RW is seen in the pink and light blue curves. The generation of the Mach stem is seen in the second red and second green lines of Fig. 15c, accompanied with the thermal fronts shown in Fig. 15d. This would be another region where ignition is expected to be observed in reactive simulations.

Maximum temperature in 2D vs 3D
To demonstrate the importance of three-dimensional simulations we perform the equivalent two-dimensional simulation of the configuration presented in Sec. 4 and compare directly the temperature fields obtained from the two set ups. In Fig. 16 we compare the maximum temperature field observed in the computational domain over time in 2D and 3D simulations.
It is observed that in the three-dimensional simulation the maximum temperature is seen upon collapse. This agrees with the results observed in the detailed analysis performed in Sec. 4 where the maximum temperature is observed in the BHS. The maximum temperature is seen to decrease after this peak, which also agrees with the earlier observation that the temperature reduces as the front and back hot spots grow. A new peak in the maximum temperature is observed at late stages and this corresponds to the temperature increase due to superposition of waves along the centreline of the cavity. In the two-dimensional simulation, the collapse of the cavity is observed slightly later than in the three-dimensional scenario, leading to a lower temperature rise than in the three-dimensional case. However, the decrease of this temperature after the initial collapse peak is almost identical between the two cases. Of great significance is the fact that the maximum temperature in the twodimensional simulation is observed in the late stages of the collapse.
The maximum temperature results in the two-and three-dimensional configurations follow the same trend, although the first peak temperature is under-predicted and the second peak temperature is over-predicted in the two-dimensional simulations. This would be of great importance in reactive simulation scenarios as, depending on the details of the reaction rate law, it could give the difference of ignition in 3D and no ignition in 2D (at the first peak) or even false late ignition (at the second peak) in 2D. This plot also highlights the importance of the detailed study of the waves and corresponding temperature fields as presented in Sec. 4, as the maximum temperature can not be entirely accounted for without the detailed wave-pattern analysis. These are the regions identified as the most probable to observe local ignition at. The three-dimensional cavity boundary is seen as a blue surface (contour z = 0.5). On the left side is a colour plot of the total density field (ρ) on a plane through the centre of the cavity and projected backwards for illustration purposes. The black contour represents the cavity boundary (z = 0.5) and the white contours (T N M = 2000K) enclose the regions of high over-temperatures. Note that the collapse process is shown here to move from bottom to top rather than left to right. The image here can be considered as rotated by 90-degrees counter-clockwise compared to Figs. 6-10 and 12-13. This is done solely for illustration purposes.

Conclusions
In this work we perform resolved numerical simulations of cavity collapse in liquid nitromethane using a multi-phase formulation, which can accurately recover temperatures in the vicinity of the cavity. Considerable care is taken regarding the form of equations and numerical algorithm to prevent spurious numerical oscillations in the temperature field. The mathematical formulation and its implementation are validated against results of shock-bubble collapse in water studied extensively in the literature, demonstrating a good match between our work and previous studies. The suitability of the model (evolution equations and equations of state) to capture realistic temperatures in the explosive is demonstrated in a second validation problem. The shock temperature of nitromethane at various shock pressures is compared between our work and a combination of analytic and experimental results showing very good agreement.
Following the validation, we study in detail the shock-induced cavity collapse in inert liquid nitromethane in three dimensions. We follow the interaction of the incident shock wave with the cavity and the generation of waves during the collapse. For ease of interpretation we split the collapse phenomenon in two regimes; early stages and late stages. We follow the events leading to the generation of locally high temperatures and identify the effect each wave has on the nitromethane temperature field, something that would be very difficult to achieve this experimentally. Moreover, numerical studies on this are very limited, as several challenges have to be overcome for this to be possible, as mentioned in the introduction of his work. We also define as over-temperature the difference between the local nitromethane temperature and the post-shock nitromethane temperature. This concept identifies potential regions that could lead to ignition due to cavity collapse that would not ignite in a neat material shocked at by same shock loading. We study the temperature field using plots on two-dimensional slices of the domain, one-dimensional, constant latitude lines on these slices as well as two-and three-dimensional temperature and over-temperature contours.
Commonly in the literature it is suggested that the locus of highest temperature is in front of the cavity, due to the collapse. We demonstrate that upon collapse two shock wave components are generated that lead to two high-temperature loci. The downstream-travelling shock wave (FCSW) leads to the generation of a high temperature region in front of the cavity (FHS) and the upstream-travelling shock wave (BCSW) leads to the generation of a high temperature region behind the cavity (BHS). We identify that, in fact, the highest temperature during the collapse is observed in the BHS behind the cavity. This is attributed to the fact that this region is shocked twice; once by the ISW before the collapse and one by the BCSW after the collapse. In contrary, the FHS region is only shocked once, by the FCSW having similar strength as the BCSW. During the early stages of the collapse, temperatures more than 2.5 times of the post-shock temperature are observed in the BHS. In the FHS temperatures of 1.6 times the post shock temperatures are observed.
A lot of attention is often given to the shock waves emanating from the lobes due to multiple wave and material interface interactions. It is identified here, though that these waves do not contribute much to the temperature elevation, leading to temperature increases of 0-100K.
Another region identified as a high-temperature locus with the potential to sustain a hot spot in reactive simulations is the Mach stem region generated at late stages after the collapse. This region is neglected in the literature. However, we demonstrate that the temperatures observed in the MSHS are higher than the temperatures in the FHS and comparable to the temperatures in the  BHS. In particular, temperatures of more than two times the post-shock temperatures are seen in the MSHS.
A final region with positive over-temperatures is identified far upstream, denoted as rear overtemperature region. Initially this region attains temperatures higher than the post-shock temperature due to the jet flow during the collapse. At late times after the col-lapse, however, this is also the region where the waves emanating from the lobes are superimposed. For the duration of our simulations this superposition leads to temperature increase of ∼450K.
The importance of performing three-dimensional simulations compared to two-dimensional ones is demonstrated in the final section. We study the maximum temperature observed over time in the computational domain and we conclude that the twodimensional scenario leads to could lead to underpredicting the initiation time or predict false late-time ignition.
The identified positive over-temperature regions will be studied in the second part of this work where threedimensional reactive simulations will be used to clarify which of these in fact lead to ignition. Having studied the hydrodynamic effects of the collapse on the temperature field, in Part II we study the additional effect of chemical reactions in the timescales and lengthscales of the configurations in consideration.

Appendix A. Convergence study
To establish the framework for the simulations of this work, grid independence of one dimensional solutions is assessed for a typical test case of an air cavity collapsing in nitromethane. The cavity has radius 0.08 mm and its centre is initially located at x = 0.25 mm, in a domain of size 0.5 mm. The initial conditions are as in Table 3. Nine different resolutions are considered, starting from 200 cells (64 cells across the cavity) and doubling them each time until reaching the finest resolution of 512000 cells (16384 cells across the cavity). In Fig. A.17, the effect of grid resolution on the distribution of the maximum temperature in the domain observed during the collapse is illustrated. For the 200 and 400-cell meshes, the time of the initial temperature increase is not converged. For the 800 and 1600-cell meshes, the time of initial increase is close to the converged value, but the final temperature plateau reached is slightly over-estimated. The resolution of 6400 cells is in very good agreement with the highest resolutions of 12800 and 51200 cells, for the time of initial increase, the final temperature plateau and the finer features of the graph. The resolution of 3200 cells is in sufficient agreement as well.
The two-dimensional analogue of the onedimensional cavity collapse simulation used in the convergence study described above is performed, with x-resolutions of 1600 and 3200 cells. Only half of the cavity is simulated, with the centreline of the cavity aligned with the bottom domain boundary in a domain of size 0.5 mm× 0.2 mm. The bottom boundary is reflective, exploiting the symmetry of the collapse process about the cavity centreline. In this way the effect of wave patterns generated during the collapse of the omitted cavity half are taken into account.
The y-resolution is adjusted accordingly to 640 and 1280 cells respectively, to obtain square cells. In Fig.  A.18 the waves generated during various stages of the two-dimensional simulations are compared in a mock-schlieren plot between the two grid resolutions. The upper half of each figure is generated from the 3200 × 1280 simulation while the lower half is generated from the 1600 × 640 simulation (reflected about its centreline).
The comparison shows that all the waves captured in the finer-grid simulation are also captured in the coarsergrid simulation. A comparison between the corresponding 3200×1280 and 6400×2560 simulations has shown that the wave pattern captured by these two simulations is the same. The same comparison between a 1600×800 and a 800 × 320 simulations revealed that the 800 × 320 simulation is less effective at capturing the small-scale waves in the lobes and in the ambient nitromethane during the late stages of the collapse.
Since the difference in the one-dimensional maximum temperature between resolutions 1600 and 3200 is not significant and the comparison between the twodimensional equivalent simulations shows no remarkable difference in the wave pattern, the choice of grid resolution is made to minimise computational cost. A three-dimensional simulation of resolution 3200 × 1280 × 1280 would require much more than 256 GB RAM (∼ 512 GB) and significantly more computational time than a simulation of resolution 1600 × 640 × 640. The latter requires ∼ 96 GB of RAM, a value that can be reduced for the early stages of the collapse if adaptive mesh refinement is used effectively. In general, without taking into account the effect of AMR, the memory required for a simulation scales like the cube of the resolution. As during the late stages of the collapse the wave pattern is complex, AMR does not reduce this scaling significantly. Hence, the choice of the 1600-cell resolution (∆x =0.3125 µm) is deemed most appropriate for the purpose of this work.