Analysis of the early stages of liquid-water-drop explosion by numerical simulation

The early stages of the shock-driven explosion of liquid-water microdrops are studied numerically with a high-resolution discretization of the axisymmetric Euler equations. A level-set based conservative interface-interaction method is extended to allow for phase transition. The method is applied to a conﬁguration that has been investigated in recent experiments [Stan et al. , Nat. Phys. 12 , 966 (2016); Stan et al. , J. Phys. Chem. Lett. 7 , 2055 (2016)]. The presented results show that the numerical model predicts the initial stages of the violent liquid-drop explosion dynamics accurately. Our results indicate that the deformation of the cylindrical vapor cavity within the droplet is not induced by a torus-shaped negative-pressure wave as was implied from the experimental data. We rather ﬁnd that this torus-shaped wave is a shock wave, and that the observed vapor-cavity deformation is caused by interaction with a negative-pressure region preceding the torus shock. The simulation results show deviations from experimental results at later stages when the drop deformation is dominated by off-center cavitation, i.e., by effects that require extension of the underlying model to take into account generalized nucleation and recondensation.


I. INTRODUCTION
Breakup of liquid droplets plays an important role for a wide range of technological applications, including coating, fuel combustion, or nanoparticle production. During secondary atomization breakup is caused by strong velocity gradients due to ambient gas cross flow, e.g., resulting from a shock wave hitting the droplet. This aerobreakup phenomenon is subject of intense experimental and numerical investigations [1][2][3][4][5].
High-speed liquid-droplet breakup can also be induced by laser irradiation. The fundamental process of the resulting breakup mechanism is still an open problem. In an opaque liquid most of the laser energy is absorbed by a thin layer at the irradiated drop surface. This leads to local evaporation and mass emission. The remaining liquid is accelerated and deformed by a recoil pressure which eventually can cause full droplet breakup [6,7]. For a transparent liquid an optical laser pulse induces a different breakup mechanism. The laser focuses due to optical effects such that the local intensity can be sufficiently high to cause plasma formation [8]. As the plasma inside the drop absorbs the remaining laser energy, explosive breakup can occur, which has been studied experimentally [9,10] and numerically [11,12]. In recent years the availability of high-energy x-ray facilities gave rise to  [13,14]. new experiments. X-ray free-electron laser (XFEL) pulses lead to hydrodynamically simpler but more violent breakup dynamics than optical lasers [13]. To our knowledge, no detailed numerical studies on the underlying droplet-breakup mechanism have been performed yet.
In two recent publications, Stan et al. [13,14] investigated high-speed drop breakup due to energy deposition by a high-energy XFEL. A femtosecond pulse initiates explosion of a water microdrop by homogeneous deposition of energy in the liquid along its trajectory on a time scale much smaller than that of any possible energy-redistribution effects inside the drop. A cylindrical water filament with the diameter of the laser beam is instantaneously heated to a very high and uniformly distributed internal energy density, equivalent to a pressure on the order of 100 GPa. This situation is depicted in Fig. 1. The high pressure in the filament initiates a radially spreading shock wave and a focusing expansion wave at the interface between filament and surrounding liquid which eventually cause the drop to explode.
In the present work we employ a high-resolution numerical model to simulate shock-driven water-drop explosions for the experimental configuration of Stan et al. [13,14]. A finite-volume method using a level-set based conservative interface-interaction method [15] with a fifth-order WENO scheme [16] is extended for phase transitions. For this purpose, a generalized Riemann problem including evaporation across the phase interface is formulated and incorporated into the interface-interaction method. Furthermore, we propose a surrogate model to initiate the rapid formation of the core vapor cavity without need for a nucleation and growth model. We perform oneand two-dimensional axisymmetric simulations to validate the numerical model and to investigate in detail the early stages of droplet explosion, before the onset of significant deviations from axisymmetry.
In Sec. II we introduce the basic equations governing the flow problem, and in Sec. III we present the numerical model. The simulation setup and numerical results for the one-dimensional model problem and the axisymmetric drop simulation are given in Sec. IV together with a comparison of our results to the experiments of Stan et al. [13,14]. Conclusions are given in Sec. V.

