The evolution of the temperature field during cavity collapse in liquid nitromethane. Part II: Reactive case
TThe evolution of the temperature field during cavity collapse in liquidnitromethane. Part II: Reactive case
L.Michael a , N. Nikiforakis a a Laboratory for Scientific Computing, Cavendish Laboratory, Department of Physics, University of Cambridge, UK
Abstract
This work is concerned with the e ff ect of cavity collapse in non-ideal explosives as a means of controlling theirsensitivity. The main objective is to understand the origin of localised temperature peaks (hot spots) which playa leading order role at the early stages of ignition. To this end we perform two- and three-dimensional numericalsimulations of shock induced single gas-cavity collapse in liquid nitromethane. Ignition is the result of a complexinterplay between fluid dynamics and exothermic chemical reaction. In the first part of this work we focused on thehydrodynamic e ff ects in the collapse process by switching o ff the reaction terms in the mathematical formulation. Inthis part, we reinstate the reactive terms and study the collapse of the cavity in the presence of chemical reactions. Byusing a multi-phase formulation which overcomes current challenges of cavity collapse modelling in reactive mediawe account for the large density di ff erence across the material interface without generating spurious temperature peaksthus allowing the use of a temperature-based reaction rate law. The mathematical and physical models are validatedagainst experimental and analytic data. In Part I, we demonstrated that, compared to experiments, the generatedhot spots have a more complex topological structure and additional hot spots arise in regions away from the cavitycentreline. Here, we extend this by identifying which of the previously-determined high-temperature regions in factlead to ignition and comment on the reactive strength and reaction growth rate in the distinct hot spots. We demonstrateand quantify the sensitisation of nitromethane by the collapse of the isolated cavity by comparing the ignition timesof nitrometane due to cavity collapse and the ignition time of the neat material. The ignition in both the centrelinehot spots and the hot spots generated by Mach stems occurs in less than half the ignition time of the neat material.We compare two- and three-dimensional simulations to examine the change in topology, temperatures and reactivestrength of the hot spots by the third dimension. It is apparent that belated ignition times can be avoided by the useof three-dimensional simulations. The e ff ect of the chemical reactions on the topology and strength of the hot spotsin the timescales considered is also studied, in a comparison between inert and reactive simulations where maximumtemperature fields and their growth rates are examined. Keywords: condensed phase explosives, cavity collapse, temperature, nitromethane, hot spots, ignition
1. Introduction
This work is motivated by the necessity for optimis-ing the performance of non-ideal, inhomogeneous ex-plosives such as those used in mining. In order to in-crease their sensitivity and to control their performance,cavities are introduced in the bulk of the explosive, of-ten by means of glass micro-balloons. When a precur-sor shock wave passes through the explosive, these cav-ities collapse, generating regions of locally high pres-sure and temperature, which are commonly referred to
Email addresses: [email protected] (L.Michael), [email protected] (N. Nikiforakis) as hot spots. These lead to multiple local ignition sites,which cumulatively result in a shorter time to ignitionthan that of the neat material. Parameters such as thenumber, size, shape and distribution of the cavities af-fect the degree of sensitisation of the explosive. Under-standing the correlation between these parameters andthe reduced ignition time will allow better control of thebehaviour of the explosive.To this end, the collapse of cavities has been exten-sively studied in the past by means of experiment andnumerical simulation. An extensive literature review onprevious studies in inert materials is given in Part I ofthis work [1]. Hence here we indicatively only mention a r X i v : . [ phy s i c s . c o m p - ph ] O c t ome of the studies for cavities collapsing in reactiveliquid, gelatinous and solid explosives. To determine thee ff ect of the collapse on the hot spot generation and ex-plosive ignition, cavity collapse experiments in reactivematerials were performed by Bourne and Field [2, 3].The principal ignition mechanism was determined to bethe hydrodynamic heating due to jet impact. Conduc-tion from the highly-compressed, heated gas inside thecavity does not occur in the short time-scales governingthe particular shock wave-cavity configurations studied.Limitation of computational power restricted earlynumerical work on cavity collapse in reactive materi-als. Amongst the first studies on cavity collapse in areactive medium was the work by Bourne and Milne[4, 5] who observed numerically the same loci of hotspots as in experiments. More recently the authors pre-sented limited results on simulating cavity collapse inreacting nitromethane [6]. Kapila et al. [7] presented thecollapse process and the detonation generation in an ex-emplary explosive and studied the e ff ect of cavity shapein the detonation generation using a pressure-dependentreaction rate. In an elastoplastic framework, Tran andUdaykumar [8] studied the response of reactive HMXwith micron-sized cavities, while Rai et al. [9, 10] con-sidered the resolution required for reactive cavity col-lapse simulations and the sensitivity behaviour of elon-gated cavities in HMX.Despite technological advances, performing com-plete numerical simulations of ignition due to shock-induced cavity collapse still poses many challenges.These include the use of complex equations ofstate to describe the explosive materials, maintainingoscillation-free pressure, velocity and temperature fieldsacross the cavities material boundaries upon their inter-action with shock waves, sustaining (at least) 1000:1density di ff erence across these boundaries, retrievingphysically accurate temperature fields in the explosivematrix and numerically modelling the ignition of thematerial as a temperature-driven phenomenon. More-over, the computational power needed for accurately re-solving the complete phenomenon in three-dimensionsis still large. As the explosive initiation is a temperature-driven phenomenon, the challenges regarding the tem-perature field are of critical importance. These chal-lenges are described in detail in art I of this work.A complete physical simulation of the initiation ofa condensed phase explosive due to cavity collapsehas several requirements; a three-dimensional frame-work, realistic material models (equations of state),oscillation-free material interfaces, the ability to re-cover accurate and oscillation-free temperature fields.A temperature-dependent reaction rate law is also de- sirable, to mathematically describe ignition to be drivenby the heating of the material. Moreover, each com-ponent used should be validated alone and in combina-tion with all the components composing the numericalframework.Simulations presented in the literature satisfy somebut not all of these requirements. In this work we exer-cise the mathematical model (MiNi16) proposed by theauthors in previous work [11] to overcome the di ffi cul-ties in numerically simulating the cavity collapse andmove towards a complete simulation of explosive initi-ation due to cavity collapse. In Part I of this work [1]we simulated the three-dimensional collapse of isolatedair cavities in nitromethane using a validated equationof state (Cochran-Chan). We looked in depth how thehydrodynamical e ff ects in absence of reaction (e.g. gen-eration and propagation of waves) lead to local temper-ature elevation. Such regions were identified as candi-dates for critical hot spots in a reactive simulation andthe necessity for three-dimensional (as opposed to two-dimensional) simulation was justified. In this secondpart of the work, we take advantage of oscillation-freeand reliable temperature fields that can be recovered byusing the MiNi16 model and extend the work of Part I toreactive scenarios. We perform three-dimensional simu-lations of the collapse of isolated air cavities in reacting,liquid nitromethane, using equations of state in Mie-Gr¨uneisen form and a temperature-dependent reaction-rate law. We study in detail the ignition process and welink the evolution of the reaction-progress variable tothe temperature elevations and the wave-pattern gener-ated during the collapse process. We identify the react-ing hot spots and study their relative reactive strengthand reaction growth rates. The e ff ect of the cavity col-lapse on shortening the time to ignition is illustratedexplicitly by comparing the ignition due to cavity col-lapse against the ignition of the neat material. We alsodemonstrate the necessity for three-dimensional simu-lations (compared to 2D) by looking at the percentageof burnt material over time in the two scenarios but alsolooking at the evolution of waves and temperature fields.Moreover, we compare inert and reactive simulations toexamine the added e ff ect of the reactions on the tem-perature fields and the topology of the hot spots at thetimescales considered.The rest of the paper is outlined as follows: the nextsection presents the underlying mathematical formula-tion in terms of the governing partial di ff erential equa-tions, the equations of state that close the system andthe form and calibration of the reaction rate law for ni-tromethane combustion. A section on validation fol-lows, where we compare numerical results against the-2retical and experimental temperatures in shocked ni-tromethane and the CJ and von Neumann values forsteady state detonation. The ignition regime is vali-dated in a similar way, by comparing numerical and ex-perimental times-to-ignition for various input pressureswith the use of an ignition Pop-plot. In the results sec-tion, we consider the collapse of a gas cavity in liquidnitromethane, follow the events leading to the genera-tion of local temperatures which are more than threetimes the post-shock temperature of the neat materialand compare and analyse the di ff erence between the 2Dand 3D simulations as well as inert and reactive simula-tions in the context mentioned earlier.
2. Mathematical and physical model
A formulation which rectifies the issues commonlypresented in shock-bubble simulations was proposed byMichael and Nikiforakis [11]. This formulation consid-ers the cavity as an inert phase ( phase 1 ) and the sur-rounding material as a reacting phase ( phase 2 ) com-posed of two materials; the reactant nitromethane as ma-terial α and the gaseous products of reaction as material β . Mixing rules are employed to determine the proper-ties of phase 2 from the properties of the two materi-als. Mixing rules are also in e ff ect between phase 1 andphase 2 across material interfaces where the di ff usionzones lie. In this work, we neglect the e ff ect of reactionproducts and thus use a reduced form of the formulation.Consider the gas inside the cavity to be phase 1 andthe liquid nitromethane around the cavity to be phase 2 .Then, the governing equations for this system take theform: ∂ z ρ ∂ t + ∇ · ( z ρ u ) = ,∂ z ρ ∂ t + ∇ · ( z ρ u ) = ,∂∂ t ( ρ u k ) + ∇ · ( ρ u k u ) + ∂ p ∂ x k = , (1) ∂∂ t ( ρ E ) + ∇ · [( ρ E + p ) u ] = ,∂ z ∂ t + u · ∇ z = ,∂ z ρ λ∂ t + ∇ · ( z ρ u λ ) = z ρ K , where for i = , ρ i are the densities for the air andnitromethane, z i are their corresponding volume frac-tions ( z + z = ρ is the total density given by ρ = z ρ + z ρ , u is the velocity vector and p is thetotal pressure. The total specific energy is given by E = u + e , where e is the total specific internal en-ergy. Also, λ is a reaction progress variable and K rep-resents the reaction source terms, to be defined later.The mixture rule for the total internal energy is given by ρ e = ρ z e + ρ z e , where e i for i = , ξ = γ − , where γ is the total adiabatic index, is alsorequired and in this case is given by ξ = z ξ + z ξ . Thesound speed for the total mixture is given by ξ c = y ξ c + y ξ c , (2)where y i is the mass fraction of phase i , given by y i = ρ i z i ρ , for i = , Q is included in thereference energy function of the reactant, representingthe heat of detonation released upon reaction, such that e ref = e ref + Q . Alternatively, this term can be in incor-porated as a source term to the energy equation. In order to model the reactions in nitromethane, asingle-step, temperature-dependent Arrhenius reactionrate law is used, of the form K = d λ dt = − λ Ce − T A / T NM , (3)where C is a constant pre-exponential factor and T A isthe activation temperature of the material.Many sets of the reaction rate parameters are avail-able in the literature for liquid nitromethane. For ex-ample, ( C , T A ) = (2 . × s − , 11 500 K) is used byMeniko ff and Shaw [13], (6 . × s − ,
14 400 K) byTarver and Urtiew [14] and (1 . × s − , 20 110 K)by Ripley et al. [15]. In this work we adopt the pre-exponential factor proposed by Meniko ff and Shaw [13]and adjust the activation temperature to T A =
11 350 Kto match the experimentally calculated overtake timeand the shape of the velocity versus distance graph ofthe shock-induced ignition experiment in She ffi eld et al.[16] (neat nitromethane, shocked at 9 . Q , we follow the approach described by Ari-enti [17]. This involves varying the parameter Q to3atch as closely as possible the experimental value ofthe pressure at the CJ point of 12 . p - v plane. In general, varying the value of Q shifts the re-active Hugoniots upwards. When the Rayleigh line be-comes tangent to the reactive (pseudo-product) Hugo-niot, the appropriate value for Q is determined. Here, Q = . × J kg − .We note that the set of calibrated parameters givenabove is valid for steady-state detonation propagationof neat nitromethane and no further adjustment is nec-essary for ignition case studies, as the results in the val-idation section indicate.
3. Validation
System (1) is integrated numerically using a high res-olution shock-capturing numerical scheme, namely theMUSCL-Hancock finite volume method with an under-lying HLLC Riemann solver. Hierarchical, structured,adaptive mesh refinement (AMR) is used to dynami-cally increase the resolution locally [19]. The sourceterms that describe chemical reactions are integratedwith a high-order Runge-Kutta scheme.The aim of this section is to validate the resultingcode and also to assess whether the combination of pa-rameters described in the previous section can matchthe experimentally determined von Neumann spike, theCJ values and the ignition pop-plot, without any fur-ther user adjustment. The inert formulation and physi-cal model were validated in Part I of this work, showingnon-oscillatory hydrodynamic fields and recovery of re-alistic temperature fields in shocked nitromethane.
Having established a reasonable post-shock tempera-ture, we now assess whether the code predicts the cor-rect CJ and von Neumann peak values. To this end, aone-dimensional slab of liquid nitromethane is shockedto 23 GPa, a value close to the von Neumann pressure.Following an initial unsteady evolution, the simu-lation settles to a steady detonation wave, yieldinga constant value of the von Neumann spike pressureof 22 . . For the final part of this assessment, it is worth re-calling that we are mainly interested in the early igni-tion stages of the shock-to-detonation process. It is not unusual to find that, depending on the mathematical for-mulation of the model, the reaction rate parameters usedfor detonation modelling are not suitable for ignitionstudies. This is attributed to the fact that the parametersare adjusted to fit post-shock temperatures and steadydetonation values. To assess whether the current set-up can be employed for arbitrary studies without anyfurther adjustment, we compare our numerical resultsagainst experimental Pop-plot data that show inductiontime (i.e. time to ignition) versus input pressure. The ex-
Figure 1: Ignition time vs input pressure Pop-plot. The black filledcircles represent experimental data and the open squares numericallycalculated data. act time when ignition occurs is problematic, since ide-ally, one would use the same definition of time of igni-tion for both numerical and experimental results. In nu-merical simulations, the time of ignition can be definedas the time when a specific fraction of the explosive ma-terial has reacted. The selection of the appropriate per-centage should be based on experiments. However, ex-perimentally it is rather di ffi cult, where possible at all,to measure the chemical species, and ignition is usuallymeasured using luminosity, something that is not trivialto accurately and consistently retrieve from the simula-tions. Thus, we arbitrarily define ignition to be the timewhen a small percentage (namely 10%) of the explosivematerial has reacted and compare our numerical resultsagainst experimental data taken from Berke et al. [21],Hardesty [22] and Chaiken [23]. This comparison ispresented in Fig. 1 and shows that the numerical resultsfall well within the range of the experimental ignitiondata.
4. Ignition of neat nitromethane
In condensed phase explosives, initiation can beachieved by the passage of a shock wave through the4uiescent explosive, raising the temperature and pres-sure and triggering the start of reaction.To study the initiation process of nitromethane, sim-ulations of the shock-induced ignition of neat, liquid ni-tromethane are performed.Since the ignition process of this application is purelyone-dimensional, we consider a domain of dimensions[0 mm , . ff ective reso-lution of 2560 cells ( dx = . µ m). This resolutionresolves well the steady state reaction zone of liquid ni-tromethane which is calculated from experiment [20] tohave width of ∼ µ m. The initial conditions for thistest are given in Table 1.When a shock wave is set up numerically as an ini-tial condition, a start up error is generated due to thesymmetric Riemann problem [24]. This error has theform of a small well in the density distribution behindthe shock wave, which translates into a small hill inthe temperature field. Since the reaction rate we useto model the reactions in liquid nitromethane dependsexponentially on the temperature, even a disturbance ofa small magnitude in this field (here ∼
20 kg m − ) wouldrapidly grow, generating a spurious hot-spot. The hot-spot would, in turn, lead to ignition earlier than it wouldbe expected if the density field was clean. In order toremove the disturbance by extrapolating from the non-disturbed shocked state, or by cutting the domain at apoint after the start up error, the disturbance has to besu ffi ciently formed. Thus, before treating the error, theshock wave is allowed to travel some small but signifi-cant distance from its initial position. If reactions wereturned on during this travel time, the numerical hot spotwould a ff ect the state behind the shock wave. To over-come this, an inert shock wave is allowed to travel thedistance required for the start up error to be adequatelyformed and then the part of the domain that is more than5 cells behind the shock wave is cut o ff . After the cuto ff ,the simulation is restarted with the reactions turned on.The cut of the domain is performed at time 0 . µ sand reaction is allowed to start at time 0 . µ s. Allthe times referred hereafter are relative to the reactionstart time. Also, all the positions are relative to the cut-o ff position ( x = .
123 mm), as the domain is consid-ered to be repositioned at x = ff .The evolution of the reaction progress variable ( λ ),pressure ( p ) and nitromethane temperature ( T NM ) is il-lustrated in Fig. 2. As can be seen in the early stages ofFigs. 2(a), (c) and (e), the fuel that is closer to the inci-dent shock wave is shocked and heated first. Hence, ithas more time to react than the fuel that is further awayfrom the shock. As a result, temperature and mass- fraction gradients are generated. During this process,explosive material is burnt and the reaction progressvariable starts to deviate away from 1. By defining asthe ignition time the time when λ = .
9, we observe thathere ignition occurs at stage 7 of Fig. 2(a) [ λ -plot], cor-responding to t ign = . µ s. At the end of this slowly-evolving induction phase, the fluid cannot sustain thehigh pressure and temperature behind the shock wave,resulting to the generation of a signal, as seen duringstages 13–14 of Fig. 2(c) [pressure plot]. At this time( t ∼ . µ s), thermal runaway is considered to occur.A rapid reaction stage follows the ignition stage, usuallycalled the transition to detonation phase, during whichthe generated pulse is growing (Fig. 2(d)). At stage 23,the reaction wave overtakes the leading shock wave de-picted in Fig. 2(d). The overtake is accompanied by arapid increase of pressure, temperature and reaction (de-crease of λ ). Thereafter, the detonation structure settlesdown towards a steady-state solution.
5. The collapse of a single cavity in reactive liquidnitromethane
In this section, we consider an isolated gas-filled cav-ity of radius 0 .
08 mm collapsing in reacting liquid ni-tromethane in the domain spanning [0 mm , . × [0 .
75 mm , .
54 mm]. The initial conditions in theshocked, pre-shocked and cavity regions are given inTable 2. The domain is longer than in the inert simula-tions to allow for su ffi cient reaction to take place. It isalso taller, to avoid the minor reflection of waves fromthe top boundary (even with transmissive boundaries)that could a ff ect the ignition process. We adopt the sameabbreviations as in Part I.The evolution of the reaction progress variable λ (left), pressure, p (middle) and nitromethane tempera-ture, T NM (right) is illustrated at selected times in Fig.3. The incident shock wave (ISW) travels within the ni-tromethane, compressing the material to 10 . ff ect of eachwave on the temperature field. The wave patterns inthis scenario are the same as those in Part I and thus wewill not repeat the details of their generation. We willonly describe the additional e ff ects observed due to thepresence of reactions.At the initial stages of the collapse (see Fig. 3(a,b))and away from the cavity (where the rarefaction wave(RW) from the collapse has not yet arrived), the re-action progress variable evolves as it would in a 1D5 able 1: Initial conditions for the initiation of neat nitromethane. Region ρ [kg m − ] ρ [kg m − ] u [m s − ] p [Pa] λ z shocked nitromethane 1934 1934 2000 10 . × − ambient nitromethane 1134 1134 0 10 − (a) (b)(c) (d)(e) (f) Figure 2: One-dimensional evolution of the ignition and transition to detonation of neat nitromethane under the influence of a10 .
98 GPa shock-waveat selected times. On the left, stages 0 −
14 are illustrated and on the right stages 15 −
29. The top row illustrates the reaction progress variable ( λ ),the middle row the pressure ( p ) and the bottom row the nitromethane temperature ( T NM ). neat nitromethane experiment. The expansion gener-ated by the RW leads to a lowering of the tempera-ture within the jet and this a ff ects the shape of the re-action zone. As a result, before the cavity collapses,the highest temperature in the nitromethane is in theuniform, unperturbed region behind the ISW and onlyminimal reaction is observed. After the cavity col- lapse (at t = . µ s) we observe for the first timein the collapse process, temperatures that are abovethe post-shock temperature. Specifically, the back col-lapse shock wave (BCSW) generates temperatures inthe range 1263 K − λ ≈ .
97. Thefront collapse shock wave (FCSW) generates tempera-6 able 2: Initial conditions for the initiation of nitromethane by cavity collapse.
Region ρ [kg m − ] ρ [kg m − ] u [m s − ] v [m s − ] p [Pa] λ z shocked nitromethane 2.388 1934 2000 0 10 . × − ambient nitromethane 1.2 1134 0 0 10 − air cavity 1.2 1134 0 0 10 − − Figure 3: The ignition process in the vicinity of a cavity embedded in liquid nitromethane following the passage of a 10 .
98 GPa shock wave.The evolution of mass fraction λ (left), pressure p (middle) and nitromethane temperature T NM (right) at times t = (a) 0 . µ s, (b) 0 . µ s, (c)0 . µ s and (d) 0 . µ s is presented. The horizontal and vertical axes are space ordinates in µ m. tures within the range 1263 K − λ ≈ . t ign = . µ s in the back hot spot (BHS).The maximum temperature in the front hot spot (FHS)reaches at this point 2068 K and 3089 K at the BHS.For this set up, the ISW proved slower than the ni- tromethane jet. Hence, at the time of collapse, the ISWis still traversing around the lobes. As a result of pass-ing over the upper part of the upper lobe (and equiva-lently the lower part of the bottom lobe) many transmis-sion / reflection processes take place. These processes re-sult in a coalescence of waves (LSW) that are transmit-ted from the lobes into the pre-shocked (by the ISW)nitromethane residing around the lobes (Fig. 3(d)). Thisleads to the generation of temperatures higher than thepost-shock temperature in the regions around the lobes.However, this temperature rise is not su ffi cient to gener-7te a new ignition site in these timescales.The superposition of the ISW and the FCSW gener-ates the Mach stem hot spot (MSHS) as described inPart I (Fig. 3(d)), which in the reactive scenario enclosestemperatures of the range 1263K-1493K. At this point,the highest overall temperature is still observed in theBHS (2773 K). As the Mach stem region grows, thetemperature it encloses increases rapidly, reaching val-ues of the order of 2700 K.Ignition is observed in theMSHS at t = . µ s.To the rear of the cavity, the waves emanating fromthe top and bottom lobe are superposed. Parts of themare superposed with the BCSW as well (see Fig. 3(d)).This results in regions that have been shocked twice oreven three times and hence leads to temperatures of theorder of 1600 K.In the final stages of the collapse process, the remainsof the cavity are advected downstream and more burn-ing is observed in the reaction sites, with λ reaching avalue of 0 . .
265 in theMSHS by the end of the simulation.
Figure 4: Percentage of explosive burnt (top) and rate of burning (bot-tom) in the centreline hot spots and Mach stem hot spots from the timeof their birth.
In Fig. 4 (top) the percentage of burning given as themaximum value of 100 × (1 − λ ) in the centreline hotspots (BHS and FHS) and in the MSHS is shown overtime. As time 0 we denote the time of birth of the hotspots which is the ignition time in each location iden-tified earlier. It can be seen that the burning in thecentreline hot spots (CHS) follows a √ x graph whilethe MSHS burning increases roughly linearly. In Fig.4 (bottom) the rate of burning in the two hot spots ispresented. It can be seen that the reactions in the CHSgrow faster initially than the reactions in the MSHS al-though this is reverted at later times and the reactions inthe MSHS grow faster than the reaction in the CHS. It is informative to consider the evolution of the flowfield and its e ff ect on the temperature along lines ofconstant latitude. First, we discuss events along y = .
21 mm. This provides insight to the hot spot genera-tion at the rear of the cavity (BHS), the minimal reactionat the front of the cavity (FHS) and the temperature dis-tribution close to the centreline of the cavity. The evo-lution of the reaction rate variable λ , pressure and ni-tromethane temperature along this line is shown in Fig.5. At the early stages, a slight increase in the tempera-ture field by the passage of the ISW is seen. Upon theinteraction of the ISW with the cavity, a downstream-travelling air shock is generated and an upstream travel-ling RW. The e ff ect of the RW is a decrease in pressureand temperature. As it propagates, its e ff ect becomesmore profound, manifested as a further local decrease inthe pressure and temperature fields and rate of reaction.The evolving crest-like feature (following the local de-crease) seen in the pressure and temperature plots (Fig.5(a)), is due to the formation of the nitromethane jet.As the cavity collapses, the lineout crosses bothshock waves (FCSW and BCSW), seen as two high-pressure and high-temperature fronts moving awayfrom each other, labelled as 1 , , , λ -plot in Fig. 5(b)). The reaction in theFHS is observed to be comparatively slow. Within theBHS and along this lineout, ignition is seen to occur at t = . µ s.Waves emanate from the lobes after the collapse ofthe cavity and are are labelled 5 as they intersect thisline. By stage 26 they are superposed with the BCSWalong this lineout. Neither the downward waves northeir superposition with the BCSW result in a consid-erable increase in the temperature field. In fact, thehighest temperature still occurs in the BHS. At stage27, the upward wave from the lower lobe (not mod-elled explicitly but with its behaviour implicitly takeninto account by the reflective lower boundary) reachesthe line y = .
21 mm. Part of this wave is superposed8 igure 5: Lineouts along y = .
21 mm from figure 3. The mass fraction (top row), the pressure (middle row) and the nitromethane temperature(bottom row) are shown at stages (a) 1–18, (b) 19–37, (c) 38–44, corresponding to times (a) 0 .
002 28 µ s–0 . µ s, (b) 0 . µ s–0 . µ s and(c) 0 . µ s–0 . µ s. with both the BCSW and the wave entities emanatingfrom the upper lobe and part of it is superposed withthe waves emanating from the upper lobe only. Thesesuperpositions increase the pressure locally but do nothave a significant e ff ect on the temperature. As theseshock waves (labelled 6 ) move away from the cavity,the temperature and pressure along the line y = .
21 mmdo not increase further, but the reaction in the BHS stillincreases.Focusing on the reaction progress variable field illus-trated in the λ − plot of Fig. 5(c), it is observed that, atall late stages, the largest amount of reaction occurs inthe BHS. Some reaction is also seen in the FHS and inthe region that is traversed by the BCSW and the wavesemanating from the lobes. After the ignition at stages21–22, the fuel continues burning until it reaches thethreshold of λ = .
01, where no more fuel is consideredto be available for burning.After stage 35, the advection of the rear remainingparts of the cavity reach the lineout, leading to the lowpressure and temperature parts of the lineouts and the λ = y = .
29 mm are used to give insightabout the MSHS and in general about the temperature distribution above the cavity. The evolution of λ , pres-sure and nitromethane temperature on this line is illus-trated in Fig. 6.The e ff ect of the rarefaction wave is seen in all threefields from stages 8-9 in this case. At stage 26 the ef-fect of the S , , wave entity as described in Part I ofthis work is seen in the plots of Fig. 6(b). The FCSWovertakes the ISW at stage 28 This results in the in-crease of the temperature and a subsequent increase ofthe reaction. As the line y = .
29 mm goes though alarge part of the MSHS, the elevated temperature is notfound in an isolated peak, but is distributed within a re-gion behind the moving Mach stem front (Fig. 6(b), la-belled as 7 ). This elevated temperature region resultsin the increase of reaction, as illustrated by the λ plotin Fig. 6(b). From stages 44 onward, the temperatureis first seen to increase in a hill-like form (labelled as8 ), then decreases taking a well-like form ( 9 ) andthen increases again taking the form of a hill, but witha shallower gradient than before ( 10 ). This is becauseat these late stages, the Mach stem feature has attaineda horn-like shape and it is advected upwards. As a re-sult, the line y = .
29 mm intersects the stem, then goesthrough a region below the stem and then intersects the9 igure 6: Lineouts along y = .
29 mm from figure 3. The mass fraction (top row), the pressure (middle row) and the nitromethane temperature(bottom row) are illustrated at stages (a) 1–18, (b) 19–37, (c) 38–44, corresponding to times (a) 0 . µ s–0 . µ s, (b) 0 . µ s–0 . µ s and(c) 0 . µ s–0 . µ s. stem again. At the location of the second intersection,the temperature is elevated compared to the regions out-side the stem, but it is lower than at the location of thefirst intersection. This translates as a well-hill-well fea-ture in the λ − plot, indicating more reaction along thepart of the line y = .
29 mm that intersects the stem thefirst time, far less reaction outside the stem and a mod-erate amount of reaction at the second intersection.Another horizontal lineout, at y = .
35 mm, is usedto give insight about the generation of the MSHS andtemperature distribution and initiation around the cav-ity. The evolution of the λ -field, the pressure and the ni-tromethane temperature along the line are given in Fig.7. On this line, the e ff ect of all the waves emanatingfrom the cavity collapse process are seen later than onthe previous lineouts considered.In the early stages of Fig. 7(a), a slight increase in thepressure and a more noticeable increase in T NM due tothe passage of the shock wave is seen (labelled as 11in the temperature plot). This is accompanied by thestart of the reaction behind the shock wave, seen in the λ − plot. The slight increase of the post-shock pressure,the formation of a temperature gradient and the shapeof the λ − plot are as observed in the shock-induced ig- nition of neat nitromethane, because no e ff ects from thecollapse of the cavity have reached the line y = .
35 mmyet.By stage 14, the rarefaction wave (RW) has reached y = .
35 mm (as opposed to stages 8-9 on y = .
29 mm)and its e ff ect is seen as a descent in the pressure andtemperature plot (labelled as 12 in the pressure plot).The e ff ect of the RW is increased as the wave propa-gates upwards, seen as growth of the dip in pressureand temperature (Fig. 7(b)). In the λ − plot, the shape ofthe lower part of a graph is a ff ected, presenting a slightincrease of the gradient of the straight-line part of thegraph. This indicates that less reaction is taking placewhen the rarefaction wave is present compared to thereaction that would have taken place in the absence ofthe rarefaction wave.At stage 34, the entity of waves S , , is seen tohave crossed the line y = .
35 mm (as opposed to stage26 for y = .
29 mm). Upon reaching y = .
35 mm, thewave increases the pressure and temperature along thepart of the line that intersects it. As this wave is cir-cular (spherical in 3D) and propagates outwards fromthe cavity, the area cut by the line increases with time.As a result, in the one-dimensional pressure plot of Fig.10 igure 7: Lineouts along y = .
35 mm from figure 3. The mass fraction (top row), the pressure (middle row) and the nitromethane temperature(bottom row) are illustrated at stages (a) 1–18, (b) 19–37, (c) 38–44, corresponding to times (a) 0 .
002 28 µ s–0 . µ s, (b) 0 . µ s–0 . µ s and(c) 0 . µ s–0 . µ s. ff ect of thewave entity S , , on the one-dimensional tempera-ture field is seen in Fig. 7(b) as an increase in T NM , inthe form of two temperature fronts. Its e ff ect on the re-action progress variable is seen in 7(b). It appears as adip in λ attributed to the backward moving part of thewave, as a sudden but not very distinct drop in λ at-tributed to the front part of the wave and as an overalldecrease of the gradient of the straight-line part of thegraph. All these features of the λ − plot indicate rapidincrease of reaction as soon as the new shock waves re-shock the area.Stage 38 (Fig. 7(c)) corresponds to the time before theFCSW overtakes the ISW on line y = .
35 mm (as op-posed to stage 28 for y = .
29 mm) and stage 39 to thetime immediately after the overtake (labelled as 15in the pressure plot). After the overtake, the lineout Overtake in this context refers to the overtake of the ISW by theCSW and not the overtake of the ISW by a detonation wave. goes through the high-temperature Mach stem region,leading to the generation of the temperature peaks inthe temperature plot (labelled as 7 ). This leads to anincrease of the reaction, as indicated by the sudden de-crease of the reaction progress variable in the λ − plot.The temperature of the hot spot in the Mach stem con-tinues to gradually increase. This translates into thetemperature peaks of Fig. 7(c), resulting in the increaseof the reaction and the sudden drop in λ in Fig. 7(c).
6. Comparing the three-dimensional and two-dimensional, cavity collapse in reacting ni-tromethane
In this section, the three-dimensional (3D) collapse ofa cavity in reacting nitromethane is presented and com-pared to the two-dimensional (2D) equivalent simula-tion presented in the previous section. Selected stages ofthe collapse are shown in Fig. 8. The three-dimensional z = . λ = . λ field throughthe centre of the cavity and project it on the left half of11 igure 8: Snapshots of the three-dimensional collapse of a cavity in reacting liquid nitromethane. The blue, three-dimensional z = . λ = . λ and nitromethane temperature fields. The slices are projected on the left half( λ field) and right half (nitromethane temperature field) of each figure. Note that the collapse process is shown here to move from bottom to toprather than left to right. The images here can be considered as rotated by 90-degrees counter-clockwise compared to Fig. 3. This is done solely forillustration purposes. each figure. This again illustrates the reaction regions.Similarly, the projection of the nitromethane tempera-ture field on the same plane is seen on the right half ofeach figure. It is therefore presented how the genera-tion of locally high temperatures leads to the generationof hot-spots. The highest temperatures are located infront and behind the point of collapse of the cavity (FHSand BHS) as well as in the Mach stem generated at latestages of the collapse due to wave superposition. Theselead to the generation of bell-shaped hot spots, mainlybehind the point of collapse (BHS) and a torus-shapedhot spot corresponding to the three-dimensional Machstem. The comparison between the 3D and 2D collapseprocess is performed by using planar pseudocolour plotsof the temperature field (Fig. 9) and the evolution ofthe temperature and λ fields along lines of constant lat-itude, specifically along the centreline of the cavity and y = µ m (Figs. 10,11).The three-dimensional nature of the flow field aroundthe cavity results to a faster jet compared to the two-dimensional equivalent scenario. This is seen in Fig.9(a)-(b) and Fig. 10(a) and e ff ectively results in a fasternitromethane jet (1.4 times faster) in the 3D case thanin the 2D case and the earlier observation of all subse-quent features of the collapse phenomenon. The cavitycollapses, as a result, faster by ∼ . µ s in the 3D case(Fig. 9(c)-(d)) and the generated collapse shock waveslead to a quicker ignition in the front and back hot spots (Fig. 10(b)). The temperatures achieved upon collapserange between 400K and 1500K higher in the 3D sce-nario than the 2D. The percentage of explosive burnt inthe centreline hot-spots in the centreline hot spot regionis shown as a function of time in Fig. 12. It can be seenthat the higher temperature attained upon the 3D col-lapse generates a higher immediate burn of the materialin 3D (59% burnt) compared to 2D (11% burnt) and afaster overall ignition and burning procedure. The Machstem and hence the MSHS are also generated quickerthan in the 2D case (Fig. 9(e)-(f) and Fig. 11(b)). TheMach stem and the MSHS grow faster in the 3D caseas well, as seen in Fig. 9(g)-(h) and Fig. 11(c)-(d). Ingeneral, in the 3D case, higher temperatures are alsoachieved compared to the 2D case, as it is evident inFig. 10(c)-(d) and Fig. 11(c)-(d) leading to faster imme-diate ignition.
7. Comparing the cavity collapse in inert and react-ing nitromethane
In Part I of this work the collapse of a cavity in non-reacting nitromethane was studied and this work ex-tended that to the equivalent reacting medium. In thissection the temperature field generated through the col-lapse process in the inert and reacting media is com-pared, to examine the e ff ect of the reactions on the tem-perature field and the hot spot topology. The compar-12 a) 2D t = . µ s (b) 3D t = . µ s(c) 2D t = . µ s (d) 3D t = . µ s(e) 2D t = . µ s (f) 3D t = . µ s(g) 2D t = . µ s (h) 3D t = . µ s Figure 9: Comparison of the temperature field and collapse times between a two-dimensional and a three-dimensional of the cavity collapse undera 10 . ison uses the two-dimensional scenarios for simplicity.The number of hot spots generated in the two cases isidentical and the location is roughly the same; any dif-ference observed is of the order of 3 −
4% of the bubbleradius, and mainly in the BHS and FHS. The signifi-cant di ff erence in the two collapse cases is seen in thetemperature field. At the early stages of the simulation, before the collapse takes place, some initial reaction isobserved behind the shock wave, which is translatedto some increase in the temperature field in the reac-tive case, compared to the inert case. The di ff erence inthe temperature fields of the two cases is small, rangingfrom 0 to 50K. Once the cavity collapses, the BCSWand FCSW generate significant reaction in the BHS and13 a) t = . µ s (b) t = . µ s(c) t = . µ s (d) t = . µ s Figure 10: Comparison of the temperature field and reaction along y = .
21 mm in the 2D and 3D configurations. (a) t = . µ s (b) t = . µ s(c) t = . µ s (d) t = . µ s Figure 11: Comparison of the temperature field and reaction along y = .
29 mm in the 2D and 3D configurations.
FHS leading to an even larger increase of temperature inthe reactive case and a di ff erence of 65K from the inertcase. This di ff erence grows as the reaction progressesand the hot spots grow in size and strength. This can be seen in Fig. 13, where the maximum temperature inthe centreline hot spots (CHS) i.e. BHS and FHS andin the MSHS is shown. The general trend is that af-ter the collapse, the maximum temperature observed in14 igure 12: Percentage of explosive burnt over time in the three-dimensional and two-dimensional scenarios.Figure 13: Maximum nitromethane temperature observed in the cen-treline hot spots (BHS and FHS, labelled collectively as CHS) plottedusing circles and in the MSHS plotted using squares. The red colordenotes results from the reactive setup and the blue color from theinert set up. the simulation in the reactive case increases, whereas itdecreases in the inert scenario. This is expected as thereaction supports the shock waves and vice versa. Inthe reactive case, the maximum temperature observedcan be higher than in the inert scenario by as much as1500K. For this configuration and the timescales stud-ied the temperatures in the MSHS are always lower thanthe CHS. Moreover, the di ff erence between the CHSand MSHS is larger in the reactive case than in the inertcase.The comparison of the rate of increase of the tem-perature in the CHS and MSHS between the inert andreactive scenarios is presented in Fig. 14. It is clear thatthe rate of increase of the maximum temperature in bothhot spots is positive in the reactive case, i.e. the maxi-mum temperature is always increasing, whereas this isnot true for the inert scenario. The first peak in both locicorresponds to the birth of the hot spots. In the CHSa second, late peak is observed which corresponds tothe superposition of the waves emanating from the twolobes along the centreline of the cavity. In the MSHS (a) CHS(b) MSHS Figure 14: Comparison of the rate of increase of maximum tempera-ture (a) in the centreline hot spots (CHS) and (b) in the Mach stem hotspot (MSHS) in the inert and reactive scenarios. the increase of the maximum temperature relies purelyon the waves generating the Mach stem. In Fig. 15 weobserve that the rate of increase in the MSHS in the inertcase is always higher than the CHS (except at the birthof the hot spot). This is, however not true for the reac-tive case. Comparing information from Figs. 13-15, wecan conclude that in maximum temperature terms andin this configuration, the CHS are stronger in absolutevalue than the MSHS, although the growth rate of reac-tion in the MSHS was seen to be higher than the CHSat late stages of the collapse in Fig. 4. In the inert case,however, the rate of temperature increase in the MSHSis higher than in the CHS (except at the bubble collapsetime) so the MSHS could, in longer timescales lead tohigher temperatures than the CHS. This could also betrue in the reactive case, though more detailed, late timesimulations would be needed to determine that.Using lineouts, we look at the temperatures in thedi ff erent hot spots. The local temperature di ff erence ishigher in the late stages after the collapse so we omitillustrations of the early times. On y = .
210 mm at t = . µ s (Fig. 16(a)), the temperature di ff erence inthe BHS is more significant than the FHS. This con-15 a) inert(b) reactive Figure 15: Comparison of the rate of increase of maximum tempera-ture in the two hot spots (a) in the inert and (b) the reactive scenarios. tinues to be true at later times, until the two hot spotsmerge ( t = . µ s) in Fig. 16(b). The temperaturedi ff erence in the MSHS is also significant, as seen on y = .
29 mm in Fig. 17(a) on y = .
35 mm in Fig. 17(b)for t = . µ s). The width of the MSHS is also largerin the reactive case.
8. Conclusions
In this work we perform resolved numerical sim-ulations of cavity collapse in liquid nitromethane us-ing a multi-phase formulation, which can recover re-liable temperatures in the vicinity of the cavity. Con-siderable care is taken regarding the form of equationsand numerical algorithm to eliminate spurious numeri-cal oscillations in the temperature field. The model isvalidated against experimental data. Specifically, wedemonstrate that the deduced CJ and von Neumann val-ues match the values found in the literature and alsothat the experimentally determined ignition Pop-plotdata are matched. Following the validation, we studythe shock-induced cavity collapse in reacting liquid ni-tromethane, in two and three dimensions and follow the events leading to the generation of local temperaturesand initiation of the explosive.Working towards elucidating the relative contributionof fluid dynamics and chemical reaction, we examinedin Part I of this work the details of the hydrodynamic ef-fects that lead to local temperature elevations and iden-tified a more complex hot spot topology than previ-ously described in the literature. In this second partof the work we demonstrate the e ff ect of chemical re-actions that acts additively yo the fluid dynamics. Weidentify which high temperature regions lead to reactivehot spots and observe much higher temperatures (40K-1500K) than in Part I.We examine additionally the ignition of nitromethanein the absence of cavities and compare the ignition timesfor the neat and single-cavity material. We observe thatthe initiation of nitromethane in the presence of an iso-lated collapsing cavity is reduced to less than one third(in 2D simulations or less than a quarter in 3D simula-tions) of the required time for igniting the neat material.This quantifies the sensitisation character of the cavitiesin this configuration.It is observed that the highest nitromethane temper-atures still occurs in the back hot spot (BHS) as alsoobserved in Part I and this leads to the first ignition sitethat resides along the centreline of the cavity. Moreover,the temperatures in the Mach stem hot spot (MSHS) areproven high enough that ignition occurs in this hot spotas well, at a time less than half (both for 2D and 3D sim-ulations) the time required for the ignition of the neatmaterial. By studying the maximum percentage burnt inthe two centreline hot spots (CHS) and Mach stem hotspot (MSHS) we observed that the burning in the CHSis, at these timescales, always higher than in the MSHS.However, the maximum burning of the CHS grows as √ x while the maximum burning in the MSHS growslinearly. As a result, the growth rate of the maximumburning in the CHS is higher than in the MSHS at firstbut this is reversed at late times.By comparing two- and three-dimensional simula-tions we identify the change in topology of the hot spotsdue to the third dimension. The faster jet (1.4 timesfaster) in the 3D case results in an earlier collapse ofthe cavity ( ∼ . µ s) and subsequent hot spot genera-tion compared to 2D. The ignition in the BHS in the 3Dcase occurs between 3 . µ s and 4 µ s and in the MSHSbetween 5 . µ s and 6 µ s. E ff ectively the ignition is ob-served earlier in the 3D case by 0 . µ s, which is theamount of time by which the jet impact occurs earlierin 3D compared to 2D. In the 3D scenario the temper-atures can reach values of more than three times higherthan the post-shock temperatures and in the 2D scenario16 a) t = . µ s (b) t = . µ s Figure 16: Comparison of the temperature field and reaction along y = .
21 mm in the inert and reactive configurations (a) y = .
29 mm (b) y = .
35 mm
Figure 17: Comparison of the temperature field and reaction along (a) y = .
29 mm and (b) y = .
35 mm at t = . µ s between the inert andreactive configurations. more than twice the post-shock temperature. This leadsto a higher percentage of maximum immediate burn ofthe material upon collapse in 3D (59%) in 3D comparedto 2D (11%). The growth of the burning follows a simi-lar trend, however, in the two cases.By comparing inert and reacting simulations we con-clude that the e ff ect of the reaction on the topology ofthe hot spots is negligible, whereas a large, additive ef-fect on the temperature field is observed. We examinethe maximum temperature in the centreline hot spots(CHS) and the MSHS in both inert and reactive sce-narios. We demonstrate that after the birth of the hotspot the maximum temperature in the reactive case isincreasing, as expected since the shock waves supportthe reactions and vice versa. In contrary, the maximumtemperature in the inert case is decreasing ad the shockwaves are not supported. An interesting feature ob-served is the superposition of waves emanating from thetwo lobes along the cavity centreline, leading to an ad-ditional short-lived maximum temperature peak in bothcases.The maximum temperatures describing the relative‘strength’ of the CHS and MSHS is studied as well. In the timescales considered, both in the inert and reac-tive scenarios, the CHS always exhibits higher temper-atures than the MSHS. It is interesting to note that thetemperatures of the MSHS in the reactive scenario arecomparable to the temperatures in the CHS in the inertsimulations.The maximum temperature growth of the hot spotsis also studied. Comparing the inert and reactive sim-ulations, we observe that the rate of increase of max-imum nitromethane temperature is positive for the re-active configuration but not in the inert case, for boththe CHS and the MSHS. Comparing the rate of increaseof maximum nitromethane temperature between the twohot spots it is observed that the rate of increase of theMSHS in the inert case is always higher than the CHS(except at the birth moment of the CHS i.e. upon thecollapse of the cavity). In the reactive case this is how-ever, not true. This is likely to have implications in con-figurations with multiple cavities collapsing in reactivemedia.17 cknowledgements This project was kindly funded by ORICA. The au-thors also benefited from conversations with A. Minch-inton, S.K. Chan and I.J. Kirby.
References [1] L. Michael, N. Nikiforakis, Shock Waves (2017)[2] N. Bourne, J. Field, Proceedings of the Royal Society of Lon-don. Series A: Mathematical and Physical Sciences (1894),423 (1991)[3] N. Bourne, J. Field, Philosophical Transactions of the Royal So-ciety of London. Series A: Mathematical, Physical and Engi-neering Sciences (1751), 295 (1999)[4] N. Bourne, A. Milne, in
Proceedings of the Twelfth Symposium(International) on Detonation (2002)[5] N. Bourne, A. Milne, Proceedings of the Royal Society of Lon-don. Series A: Mathematical, Physical and Engineering Sci-ences (2036), 1851 (2003)[6] L. Michael, N. Nikiforakis, in (2014), pp. 60–70[7] A. Kapila, D. Schwendeman, J. Gambino, W. Henshaw, ShockWaves (6), 545 (2015)[8] L. Tran, H. Udaykumar, Journal of propulsion and power (5),959 (2006)[9] N.K. Rai, M.J. Schmidt, H. Udaykumar, Physical Review Fluids (4), 043202 (2017)[10] N.K. Rai, M.J. Schmidt, H. Udaykumar, Physical Review Fluids (4), 043201 (2017)[11] L. Michael, N. Nikiforakis, Journal of Computational Physics , 193 (2016)[12] N. Nikiforakis, J. Clarke, Mathematical and Computer Mod-elling (8), 149 (1996)[13] R. Meniko ff , M. Shaw, Combustion and Flame (2011)[14] C. Tarver, P. Urtiew, Journal of Energetic Materials (4), 299(2010)[15] R. Ripley, F. Zhang, F. Lien, in (2006)[16] A. Stephen, A. She ffi eld, D. Datterbaum, R. Engelke, R. Alcon,B. Crouzet, D. Robbins, D. Stahl, R. Gustavsen, in (2006)[17] M. Arienti, E. Morano, J. Shepherd, Shock and detonation mod-eling with the Mie-Gr¨uneisen equation of state. Tech. rep.,Graduate Aeronautical Laboratories, California Insitute of Tech-nology (2004)[18] D. Dattelbaum, S. She ffi eld, D. Stahl, A. Dattelbaum, M. Elert,M. Furnish, W. Anderson, W. Proud, W. Butler, in Aip Confer-ence Proceedings , vol. 11 (2009), vol. 11[19] K. Bates, N. Nikiforakis, D. Holder, Physics of Fluids ,036101 (2007)[20] S. She ffi eld, R. Engelke, R. Alcon, R. Gustavsen, D. Robins,D. Stahl, H. Stacy, M. Whitehead, in (2002)[21] J. Berke, R. Shaw, D. Tegg, L. Seely, in (1970), p. 168[22] D. Hardesty, Combustion and flame , 229 (1976)[23] R. Chaiken, in Proceedings of the High Dynamic PressureSymposium, Commissariat d’Energie Atomique, Paris, France (1978)[24] R. LeVeque, Computational Methods for Astrophysical FluidFlow, RJ Leveque, D. Mihalas, E. Dorfi, and E. Mueller, 27thSaas-Fee Advanced Course Lecture Notes, Edited by O. Steinerand A. Gautschy, Springer-Verlag (1998)(1978)[24] R. LeVeque, Computational Methods for Astrophysical FluidFlow, RJ Leveque, D. Mihalas, E. Dorfi, and E. Mueller, 27thSaas-Fee Advanced Course Lecture Notes, Edited by O. Steinerand A. Gautschy, Springer-Verlag (1998)