On the theory and numerical modeling of spall fracture in pure liquids
aa r X i v : . [ phy s i c s . p l a s m - ph ] F e b On the theory and numerical modeling of spall fracture in pure liquids
Mikhail M. Basko a) Keldysh Institute of Applied Mathematics, Miusskaya square 4, 125047 Moscow,Russia (Dated: 3 February 2021)
A novel phase-flip model is proposed for thermodynamically consistent and computationally efficient descrip-tion of spallation and cavitation in pure liquids within the framework of ideal hydrodynamics. Aiming atultra-fast dynamic loads, the spall failure of a liquid under tension is approximated as an instantaneousdecomposition of metastable states upon reaching the spinodal stability limit of an appropriate two-phaseliquid-gas equation of state. The spall energy dissipation occurs as entropy jumps in two types of discontinu-ous solutions, namely, in hypersonic spall fronts and in pull-back compression shocks. Practical application ofthe proposed model is illustrated with numerical simulations and a detailed analysis of a particular problemof symmetric plate impact. The numerical results are found to be in good agreement with the previously pub-lished molecular-dynamics simulations. Also, new approximate nonlinear formulae are derived for evaluationof the strain rate, fractured mass, and spall strength in terms of the observed variation of the free-surfacevelocity. The new formula for the spall strength clarifies complex interplay of the three first-order nonlinearcorrection terms and establishes a universal value of the correction factor for attenuation of the spall pulse,which in the limit of weak initial loads is independent of the equation of state.Keywords: spallation and cavitation in liquids, fluid dynamics with phase transitions, metastable liquids,spall strength, post-acoustic approximation
I. INTRODUCTION
For many decades, investigation of spall failure inmatter under dynamic tension was conducted mainlyfor solids . In recent years, measurements of thespall strength in liquid metals, loaded by the tradi-tional method of impactor plates, have been reported .Also, many new experiments have been performed wherespall failure occurs in liquids, often in molten metals, atvery high loading rates under the action of nanosecondand sub-picosecond laser pulses , or very short X-raybursts . Adequate modeling of dynamic spall fracture inliquids is a challenging problem, especially when the ge-ometry of the experiment is not one-dimensional and/orspall fracture is only part of a more complex problemas, for example, by modeling liquid tin targets for EUVlithography sources .A well-established approach is to use the first-principlemolecular dynamics (MD) simulations , which pro-vide an important insight and allow calculation of theultimate spall strength in liquids at high strain ratesbut are limited to sub-micron targets, are computation-ally very costly and, for that reason, not quite practi-cal. A much simpler and computationally inexpensivealternative is to simply adopt the fully equilibrium (EQ)equation of state (EOS), obtained by applying Maxwell’srule to the primitive EOS of the van-der-Waals (vdW)type . However, because the EQ EOS allows nonegative pressures, it cannot describe tensile states, pull-back shocks, and only approximately reproduces massdistribution in the regions of spallation and cavitation . a) Another often used model is based on combining anappropriate empirical EOS that admits negative pres-sures with an ad hoc limit on the tensile stresses: onceand where the spall threshold is exceeded, an expand-ing void cavity with zero-pressure boundaries is intro-duced into the simulation. The principal flaw of thisapproach is its inability to adequately account for theliquid-vapor phase transition, especially when spall fail-ure occurs practically simultaneously over a large quasi-uniform mass of metastable liquid, converting it into afinely dispersed liquid-vapor mixture.In the present work, we investigate a novel phase-flip(PF) model based on the premise that all the informa-tion, needed for modeling cavitation and spallation inpure liquids, is provided by a two-phase liquid-gas vdW-type EOS, and no other material properties such as sur-face tension, viscosity, etc. are to be invoked. This isa clear aspect where our model differs from the energy-based scheme by Grady , preserving nonetheless thecapability to adequately account for the energy dissipa-tion by spallation. Aiming at description of ultra-fastdynamic processes, we assume the spall strength to at-tain its maximum possible value, i.e. to be defined bythe thermodynamic stability limit of the vdW-type EOSalong the spinodal curve. The spall fracture is then natu-rally associated with the spinodal decomposition, which,in line with the adopted minimalistic approach, is as-sumed to take place instantly. To justify the latter as-sumption, we invoke the theory of homogeneous bubblenucleation in metastable liquids.The primary motivation for proffering the PF model isits potential to be easily and computationally cheaplyintegrated into large multi-dimensional fluid dynamicscodes that are already loaded with many other complexphysical processes like spectral radiation transfer, laserenergy deposition, etc. when need arises to describe tar-get fragmentation by spallation and cavitation under fastand intense initial loads. A notable example is the task toadequately describe the initial shattering of small liquid-tin drops by picosecond laser pulses in EUV-lithographyapplications , where MD codes cannot be directly ap-plied, and where the PF model promises to become amajor improvement over the more primitive models em-ployed in Refs. 9 and 10.The full formulation and justification of the PF modelis given in Sect. II, where we start by introducing the em-ployed vdW-type EOS, describe and analyze the adoptedspall criterion, and explore the kinematics of hypersonicspall fronts. To illustrate how the proposed model couldbe applied in practice, we perform a detailed analysis ofthe spallation process in a planar collision between identi-cal flyer and target plates. Although real experiments areusually done for asymmetric configurations with target-to-flyer thickness ratios about 3–5 , we choose the theo-retically simplest symmetric case to display the key prop-erties of our model, and to compare it with the avail-able MD simulations. Our numerical results, presentedin Sect. III, indicate that the PF model can reproducethe key aspects of the spall dynamics no less accuratelythan the MD simulations — provided that an adequatetwo-phase EOS is available. The difference emerges onlyin the microstructure of fracture zones, stipulated by sta-tistical atomic-scale properties of an MD model, and de-termined by small-scale perturbations in the PF model.Interpretation of plate impact experiments for deter-mination of the spall strength σ sp relies on the theoreticalrelationship between σ sp and the observed variation ∆ u fs of the free-surface velocity (the “velocity pullback”), forwhich usually the linear acoustic approximation is used.In Sect. IV we develop a new variant of the post-acousticapproximation and derive a new formula relating σ sp to∆ u fs in the case of a symmetric plate impact. The newformula is significantly more accurate than the acousticone, and provides an important insight into the complexinterplay of the first-order nonlinear correction terms.New post-acoustic formulae are also obtained for thestrain rate and the fractured mass fraction. II. THE PHASE-FLIP MODEL OF SPALL FRACTUREIN LIQUIDSA. The employed two-phase EOS
Conceptually, the thermodynamics of a liquid-gasphase transition is adequately represented by the clas-sical van der Waals EOS. In this work we make use ofone of the simplest generalizations of this EOS (termed GWEOS) in the form p ( v, θ ) = αθv − κ − − κv n , (1) e ( v, θ ) = αc V θ − κ ( κ − v − n , (2) s ( v, θ ) = α (cid:2) c V (1 + ln θ ) + ln( v − κ − ) (cid:3) , (3)where κ = ( n + 1) / ( n − α = κ − κ − , and the expo-nent n > c V > p , thespecific volume v ≡ ρ − and the temperature θ are nor-malized to the corresponding values P cr , V cr = ρ − cr and T cr at the liquid-gas critical point; the mass-specific in-ternal energy e is in units of P cr V cr , the entropy s in unitsof P cr V cr /T cr . The dimensional quantities P cr , V cr , T cr make up additional three free parameters of GWEOS;for more details see Ref. 25. The main motivation for us-ing GWEOS in this work is its mathematical simplicity,which allows the in-line use of its both the metastable(MS) and the equilibrium (EQ) branches in hydrody-namic codes with the rounding-error accuracy — thusexcluding numerical errors due to imperfect EOS.In numerical examples, discussed below, the few freeparameters of GWEOS were chosen such as to providethe best fit to the properties of water near normal condi-tions. The well known critical parameters of water allowconfident normalization of its other measured quantitiesto the system of reduced variables. For simplicity, weassume the normal state ( v , θ ) to have zero pressure p ( v , θ ) = 0. The values of n , c V and θ are adjustedto fit as best as possible the experimental values of thenormal specific volume v , the isentropic sound speed c ,the vaporization enthalpy h vap , and the specific heat αc V at T = 20 ◦ C. Table I shows how good a compromise hasbeen achieved with the values n = 1 . , c V = 8 . , θ = 0 . , (4)used throughout this paper. TABLE I. Comparison between the normal properties ofwater and those produced by GWEOS with parameters (4).For all quantities the reduced values are given θ ρ c h vap αc V b GWEOS 0.450 3.9073 5.9939 21.5 39.26 3.17H O 0.453 3.10 5.66 35.8 39.26 2.0
As is shown in Fig. 1, Eq. (1) exhibits all the standardfeatures of the vdW EOS. The phase coexistence regionlies between the binodal curve bi and the p = 0 line.The EQ EOS in this region, denoted p EQ ( v, θ ), e EQ ( v, θ ), s EQ ( v, θ ), and the position of the binodal itself are de-termined by applying the Maxwell rule ( §§ θ = 0 .
93 isotherm in Fig. 1. The spin-odal curve sp , defined by the condition ∂p ( v, θ ) /∂v = 0,lies under the binodal and delimits the region of ab-solute thermodynamic instability; hence, only the EQstates with p >
0, obtained by applying the Maxwellrule, are possible below the spinodal. In the intermediatemetastable region between the binodal and the spinodalboth the MS and the EQ branches of EOS are thermody-namically admissible and compatible with the equationsof fluid dynamics. The latter is ensured by the positive-ness of the square of the isentropic sound speed c = (cid:18) ∂p∂ρ (cid:19) s > p < v <
1. Note that, byassuming p ( v , θ ) = 0, we put the normal state ( v , θ )into the category of metastable ones. GFBO SO h sp biP E spp CP E Q i s en t r ope MS isentropeEQ isotherm v FIG. 1. Thermodynamic ( v, p ) plane of a fluid with a liquid-gas phase transition. Shown are the binodal bi and the spin-odal sp curves with the critical point CP at ( v, p ) = (1 , P BO S shows an MS isentrope, along which therelease-wave states would evolve after being brought to P along the Hugoniot OP (dashed) starting from the normalstate O . Dash-dotted curve h sp (orange) is the locus of iso-choric post-flip states E as the initial state S slides along thespinodal sp ; dash-dotted curve F EG (red) is the EQ isentropepassing through the chosen state E . The behavior of the MS(solid green) and EQ (dashed green) isotherms is illustratedfor θ = 0 . A non-trivial point is the interpretation of the EQstates from the viewpoint of matter structure. Whilethere is no ambiguity about the metastable liquids, whichcan always be assumed homogeneous, the EQ states arenormally believed to be necessarily heterogenous ( § §
64 in ]. Below we apply the term “fog”to such EQ states and treat them as perfectly homoge-neous in the context of spallation and cavitation in pureliquids. B. Spall criterion and the phase-flip approximation
Spall failure inside (i.e. away from boundaries) a vol-ume of pure liquid occurs in the process of sponta-neous decay of its metastable state due to thermody-namic fluctuations . This implies that, other factorsbeing equal, the spall strength σ sp must increase withthe increasing strain rate, which is confirmed by manyexperiments and the MD simulations . Hence, thetheoretical maximum of σ sp — the ultimate strength —must be attained at the spinodal σ sp = − p sp ( v ) = n − v ( n + 1) v n +1 , n − n + 1 < v < nn + 1 , (6)where the metastability decay rate is maximum .In this work we adopt a simple and universal criterionthat a given fluid element undergoes spall failure when-ever its thermodynamic trajectory reaches the spinodalfrom the side of metastable liquid. The results obtainedunder this assumption should be in the first place applica-ble to tensile loads with the highest possible strain rates,like, for example, ˙ ε ≃ –10 s − observed in laser ex-periments. In reality, of course, the tensile strength canonly approach the spinodal limit from below , but howclosely — about 0.3–0.5 or up to 0.8–0.9 of the lim-iting value — is still debatable. Given this uncertaintyand the uncertainty in the spinodal position for real ma-terials, we do not try to soften this criterion or make itdependent on ˙ ε , the more so that it would not affect themain results of this work.Having set the pressure threshold for spallation, weneed to add further assumptions about the kinetics andmorphology of the ensuing fracture. Because the rate ofspontaneous bubble nucleation very rapidly increases bymany orders of magnitude as the MS liquid approachesthe spinodal , we can make a simplifying assumptionthat spall failure occurs as an instantaneous and irre-versible MS → EQ phase jump — the phase flip . Theextremely short timescale, on the order of a few picosec-onds, of spinodal decomposition in real liquids justifiesthis approximation for strain rates up to ˙ ε ≃ s − .Because the rate of density variation ˙ ε is controlled byhydrodynamics (i.e. by fluid inertia and pressure gradi-ents), the phase flip must occur at a constant ρ = v − whenever ˙ ε is finite. The only possibility where it couldbe accompanied by a density jump, would be inside ararefaction shock , which can be excluded in the con-text of spallation because relaxation of tension meansgrowth of pressure. In addition, energy conservation re-quires the local value of the specific internal energy e tobe preserved as well, which implies that the phase flip isalways accompanied by a jump-like increase in pressure,temperature, and entropy.In Fig. 1 the trajectory P BO SE illustrates possibleevolution of the thermodynamic state of a fluid element,undergoing spall failure: expansion along the MS isen-trope P BO S brings it to a negative-pressure state S onthe spinodal, followed by an instantaneous jump to theEQ “fog” state E at a higher entropy; its subsequent evo-lution is governed by the EQ EOS branch so long as itremains under the binodal. If the tracked mass of “fog”emerges from under the binodal (say, by recompressionalong the EQ isentrope EF ), it may later return to thetwo-phase region again as metastable liquid and undergoa secondary spall failure. Thus, the fluid elements in ourscheme can exhibit a hysteresis-like thermodynamic be-havior that is fully consistent with the first and secondlaws of thermodynamics.Because metastable liquid can border vacuum at a fi-nite density, one might conjecture that it would ruptureby forming an expanding vacuum cavity, originating fromthe inception point where the spall threshold is first at-tained. To assess the feasibility of such scenario, we haveto consider the compression shock pulse that would nor-mally be launched into the surrounding stretched liquidby pressure relaxation in the fractured material (oftencalled the “spall pulse”). Because this shock relaxes thetension in metastable liquid and pulls its thermodynamicstate back from the spall threshold, we call it the pull-back (pb) shock . For an expanding vacuum cavity toform, the pb shock would have to start immediately fromthe inception point. However, a closer examination inthe next section reveals a different general pattern: oncethe spall threshold is reached at the inception point, theregion of spall failure very rapidly (hypersonically) ex-tends over some finite volume, bounded by a surface fromwhere the pb shock does actually start, and which maybe called the “spall horizon” of a given inception point.Due to the continued stretching, sustained by inertia ofthe fluid motion, no fluid element inside the spall horizoncan be reached by the pb shock before it undergoes spallfailure. In other words, the spall fracture in our modeldevelops as a rapid phase transition between the MS andEQ EOS branches, occurring practically simultaneouslythroughout a certain finite mass. Since an EQ state canborder vacuum only at zero density, there is no need topostulate any void openings, and the spalling (cavitating)liquid can be treated as everywhere continuous.At first glance, the above spall criterion appears similarto that proposed earlier by Grady , who also assumed an instantaneous MS → EQ phase jump in the frameworkof a vdW-like EOS. The essential point of difference is,however, that Grady postulated the post-jump EQ state E to have the same entropy as the pre-jump MS state S . As the EQ isentrope, branching from the same point B on the binodal, has always lower values of the inter-nal energy e EQ ( v, s ) < e ( v, s ) than the MS isentrope atthe same density, Grady interpreted the excess energy∆ e NE = e ( v, s ) − e EQ ( v, s ) as the interfacial energy (ex-traneous to EOS) of the heterogeneous EQ state, consist-ing of liquid fragments with a characteristic size l fr ; hethen used ∆ e NE and l fr to evaluate the spall strengthby invoking the surface tension and viscosity. Our model,in contrast, is based on the assumption that l fr is verysmall and cannot be resolved hydrodynamically. Thenthe eventual interfacial energy must be automatically in-cluded into the fluid EOS, which requires the MS → EQphase jump to be iso-energetic and accompanied by theentropy increase.Another highly controversial point in Grady’s argu-mentation is the postulated linear dependence l fr = 2 ct w (7)of the typical fragment size l fr on the “waiting” time t w ∝ ˙ ε − elapsed after crossing the binodal. Here we ar-gue that such a linear relation would be in stark disagree-ment with the theory of homogeneous nucleation ,which predicts that the new phase in metastable liquidsunder tension − p > r cr = 2 α st p sat − p ≈ − α st p , (8)where α st is the surface tension, and p sat = p sat ( T ) > T for which p is calculated . Because the probabilityof spontaneous creation of such bubbles, proportional toexp( − G ) where G = 16 πα st p sat − p ) T = 4 πα st r cr T ≫ p and T ,the spall fracture of a uniformly stretched liquid occursas a sudden emergence of myriad critical vapor bubbleswith radii of a few nanometers , which practically in-stantaneously (i.e. with virtually no time for subsequentgrowth) begin to coalesce, breaking the molecular bondsthat sustain tension. Within this picture, it is expectedthat l fr ≈ r cr and t w ∝ exp( G ), which, together withEq. (9), implies a highly nonlinear relation between l fr and t w ∝ ˙ ε − , and significantly smaller values of l fr thanpredicted by Eq. (7). Therefore, when the goal is to ex-plore the dependence of σ sp on the strain rate ˙ ε , thesimple relationship (7) must be replaced by a more ad-equate model based on the nucleation theory — as, forexample, was proposed in .In summary, unlike Grady’s scheme, our model de-scribes energy dissipation due to spall fracture withoutinvolving the surface tension and viscosity, simply as theentropy increase in the phase-flip transitions and behindthe pull-back shock fronts. In this respect it is similar tothe viscous-free ideal hydrodynamics with dissipation dueonly to discontinuous shocks. The weak dependence ofthe spall strength on the strain rate is ignored by adopt-ing the theoretically maximum value of σ sp allowed bythermodynamics. C. Kinematics of spall fronts and the size of fracture zone
Consider a flow of liquid with a velocity field u ( t, x ),where a certain region is stretched to develop negativepressures. The local rate of stretching is characterizedby the strain rate , defined as˙ ε def = − d ln ρdt = d ln vdt = ∇ · u , (10)where d/dt = ∂/∂t + u · ∇ is the material time derivative,and the use is made of the fluid continuity equation ∂ρ∂t + ∇ · ( ρ u ) = 0 . (11)Next, we assume for simplicity that the flow is isen-tropic with a fixed value of entropy s throughout theconsidered volume. Then the pressure becomes a knownmonotonic function p ( v ) = p ( v, s ) of a single variable v with dp/dv < P BO S in Fig. 1); the spall strength σ sp = − p ( v sp ) andthe spall density ρ sp = v − sp take on universal values thatare determined by the intersection of the isentrope p ( v )with the liquid branch of the spinodal (point S in Fig. 1). xt > t sp x = u spf t x+ xCA Bp sp p x t sp t + tfog bubble FIG. 2. Schematics of the spall front propagation. As thelocal pressure minimum is forced to drop ever deeper belowthe spall threshold p sp , the position of the spall front shiftsfrom the inception point A at t sp to points B and C at latertimes t and t + ∆ t Spall failure is initiated at time t sp when the absoluteminimum in the spatial pressure distribution p ( t sp , x ) at-tains the threshold value p sp = − σ sp at some inceptionpoint A , as is illustrated in Fig. 2. Generally, when p ( t, x )and p ( v ) are differentiable at A , one has ∇ p = ∇ ρ = 0at this point. At later times t > t sp the condition p ( t, x ) = p sp will be satisfied on a certain surface aroundthe inception point, which we call the spall front andwhich encompasses a fog bubble. As suggested by theterm, all the liquid inside the fog bubble has been, uponpassing through the spall front, instantaneously trans-formed into an EQ “fog” state at a positive pressure.The condition ∇ p = 0 by inception implies that the ex-pansion of the fog bubble starts infinitely fast, but grad-ually slows down as |∇ p | increases. Because p is a single-valued monotonically growing function of ρ , the normal(i.e. along the colinear ∇ p and ∇ ρ vectors) kinematicvelocity of the spall front u spf = u · ∇ ρ |∇ ρ | + ˙ ε |∇ ln ρ | (12)is readily calculated by differentiating the equation ρ ( t, x ) = constant , as is illustrated in Fig. 2. Note that inthe comoving reference frame, where u = 0, the front ve-locity u spf is given by the second term on the right-handside of Eq. (12).The merely kinematic propagation of the spall frontwith velocity (12) continues until it is overtaken by the pbshock, generated by the pressure jump ∆ p pb = p EQ,sp + σ sp behind the spall front; here p EQ,sp > E in Fig. 1),calculated from the EQ EOS branch for the same valuesof v sp and the internal energy e sp as in the pre-flip MSstate ( v sp , p sp ) (point S in Fig. 1). The pb-shock frontis locally coplanar with the spall front and propagates inthe same direction with the normal velocity u pbf = u · ∇ ρ |∇ ρ | + D pb , (13)where D pb = D pb ( s ) is the velocity of the shock front,driven by the pressure jump ∆ p pb , relative to the liquidin the pre-shock state ( v sp , p sp ); it is, of course, calculatedby using the MS EOS branch. Similar to v sp , σ sp , and∆ p pb , the shock velocity D pb is uniquely determined bythe entropy value s .Once the condition u pbf ≥ u spf is met and the pbshock breaks out before the spall front, it recompressesand heats up the stretched liquid, preventing it fromreaching the spall threshold. As a result, the expansionof the fog bubble in Lagrangian coordinates, i.e. relativeto the fluid particles, comes to a halt. In other words,the mass of fractured liquid continues to grow so long asthe effective spall Mach number , defined as M sp def = (cid:18) u spf u pbf (cid:19) u =0 = ˙ εD pb |∇ ln ρ | = ˙ ερc D pb |∇ p | , (14)stays above unity, M sp >
1. Note that M sp is de-fined in the comoving frame, and relative not to the localsound speed c but to the pb-shock velocity D pb > c ;hence, we use the term hypersonic to a spall front with M sp >
1. The condition M sp = 1, when fulfilled on thespall-front surface, delimits the size and mass of the fogbubble. The bubble mass is always finite when ˙ ε > ∇ ρ = ∇ p = 0 at the inception point. It may become zeroin some singular cases where either ˙ ε changes sign or ∇ p is discontinuous at inception. In practice, a fair estimateof the fractured mass can often be made by analyzingthe spatial distribution of M sp ( t sp , x ) at the moment ofinception t sp . The above arguments and formulae can bereadily generalized to non-isentropic flows, where M sp becomes a function of gradients of two independent ther-modynamic variables. III. NUMERICAL RESULTS FOR A SHOCK-LOADEDPLANAR TARGETA. Formulation of the problem and the initial state
In a common experiment for determination of the spallstrength, two planar plates are collided to launch the ini-tial load shocks in two opposite directions from the im-pact plane . In this work we focus our attention on thesimplest version of such an experiment, where both plateshave the same composition and dimensions. Then theproblem is one-dimensional (1D) and symmetrical withrespect to the impact plane, which we place at the ori-gin x = 0 of the center-of-mass coordinate system andconsider only one (right) plate at x ≥
0, which we referto as the target. We set t = 0 at the moment when theload shock breaks out at the outer edge x = h , leav-ing behind the shocked material in a compressed stateat zero initial velocity. Thus defined initial state ( v , p )(point P in Fig. 1) lies on the Hugoniot OP , originatingfrom the normal state ( v , p ) at point O in Fig. 1. Theintensity of the initial load is characterized by the parti-cle (piston) velocity u p , measured relative to the fluid infront of the loading shock. The latter means that, beforethe impact, the flyer and the target plates had, respec-tively, the velocities + u p and − u p in the center-of-massframe; their density and thickness were ρ and h . Thevalue of u p uniquely determines the entropy s through-out the subsequent flow, the spall strength σ sp , as well asthe pressure amplitude ∆ p pb and the front velocity D pb of the pb shock. Table II lists the key flow parametersthat can be calculated directly from EOS for two chosenreference cases of the shock load.Before specifying the initial target thickness h (or itsthickness h before the shock load), we remark the fol-lowing. One easily verifies that the above formulatedproblem is invariant with respect to rescaling of all thelengths and times by one and the same arbitrary factor.Therefore, having obtained the full solution for an initialtarget thickness h , we can rescale it to any other thick-ness h ′ by multiplying times and lengths by the factor h ′ /h . In view of such scalability, all the simulations dis- cussed below were done for h = 1. Recalling that v , p ,and u are, respectively, everywhere in units of V cr , P cr ,and ( P cr V cr ) / , the latter is equivalent to setting theunits of length and time to h and h ( P cr V cr ) − / . B. Simulation results for u p = 1 . and u p = 4 . Given an adequate two-phase EOS, implementation ofour spall criterion into a standard Lagrangian hydrocodeis straightforward. Here we present numerical results, ob-tained with the 1D DEIRA code that has been previouslydeveloped and extensively used for simulations of inertialconfinement fusion targets . Figures 3 and 5 displaythe target density evolution on the ( x, t ) plane calculated,respectively, for u p = 1 .
25 and u p = 4 .
0. For a detaileddiscussion, we focus on the u p = 1 .
25 case, which liesmoderately above the spall threshold of u p,sp ≈ . u p = 4 . x = 1,arrives at and is reflected from the symmetry plane x = 0,a tension region with ρ < .
90 develops near the targetcenter at t & .
12. Further on, as the pressure in thisregion falls to p = − σ sp = − .
15, all the liquid massat 0 < x . .
19 undergoes spall fragmentation within ashort time interval of t sp = 0 . < t < . t = 0 . t = 0 . FIG. 3. Color density map on the space-time diagram ascalculated with the 1D DEIRA code for u p = 1 . More insight into the spall dynamics can be gained byexamining the profiles in Figs. 4 and 6 along the dimen-sionless Lagrangian coordinate¯ m = mm , m = x Z ρ ( t, x ′ ) dx ′ , m = ρ h = ρ h , (15) TABLE II. Shock-load and subsequent flow parameters that can be calculated directly from EOS for two selected values of thepost-shock particle velocity u p . For detailed explanations see the text u p u ρ = v − c ρ c b z p EQ,sp D pb,sp c sp σ sp m k m. strain rate u p = 1.25, t = 0.21818 spall Mach number M sp -12-11-10-9-8-7-6 p pressure p FIG. 4. Profiles of the strain rate ˙ ε , the spall Mach number M sp (left ordinates), and of the pressure p (right ordinates)versus the fractional target mass ¯ m just before the spall failurefor u p = 1 .
25. Only the inner half of the target mass is shown plotted for the time moment immediately preceding t sp .(Note that these profiles are invariant with respect torescaling of the target thickness.) In Fig. 4 the spallfailure occurs at t sp = 0 . m = 0, wherethe pressure is minimum and M sp = ∞ . The outer limitof the spall zone is clearly marked by a precipitous dropof the effective spall Mach number M sp , which crossesthe M sp = 1 line at ¯ m k = 0 . t = 0 . m sp = 0 . m k = 0 .
572 and ¯ m sp = 0 . FIG. 5. Same as Fig. 3 but for u p = 4 . . strain rate m k mu p = 4.0, t = 0.14984 spall Mach number M sp -10-8-6-4-20 p pressure p FIG. 6. Same as Fig. 4 but for u p = 4 . Two more observations to be made from Figs. 4 and6 are (i) a virtually flat pressure profile across the incip-ient spall zone 0 < ¯ m < ¯ m k , and (ii) a jump by abouta factor 2 in the value of the strain rate ˙ ε at ¯ m = ¯ m k .Here we must recall that the exact solution to our prob-lem contains weak discontinuities (see problem 1 to § ), i.e. discontinuities of ∂p/∂x , ∂ρ/∂x , ∂u/∂x that aresmeared by the artificial viscosity in numerical profiles ofFigs. 3–6. In Figs. 4 and 6 the kink point ¯ m k in the pres-sure profile is just such a discontinuity. As a consequence,the exact profiles of M sp and ˙ ε would both have jumpsat ¯ m = ¯ m k . In our simulations the threshold value of M sp = 1 is bracketed with a wide margin by its pre- andpost-jump values. Due to nonlinearity of the problem,the pressure profile at ¯ m < ¯ m k is not exactly flat: it ismonotonic with a minimum at ¯ m = 0 for u p = 1 .
25, andnon-monotonic with a minimum at the kink point ¯ m k for u p = 4 .
0. That is, the exact position of spall initiationmay differ depending on the EOS details and the loadintensity. The narrow peak of M sp just before ¯ m = ¯ m k in Fig. 6 is caused by numerical smoothing of this cornerpoint in the p and u profiles, and would not occur in theexact solution. C. On morphology of the spall zone
The space-time image of Fig. 3 offers a convenientopportunity to make a comparison with the publishedMD simulation of a similar problem for liquid copper .When rescaled to the critical point of copper , the pa-rameters of our u p = 1 .
25 run are close to the MD caseshown in Fig. 1 of Ref. 12. The overall spall dynamics, thefinal size, and the qualitative structure of the spall zonelook also quite similar. The spall strength σ sp ≈ P cr ,calculated in Ref. 12, practically coincides with the spin-odal limit given in Table II. Significant differences emergeonly when the structure of fractured liquid is examinedin more detail.In the results of Ref. 12, the spall zone has a width ofsome 600 interatomic distances and the fractured masssplits into three liquid drops separated by voids. Gener-ally, one would expect the MD approach to produce somestatistical distribution of liquid fragments, whose typicalsizes would be determined by the atomic length scale, in-herent in any MD model, and by statistical properties ofthe microscopic initial states used in the MD code runs.The present model, in contrast, is expected to producea fully deterministic pattern of spall fragmentation be-cause it is based on the equations of ideal hydrodynam-ics and contains no length scale to stipulate a specificfragment size. According to the arguments of Sect. II C,the entire spall zone in the exact solution to our prob-lem should look like a single fog bubble without any re-maining liquid drops because the condition M sp ≫ ρ , p , and u but those of their spatial derivatives are kept sufficientlysmall — hence the strong sensitivity of the fragmenta-tion morphology to small perturbations with high wavenumbers.As an illustrative example, Fig. 7 displays the samecase as in Fig. 3 but simulated with a perturbed initialdensity ρ (0 , x ) = ρ [1 − .
001 cos(46 πx )] (16)under the unaltered pressure p (0 , x ) = p . One sees thateven as small as a 0.1% density perturbation results in adramatic change of the spall pattern: instead of a singlefog bubble we have three liquid drops separated by fourfog bubbles. Once well resolved by the mesh (43 cellsper each drop in Fig. 7), the mass and other dynamicproperties of these three drops remain stable with respect to further mesh refinement. FIG. 7. Same as Fig. 3 but calculated with the initial densityperturbation (16)
In summary, being fully deterministic for given initialand boundary conditions, our model should predict real-istic morphology of the spall zone when realistic smallperturbations are included into the problem formula-tion. By numerical implementation, a special care maybe needed to suppress spurious short-wavelength pertur-bations due to finite differencing.
IV. FIRST POST-ACOUSTIC APPROXIMATION FOR ASYMMETRIC PLANAR SPALL
In most plate impact experiments, the spall strength σ th is inferred from the time history of the free surfacevelocity u fs ( t ) measured at the outer target edge . Fig-ure 8 shows the plots of u fs ( t ), calculated for two casesof the initial load in the above stated problem, one be-low ( u p = 0 .
8) and the other above ( u p = 1 .
25) the spallthreshold. Given the density ρ and the sound speed c of the normal state, σ th is most often evaluated in thelinear acoustic approximation from the formula σ th,ac = 12 ρ c ∆ u fs , (17)where ∆ u fs = u fs,max − u fs,min (18)is the amplitude of the first dip in the u fs ( t ) profile. How-ever, contrary to what might be expected, formula (17)never becomes asymptotically accurate even in the limitof an infinitely weak load because of the pb-shock atten-uation: below we show its error to be 40% for the con-sidered case of symmetric plate collision. The situationbecomes even worse when the acoustic formulae ˙ ε sp,ac = ∆ u fs c ( t pb − t ) , ¯ m sp,ac = 1 − c t pb h (19)for the spall strain rate ˙ ε sp and the fractured mass frac-tion ¯ m sp are invoked: the error in ˙ ε sp is often a factor 3–5,while absurd negative values are easily obtained for ¯ m sp .In this section we go to the next order in the nonlinearterms and derive new formulae in what we call the firstpost-acoustic (PA) approximation, which provide signif-icantly better accuracy than Eqs. (17) and (19) for therespective quantities. t pb u fs,min u p =0.80 u p =1.25 u fs tu fs,max = u t FIG. 8. Time dependence of the free-surface velocity u fs ( t ) atthe outer target edge ¯ m = 1, calculated with the DEIRA codefor two cases of the initial shock load: one below ( u p = 0 . u p = 1 .
25, black solid) thespall threshold
A. Characteristics chart
As the loading shock arrives at the outer target edge¯ m = 1 at t = 0, a rarefaction wave R1 is launched in thenegative direction towards the symmetry plane ¯ m = 0, asis schematically shown in Fig. 9a. The fluid before the R1front is at rest in the shock-compressed state ρ = v − , p = p ( ρ , s ) (point P in Fig. 1); the fluid behind itmoves with a constant positive velocity u = p Z dpρc = ρ Z ρ c d ln ρ (20)and is in the thermodynamic state ρ , p ≡ p ( ρ , s ) =0 (point O in Fig. 1). Here and below we assume that theMS EOS (1)-(3) is cast in the form p = p ( ρ, s ). So longas no new shocks emerge and the entropy remains at itspost-load-shock value s = s , the pressure p = p ( ρ, s ) = p ( ρ ) and the sound speed c = ( dp/dρ ) / = c ( ρ ) arefunctions of ρ only.The key to further analysis is the behavior of charac-teristics for the two governing hydrodynamic equations,whose Lagrangian form in our case is ∂v∂t − ∂u∂m = 0 , (21) ∂u∂t + ∂p∂m = 0 , (22) R1 tailR1 head p u mR1 wave p,u (a) u pb p pb u mR2 wave p,u (b) pb shock R2 head FIG. 9. Schematic velocity and pressure profiles in the left-propagating R1 (a) and right-propagating R2 (b) rarefactionwaves. The right-propagating R2 wave is “propped up” frombehind by the pull-back shock where the mass coordinate m and its normalized counter-part ¯ m are defined in Eq. (15). In terms of the normalizedacoustic impedance z = z ( ρ ) def = ρcρ c , z = ρ c ρ c , (23)and the normalized time¯ t def = ρ c ρ h t = c th z , (24)the two families of C ± characteristics on the ( ¯ m, ¯ t ) planeare given by d ¯ md ¯ t = ( + z, C + : dJ + ≡ du + dp/ρc = 0 , − z, C − : dJ − ≡ du − dp/ρc = 0 , (25)where J ± are the two respective Riemann invariants .Figure 10 shows the pattern of C ± characteristics, rep-resenting the solution to our problem in the normal-ized ( ¯ m, ¯ t ) plane. The R1 wave starts at point 1 with( ¯ m, ¯ t ) = (1 ,
0) as a centered C − fan between the lead-ing C − h and the trailing C − t lines. As each of these C − characteristics reaches the symmetry plane ¯ m = 0, it isreflected and continues as a C + characteristic, i.e. C + h is the reflected continuation of C − h . Where the R1 fanborders a region of constant flow, it is a simple wave withall the C − characteristics being straight. In Fig. 10 theR1 simple wave (shaded cyan) is confined to a partiallycurvilinear triangle bounded by the contour 1-¯ t - M -1,inside which J + is constant. In the left reflection zoneR12 (shaded grey), bounded by the contour ¯ t - M -¯ t -¯ t ,the flow is not a simple wave, and all the C ± character-istics are curved.As the R1 wave is reflected from the left edge ¯ m = 0,it becomes the R2 simple wave (shown schematically inFig. 9b) beyond the R12 reflection zone. If the pb shockis launched from point M sp , lying on the C − t characteris-tic, the R2 simple-wave zone (shaded cyan in Fig. 10) isbounded by the contour M - M sp - M -¯ t - M , where M sp - M is the curved trajectory of the pb shock and ¯ t - M is a curved segment of the C − h characteristic, being the0 m o2 t t o2 C +t2 C +o2 C -h3 C +h2 C -t1 z = z z = z z = m m m sp t t pb t t m t sp t t M sp M o2 M M R12 O1R23R1 C -h1 R2 z = pb s ho ck FIG. 10. The pattern of characteristics in the normalized( ¯ m, ¯ t ) plane, depicting interaction of the release waves initi-ated by the shock load with u p = 1 .
25. Cyan-shaded areasshow the R1 and R2 simple wave zones. Grey-shaded areasare the R12 and R23 reflection zones. Thick red curve M sp - M -¯ t pb is the trajectory of the pull-back shock continuation of C + h upon its reflection off the free sur-face ¯ m = 1. The curvilinear triangle ¯ t - M -¯ t pb (shadedgrey), where ¯ t pb is the time of the pb-shock arrival at¯ m = 1, encompasses the right reflection zone R23. Inthe R2 zone, all the C + characteristics are straight —like the C +2 o characteristic passing through M . In theabsence of spall, the R1 trailing characteristic C − t is re-flected from the ¯ m = 0 edge at ¯ t and continues asthe C + t straight line, which replaces the pb shock as thedownstream boundary of the R2 simple wave.The O1 zone inside the rectilinear isosceles triangle 1- M -¯ t is the region of constant flow with u = u , ρ = ρ , c = c , p = 0. The times t and t pb in Fig. 8, which markthe beginning and the end of the depression in the u fs ( t )curve, are the non-normalized counterparts of ¯ t and ¯ t pb in Fig. 10. Note that all the head and tail characteristics,plotted in Fig. 10, are also the trajectories of weak flowdiscontinuities. B. Release isentrope in the PA approximation
To construct the desired PA approximation, we needcertain information about the isentrope p ( ρ ) (curve P O S in Fig. 1), passing through the zero-pressure state O at ρ = ρ . We begin by assuming that, in addition to ρ and c = c ( ρ ), one knows from shock experiments the coefficient b in the widely used linear approximation D = c + bu ′ (26)to the compression Hugoniot with the initial state at O ;here D and u ′ are, respectively, the shock-front and thepost-shock particle velocities relative to the fluid in zoneO1, which itself moves with velocity u in our coordinatesystem. Having applied the mass and momentum balancerelations across the shock front, we can convert (26) intothe parametric representation of the Hugoniot p H ( u ′ ) = ρ u ′ ( c + bu ′ ) , (27) ρ H ( u ′ ) = ρ c + bu ′ c + ( b − u ′ . (28)If we consider now an arbitrary isentropic simple wave,bordering zone O1 in Fig. 10 and propagating in the posi-tive direction, we can similarly parametrize the isentrope p ( ρ ) in terms of the particle velocity increment u ′ by in-tegrating the equation du ′ = dp/ ( ρc ) = c d ln ρ, (29)which expresses the fact that dJ − = 0 throughout boththe O1 zone and the bordering simple-wave regions. Indoing this, we can simply put p ( u ′ ) = p H ( u ′ ) becausethe Hugoniot and the Poisson adiabats have a second-order contact at the common origin point, which implies dp/du ′ = dp H /du ′ and d p/du ′ = d p H /du ′ at u ′ = 0.Having applied this argument to our R1 and R2 simplerarefaction waves, we obtain the following PA approxi-mation to the required isentrope p (¯ u ′ ) = ρ c (cid:0) ¯ u ′ + b ¯ u ′ (cid:1) , (30) z (¯ u ′ ) = 1 + 2 b ¯ u ′ , (31)where¯ u ′ = u ′ c , u ′ = ( u − u > , R1 wave ,u − u < , R2 wave . (32)Below, we refer to the second terms on the right-handsides of Eqs. (30) and (31) as the pu-nonlinearity correc-tions . They characterize the nonlinearity of the pertinentEOS and are proportional to the value of the fundamen-tal gasdynamic derivative Γ = c v (cid:18) ∂ v∂p (cid:19) s = v c (cid:18) ∂ p∂v (cid:19) s = 2 b (33)in the state where ¯ u ′ = 0. The accuracy of the PA ap-proximation to the GWEOS isentrope, passing throughthe normal state with parameters from Table I, is illus-trated in Fig. 11, where the essentially positive dimen-sionless ratio p/ ( ρ c u ′ ) is plotted versus ¯ u ′ . The errorin p stays below 5% for | ¯ u ′ | . .
15, but exceeds 30% for¯ u ′ & . -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6012345 p / ( c u ’ ) u’/c acoustic approximation PA approximation exact Hugoniot exact Poisson FIG. 11. PA approximation to the release isentrope (dash-dotted), passing through the normal state ρ , p = 0, versusparticle velocity increment u ′ . Shown also are the exact Hugo-niot (solid gray) and Poisson (solid blue) adiabats calculatedfrom Eqs. (1)–(3) If Eqs. (30), (31) are deemed to be exact along a con-sidered isentrope, they yield ρ (¯ u ′ ) = ρ (cid:2) − (2 b ) − ln(1 + 2 b ¯ u ′ ) (cid:3) − = ρ (cid:2) u ′ − ( b − u ′ + . . . (cid:3) , (34) c (¯ u ′ ) = c (1 + 2 b ¯ u ′ ) (cid:2) − (2 b ) − ln(1 + 2 b ¯ u ′ ) (cid:3) − = c (cid:2) b − u ′ − b ¯ u ′ + . . . (cid:3) . (35)Because neither ρ nor z can be negative, the applicabilityof the above formulae is in any case limited to the range − / b < ¯ u ′ < (cid:0) e b − (cid:1) / b, (36)which is usually sufficient for application to spall exper-iments. In fact, if separate values of ρ ( u ′ ) or c ( u ′ ) areneeded, one can justifiably restrict oneself to the zero-and first-order terms in powers of ¯ u ′ . C. Evaluation of the strain rate
Introduction of nonlinear terms opens up a possibil-ity to obtain an a priori estimate of the strain rate atspallation without invoking the u fs ( t ) curve — the pos-sibility not available in the acoustic approximation. Ourestimate is based on the assumption of time-independentlinear velocity profile ∂u ( t, m ) ∂m ≈ ∂u ( t m , m ) ∂m ≈ u m m (37)throughout the reflection zone R12 — the LU approxi-mation; here the mid-reflection time t m is defined as themoment when the characteristics C − t and C + h in Fig. 10intersect at point M , and the zone O1 of constant flowbehind R1 has its maximum extension m − m m along the m coordinate. Figure 12 confirms that this approxi-mation remains quite accurate even in the highly nonlin-ear case of u p = 4 . u p =4.0, t=t u p =4.0, t=t sp m u / u FIG. 12. Normalized velocity profiles at the moment t m ofR1 reflection and at the moment t sp of spallation, calculatedwith the DEIRA code for the two reference cases of u p = 1 . u p = 4 . The LU approximation (37) can be supported by theobservation that the piece-wise linear profile of u , dis-played in Fig. 9a, is in fact the exact solution of thenonlinear Euler equation for u ( t, m ), obtained by substi-tuting Eq. (30) into Eq. (22); in this solution the width∆ m = ( ρ c − ρ c ) t of the R1 wave grows linearlyin time. When the reflection of such a linear in m rar-efaction pulse from the m = 0 boundary is treated inthe acoustic limit of c − c ≪ c , both approximateequalities in (37) become exact. Also, the fluid state atmid-reflection is then everywhere uniform with p = 0and ρ = ρ . Note that the slope (37) remains finite asboth u and m m become infinitely small in the limit of u p → c → c .To evaluate the mid-reflection coordinate ¯ m m = m m /m , we use the obvious fact that the integral H d ¯ t =0 along any closed contour in the ( ¯ m, ¯ t ) plane. Then wecombine this circulation rule with Eq. (25) and apply it tothe perimeter of the R1 simple-wave zone. Whereas inte-gration along straight characteristic segments, on whichall flow parameters remain constant, is trivial, for anycurved segment — like ¯ t - M in Fig. 10 — we assumethe acoustic impedance to be approximately constantand equal to its arithmetic mean z = ( z a + z b ) / z a and z b . Thus, the condition ofzero H d ¯ t = H ( d ¯ m/d ¯ t ) − d ¯ m along the triangular contour1-¯ t - M -1 yields the equation z − + 2(1 + z ) − ¯ m m − (1 − ¯ m m ) = 0 , (38)from which we get¯ m m = z − z − z + 3 = 2 β (1 + β )(2 + β )(1 + 2 β ) , (39)2where, according to (31), z = 1 + 2 β , β = bu /c . (40)Note that the acoustic limit corresponds to β ≪ ∂u/∂m known from Eq. (37), the strain rate(10) can be calculated from Eq. (21) as˙ ε = ρ ∂v∂t = ρ ∂u∂m . (41)Before finalizing our estimate of ˙ ε , we recall that the spa-tial derivatives of u and p are discontinuous on the C − t characteristic — as is manifested by the kinks at ¯ m =¯ m sp on the velocity profiles in Fig. 12. Hence, we haveto choose between the two different values ( ∂u/∂m ) ± obtained on the two sides m = m sp ± m < m sp whichundergoes spall fracture, we make the logical choice of( ∂u/∂m ) − = u /m m . In contrast, the often used acous-tic formula (19) represents the value ( ∂u/∂m ) + in the un-fractured material, which is exactly one half of ( ∂u/∂m ) − in the acoustic limit. Thus, having combined Eqs. (37)–(41), we obtain the following PA estimates for the strainrates at mid-reflection,˙ ε m = ρ u m m = c h z (2 + β )(1 + 2 β )2 b (1 + β ) , (42)and at spallation˙ ε sp = ρ sp ρ ˙ ε m = (cid:0) − (cid:12)(cid:12) ¯ u ′ sp (cid:12)(cid:12)(cid:1) ˙ ε m . (43)The small correction (cid:12)(cid:12) ¯ u ′ sp (cid:12)(cid:12) ≃ σ th /ρ c , accounting forthe fluid expansion between the t m and t sp moments, istypically a few percent even for the highest theoreticalvalues of σ th and can usually be neglected; the PAformula for | ¯ u ′ sp | is given in the next section.If time t is known from the measurements of the freesurface velocity u fs ( t ), one readily obtains an alternativea posteriori LU estimate¯ m m = 1 − c t h z , (44)˙ ε m = c h u /c z − c t / h , (45)which is a direct consequence of the LU approximation(37) and of the obvious fact that t m = t /
2. Notethat, having the dimension of t − , the strain rate ˙ ε m in Eqs. (42) and (45) scales in inverse proportion to theinitial target thickness h .Regarding practical application of Eqs. (42) and (43),we will distinguish two subcases of the PA approximationdepending on the available thermodynamic informationfor the material in question, namely, a cruder PA0 versionand a more accurate PA1 variant:PA0: ρ , c , b → known, ρ , c , b, u → replaced by ρ , c , b , u p ;PA1: ρ , c , ρ , c , b, u → all known. In the PA0 version, only the normal-state values of ρ , c , b are assumed to be known, while for the PA1 caseone needs also the parameters of the p = 0 state O on the release isentrope behind the initial load shock.Analogously, the LU approximation is subdivided into acruder LU0 and a more accurate LU1 versions. TABLE III. Comparison of the strain rates, evaluated in dif-ferent versions of the PA and LU approximations, with thoseobtained in numerical simulations with the DEIRA codeDEIRA PA1 PA0 LU1 LU0 ac u p = 1 .
25, ∆ u fs = 1 . m m h /c ) ˙ ε m h /c ) ˙ ε sp u p = 4 .
0, ∆ u fs = 1 . m m h /c ) ˙ ε m h /c ) ˙ ε sp Table III compares the values of ¯ m m , ˙ ε m , and ˙ ε sp ,evaluated in different versions of the PA and LU approx-imations, with those obtained in the numerical DEIRAruns for the two reference cases of the shock load; the ac column lists the values given by the acoustic formula (19);the strain rate is given in terms of the scale-invariantdimensionless product ( c /h ) ˙ ε . All the required EOSparameters are taken from Tables I and II; the valuesof ∆ u fs , used to calculate ˙ ε sp,ac and ¯ u ′ sp , as well as thetimes t are inferred from the DEIRA produced plots of u fs ( t ) as shown in Fig. 8 for the case of u p = 1 .
25. Thenumbers in Table III demonstrate that the PA formulae(42), (42) provide a robust estimate of the strain rate atspallation with an accuracy of about 15–30%, which isnot sensitive to the loading shock strength. The conven-tional acoustic formula (19), on the contrary, tends tounderestimate ˙ ε sp by a significant factor of about 3–5, aswas already noticed in MD simulations . If one canreliably measure the arrival time t of the R2 rarefactionwave, Eqs. (43)–(45) yield even more accurate values of¯ m m and ˙ ε sp , which confirms the high accuracy of theLU approximation (37). D. Evaluation of the spall strength
Similar to the acoustic formula (17), our PA estimate ofthe spall strength is based on the measured free-surfacevelocity pullback ∆ u fs . We begin by arguing that theouter boundary ¯ m sp of the spall zone, represented by thepb-shock starting point M sp in Fig. 10, should generallylie on the R1 trailing characteristic C − t , where ∂p/∂m and ∂u/∂m have a jump. Firstly, we note that every-where above and to the right of the M -¯ t segment ofthe C − t curve, where the pressure is negative in the R23simple-wave zone, we have M sp = c/D pb < . (46)This inequality, obtained by substituting Eq. (41) andthe relation dp = ± ρc du into the definition (14) of theeffective spall Mach number M sp , implies that the pbshock must originate either from the C − t characteristic orfrom the interior of the R12 reflection zone. On the otherhand, we know that in the acoustic limit β ≪ ∂p/∂m → M sp → ∞ everywhere inside R12. This proves thatat least for not too large values of u p and β the pb shockcan only start from the C − t characteristic. In stronglynonlinear cases with β &
1, this fact is confirmed bynumerical simulations (at least for the present EOS), asis shown in Fig. 6.With M sp lying on the C − t characteristic above the M point, we can use the analytic properties of the R2simple wave, represented by Eqs. (30)–(35), to relate σ sp to ∆ u fs . The simple acoustic formula (17) is based onthe premise that the negative-pressure spall pulse (seeFig. 9b), created at M sp with the initial velocity decre-ment u ′ sp ≡ c ¯ u ′ sp def = u ( t sp , m sp ) − u < , (47)propagates without attenuation down to the free surface,where the well-known ( §
11, Ch. XI) velocity doublingrule ∆ u fs = − u ′ sp applies; then (17) follows directlyfrom Eq. (30) in the limit of b | ¯ u ′ | ≪ C − h at point M and apply the doubling rule∆ u fs = − u ′ (48)to the velocity decrement u ′ ≡ c ¯ u ′ = u ( t m , m m ) − u < . (49)at that point, which has coordinates ( ¯ m m , ¯ t m ) inFig. 10. Having introduced the pb-attenuation correction δ pb def = u ′ sp /u ′ − σ sp = − p (¯ u ′ sp ) = ρ c | ¯ u ′ sp | (cid:0) − b | ¯ u ′ sp | (cid:1) = 12 ρ c ∆ u fs z − (1 + δ pb ) [1 − β (1 + δ pb )] , (51)where β = b | ¯ u ′ | = b ∆ u fs c , | ¯ u ′ sp | = (1 + δ pb ) ∆ u fs c . (52) The attenuation correction δ pb = 2(1 − β ) ω + [ ω − β (1 − β )(3 µ + 1)] / (53)is found as the solution of the quadratic equation (A14),where µ = 1 + β β − β − β , ω = 12 + µ (2 − β ) . (54)Note that formula (51) may be considered as a general-ization of the simple extrapolation of the Hugoniot intothe negative-pressure region, advocated in Ref. 29, whichtakes a more consistent account of the first-order nonlin-ear terms. For the fractional mass of the spall zone weobtain¯ m sp = ¯ m m (cid:18) − | u ′ sp | u (cid:19) = ¯ m m (cid:20) − (1 + δ pb )∆ u fs u (cid:21) , (55)where ¯ m m is given by Eq. (39). The derivation ofEqs. (53)–(55) is detailed in Appendix A.Although the effect of the pb-shock attenuation is fairlywell known and has been discussed in literature , theexplicit formula (53) for the corresponding correction δ pb is an important new result of this paper. For the par-ticular considered problem of symmetric plate collision,this correction is a function of two dimensionless param-eters β and β representing (in relative terms), respec-tively, the intensity of the initial load and the amplitudeof the spall pulse. Figure 13, which displays the behav-ior of δ pb ( β , β ) in the allowed range of 0 ≤ β ≤ . β . . . ≤ δ pb < .
56 for anyvalue of β >
0, i.e. for arbitrary velocity of the impactorplate. This observation is important for applications be-cause in real experiments the limit of β ≈ . pb 2 FIG. 13. The pb-attenuation correction δ pb = δ pb ( β , β ) asa function of β for three selected values of β β , β → δ ≡ δ pb (0 ,
0) = 25 (56)that is not much lower than the upper limit of δ pb,max ≈ .
56. The latter is explained by the fact that for β ∼ β ≪ ( § δ turns out to be a universal constant which dependsneither on Γ nor on target parameters. Its value can becalculated independently and more rigorously by directintegration of the equation (A1) for the pb-shock trajec-tory in the asymptotic limit of β , β →
0, which yields δ = √ − δ = 0 . ρ c = z − ρ c , the pb-attenuation cor-rection δ pb , and the negative pu-nonlinearity correction − β (1 + δ pb ). The relative importance of these correctionterms can be judged from the data in Tables II and IV.The negative impedance correction z − − u p . The positive pb-attenuation correction δ pb is alwayssignificant and must never be ignored, at least not in theexperiments where the flyer and the target plates havecomparable masses.The relative impact of the pu-nonlinearity term de-pends on the value of β , which is insensitive to u p andreaches its maximum for the maximum possible spallstrength. But even then the value β ≈ .
3, obtainedin our simulations and quoted in Table IV, appears asrather an absolute upper limit; more typical for liquidmetals at maximum possible tension would be β . . β ≈ .
12 obtained in the MD sim-ulations for liquid copper . The latter implies that thepu-nonlinearity should be taken into account near thetheoretical limit for σ th , but can be ignored in exper-iments where the spall failure occurs at tensions 3–10times below this limit.A characteristic feature of the formula (51) is that thecontributions of its three PA correction terms tend tocompensate for one another. For large u p , the decreas-ing impedance ρ c of the zero-pressure release statestrongly suppresses the impact of the 1 + δ pb factor; for β ≈ . , and why in TABLE IV. Comparison of the PA and acoustic (ac) approx-imations with the exact (EOS or DEIRA) results for the spallstrength σ th and the spalled mass fraction ¯ m sp , calculated forthe two reference values of u p . The values of ∆ u fs are takenfrom the u fs ( t ) curves produced by the DEIRA runsEOS/DEIRA PA1 PA0 ac u p = 1 .
25, ∆ u fs = 1 . β – 0.6442 0.6614 – β – 0.2941 0.2939 – δ pb – 0.4370 0.4378 – σ sp m sp u p = 4 .
0, ∆ u fs = 1 . β – 2.031 2.117 – β – 0.2879 0.2812 – δ pb – 0.4740 0.4800 – σ sp m sp practice the simple formula (17) may frequently be moreaccurate than more complex expressions with only a par-tial account for the nonlinearity effects. On the whole,the data in Table IV confirm that Eqs. (51) and (55)provide a significant improvement over the traditionalacoustic formulae for interpretation of spall experimentswith symmetric initial loading.Our final remark is on the applicability of the acousticlimit. Naively, one might expect that in the limit of smallamplitudes of both the loading shock and the spall pulse,i.e. in the limit of β , β → z →
1, the acoustic for-mula (17) should be asymptotically exact. Our presentanalysis demonstrates that, whenever a pull-back shockis present, this is not the case. Instead, for symmetricallyloaded impact plates, the correct small-amplitude limitis given by the formula σ sp = 12 (1 + δ ) ρ c ∆ u fs = 1 √ ρ c ∆ u fs , (57)which exceeds the acoustic result by 41%, and whichwould be suitable for interpretation of experiments where β . . β . . V. CONCLUSIONS
The paper presents the results of detailed analysis ofthe phase-flip (PF) model of spallation in pure liquids,where spall failure is treated as an infinitely fast irre-versible liquid-gas phase transition. The spall strength isset equal to its theoretical maximum determined by thethermodynamic stability limit, i.e. by the negative pres-sure on the liquid branch of the spinodal in a two-phaseEOS. The main advantage of the proposed model is that5it can be easily and in a thermodynamically fully consis-tent way incorporated into the framework of ideal hydro-dynamics. Unlike in the approach proposed in Ref. 19and 20, energy dissipation due to spall fracture is de-scribed without invoking the concepts of surface tensionand/or viscosity, simply as the entropy jumps in the PFtransitions and in shocks generated by the PF pressurejumps.The process of spall fracture in the proposed model de-velops as a sudden appearance of one or more “fog” bub-bles behind hypersonic spall fronts (representing a newtype of discontinuity in ideal hydrodynamics) that arelaunched from one or more isolated inception points at lo-cal pressure minima. In our context, the term “fog” refersto the equilibrium matter state, to which a metastableliquid relaxes in the process of spinodal decomposition.The volume and mass of each fog bubble continue to growuntil the spall front, precipitating the PF transition, isovertaken by the pull-back compression shock.To illustrate how the PF model could be used in prac-tice, the results of numerical modeling of a symmetricplate impact experiment are presented. In most aspects,they look quite similar to those found in the MD simu-lation of an analogous problem — but obtained at im-measurably lower computational cost. The principal dif-ference between the two approaches manifests itself onlyin the microstructure of the fracture zone: while in theMD case it is stipulated by the statistics and the inherentatomic scale, in the fully deterministic PF approach thesizes and distribution of liquid fragments are determinedby small flow perturbations and inhomogeneities of theinitial state. When equipped with an adequate two-phaseEOS, the PF model has a potential to become a power-ful tool for modeling experiments where spallation andcavitation in liquids is an issue.Using the PF numerical results as a reference solution,we derive new post-acoustic formulae for evaluation ofthe strain rate, the fractured mass fraction, and the spallstrength in plate impact experiments where these quan-tities are to be inferred from the measured free-surfacevelocity. Though valid for the particular case of a sym-metric plate impact, the new approximate formula (51)for σ th reveals a non-trivial interplay of the three nonlin-ear correction terms, which in reality often compensatefor one another. A notable new result is that the correc-tion factor for the pb-shock attenuation never becomesclose to unity, and in the limit u p ≪ c of weak initialloads approaches a universal value 1 + δ = √ τ pr re-mains finite and limited to τ pr & ∼ u spf τ pr , where u spf is the spall- front velocity in the comoving system. Consequently,the proposed model becomes applicable on time scales t ≫ τ pr and length scales l ≫ c τ pr &
10 nm, where c is the sound velocity in the relevant zero-pressuremetastable state. Also, the length scale of emerging frac-ture zones must exceed the typical size of critical va-por bubbles, which again sets a limit of l ≫
10 nm. Ifthe strain-rate dependence of the spall strength is to beaccounted for, the present assumption of the spinodal-limited spall strength can be relaxed along the lines pro-posed in Refs. 35, , and in other similar works. ACKNOWLEDGMENTS
This work was funded by the Keldysh Institute of Ap-plied Mathematics of the Russian Academy of Sciencesunder the state contract with the Ministry of Science andHigher Education of the Russian Federation.
Appendix A: Derivation of formula (53) for thepb-attenuation correction
We begin by ascertaining that the trajectory ¯ m pb (¯ t ) ofthe pb shock, as it propagates from point M sp to point M in Fig. 10 against the background of the R2 simplewave where all flow characteristics are functions of a sin-gle parameter u ′ < d ¯ m pb d ¯ t ≡ ζ (¯ u ′ ) = 1 − b | ¯ u ′ | . (A1)Indeed, the trajectory m ( t ) of any shock front, travelingwith a speed D if relative to the fluid in an initial state( v i , p i ), obeys the equation ( §
16, Ch. I) dmdt = ρ i D if = (cid:18) p f − p i v i − v f (cid:19) / , (A2)where ( v f , p f ) is the final post-shock state. Importanthere is that the shock impedance (A2) remains invari-ant by interchanging the initial and final states. Inour context, where the compression Hugoniot for a fi-nal state p f = p H ( u ′ ) > u ′ > p i = p H (0) = 0 is given by Eqs. (27), (28), the normalizedshock impedance takes the form d ¯ md ¯ t = 1 + b ¯ u ′ . (A3)If we consider now the pb compression shock, whose ini-tial pressure p i = p H ( u ′ ) < u ′ < p f = p H (0) = 0, we conclude that, due to the i ⇄ f in-terchange symmetry, its impedance must be equal to thatof the rarefaction shock (perhaps non-physical but math-ematically conceivable) represented by the extension of6Eqs. (27), (28), and (A3) into the u ′ < ζ ,the front speed of the pb shock D pb = c − ( b − | u ′ | (A4)is not equal to the extension of Eq. (26) to negative u ′ .The fact that ρ H ( u ′ ) in Eq. (28) is not exactly equal to ρ (¯ u ′ ) from Eq. (34) has no practical significance because,when expanded in powers of ¯ u ′ , the two formulae begin todiverge in O (¯ u ′ ) terms only, which can be safely ignoredfor | ¯ u ′ | ≤ | ¯ u ′ sp | . . p EQ,sp (see Sect. II C), can practically always be ne-glected in comparison with σ sp , as is seen from Table II.Next, we introduce an auxiliary characteristic C + o intothe chart in Fig. 10, which belongs to the C + fan of theR2 wave and crosses the right-reflected characteristic C − h at the same point M as the pb shock. Because for u ′ < ζ , C + o is straight and lies below the pb-shocktrajectory everywhere between M and its intersectionwith the C − t characteristic at point M o with coordinates( ¯ m o , ¯ t o ), which lie in the range¯ m sp < ¯ m o < ¯ m m , ¯ t m < ¯ t o < ¯ t sp . (A5)When applied to points M o and M sp , the LU approxi-mation (37) takes the form u m m = u − | u ′ | m o = u − | u ′ sp | m sp , (A6)which relates ¯ m o and ¯ m sp to the known quantity ¯ m m ,leading, in particular, to Eq. (55).Now, the two unknown quantities δ pb and ¯ m m can becalculated from the system of two equations1 − ¯ m m + 1 − ¯ m m (1 + z ) / − ¯ m m − ¯ m o z − ¯ m m − ¯ m o (1 + z ) / , (A7)¯ m m − ¯ m o z − ¯ m m − ¯ m sp ( ζ + ζ sp ) / − ¯ m o − ¯ m sp ( z + z sp ) / , (A8)obtained, respectively, by applying the H d ¯ t = H z − d ¯ t =0 circulation rule to the two closed contours M -¯ t - M - M o - M and M o - M - M sp - M o in the ( ¯ m, ¯ t ) plane, andby using the mean impedance values along the curvilin-ear characteristic segments and the pb-shock trajectory.Here z = 1 − β , ζ = 1 − β , (A9)are the normalized acoustic and pb-shock impedances atpoint M , and z sp = 1 − β (1 + δ pb ) , ζ sp = 1 − β (1 + δ pb ) (A10) are the same impedances at point M sp ; β is defined inEq. (52). Equations (A7) and (A8) are easily transformedto ¯ m m − ¯ m o = (1 − ¯ m m ) z (3 + z )1 + 3 z , (A11)¯ m m − ¯ m o ¯ m o − ¯ m sp = 2 z ( z + z sp + ζ + ζ sp )( z + z sp )( ζ + ζ sp − z ) , (A12)while from Eq. (A6) we get¯ m o − ¯ m sp = ¯ m m ( β /β ) δ pb . (A13)Having substituted Eqs. (A9)–(A11) together withEqs. (A13) and (39) into Eq. (A12), we obtain thequadratic equation(3 µ + 1) β δ pb − ωδ pb + 2(1 − β ) = 0 (A14)for the attenuation correction δ pb , where µ and ω aredefined in Eq. (54). The basic equations (A7), (A8) andthe ensuing formulae are physically meaningful only ifthe impedance 1 ≥ z ≥
0, i.e. under the condition0 ≤ β ≤ . (A15)It is not difficult to prove that the above condition guar-antees positiveness of the discriminant in Eq. (A14). T. Antoun, L. Seaman, D. R. Curran, G. I. Kanel, S. V. Ra-zorenov, and A. V. Utkin,
Spall Fracture , Shock Wave and HighPressure Phenomena (Springer-Verlag New York, 2003). G. I. Kanel, V. E. Fortov, and S. V. Razorenov, “Shock waves incondensed-state physics,” Physics-Uspekhi , 771–791 (2007). G. I. Kanel, A. S. Savinykh, G. V. Garkushin, andS. V. Razorenov, “Dynamic strength of tin and lead melts,”JETP Letters , 548–551 (2015). E. B. Zaretsky, “Experimental determination of thedynamic tensile strength of liquid Sn, Pb, and Zn,”Journal of Applied Physics , 025902 (2016). T. de Ress´eguier, L. Signor, A. Dragon, M. Boustie,G. Roy, and F. Llorca, “Experimental investi-gation of liquid spall in laser shock-loaded tin,”Journal of Applied Physics , 013506 (2007). M. B. Agranat, S. I. Anisimov, S. I. Ashitkov, V. V. Zhakhovskii,N. A. Inogamov, P. S. Komarov, A. V. Ovchinnikov, V. E. Fortov,V. A. Khokhlov, and V. V. Shepelev, “Strength properties of analuminum melt at extremely high tension rates under the actionof femtosecond laser pulses,” JETP Letters , 471–477 (2010). M. S. Krivokorytov, A. Y. Vinokhodov, Y. V. Sidelnikov, V. M.Krivtsun, V. O. Kompanets, A. A. Lash, K. N. Koshelev, andV. V. Medvedev, “Cavitation and spallation in liquid metaldroplets produced by subpicosecond pulsed laser radiation,”Phys. Rev. E , 031101 (2017). C. A. Stan, P. R. Willmott, H. A. Stone, J. E. Koglin, M. Liang,A. L. Aquila, J. S. Robinson, K. L. Gumerlock, G. Blaj, R. G.Sierra, S. Boutet, S. A. H. Guillet, R. H. Curtis, S. L. Vetter,H. Loos, J. L. Turner, and F.-J. Decker, “Negative pressures andspallation in water drops subjected to nanosecond shock waves,”J. Phys. Chem. Lett. , 2055–2062 (2016). M. M. Basko, M. S. Krivokorytov, A. Y. Vinokhodov, Y. V.Sidelnikov, V. M. Krivtsun, V. V. Medvedev, D. A. Kim, V. O.Kompanets, A. A. Lash, and K. N. Koshelev, “Fragmentationdynamics of liquid-metal droplets under ultra-short laser pulses,”Laser Physics Letters , 036001 (2017). S. Y. Grigoryev, B. V. Lakatosh, M. S. Krivokorytov, V. V.Zhakhovsky, S. A. Dyachkov, D. K. Ilnitsky, K. P. Migdal, N. A.Inogamov, A. Y. Vinokhodov, V. O. Kompanets, Y. V. Sidel-nikov, V. M. Krivtsun, K. N. Koshelev, and V. V. Medvedev,“Expansion and fragmentation of a liquid-metal droplet by ashort laser pulse,” Phys. Rev. Applied , 064009 (2018). V. V. Pisarev, A. Y. Kuksin, G. E. Norman, V. V. Stegailov, andA. V. Yanilkin, “Microscopic theory and kinetic model of spallin liquids,” AIP Conference Proceedings , 801–804 (2009). Y. Cai, H. A. Wu, and S. N. Luo, “Spall strengthof liquid copper and accuracy of the acoustic method,”Journal of Applied Physics , 105901 (2017). A. E. Mayer and P. N. Mayer, “Strain rate dependenceof spall strength for solid and molten lead and tin,”International Journal of Fracture , 171–195 (2020). S. I. Anisimov, N. A. Inogamov, A. M. Oparin,B. Rethfeld, T. Yabe, M. Ogawa, and V. E. For-tov, “Pulsed laser evaporation: equation-of-state effects,”Applied Physics A , 617–620 (1999). J. P. Colombier, P. Combis, F. Bonneau, R. L. Harzic, andE. Audouard, “Hydrodynamic simulations of metal ablation byfemtosecond laser irradiation,” Phys. Rev. B , 165406 (2005). N. Zhao, A. Mentrelli, T. Ruggeri, and M. Sugiyama, “Admis-sible shock waves and shock-induced phase transitions in a vander Waals fluid,” Phys. Fluids , 086101 (2011). M. M. Basko, “Centered rarefaction wave with a liquid-gas phasetransition in the approximation of “phase-flip” hydrodynamics,”Physics of Fluids , 123306 (2018). J. M. Boteler and G. T. Sutherland, “Tensile fail-ure of water due to shock wave interactions,”Journal of Applied Physics , 6919–6924 (2004). D. E. Grady, “The spall strength of condensed matter,”Journal of the Mechanics and Physics of Solids , 353–384 (1988). D. E. Grady, “Spall and fragmenta-tion in high-temperature metals,” in
L. Davison and D.E. Grady and M. Shahinpoor (eds) High Pressure Shock Compression of Solids II: Dynamic Fracture and Fragmentation (Springer, New York, 1996) pp. 219–236. V. Bakshi,
EUV Sources for Lithography , Press Monographs(SPIE Press, 2006). H. Mizoguchi, H. Nakarai, T. Abe, K. M. Nowak, Y. Kawa-suji, H. Tanaka, Y. Watanabe, T. Hori, T. Kodama,Y. Shiraishi, T. Yanagida, G. Soumagne, T. Yamada,T. Yamazaki, S. Okazaki, and T. Saitou, “Perfor-mance of one hundred watt HVM LPP-EUV source,”Proc. SPIE , 94220C–94220C–13 (2015). M. M. Martynyuk, “Generalized van der Waals equation of statefor liquids and gases,” Zh. Fiz. Khim. , 1716–1717 (1991). M. M. Martynyuk, “Transition of liquid metals intovapor in the process of pulse heating by current,”International Journal of Thermophysics , 457–470 (1993). , 28 p. (2018). W. M. Haynes,
CRC Handbook of Chemistry and Physics, 96thEdition , 100 Key Points (CRC Press, 2015). L. D. Landau and E. M. Lifshitz,
Statistical Physics , 3rd ed.(Butterworth Heinemann, 1996). L. D. Landau and E. M. Lifshitz,
Fluid Mechanics , Course oftheoretical physics (Pergamon Press, 1987). G. I. Kanel, “Spall fracture: methodologi-cal aspects, mechanisms and governing factors,”International Journal of Fracture , 173–191 (2010). A. R. Imre, A. Drozd-Rzoska, T. Kraska, S. J. Rzoska, and K. W.Wojciechowski, “Spinodal strength of liquids, solids and glasses,”Journal of Physics: Condensed Matter , 244104 (2008). V. L. Malyshev, D. F. Marin, E. F. Moiseeva, N. A.Gumerov, and I. S. Akhatov, “Study of the tensilestrength of a liquid by molecular dynamics methods,”High Temperature , 406–412 (2015). M. Blander and J. L. Katz, “Bubble nucleation in liquids,”AIChE Journal , 833–848 (1975). V. P. Skripov and A. V. Skripov, “Spinodal de-composition (phase transitions via unstable states),”Sov. Phys. Usp. , 389 (1979). V. P. Skripov,
Metastable Liquids (Wiley, New York, 1974). E. Dekel, S. Eliezer, Z. Henis, E. Moshe, A. Ludmirsky, and I. B.Goldberg, “Spallation model for the high strain rates range,”Journal of Applied Physics , 4851–4858 (1998). S. Faik, M. M. Basko, A. Tauschwitz, I. Iosilevskiy, andJ. A. Maruhn, “Dynamics of volumetrically heated mat-ter passing through the liquid-vapor metastable states,”High Energy Density Phys. , 349–359 (2012). M. M. Basko, “DEIRA: A 1-D 3-T hydrodynamic code for simulating ICF targets driven by fast ion beams. Version 4,”(2001), (accessed 27 September 2020; unpublished). M. M. Basko, “Spark and volume ignition of DT and D2 micro-spheres,” Nuclear Fusion , 2443–2452 (1990). V. Fortov, I. Iakubov, and A. Khrapak,
Physics of StronglyCoupled Plasma , 1st ed., International series of monographs onphysics (Clarendon Press, Oxford, 2006). P. A. Thompson, “A fundamental derivative in gasdynamics,”Phys. Fluids (1958-1988) , 1843–1849 (1971). A. E. Mayer and P. N. Mayer, “Continuum model oftensile fracture of metal melts and its application to aproblem of high-current electron irradiation of metals,”Journal of Applied Physics , 035903 (2015). S.-N. Luo, Q. An, T. C. Germann, and L.-B. Han, “Shock-induced spall in solid and liquid cu at extreme strain rates,”Journal of Applied Physics , 013502 (2009). Y. B. Zel’dovich and Y. P. Raizer,