II. GOVERNING EQUATIONS
For a perfectly spherical droplet and an XFEL pulse passing through its axis, early stages of the flow evolution can be assumed to be axisymmetric as indicated in Fig. 1. For the dominating rapid inertial processes considered in this paper, capillary and viscous effects can be neglected. To reduce computational cost as compared to the fully three-dimensional treatment, we employ the axisymmetric compressible Euler equations, where r and z denote the radial and axial direction, respectively. U is the vector of conserved variables and F and H are the corresponding radial and axial flux functions, respectively. The transformation to cylindrical coordinates generates the source vector S [17]. These vectors are 044003-2 defined as where ρ is the mass density and p the pressure of the fluid. The velocities in r and z direction are denoted by u and w, respectively. The total energy density E is calculated from the specific internal energy e and the kinetic energy as To close the Euler equations, an equation of state (EOS) is used that relates pressure to density and internal energy. We employ the stiffened-gas EOS [18], where γ is a model constant and p stiff accounts for a pre-compression of the fluid. For liquid water we use γ = 4.4 and p stiff = 0.6 GPa as suggested by Hawker and Ventikos [19], who fitted the parameters to experimental shock data [20] and isotherms [21] for shock problems with pressures up to 100 GPa. The EOS parameters have been calibrated for the positive pressure range, accordingly quantitative deviations from the correct pressure values, e.g., [22], occur in the negative range. In the vapor phase we use γ = 1.33 and p stiff = 0. For the latter, Eq. (4) degenerates to the ideal-gas EOS with the ratio of specific heats γ . Another version of the ideal-gas EOS, relates pressure and density to the vapor temperature with the specific gas constant R = 461.52 J/(kgK). For the liquid water phase, Hawker and Ventikos [19] proposed an approximation of the temperature as Here,R = 1700 J(kg K) is a pseudogas constant and α = −1.8 × 10 −6 J m 10.2 /(kg 4.4 ) denotes a thermal energy factor. Finally, the speed of sound c for the stiffened-gas EOS is given by Since the EOS does not predict cavitation, arbitrarily small densities in the liquid phase are allowed. According to Eq. (4) pressure can reach −γ p stiff . However, Eq. (7) limits the maximum negative pressure to −p stiff . The chosen material model neglects supercritical states and states of ionization and dissociation, which may arise within the filament for a very short period of time at the very beginning of the experiment. We believe that these simplifications imply only small quantitative errors while the overall qualitative behavior is not significantly affected, which is confirmed by our results shown below.

III. METHODS
The cylindrical computational domain is represented by a two-dimensional Cartesian grid with grid spacing Δr and Δz in radial and axial direction, respectively. To obtain a finite-volume discretization of the Euler equations on the grid, the left-hand side of Eq. (1) is integrated over where ∂ i,k denotes the boundary of the cell and n r and n z are the r and z component of the outward-pointing normal vector n, respectively. Using an explicit first-order discretization in time, Eq. (8) results in an evolution equation for U i,k : F n n r + H n n z dr dz.
Here, U i,k denotes the conserved variables averaged over the cell volume V i,k (t ). Higher-order Runge-Kutta schemes are obtained by considering Eq. (9) as building block for the sub-steps. A multiresolution (MR) method is employed to improve computational efficiency. It creates a hierarchy of nested Cartesian grids that represent the flow problem with different resolutions depending on the locally most-demanding scales rather than the globally smallest scales. Discontinuities such as shocks and phase boundaries are always captured on the finest resolution level. For further information on the MR method in general and its present implementation the reader is referred to Harten [23] and Han et al. [24], respectively.

A. Conservative interface-interaction method
At any time t the computational domain is filled with two different fluids separated into two phases by a sharp interface (t ). Cut cells contain segments of the interface and have two fractions each containing one of the fluids. The discretization of a cut cell is illustrated in Fig. 2(a). The volume fraction occupied by the first fluid is given by 0 < α i,k (t ) < 1. Consequently, for the second fluid α i,k (t ) = 1 − α i,k (t ). The lengths of the segments in contact with the first fluid can be expressed as A i±1/2,k (t )Δz and A i,k±1/2 (t )Δr, i.e., by the so-called cell-face apertures A · , · (t ) ∈ [0; 1]. As before, A · , · (t ) = 1 − A · , · (t ) holds for the second fluid, see Fig. 2(a). These segments, however, do not form a closed surface around the volume occupied by each fluid. The segment of the interface within the cut cell, Δ , closes the surface. Its area can be approximated with the cell-face apertures as Applying the cut-cell relations to the finite-volume discretization Eq. (9) gives the discretized evolution of the conserved variables in a cut cell Note that this discretized form holds for each fluid individually with its respective volume fraction and apertures. The last term on the right-hand side expresses the interaction of the two fluids across the interface segment Δ and originates from the closure of the surface integral. X is called interface-exchange term and is determined by an interface-interaction model, which will be discussed later. For the single fluid in a noncut cell, the volume fraction and cell-face apertures are unity and there is no interface exchange. The finite-volume discretization for single-fluid cells thus reduces to

B. Level-set function
The level-set function φ(r, z, t ) describes the signed distance from the phase interface, thus the interface between the two fluids is given by the zero level set of this function, The motion of the interface is tracked via advection of the level-set field by whereû andŵ are the advection velocities in r and z direction, respectively. To maintain the signeddistance property, the level-set function is reinitialized after each advection time step [15,25]. For efficiency, the level-set field φ is tracked only in a narrow band of cells along the interface. As indicated in Fig. 2(b), we consider φ within cut cells and a band of four cells on each side of the interface. In cut cells the interface velocity is determined from an interface Riemann problem [26,27]. In narrow-band cells we set the level-set advection velocity by interpolating the interface velocity from the closest cut cells.

C. Interface interaction
In cut cells, the interaction between both fluids is determined to calculate the interface velocity and the interface-exchange term X required in Eq. (11). The interaction is assumed to be onedimensional and normal to the interface. With respect to the interface-normal vector N, which points from the negative to the positive level-set region, Fig. 2(a), we denote the negative and positive fluid as left and right fluid, respectively. The general form of the exchange term for each fluid is given where + and − apply for the left and right fluid, respectively. The first term on the right-hand side is due to pressure forces acting on the interface and the second term corresponds to mass transfer. The state variables p * L/R , u * L/R , e * L/R and the mass flow rate m are the results of an interfacenormal Riemann problem. Additionally, the interface velocity u I is calculated in the cut cell as advection velocity for the level set. It should be noted that the interactions are formulated normal to the interface. Thus, the fluid velocities are projected onto normal direction and back with the vector If no mass transfer across the phase interface is considered, the solution of the interface Riemann problem is similar to the single-phase solution. The phase interface corresponds to the contact discontinuity, and on both sides either a shock or rarefaction wave propagates away from the interface [17]. If phase transition is possible, however, mass transfer across the interface occurs and needs to be taken into account in the Riemann solver.
In this work we consider evaporation only and develop an extended Riemann solver which employs a transition model proposed by Ytrehus and Østmo [28]. The transition process takes place in a finite vapor layer, called the Knudsen layer, adjacent to the surface of the condensed phase with a thickness on the order of the molecular mean free path [28,29]. Considering vapor to be an ideal gas, Ytrehus and Østmo [28] used kinetic theory and molecular conservation laws to find universal jump conditions relating the saturation vapor state variables ( · ) sat at the liquid temperature T l to the actual vapor state at the edge of the layer ( · ) v . The resulting relations are given in Appendix A.
We approximate the saturation pressure, which is essential for finding the Knudsen-layer solution, by the Clausius-Clapeyron relation as a function of temperature, where L v is the specific latent heat of vaporization, and ( · ) ref denotes a known reference saturation state. We apply the critical point, T crit = 647.096 K and p crit = 22.064 MPa [21], as reference state and use L v = 2.26 kJ/kg [30].
Since the Knudsen layer has molecular length scale, it can be considered as a discontinuity on the continuum scale called Knudsen front. Our modified Riemann solution scheme including this additional discontinuity and its numerical solution by an iterative Newton-Raphson method is detailed in Appendix B.

D. Cavity-formation model
Originally, the level-set approach [31] was used to model the evolution of two-phase flow starting from an initial phase interface, e.g., Refs. [25,[32][33][34]. The development of very small scales, i.e., filaments, bubbles or droplets, close to the interface requires high computational effort. To overcome this problem, these scales are detected and treated in a special way. Depending on the model they are either deleted, converted to Lagrangian particles [35] or transferred to the other phase [36,37]. Additionally, a subgrid representation of the level-set field might be used [38,39]. In contrast to these approaches, we need to capture the development of small-scale two-phase regimes remote from the existing phase interface. Instead of resolving these small scales exactly, we introduce a model motivated from the considered flow problem and utilizing the level-set concept.
Due to cylindrical focusing, the strength of the inwards-running rarefaction wave, see Fig. 1, increases with decreasing radius. As a consequence, tension tends to its theoretical maximum at the The Knudsen-front model allows evaporation only across the liquid-vapor interface, i.e., no new vapor regions can emerge spontaneously within the bulk liquid phase. To avoid the need for a full nucleation model we propose a surrogate model for the formation of the core vapor cavity. It is formulated as an interface interaction and is linked to a pre-initialized pseudointerface on the axis, where cavity formation is to be expected for the considered problem. We assume the onset of cavity formation as soon as pressure drops below zero at the pseudointerface.
Numerically, the initial pseudointerface is located in the cavity initialization zone, which is formed by the second and third cell layer next to the axis, Fig. 3. It may not be located within the first layer as this would require the definition of separate numerical level-set boundary conditions. Thus, a small area of initially inactive vapor between the boundary and the interface exists. As long as cavity formation has not started, the interface-exchange term X at the pseudointerface models the symmetry boundary condition for the liquid phase which is not in contact with the actual boundary, where u i−1/2,k is linearly interpolated from the cell velocity and the prescribed zero velocity at r = 0. The full algorithm for cut-cell treatment is outlined in Appendix C. It requires a modified Knudsen-front model to calculate a solution that solely depends on the current liquid state for cut cells within the cavity initialization zone. To find this solution, we assume that the liquid is compressed to its corresponding saturation pressure, p L = p sat (T L ), and that evaporation is critical, M v = 1; see Appendix A. Outside of the cavity initialization zone, i.e., after the pseudointerface has turned into a real physical liquid-vapor interface, it is advected according to the Knudsen-front interface-interaction method.

E. Numerical schemes
The intercell fluxes F i±1/2,k and H i,k±1/2 are approximated with a fifth-order weighted essentially non-oscillatory (WENO) scheme [16] using characteristic decomposition with local Lax-Friedrichs flux splitting [40] based on Roe averages [32,41]. The states of both fluids are extended across the interface to fill so-called ghost cells [27]. This extension allows the use of the same high-order WENO stencil in the vicinity of the phase interface. Both the flow variables and the level-set function are propagated in time using a strongly stable explicit second-order Runge-Kutta scheme [42]. The time-step size Δt is chosen according to the Courant-Friedrichs-Lewy (CFL) condition whereû andŵ are nonzero in narrow-band cells only and C CFL = 0.6. As the CFL condition might be violated in cut cells with small volume fractions, a mixing procedure [15,43] is employed for stabilization.

IV. NUMERICAL SIMULATIONS
The water drop initially is perfectly spherical and centered on the laser beam. Figure 4 shows the two-dimensional computational domain. Note that we impose rotational and axial symmetry with respect to the z and r axis, respectively. Besides the symmetry conditions, zero-gradient boundary conditions are applied at the remaining far-field boundaries. The cavity-formation model is activated at the drop axis. The extent of the domain is three times the drop radius in both radial and axial direction. This ensures that boundary effects do not distort the flow field in the vicinity of the drop during the runtime of our simulations. Stan et al. [13,14] used drop diameters ranging from 32 μm to 71 μm. For the present simulations, a diameter of 40 μm is chosen. The radius of the heated filament is set as r f = 0.5 μm in accordance with the experimental parameters.
The computational domain is discretized with a Cartesian grid with equal spacing in radial and axial direction. As a MR method is applied, the actual resolution is not uniform across the domain and changes with time. The theoretical maximum resolution is 2048 cells per dimension at the highest MR level and thus the smallest possible grid size is Δr min = Δz min ≈ 2.93 × 10 −8 m. The chosen resolution was checked for convergence by comparison of the results with a higher resolution. At the given mesh size, all relevant flow structures were already present and well resolved.
The level-set function is initialized to represent the liquid drop. For the pseudointerface of the cavity-formation model, an additional zero level is initialized in the second cell layer at r = 4 × 10 −8 m for 0 z R. Although slightly larger r do not significantly affect cavity dynamics, it should be chosen as small as possible to minimize modeling errors. 044003-8 Since we employ a two-phase approach, we assume the ambient gas to be atmospheric vapor instead of pure air. We have performed simulations with different initial filament pressures between 50 and 400 GPa. As we observed only differences in the time scaling but not of the overall wave and drop-breakup dynamics, however, we restrict the present study to initial filament pressures of p f = 100 GPa. The initial conditions of the presented simulations are summarized in Table I.

A. Model analysis
In the following a one-dimensional model problem is considered to demonstrate the effects of the cavity-formation model. As depicted in Fig. 4 the one-dimensional domain is equivalent to the radial axis of the two-dimensional version. It represents an infinite water cylinder with a heated filament at the axis. The model problem was simulated both with and without cavity-formation model to analyze its influence.
In Fig. 5(a) the radial pressure profile is plotted for both simulations at time t = 8.1 ns, which is approximately 0.05 ns before the pressure wave hits the drop interface for the first time. The shock front is identical for both cases since it propagates prior to cavity formation. Without cavityformation model, an expansion in the post-shock fluid region decreases the pressure in large parts of the drop below zero. In contrast, the pressure in the bulk fluid is mostly positive in case of cavity formation. This shows that the cavity formation triggered by negative pressure in the liquid counteracts the tension and results in positive pressure in the remaining nonevaporated water.
The evolution of the phase interface locations is presented in Fig. 6(a). The outer interface is accelerated away from the axis as the shock wave impinges. The shock is reflected as negativepressure wave propagating towards the axis; see Fig. 5(b). As soon as the wave is fully developed and detached from the the phase interface, the latter is affected by the pressure in the bulk fluid. Therefore, the interface moves inwards again in the case without vapor cavity whereas its position remains almost constant in the other case. Figure 6(a) shows that the inner cavity starts to grow shortly after the simulation is started. Each time the negative pressure front hits the cavity interface it is accelerated towards the liquid side, e.g., at roughly 20 ns. Upper and lower plot of Fig. 6(a) illustrate that the pressure wave travels between the inner and the outer interface and pushes the latter more and more away from the axis. In contrast, if no vapor cavity is considered, the interface oscillates between its initial position and a small displacement such that there is no residual increase of the cylinder radius; see the dashed curve in Fig. 6(a). From these findings we can conclude that the formation of an inner cavity is essential to allow an expansion of the liquid cylinder away from its axis. In that case the liquid forms a hollow cylinder which keeps expanding as illustrated by the long-term evolution depicted in Fig. 6(b).

B. Simulation results
A main difference between the full axisymmetric simulation and the simplified one-dimensional model problem is the radial spreading of the spherical wave dynamics causing focusing and attenuation. The following results show that this leads to a significantly different evolution.
Note, in all following contour plots we have scaled up the density gradient in the vapor phase by a factor of 50 to facilitate legibility. Furthermore, the color maps of pressure do not reflect the actual minimum and maximum values but are chosen to best illustrate all relevant features.
At the beginning of the simulation the ends of the heated filament are in contact with the drop surface. There the interface is strongly accelerated in axial direction and the water reaches negative pressures much faster than in the remaining filament. Figure 7 shows a close-up view of that region. The main shock emerges from the edge of the filament and spreads radially. A secondary shock wave forms at the expanding tip of the filament and propagates in axial direction towards the equatorial plane; see Fig. 8. This shock wave is termed axial shock in the following.
While the main shock wave propagates in radial direction, it impinges obliquely on the drop interface and deflects it outwards with a kink at the impingement point [44]. First, the main shock is nearly perpendicular to the drop surface and an anomalous reflection (AR) occurs [45]. It causes the shock to bend in the vicinity of the interface; see Fig. 7. With increasing distance from the axis the angle of incidence decreases and the shock remains straight during reflection, Fig. 8(a). Due to the interaction with the interface, the main shock is reflected as a Prandtl-Meyer expansion fan (PME) [46]. Another expansion wave connects the PME with the deflected phase interface. Figure  8  shows the drop shortly after the main shock wave has reached the drop apex. The curved negative-pressure wave resulting from the refraction collapses towards the equatorial plane, which increases the tension. This is termed spherical focusing in contrast to cylindrical focusing. According to theoretical considerations based on optics by Stan et al. [14] the focal point is located at r = R 2 , which is indicated in Fig. 8(c). When the axial shocks of both ends of the drop collide in the equatorial plane (at around t = 16.1 ns), they coalesce into a stronger shock wave propagating in opposite direction; see Figs. 9(a) and 9(b). This wave has a toroidal shape and is termed torus shock in the following. Figure 10 shows the locations of the inner and outer interface along the radial axis. The kink at t = 17.6 ns in the right plot indicates that the torus shock hits the drop surface accelerating it outwards.
In the vicinity of the radial axis a negative pressure region is present between the cavity interface and the torus shock; see Fig. 9(a). Due to this negative pressure the inner phase interface is bent towards the liquid as indicated in Fig. 9(b) (red arrow) and Fig. 10(a). As soon as the torus shock hits the cavity boundary it is accelerated towards the drop axis. This is evident from the kink in Fig.  10(a) at t = 19.8 ns. Since the interface already has small-scale perturbations prior to the impact, the shock impingement causes Richtmyer-Meshkov instabilities to develop [47,48].
While the remainder of the torus shock propagates further towards the ends of the drop it interacts with both the inner and the outer phase interface. Subsequently, several additional reflections occur. The liquid is stretched in axial direction, which is evident from comparing Figs. 9(b) and 9(c). However, due to the increasing distortion of the cavity boundary, the wave reflection patterns become increasingly chaotic, Fig. 9(c), and there is no bulk movement of the remaining drop observable anymore. Instead, instabilities grow along the entire interface. As these instabilities in reality are no longer axisymmetric, our assumption of a two-dimensional axisymmetric evolution ceases to be valid at later times.

C. Comparison with experimental results
In the experiments by Stan et al. [13,14] water drops were subjected to XFEL pulses. Shockdriven explosions of the drops occurred due to the energy deposited by the laser pulse. Figure 11 reproduces experimental results of Stan et al. [14]. In their photographs, the spherical shape of the water drop optically magnifies the interior of the drop but obscures the outer region. Since the energy-pressure relation applied through the EOS, Eq. (4), is not strictly valid throughout the experimental fluid states and the initial energy transfer to the fluid cannot be directly considered by the employed continuum model, it is not possible to match directly a specific initial pulse energy to the simulation settings. In Fig. 11(a) the temporal evolution of the drop diameter in the equatorial plane is shown. The first displacement upon shock arrival and the subsequent evolution is qualitatively in very good agreement with the results from the present simulations, compare Figs. 6(a) and 10(b). The experimental data show that the diameter appears to exhibit a transient reduction before a subsequent increase. A similar although more pronounced evolution can be seen in the simulation results. This observation is discussed in the following. The negative-pressure wave resulting from the reflection of the initial shock is apparent both in the simulation and in the experiment, Fig. 11(b) (third picture, 15.5 ns). Due to spherical focusing, the local tension in the liquid can become strong enough for it to rupture. Stan et al. [14] observed this cavitation, which they name spallation, close to the focal point as bright and dark spots, Fig. 11  The torus shock that appeared in the present simulations looks similar to a pressure wave observed in the experiment, compare Figs. 11(c) and 11(d). Stan et al. [14] found that the boundary of the cavity bends towards this wave, see blue arrow in Fig. 11(c), and thus concluded that the entire wave including the visible front is a negative-pressure wave [14]. Interestingly, the same behavior of the cavity interface is observed in our results although the wave front belongs to a shock. Note that in the experimental images the visible deformation is enlarged due to optical distortion. As discussed above, the deformation is not caused by the torus wave itself but by the negative-pressure region preceding it, which is in accordance with the situation depicted in Fig. 11(c). Moreover, we believe that a slight bending of the cavity interface in opposite direction in Fig. 11(c) (fourth picture, 53.25 ns) and in Video 2 of Ref. [14] indicates the impact of the torus shock on the cavity. Figure 12 shows the deformation of the inner cavity right before and after the torus shock impinges in a three-dimensional rendering of the drop interface. Figure 12(a) illustrates once more the similarity to Fig. 11(c), while the inflection in Fig. 12(b) resembles the last panel of Fig. 11(b).

V. CONCLUSION
We extended an existing two-phase model based on a conservative interface-interaction method with a Knudsen-front Riemann solver to allow for phase transition through rapid evaporation. Additionally, we introduced a surrogate model that allows the creation of vapor cavities at the axis of axisymmetric flow problems. This model is physically-motivated by the strong tension that can occur in liquids around the axis due to radial focusing of expansion waves.
We applied our new model to the simulation of a laser-induced shock-driven explosion of an isolated liquid drop and compare our results with observations from a recent experiment [14]. Our simulated drop dynamics agree qualitatively very well with the experimental results. The formation and growth of the cylindrical vapor cavity inside the drop cause a macroscopic deformation of the liquid drop as in the experiments. Our results show that the cavity formation is essential to obtain a considerable radial movement. Stan et al. [14] conclude from optical measurements that a reflected negative-pressure wave causes a deformation of the inner cavity. Interestingly, we observe the same deformation but related to a shock wave with similar shape as the negative-pressure wave which was suggested from the experimental observations. In our simulation, the interface deformation is actually caused by a negative-pressure region preceding the shock. At later times, the explosive expansion of the drops in the experiment is driven by additional cavitation effects. These effects  are not captured by the simulation which does not include a nucleation model. Therefore, the experimental drop dynamics at later stages are not reproduced by the simulation.
Further research will concentrate on improving the formation of vapor cavities to allow cavitation to occur at locations away from the drop axis. We expect to extend the good agreement between simulation and experiment to late stages of drop deformation. Furthermore, the axisymmetric formulation will be replaced by a fully three-dimensional model to allow for symmetry breaking at later stages. to an intermediate state (L ), from which the liquid evaporates across the Knudsen front. The initial vapor state on the right (R) is also either compressed or expanded, such that the resulting state ( * R) is in pressure and velocity equilibrium with the newly created vapor state ( * L) across the contact discontinuity.
Applying the Knudsen-layer relations and momentum conservation across the Knudsen front as well as continuity of pressure and velocity across the contact discontinuity, the Riemann solution is fully determined. We use an iterative Newton-Raphson method to obtain the exact solution, which is detailed below.

Exact iterative solver
An exact Riemann solver utilizes the fact that velocity u * and pressure p * are constant across the contact discontinuity of the Riemann solution; see Fig. 13. Thus, we can find p * as the root of the function f (p * ) = u * L (p * ) − u * R (p * ) ≡ 0, where u * L and u * R indicate the velocity u * calculated with respect to the left and right state, respectively. Using a Newton-Raphson method to find the root, the pressure p * is iterated by where n is the number of iterations and f the derivative of Eq. (B1) with respect to p * . The iteration is stopped as soon as the relative correction reaches a given threshold, with = 10 −6 . The velocity u * R is calculated from the right state by considering the change in velocity due to the right shock or expansion wave, f R (p * ), as u * R = u R + f R (p * ). (B4) The function f R and its left counterpart f L are given later. The state change from ( · ) L to ( · ) * L is more complicated due to the additional intermediate state (L ). Its velocity is determined by considering the velocity change across the leftmost wave, f L (p L ). The additional change across the Knudsen front can be calculated from mass conservation and the mass transfer rate per unit area m to yield,

044003-17
Our nested solution procedure is initialized with (p * ) 0 = p R , (p L ) 0,0 = p sat (T L ) (B10) to start the outer iteration. In case of divergence or invalid converged solution, e.g., M v > 1, our algorithm falls back to a nonevaporative solution, meaning that the Knudsen front vanishes and (L ) coincides with ( * L).

APPENDIX C: ALGORITHM FOR CAVITY FORMATION
In the following we present the core algorithm of our cavity-formation model. It requires a variable for each axial cell index k, TAG k , which initially is set to [no cavity]. The radial location of the initial pseudointerface is denoted as r 0 interface and equals 4 × 10 −8 m in our simulation setup.