# On the theory and numerical modeling of spall fracture in pure liquids

Mikhail M. Basko<EMAIL_ADDRESS>http://www.basko.net Keldysh Institute of
Applied Mathematics, Miusskaya square 4, 125047 Moscow, Russia

###### Abstract

A novel phase-flip model is proposed for thermodynamically consistent and
computationally efficient description of spallation and cavitation in pure
liquids within the framework of ideal hydrodynamics. Aiming at ultra-fast
dynamic loads, the spall failure of a liquid under tension is approximated as
an instantaneous decomposition of metastable states upon reaching the spinodal
stability limit of an appropriate two-phase liquid-gas equation of state. The
spall energy dissipation occurs as entropy jumps in two types of discontinuous
solutions, namely, in hypersonic spall fronts and in pull-back compression
shocks. Practical application of the proposed model is illustrated with
numerical simulations and a detailed analysis of a particular problem of
symmetric plate impact. The numerical results are found to be in good
agreement with the previously published molecular-dynamics simulations. Also,
new approximate nonlinear formulae are derived for evaluation of the strain
rate, fractured mass, and spall strength in terms of the observed variation of
the free-surface velocity. The new formula for the spall strength clarifies
complex interplay of the three first-order nonlinear correction 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.

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 in matter under dynamic
tension was conducted mainly for solids Antoun  _et al._ (2003); Kanel,
Fortov, and Razorenov (2007). In recent years, measurements of the spall
strength in liquid metals, loaded by the traditional method of impactor
plates, have been reported Kanel  _et al._ (2015); Zaretsky (2016). Also, many
new experiments have been performed where spall failure occurs in liquids,
often in molten metals, at very high loading rates under the action of
nanosecond and sub-picosecond laser pulses de Rességuier  _et al._ (2007);
Agranat  _et al._ (2010); Krivokorytov  _et al._ (2017), or very short X-ray
bursts Stan  _et al._ (2016). Adequate modeling of dynamic spall fracture in
liquids is a challenging problem, especially when the geometry of the
experiment is not one-dimensional Stan  _et al._ (2016) and/or spall fracture
is only part of a more complex problem as, for example, by modeling liquid tin
targets for EUV lithography sources Basko  _et al._ (2017); Grigoryev  _et
al._ (2018).

A well-established approach is to use the first-principle molecular dynamics
(MD) simulations Pisarev  _et al._ (2009); Cai, Wu, and Luo (2017); Mayer and
Mayer (2020), which provide an important insight and allow calculation of the
ultimate spall strength in liquids at high strain rates but are limited to
sub-micron targets, are computationally very costly and, for that reason, not
quite practical. A much simpler and computationally inexpensive alternative is
to simply adopt the fully equilibrium (EQ) equation of state (EOS), obtained
by applying Maxwell’s rule to the primitive EOS of the van-der-Waals (vdW)
type Anisimov  _et al._ (1999); Colombier  _et al._ (2005); Zhao  _et al._
(2011); Basko  _et al._ (2017). However, because the EQ EOS allows no negative
pressures, it cannot describe tensile states, pull-back shocks, and only
approximately reproduces mass distribution in the regions of spallation and
cavitation Basko (2018a). Another often used model Boteler and Sutherland
(2004); Grigoryev  _et al._ (2018) is based on combining an appropriate
empirical EOS that admits negative pressures with an ad hoc limit on the
tensile stresses: once and where the spall threshold is exceeded, an expanding
void cavity with zero-pressure boundaries is introduced into the simulation.
The principal flaw of this approach is its inability to adequately account for
the liquid-vapor phase transition, especially when spall failure occurs
practically simultaneously over a large quasi-uniform mass of metastable
liquid, converting it into a finely dispersed liquid-vapor mixture.

In the present work, we investigate a novel phase-flip (PF) model based on the
premise that all the information, needed for modeling cavitation and
spallation in pure liquids, is provided by a two-phase liquid-gas vdW-type
EOS, and no other material properties such as surface tension, viscosity, etc.
are to be invoked. This is a clear aspect where our model differs from the
energy-based scheme by GradyGrady (1988, 1996), preserving nonetheless the
capability to adequately account for the energy dissipation by spallation.
Aiming at description of ultra-fast dynamic processes, we assume the spall
strength to attain its maximum possible value, i.e. to be defined by the
thermodynamic stability limit of the vdW-type EOS along the spinodal curve.
The spall fracture is then naturally associated with the spinodal
decomposition, which, in line with the adopted minimalistic approach, is
assumed to take place instantly. To justify the latter assumption, we invoke
the theory of homogeneous bubble nucleation in metastable liquids.

The primary motivation for proffering the PF model is its potential to be
easily and computationally cheaply integrated into large multi-dimensional
fluid dynamics codes that are already loaded with many other complex physical
processes like spectral radiation transfer, laser energy deposition, etc. when
need arises to describe target fragmentation by spallation and cavitation
under fast and intense initial loads. A notable example is the task to
adequately describe the initial shattering of small liquid-tin drops by
picosecond laser pulses in EUV-lithography applications Bakshi (2006);
Mizoguchi  _et al._ (2015), where MD codes cannot be directly applied, and
where the PF model promises to become a major improvement over the more
primitive models employed in Refs. Basko  _et al._ , 2017; Grigoryev  _et al._
, 2018.

The full formulation and justification of the PF model is given in Sect. II,
where we start by introducing the employed vdW-type EOS, describe and analyze
the adopted spall criterion, and explore the kinematics of hypersonic spall
fronts. To illustrate how the proposed model could be applied in practice, we
perform a detailed analysis of the spallation process in a planar collision
between identical flyer and target plates. Although real experiments are
usually done for asymmetric configurations with target-to-flyer thickness
ratios about 3–5 Antoun  _et al._ (2003), we choose the theoretically simplest
symmetric case to display the key properties of our model, and to compare it
with the available MD simulations. Our numerical results, presented in Sect.
III, indicate that the PF model can reproduce the key aspects of the spall
dynamics no less accurately than the MD simulations — provided that an
adequate two-phase EOS is available. The difference emerges only in the
microstructure of fracture zones, stipulated by statistical atomic-scale
properties of an MD model, and determined by small-scale perturbations in the
PF model.

Interpretation of plate impact experiments for determination of the spall
strength $\sigma_{sp}$ relies on the theoretical relationship between
$\sigma_{sp}$ and the observed variation $\Delta u_{fs}$ of the free-surface
velocity (the “velocity pullback”), for which usually the linear acoustic
approximation is used. In Sect. IV we develop a new variant of the post-
acoustic approximation and derive a new formula relating $\sigma_{sp}$ to
$\Delta u_{fs}$ in the case of a symmetric plate impact. The new formula is
significantly more accurate than the acoustic one, and provides an important
insight into the complex interplay of the first-order nonlinear correction
terms. New post-acoustic formulae are also obtained for the strain rate and
the fractured mass fraction.

##  II The phase-flip model of spall fracture in liquids

###  II.1 The employed two-phase EOS

Conceptually, the thermodynamics of a liquid-gas phase transition is
adequately represented by the classical van der Waals EOS. In this work we
make use of one of the simplest generalizations of this EOS (termed GWEOS)
Martynyuk (1991, 1993); Basko (2018b) in the form

$\displaystyle p(v,\theta)$ $\displaystyle=$
$\displaystyle\frac{\alpha\theta}{v-\kappa^{-1}}-\frac{\kappa}{v^{n}},$ (1)
$\displaystyle e(v,\theta)$ $\displaystyle=$ $\displaystyle\alpha
c_{V}\theta-\frac{1}{2}\kappa(\kappa-1)v^{1-n},$ (2) $\displaystyle
s(v,\theta)$ $\displaystyle=$
$\displaystyle\alpha\left[c_{V}(1+\ln\theta)+\ln(v-\kappa^{-1})\right],$ (3)

where $\kappa=(n+1)/(n-1)$, $\alpha=\kappa-\kappa^{-1}$, and the exponent
$n>1$ together with a constant heat capacity $c_{V}>0$ are its two free
dimensionless parameters. The EOS (1)–(3) is cast in its reduced form, where
the pressure $p$, the specific volume $v\equiv\rho^{-1}$ and the temperature
$\theta$ are normalized to the corresponding values $P_{cr}$,
$V_{cr}=\rho_{cr}^{-1}$ and $T_{cr}$ at the liquid-gas critical point; the
mass-specific internal energy $e$ is in units of $P_{cr}V_{cr}$, the entropy
$s$ in units of $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. Basko, 2018b. The main motivation for using 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 hydrodynamic codes
with the rounding-error accuracy — thus excluding numerical errors due to
imperfect EOS.

In numerical examples, discussed below, the few free parameters of GWEOS were
chosen such as to provide the best fit to the properties of water near normal
conditions. The well known critical parameters of water allow confident
normalization of its other measured quantities to the system of reduced
variables. For simplicity, we assume the normal state $(v_{0},\theta_{0})$ to
have zero pressure $p(v_{0},\theta_{0})=0$. The values of $n$, $c_{V}$ and
$\theta_{0}$ are adjusted to fit as best as possible the experimental values
of the normal specific volume $v_{0}$, the isentropic sound speed $c_{0}$, the
vaporization enthalpy $h_{vap}$, and the specific heat $\alpha c_{V}$ at
$T=20^{\circ}$ C. Table 1 shows how good a compromise has been achieved with
the values

$n=1.5,\quad c_{V}=8.18,\quad\theta_{0}=0.45,$ (4)

used throughout this paper.

Table 1: Comparison between the normal properties of water Haynes (2015) and those produced by GWEOS with parameters (4). For all quantities the reduced values are given | $\theta_{0}$ | $\rho_{0}$ | $c_{0}$ | $h_{vap}$ | $\alpha c_{V}$ | $b_{0}$
---|---|---|---|---|---|---
GWEOS | 0.450 | 3.9073 | 5.9939 | 21.5 | 39.26 | 3.17
H2O | 0.453 | 3.10 | 5.66 | 35.8 | 39.26 | 2.0

As is shown in Fig. 1, Eq. (1) exhibits all the standard features of the vdW
EOS. The phase coexistence region lies between the _binodal_ curve $bi$ and
the $p=0$ line. The EQ EOS in this region, denoted $p_{EQ}(v,\theta)$,
$e_{EQ}(v,\theta)$, $s_{EQ}(v,\theta)$, and the position of the binodal itself
are determined by applying the Maxwell rule (Landau and Lifshitz, 1996,
§§84,85) to the MS isotherms (1), as is illustrated with the horizontal
segment of the $\theta=0.93$ isotherm in Fig. 1. The _spinodal_ curve $sp$,
defined by the condition $\partial p(v,\theta)/\partial v=0$, lies under the
binodal and delimits the region of absolute thermodynamic instability; hence,
only the EQ states with $p>0$, obtained by applying the Maxwell rule, are
possible below the spinodal. In the intermediate metastable region between the
binodal and the spinodal both the MS and the EQ branches of EOS are
thermodynamically admissible and compatible with the equations of fluid
dynamics. The latter is ensured by the positiveness of the square of the
isentropic sound speed

$c^{2}=\left(\frac{\partial p}{\partial\rho}\right)_{s}>0$ (5)

everywhere along and above the spinodal curve. Tension states with $p<0$ are
only possible for the MS states in the superheated-liquid region at $v<1$.
Note that, by assuming $p(v_{0},\theta_{0})=0$, we put the normal state
$(v_{0},\theta_{0})$ into the category of metastable ones.

Figure 1: Thermodynamic $(v,p)$ plane of a fluid with a liquid-gas phase
transition. Shown are the binodal $bi$ and the spinodal $sp$ curves with the
critical point $CP$ at $(v,p)=(1,1)$. Curve $PBO_{1}S$ shows an MS isentrope,
along which the release-wave states would evolve after being brought to $P$
along the Hugoniot $OP$ (dashed) starting from the normal state $O$. Dash-
dotted curve $h_{sp}$ (orange) is the locus of isochoric post-flip states $E$
as the initial state $S$ slides along the spinodal $sp$; dash-dotted curve
$FEG$ (red) is the EQ isentrope passing through the chosen state $E$. The
behavior of the MS (solid green) and EQ (dashed green) isotherms is
illustrated for $\theta=0.93$

A non-trivial point is the interpretation of the EQ states from the viewpoint
of matter structure. While there is no ambiguity about the metastable liquids,
which can always be assumed homogeneous, the EQ states are normally believed
to be necessarily heterogenous (Landau and Lifshitz, 1996, §83). On the other
hand, no length scale for their heterogeneity can be derived from our basic
thermodynamic relations (1), (2). Without extra elements of physics (like the
kinetics of bubble nucleation, surface tension, impurities, etc.), one can
only assume that an EQ state represents a quasi-uniform finely dispersed
liquid-gas mixture, where the characteristic sizes of either vapor bubbles or
liquid drops are always much smaller than the smallest relevant length scale
in the considered problem [see, for example, problem 1 to § 64 in Landau and
Lifshitz (1987)]. Below we apply the term “fog” to such EQ states and treat
them as perfectly homogeneous in the context of spallation and cavitation in
pure liquids.

###  II.2 Spall criterion and the phase-flip approximation

Spall failure inside (i.e. away from boundaries) a volume of pure liquid
occurs in the process of spontaneous decay of its metastable state due to
thermodynamic fluctuations Grady (1988). This implies that, other factors
being equal, the spall strength $\sigma_{sp}$ must increase with the
increasing strain rate, which is confirmed by many experiments Kanel, Fortov,
and Razorenov (2007); Kanel (2010) and the MD simulations Mayer and Mayer
(2020). Hence, the theoretical maximum of $\sigma_{sp}$ — the _ultimate
strength_ — must be attained at the spinodal

$\sigma_{sp}=-p_{sp}(v)=\frac{n-v(n+1)}{v^{n+1}},\quad\frac{n-1}{n+1}<v<\frac{n}{n+1},$
(6)

where the metastability decay rate is maximum Imre  _et al._ (2008).

In this work we adopt a simple and universal criterion that a given fluid
element undergoes spall failure whenever its thermodynamic trajectory reaches
the spinodal from the side of metastable liquid. The results obtained under
this assumption should be in the first place applicable to tensile loads with
the highest possible strain rates, like, for example, $\dot{\varepsilon}\simeq
10^{9}$–$10^{10}$ s-1 observed in laser experiments. In reality, of course,
the tensile strength can only approach the spinodal limit from below Imre  _et
al._ (2008), but how closely — about 0.3–0.5 Kanel, Fortov, and Razorenov
(2007); Agranat  _et al._ (2010) or up to 0.8–0.9 Pisarev  _et al._ (2009);
Malyshev  _et al._ (2015) of the limiting value — is still debatable. Given
this uncertainty and the uncertainty in the spinodal position for real
materials, we do not try to soften this criterion or make it dependent on
$\dot{\varepsilon}$, the more so that it would not affect the main results of
this work.

Having set the pressure threshold for spallation, we need to add further
assumptions about the kinetics and morphology of the ensuing fracture. Because
the rate of spontaneous bubble nucleation very rapidly increases by many
orders of magnitude as the MS liquid approaches the spinodal Blander and Katz
(1975); Skripov and Skripov (1979), we can make a simplifying assumption that
spall failure occurs as an instantaneous and irreversible MS $\rightarrow$ EQ
phase jump — the _phase flip_. The extremely short timescale, on the order of
a few picoseconds, of spinodal decomposition in real liquids justifies this
approximation for strain rates up to $\dot{\varepsilon}\simeq 10^{11}$ s-1.
Because the rate of density variation $\dot{\varepsilon}$ is controlled by
hydrodynamics (i.e. by fluid inertia and pressure gradients), the phase flip
must occur at a constant $\rho=v^{-1}$ whenever $\dot{\varepsilon}$ is finite.
The only possibility where it could be accompanied by a density jump, would be
inside a rarefaction shock Basko (2018a), which can be excluded in the context
of spallation because relaxation of tension means growth of pressure. In
addition, energy conservation requires the local value of the specific
internal energy $e$ to be preserved as well, which implies that the phase flip
is always accompanied by a jump-like increase in pressure, temperature, and
entropy.

In Fig. 1 the trajectory $PBO_{1}SE$ illustrates possible evolution of the
thermodynamic state of a fluid element, undergoing spall failure: expansion
along the MS isentrope $PBO_{1}S$ brings it to a negative-pressure state $S$
on the spinodal, followed by an instantaneous jump to the EQ “fog” state $E$
at a higher entropy; its subsequent evolution is governed by the EQ EOS branch
so long as it remains under the binodal. If the tracked mass of “fog” emerges
from under the binodal (say, by recompression along the EQ isentrope $EF$), it
may later return to the two-phase region again as metastable liquid and
undergo a secondary spall failure. Thus, the fluid elements in our scheme can
exhibit a hysteresis-like thermodynamic behavior that is fully consistent with
the first and second laws of thermodynamics.

Because metastable liquid can border vacuum at a finite density, one might
conjecture that it would rupture by forming an expanding vacuum cavity,
originating from the _inception point_ where the spall threshold is first
attained. To assess the feasibility of such scenario, we have to consider the
compression shock pulse that would normally be launched into the surrounding
stretched liquid by pressure relaxation in the fractured material (often
called the “spall pulse”). Because this shock relaxes the tension in
metastable liquid and pulls its thermodynamic state back from the spall
threshold, we call it the _pull-back (pb) shock_. For an expanding vacuum
cavity to form, the pb shock would have to start immediately from the
inception point. However, a closer examination in the next section reveals a
different general pattern: once the spall threshold is reached at the
inception point, the region of spall failure very rapidly (hypersonically)
extends over some finite volume, bounded by a surface from where the pb shock
does actually start, and which may be called the “spall horizon” of a given
inception point. Due to the continued stretching, sustained by inertia of the
fluid motion, no fluid element inside the spall horizon can be reached by the
pb shock before it undergoes spall failure. In other words, the spall fracture
in our model develops as a rapid phase transition between the MS and EQ EOS
branches, occurring practically simultaneously throughout a certain finite
mass. Since an EQ state can border vacuum only at zero density, there is no
need to postulate any void openings, and the spalling (cavitating) liquid can
be treated as everywhere continuous.

At first glance, the above spall criterion appears similar to that proposed
earlier by GradyGrady (1988, 1996), who also assumed an instantaneous MS
$\rightarrow$ EQ phase jump in the framework of 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 internal energy $e_{EQ}(v,s)<e(v,s)$ than the MS isentrope at
the same density, Grady interpreted the excess energy $\Delta
e_{NE}=e(v,s)-e_{EQ}(v,s)$ as the interfacial energy (extraneous to EOS) of
the heterogeneous EQ state, consisting of liquid fragments with a
characteristic size $l_{fr}$; he then used $\Delta e_{NE}$ and $l_{fr}$ to
evaluate the spall strength by invoking the surface tension and viscosity. Our
model, in contrast, is based on the assumption that $l_{fr}$ is very small and
cannot be resolved hydrodynamically. Then the eventual interfacial energy must
be automatically included into the fluid EOS, which requires the MS
$\rightarrow$ EQ phase jump to be iso-energetic and accompanied by the entropy
increase.

Another highly controversial point in Grady’s argumentation is the postulated
linear dependence

$l_{fr}=2ct_{w}$ (7)

of the typical fragment size $l_{fr}$ on the “waiting” time
$t_{w}\propto\dot{\varepsilon}^{-1}$ elapsed after crossing the binodal. Here
we argue that such a linear relation would be in stark disagreement with the
theory of homogeneous nucleation Skripov (1974); Skripov and Skripov (1979),
which predicts that the new phase in metastable liquids under tension $-p>0$
appears in the form of vapor bubbles with the critical radius

$r_{cr}=\frac{2\alpha_{st}}{p_{sat}-p}\approx-\frac{2\alpha_{st}}{p},$ (8)

where $\alpha_{st}$ is the surface tension, and $p_{sat}=p_{sat}(T)>0$ is the
saturated vapor pressure at the same temperature $T$ for which $p$ is
calculated Landau and Lifshitz (1996); Skripov (1974); Blander and Katz
(1975). Because the probability of spontaneous creation of such bubbles,
proportional to $\exp(-G)$ where

$G=\frac{16\pi\alpha_{st}^{3}}{3(p_{sat}-p)^{2}T}=\frac{4\pi\alpha_{st}r_{cr}^{2}}{3T}\gg
1$ (9)

is the Gibbs number, is a very steep function of $p$ and $T$, the spall
fracture of a uniformly stretched liquid occurs as a sudden emergence of
myriad critical vapor bubbles with radii of a few nanometers Skripov and
Skripov (1979); Pisarev  _et al._ (2009), which practically instantaneously
(i.e. with virtually no time for subsequent growth) begin to coalesce,
breaking the molecular bonds that sustain tension. Within this picture, it is
expected that $l_{fr}\approx r_{cr}$ and $t_{w}\propto\exp(G)$, which,
together with Eq. (9), implies a highly nonlinear relation between $l_{fr}$
and $t_{w}\propto\dot{\varepsilon}^{-1}$, and significantly smaller values of
$l_{fr}$ than predicted by Eq. (7). Therefore, when the goal is to explore the
dependence of $\sigma_{sp}$ on the strain rate $\dot{\varepsilon}$, the simple
relationship (7) must be replaced by a more adequate model based on the
nucleation theory — as, for example, was proposed in Dekel  _et al._ (1998);
Pisarev  _et al._ (2009); Faik  _et al._ (2012).

In summary, unlike Grady’s scheme, our model describes energy dissipation due
to spall fracture without involving the surface tension and viscosity, simply
as the entropy increase in the phase-flip transitions and behind the pull-back
shock fronts. In this respect it is similar to the viscous-free ideal
hydrodynamics with dissipation due only to discontinuous shocks. The weak
dependence of the spall strength on the strain rate is ignored by adopting the
theoretically maximum value of $\sigma_{sp}$ allowed by thermodynamics.

###  II.3 Kinematics of spall fronts and the size of fracture zone

Consider a flow of liquid with a velocity field $\mathbf{u}(t,\mathbf{x})$,
where a certain region is stretched to develop negative pressures. The local
rate of stretching is characterized by the _strain rate_ , defined as

$\dot{\varepsilon}\,\stackrel{{\scriptstyle\mbox{\scriptsize
def}}}{{=}}\,-\frac{d\ln\rho}{dt}=\frac{d\ln v}{dt}=\nabla\cdot\mathbf{u},$
(10)

where $d/dt=\partial/\partial t+\mathbf{u}\cdot\nabla$ is the material time
derivative, and the use is made of the fluid continuity equation

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mathbf{u}\right)=0.$
(11)

Next, we assume for simplicity that the flow is isentropic with a fixed value
of entropy $s$ throughout the considered volume. Then the pressure becomes a
known monotonic function $p(v)=p(v,s)$ of a single variable $v$ with $dp/dv<0$
everywhere above the spinodal (curve $PBO_{1}S$ in Fig. 1); the spall strength
$\sigma_{sp}=-p(v_{sp})$ and the spall density $\rho_{sp}=v_{sp}^{-1}$ take on
universal values that are determined by the intersection of the isentrope
$p(v)$ with the liquid branch of the spinodal (point $S$ in Fig. 1).

Figure 2: Schematics of the spall front propagation. As the local pressure
minimum is forced to drop ever deeper below the spall threshold $p_{sp}$, the
position of the spall front shifts from the inception point $A$ at $t_{sp}$ to
points $B$ and $C$ at later times $t$ and $t+\Delta t$

Spall failure is initiated at time $t_{sp}$ when the absolute minimum in the
spatial pressure distribution $p(t_{sp},\mathbf{x})$ attains the threshold
value $p_{sp}=-\sigma_{sp}$ at some inception point $A$, as is illustrated in
Fig. 2. Generally, when $p(t,\mathbf{x})$ and $p(v)$ are differentiable at
$A$, one has $\nabla p=\nabla\rho=0$ at this point. At later times $t>t_{sp}$
the condition $p(t,\mathbf{x})=p_{sp}$ will be satisfied on a certain surface
around the inception point, which we call the _spall front_ and which
encompasses a fog bubble. As suggested by the term, all the liquid inside the
fog bubble has been, upon passing through the spall front, instantaneously
transformed into an EQ “fog” state at a positive pressure. The condition
$\nabla p=0$ by inception implies that the expansion of the fog bubble starts
infinitely fast, but gradually slows down as $\left|\nabla p\right|$
increases. Because $p$ is a single-valued monotonically growing function of
$\rho$, the normal (i.e. along the colinear $\nabla p$ and $\nabla\rho$
vectors) kinematic velocity of the spall front

$u_{spf}=\frac{\mathbf{u}\cdot\nabla\rho}{\left|\nabla\rho\right|}+\frac{\dot{\varepsilon}}{\left|\nabla\ln\rho\right|}$
(12)

is readily calculated by differentiating the equation
$\rho(t,\mathbf{x})=constant$, as is illustrated in Fig. 2. Note that in the
comoving reference frame, where $\mathbf{u}=0$, the front velocity $u_{spf}$
is given by the second term on the right-hand side of Eq. (12).

The merely kinematic propagation of the spall front with velocity (12)
continues until it is overtaken by the pb shock, generated by the pressure
jump $\Delta p_{pb}=p_{EQ,sp}+\sigma_{sp}$ behind the spall front; here
$p_{EQ,sp}>0$ is the fog pressure in the post-flip EQ state (point $E$ in Fig.
1), calculated from the EQ EOS branch for the same values of $v_{sp}$ and the
internal energy $e_{sp}$ as in the pre-flip MS state $(v_{sp},p_{sp})$ (point
$S$ in Fig. 1). The pb-shock front is locally coplanar with the spall front
and propagates in the same direction with the normal velocity

$u_{pbf}=\frac{\mathbf{u}\cdot\nabla\rho}{\left|\nabla\rho\right|}+D_{pb},$
(13)

where $D_{pb}=D_{pb}(s)$ is the velocity of the shock front, driven by the
pressure jump $\Delta p_{pb}$, relative to the liquid in the pre-shock state
$(v_{sp},p_{sp})$; it is, of course, calculated by using the MS EOS branch.
Similar to $v_{sp}$, $\sigma_{sp}$, and $\Delta p_{pb}$, the shock velocity
$D_{pb}$ is uniquely determined by the entropy value $s$.

Once the condition $u_{pbf}\geq u_{spf}$ is met and the pb shock breaks out
before the spall front, it recompresses and heats up the stretched liquid,
preventing it from reaching the spall threshold. As a result, the expansion of
the fog bubble in Lagrangian coordinates, i.e. relative to the fluid
particles, comes to a halt. In other words, the mass of fractured liquid
continues to grow so long as the effective _spall Mach number_ , defined as

$\mathcal{M}_{sp}\,\stackrel{{\scriptstyle\mbox{\scriptsize
def}}}{{=}}\,\left(\frac{u_{spf}}{u_{pbf}}\right)_{\mathbf{u}=0}=\frac{\dot{\varepsilon}}{D_{pb}\left|\nabla\ln\rho\right|}=\frac{\dot{\varepsilon}\rho
c^{2}}{D_{pb}\left|\nabla p\right|},$ (14)

stays above unity, $\mathcal{M}_{sp}>1$. Note that $\mathcal{M}_{sp}$ is
defined in the comoving frame, and relative not to the local sound speed $c$
but to the pb-shock velocity $D_{pb}>c$; hence, we use the term _hypersonic_
to a spall front with $\mathcal{M}_{sp}>1$. The condition
$\mathcal{M}_{sp}=1$, when fulfilled on the spall-front surface, delimits the
size and mass of the fog bubble. The bubble mass is always finite when
$\dot{\varepsilon}>0$ and $\nabla\rho=\nabla p=0$ at the inception point. It
may become zero in some singular cases where either $\dot{\varepsilon}$
changes sign or $\nabla p$ is discontinuous at inception. In practice, a fair
estimate of the fractured mass can often be made by analyzing the spatial
distribution of $\mathcal{M}_{sp}(t_{sp},\mathbf{x})$ at the moment of
inception $t_{sp}$. The above arguments and formulae can be readily
generalized to non-isentropic flows, where $\mathcal{M}_{sp}$ becomes a
function of gradients of two independent thermodynamic variables.

##  III Numerical results for a shock-loaded planar target

###  III.1 Formulation of the problem and the initial state

In a common experiment for determination of the spall strength, two planar
plates are collided to launch the initial load shocks in two opposite
directions from the impact plane Antoun  _et al._ (2003). In this work we
focus our attention on the simplest version of such an experiment, where both
plates have the same composition and dimensions. Then the problem is one-
dimensional (1D) and symmetrical with respect to the impact plane, which we
place at the origin $x=0$ of the center-of-mass coordinate system and consider
only one (right) plate at $x\geq 0$, which we refer to as the target. We set
$t=0$ at the moment when the load shock breaks out at the outer edge
$x=h_{1}$, leaving behind the shocked material in a compressed state at zero
initial velocity. Thus defined initial state $(v_{1},p_{1})$ (point $P$ in
Fig. 1) lies on the Hugoniot $OP$, originating from the normal state
$(v_{0},p_{0})$ at point $O$ in Fig. 1. The intensity of the initial load is
characterized by the particle (piston) velocity $u_{p}$, measured relative to
the fluid in front of the loading shock. The latter means that, before the
impact, the flyer and the target plates had, respectively, the velocities
$+u_{p}$ and $-u_{p}$ in the center-of-mass frame; their density and thickness
were $\rho_{0}$ and $h_{0}$. The value of $u_{p}$ uniquely determines the
entropy $s$ throughout the subsequent flow, the spall strength $\sigma_{sp}$,
as well as the pressure amplitude $\Delta p_{pb}$ and the front velocity
$D_{pb}$ of the pb shock. Table 2 lists the key flow parameters that can be
calculated directly from EOS for two chosen reference cases of the shock load.

Table 2: Shock-load and subsequent flow parameters that can be calculated directly from EOS for two selected values of the post-shock particle velocity $u_{p}$. For detailed explanations see the text $u_{p}$ | $u_{1}$ | $\rho_{1}=v_{1}^{-1}$ | $c_{1}$ | $\rho_{01}$ | $c_{01}$ | $b$ | $z_{0}$ | $p_{EQ,sp}$ | $D_{pb,sp}$ | $c_{sp}$ | $\sigma_{sp}$
---|---|---|---|---|---|---|---|---|---|---|---
1.25 | 1.2170 | 4.4470 | 14.351 | 3.8928 | 5.9329 | 3.1405 | 1.0140 | 0.00256 | 4.144 | 1.266 | 11.149
4.0 | 3.7495 | 4.7919 | 46.732 | 3.6309 | 4.9757 | 2.6954 | 1.2963 | 0.01302 | 3.635 | 1.241 | 8.2190

Before specifying the initial target thickness $h_{1}$ (or its thickness
$h_{0}$ before the shock load), we remark the following. One easily verifies
that the above formulated problem is invariant with respect to rescaling of
all the lengths and times by one and the same arbitrary factor. Therefore,
having obtained the full solution for an initial target thickness $h_{1}$, we
can rescale it to any other thickness $h_{1}^{\prime}$ by multiplying times
and lengths by the factor $h_{1}^{\prime}/h_{1}$. In view of such scalability,
all the simulations discussed below were done for $h_{1}=1$. Recalling that
$v$, $p$, and $u$ are, respectively, everywhere in units of $V_{cr}$,
$P_{cr}$, and $\left(P_{cr}V_{cr}\right)^{1/2}$, the latter is equivalent to
setting the units of length and time to $h_{1}$ and
$h_{1}\left(P_{cr}V_{cr}\right)^{-1/2}$.

###  III.2 Simulation results for $u_{p}=1.25$ and $u_{p}=4.0$

Given an adequate two-phase EOS, implementation of our spall criterion into a
standard Lagrangian hydrocode is straightforward. Here we present numerical
results, obtained with the 1D DEIRA code that has been previously developed
and extensively used for simulations of inertial confinement fusion targets
Basko (2001, 1990). Figures 3 and 5 display the target density evolution on
the $(x,t)$ plane calculated, respectively, for $u_{p}=1.25$ and $u_{p}=4.0$.
For a detailed discussion, we focus on the $u_{p}=1.25$ case, which lies
moderately above the spall threshold of $u_{p,sp}\approx 0.82$. The high-load
case of $u_{p}=4.0$ mainly serves to make the nonlinear effects more
conspicuous

As the centered rarefaction wave, starting at $x=1$, arrives at and is
reflected from the symmetry plane $x=0$, a tension region with $\rho<3.90$
develops near the target center at $t\gtrsim 0.12$. Further on, as the
pressure in this region falls to $p=-\sigma_{sp}=-11.15$, all the liquid mass
at $0<x\lesssim 0.19$ undergoes spall fragmentation within a short time
interval of $t_{sp}=0.21819<t<0.220$, making up the spall zone that continues
to expand preserving its mass. A salient feature is the pb shock, launched
from the outer boundary of the spall zone at $t=0.220$, and giving a kick to
the free surface velocity later at $t=0.434$.

Figure 3: Color density map on the space-time diagram as calculated with the
1D DEIRA code for $u_{p}=1.25$ Figure 4: Profiles of the strain rate
$\dot{\varepsilon}$, the spall Mach number $\mathcal{M}_{sp}$ (left
ordinates), and of the pressure $p$ (right ordinates) versus the fractional
target mass $\bar{m}$ just before the spall failure for $u_{p}=1.25$. Only the
inner half of the target mass is shown

More insight into the spall dynamics can be gained by examining the profiles
in Figs. 4 and 6 along the dimensionless Lagrangian coordinate

$\bar{m}=\frac{m}{m_{0}},\quad
m=\\!\int\limits_{0}^{x}\rho(t,x^{\prime})\,dx^{\prime},\quad
m_{0}=\rho_{0}h_{0}=\rho_{1}h_{1},$ (15)

plotted for the time moment immediately preceding $t_{sp}$. (Note that these
profiles are invariant with respect to rescaling of the target thickness.) In
Fig. 4 the spall failure occurs at $t_{sp}=0.21819$ in the center $\bar{m}=0$,
where the pressure is minimum and $\mathcal{M}_{sp}=\infty$. The outer limit
of the spall zone is clearly marked by a precipitous drop of the effective
spall Mach number $\mathcal{M}_{sp}$, which crosses the $\mathcal{M}_{sp}=1$
line at $\bar{m}_{k}=0.135$; the actual fractured mass, established at a
slightly later moment $t=0.220$, has a somewhat higher value of
$\bar{m}_{sp}=0.131$. In Fig. 6 the respective numbers are $\bar{m}_{k}=0.572$
and $\bar{m}_{sp}=0.568$.

Figure 5: Same as Fig. 3 but for $u_{p}=4.0$ Figure 6: Same as Fig. 4 but for
$u_{p}=4.0$

Two more observations to be made from Figs. 4 and 6 are (i) a virtually flat
pressure profile across the incipient spall zone $0<\bar{m}<\bar{m}_{k}$, and
(ii) a jump by about a factor 2 in the value of the strain rate
$\dot{\varepsilon}$ at $\bar{m}=\bar{m}_{k}$. Here we must recall that the
exact solution to our problem contains weak discontinuities (see problem 1 to
§ 105 in Landau and Lifshitz (1987)), i.e. discontinuities of $\partial
p/\partial x$, $\partial\rho/\partial x$, $\partial u/\partial x$ that are
smeared by the artificial viscosity in numerical profiles of Figs. 3–6. In
Figs. 4 and 6 the kink point $\bar{m}_{k}$ in the pressure profile is just
such a discontinuity. As a consequence, the exact profiles of
$\mathcal{M}_{sp}$ and $\dot{\varepsilon}$ would both have jumps at
$\bar{m}=\bar{m}_{k}$. In our simulations the threshold value of
$\mathcal{M}_{sp}=1$ is bracketed with a wide margin by its pre- and post-jump
values. Due to nonlinearity of the problem, the pressure profile at
$\bar{m}<\bar{m}_{k}$ is not exactly flat: it is monotonic with a minimum at
$\bar{m}=0$ for $u_{p}=1.25$, and non-monotonic with a minimum at the kink
point $\bar{m}_{k}$ for $u_{p}=4.0$. That is, the exact position of spall
initiation may differ depending on the EOS details and the load intensity. The
narrow peak of $\mathcal{M}_{sp}$ just before $\bar{m}=\bar{m}_{k}$ in Fig. 6
is caused by numerical smoothing of this corner point in the $p$ and $u$
profiles, and would not occur in the exact solution.

###  III.3 On morphology of the spall zone

The space-time image of Fig. 3 offers a convenient opportunity to make a
comparison with the published MD simulation of a similar problem for liquid
copper Cai, Wu, and Luo (2017). When rescaled to the critical point of copper
Fortov, Iakubov, and Khrapak (2006), the parameters of our $u_{p}=1.25$ run
are close to the MD case shown in Fig. 1 of Ref. Cai, Wu, and Luo, 2017. The
overall spall dynamics, the final size, and the qualitative structure of the
spall zone look also quite similar. The spall strength $\sigma_{sp}\approx
11P_{cr}$, calculated in Ref. Cai, Wu, and Luo, 2017, practically coincides
with the spinodal limit given in Table 2. Significant differences emerge only
when the structure of fractured liquid is examined in more detail.

In the results of Ref. Cai, Wu, and Luo, 2017, the spall zone has a width of
some 600 interatomic distances and the fractured mass splits into three liquid
drops separated by voids. Generally, one would expect the MD approach to
produce some statistical distribution of liquid fragments, whose typical sizes
would be determined by the atomic length scale, inherent in any MD model, and
by statistical properties of the microscopic initial states used in the MD
code runs.

The present model, in contrast, is expected to produce a fully deterministic
pattern of spall fragmentation because it is based on the equations of ideal
hydrodynamics and contains no length scale to stipulate a specific fragment
size. According to the arguments of Sect. II.3, the entire spall zone in the
exact solution to our problem should look like a single fog bubble without any
remaining liquid drops because the condition $\mathcal{M}_{sp}\gg 1$ is
fulfilled everywhere throughout it. The many small droplets, observed in the
spall zones of our Figs. 3 and 5, are, in fact, of spurious origin due to a
short-wavelength (on the cell-size scale) numerical noise, generated by the
simple numerical scheme of the DEIRA code when phase-flip pressure jumps occur
at discrete time moments. As the mesh is refined, an ever more fine pattern of
alternating liquid and fog cells develops, and no convergence to the exact
solution is observed.

The above predicament for numerical modeling is closely related to the fact
that our base model can ensure only conditional stability of the fragmentation
morphology with respect to small perturbations. Because the spall Mach number
(14) is controlled by the gradients of principal variables, the overall
fragmentation pattern (i.e. the number, spacing and sizes of separate liquid
fragments) becomes stable with respect to small perturbations when not only
the variations of $\rho$, $p$, and $\mathbf{u}$ but those of their spatial
derivatives are kept sufficiently small — hence the strong sensitivity of the
fragmentation morphology to small perturbations with high wave numbers.

As an illustrative example, Fig. 7 displays the same case as in Fig. 3 but
simulated with a perturbed initial density

$\rho(0,x)=\rho_{1}\left[1-0.001\cos(46\pi x)\right]$ (16)

under the unaltered pressure $p(0,x)=p_{1}$. One sees that even as small as a
0.1% density perturbation results in a dramatic change of the spall pattern:
instead of a single fog bubble we have three liquid drops separated by four
fog bubbles. Once well resolved by the mesh (43 cells per each drop in Fig.
7), the mass and other dynamic properties of these three drops remain stable
with respect to further mesh refinement.

Figure 7: Same as Fig. 3 but calculated with the initial density perturbation
(16)

In summary, being fully deterministic for given initial and boundary
conditions, our model should predict realistic morphology of the spall zone
when realistic small perturbations are included into the problem formulation.
By numerical implementation, a special care may be needed to suppress spurious
short-wavelength perturbations due to finite differencing.

##  IV First post-acoustic approximation for a symmetric planar spall

In most plate impact experiments, the spall strength $\sigma_{th}$ is inferred
from the time history of the free surface velocity $u_{fs}(t)$ measured at the
outer target edge Antoun  _et al._ (2003); Kanel, Fortov, and Razorenov
(2007). Figure 8 shows the plots of $u_{fs}(t)$, calculated for two cases of
the initial load in the above stated problem, one below ($u_{p}=0.8$) and the
other above ($u_{p}=1.25$) the spall threshold. Given the density $\rho_{0}$
and the sound speed $c_{0}$ of the normal state, $\sigma_{th}$ is most often
evaluated in the linear acoustic approximation Dekel  _et al._ (1998); Antoun
_et al._ (2003) from the formula

$\sigma_{th,ac}=\frac{1}{2}\rho_{0}c_{0}\Delta u_{fs},$ (17)

where

$\Delta u_{fs}=u_{fs,max}-u_{fs,min}$ (18)

is the amplitude of the first dip in the $u_{fs}(t)$ profile. However,
contrary to what might be expected, formula (17) never becomes asymptotically
accurate even in the limit of an infinitely weak load because of the pb-shock
attenuation: below we show its error to be 40% for the considered case of
symmetric plate collision. The situation becomes even worse when the acoustic
formulae Dekel  _et al._ (1998); Cai, Wu, and Luo (2017)

$\dot{\varepsilon}_{sp,ac}=\frac{\Delta
u_{fs}}{2c_{0}(t_{pb}-t_{2})},\quad\bar{m}_{sp,ac}=1-\frac{c_{0}t_{pb}}{2h_{0}}$
(19)

for the spall strain rate $\dot{\varepsilon}_{sp}$ and the fractured mass
fraction $\bar{m}_{sp}$ are invoked: the error in $\dot{\varepsilon}_{sp}$ is
often a factor 3–5, while absurd negative values are easily obtained for
$\bar{m}_{sp}$. In this section we go to the next order in the nonlinear terms
and derive new formulae in what we call the _first post-acoustic_ (PA)
approximation, which provide significantly better accuracy than Eqs. (17) and
(19) for the respective quantities.

Figure 8: Time dependence of the free-surface velocity $u_{fs}(t)$ at the
outer target edge $\bar{m}=1$, calculated with the DEIRA code for two cases of
the initial shock load: one below ($u_{p}=0.8$, blue dashes) and the other
above ($u_{p}=1.25$, black solid) the spall threshold

###  IV.1 Characteristics chart

As the loading shock arrives at the outer target edge $\bar{m}=1$ at $t=0$, a
rarefaction wave R1 is launched in the negative direction towards the symmetry
plane $\bar{m}=0$, as is schematically shown in Fig. 9a. The fluid before the
R1 front is at rest in the shock-compressed state $\rho_{1}=v_{1}^{-1}$,
$p_{1}=p(\rho_{1},s_{1})$ (point $P$ in Fig. 1); the fluid behind it moves
with a constant positive velocity

$u_{1}=\int\limits_{0}^{p_{1}}\frac{dp}{\rho
c}=\int\limits_{\rho_{01}}^{\rho_{1}}c\,d\ln\rho$ (20)

and is in the thermodynamic state $\rho_{01}$, $p_{01}\equiv
p(\rho_{01},s_{1})=0$ (point $O_{1}$ in Fig. 1). Here and below we assume that
the MS EOS (1)-(3) is cast in the form $p=p(\rho,s)$. So long as no new shocks
emerge and the entropy remains at its post-load-shock value $s=s_{1}$, the
pressure $p=p(\rho,s_{1})=p(\rho)$ and the sound speed
$c=(dp/d\rho)^{1/2}=c(\rho)$ are functions of $\rho$ only.

Figure 9: Schematic velocity and pressure profiles in the left-propagating R1
(a) and right-propagating R2 (b) rarefaction waves. The right-propagating R2
wave is “propped up” from behind by the pull-back shock

The key to further analysis is the behavior of characteristics for the two
governing hydrodynamic equations, whose Lagrangian form in our case is

$\displaystyle\frac{\partial v}{\partial t}-\frac{\partial u}{\partial m}$
$\displaystyle=$ $\displaystyle 0,$ (21) $\displaystyle\frac{\partial
u}{\partial t}+\frac{\partial p}{\partial m}$ $\displaystyle=$ $\displaystyle
0,$ (22)

where the mass coordinate $m$ and its normalized counterpart $\bar{m}$ are
defined in Eq. (15). In terms of the normalized acoustic impedance

$z=z(\rho)\,\stackrel{{\scriptstyle\mbox{\scriptsize def}}}{{=}}\,\frac{\rho
c}{\rho_{01}c_{01}},\quad z_{0}=\frac{\rho_{0}c_{0}}{\rho_{01}c_{01}},$ (23)

and the normalized time

$\bar{t}\,\stackrel{{\scriptstyle\mbox{\scriptsize
def}}}{{=}}\,\frac{\rho_{01}c_{01}}{\rho_{0}h_{0}}\,t=\frac{c_{0}t}{h_{0}z_{0}},$
(24)

the two families of $C^{\pm}$ characteristics on the $(\bar{m},\bar{t})$ plane
are given by

$\frac{d\bar{m}}{d\bar{t}}=\left\\{\begin{array}[]{ll}+z,\;\;C^{+}:\,dJ^{+}\equiv
du+dp/\rho c=0,\\\ -z,\;\;C^{-}:\,dJ^{-}\equiv du-dp/\rho
c=0,\end{array}\right.$ (25)

where $J^{\pm}$ are the two respective Riemann invariants Landau and Lifshitz
(1987).

Figure 10 shows the pattern of $C^{\pm}$ characteristics, representing the
solution to our problem in the normalized $(\bar{m},\bar{t})$ plane. The R1
wave starts at point 1 with $(\bar{m},\bar{t})=(1,0)$ as a centered $C^{-}$
fan between the leading $C_{h1}^{-}$ and the trailing $C_{t1}^{-}$ lines. As
each of these $C^{-}$ characteristics reaches the symmetry plane $\bar{m}=0$,
it is reflected and continues as a $C^{+}$ characteristic, i.e. $C_{h2}^{+}$
is the reflected continuation of $C_{h1}^{-}$. Where the R1 fan borders a
region of constant flow, it is a _simple_ wave with all the $C^{-}$
characteristics being straight. In Fig. 10 the R1 simple wave (shaded cyan) is
confined to a partially curvilinear triangle bounded by the contour
1-$\bar{t}_{1}$-$M_{1}$-1, inside which $J^{+}$ is constant. In the left
reflection zone R12 (shaded grey), bounded by the contour
$\bar{t}_{1}$-$M_{1}$-$\bar{t}_{12}$-$\bar{t}_{1}$, the flow is not a simple
wave, and all the $C^{\pm}$ characteristics are curved.

Figure 10: The pattern of characteristics in the normalized
$(\bar{m},\bar{t})$ plane, depicting interaction of the release waves
initiated by the shock load with $u_{p}=1.25$. Cyan-shaded areas show the R1
and R2 simple wave zones. Grey-shaded areas are the R12 and R23 reflection
zones. Thick red curve $M_{sp}$-$M_{2}$-$\bar{t}_{pb}$ is the trajectory of
the pull-back shock

As the R1 wave is reflected from the left edge $\bar{m}=0$, it becomes the R2
simple wave (shown schematically in Fig. 9b) beyond the R12 reflection zone.
If the pb shock is launched from point $M_{sp}$, lying on the $C^{-}_{t1}$
characteristic, the R2 simple-wave zone (shaded cyan in Fig. 10) is bounded by
the contour $M_{1}$-$M_{sp}$-$M_{2}$-$\bar{t}_{2}$-$M_{1}$, where
$M_{sp}$-$M_{2}$ is the curved trajectory of the pb shock and
$\bar{t}_{2}$-$M_{2}$ is a curved segment of the $C_{h3}^{-}$ characteristic,
being the continuation of $C_{h2}^{+}$ upon its reflection off the free
surface $\bar{m}=1$. The curvilinear triangle
$\bar{t}_{2}$-$M_{2}$-$\bar{t}_{pb}$ (shaded grey), where $\bar{t}_{pb}$ is
the time of the pb-shock arrival at $\bar{m}=1$, encompasses the right
reflection zone R23. In the R2 zone, all the $C^{+}$ characteristics are
straight — like the $C^{+}_{2o}$ characteristic passing through $M_{2}$. In
the absence of spall, the R1 trailing characteristic $C^{-}_{t1}$ is reflected
from the $\bar{m}=0$ edge at $\bar{t}_{12}$ and continues as the $C^{+}_{t2}$
straight line, which replaces the pb shock as the downstream boundary of the
R2 simple wave.

The O1 zone inside the rectilinear isosceles triangle 1-$M_{1}$-$\bar{t}_{2}$
is the region of constant flow with $u=u_{1}$, $\rho=\rho_{01}$, $c=c_{01}$,
$p=0$. The times $t_{2}$ and $t_{pb}$ in Fig. 8, which mark the beginning and
the end of the depression in the $u_{fs}(t)$ curve, are the non-normalized
counterparts of $\bar{t}_{2}$ and $\bar{t}_{pb}$ in Fig. 10. Note that all the
head and tail characteristics, plotted in Fig. 10, are also the trajectories
of weak flow discontinuities.

###  IV.2 Release isentrope in the PA approximation

To construct the desired PA approximation, we need certain information about
the isentrope $p(\rho)$ (curve $PO_{1}S$ in Fig. 1), passing through the zero-
pressure state $O_{1}$ at $\rho=\rho_{01}$. We begin by assuming that, in
addition to $\rho_{01}$ and $c_{01}=c(\rho_{01})$, one knows from shock
experiments the coefficient $b$ in the widely used linear approximation

$D=c_{01}+bu^{\prime}$ (26)

to the compression Hugoniot with the initial state at $O_{1}$; here $D$ and
$u^{\prime}$ are, respectively, the shock-front and the post-shock particle
velocities relative to the fluid in zone O1, which itself moves with velocity
$u_{1}$ in our coordinate system. Having applied the mass and momentum balance
relations across the shock front, we can convert (26) into the parametric
representation of the Hugoniot

$\displaystyle p_{H}(u^{\prime})$ $\displaystyle=$
$\displaystyle\rho_{01}u^{\prime}(c_{01}+bu^{\prime}),$ (27)
$\displaystyle\rho_{H}(u^{\prime})$ $\displaystyle=$
$\displaystyle\rho_{01}\frac{c_{01}+bu^{\prime}}{c_{01}+(b-1)u^{\prime}}.$
(28)

If we consider now an arbitrary isentropic simple wave, bordering zone O1 in
Fig. 10 and propagating in the positive direction, we can similarly
parametrize the isentrope $p(\rho)$ in terms of the particle velocity
increment $u^{\prime}$ by integrating the equation

$du^{\prime}=dp/(\rho c)=c\,d\ln\rho,$ (29)

which expresses the fact that $dJ^{-}=0$ throughout both the O1 zone and the
bordering simple-wave regions. In doing this, we can simply put
$p(u^{\prime})=p_{H}(u^{\prime})$ because the Hugoniot and the Poisson
adiabats have a second-order contact at the common origin point, which implies
$dp/du^{\prime}=dp_{H}/du^{\prime}$ and $d^{2}p/du^{\prime
2}=d^{2}p_{H}/du^{\prime 2}$ at $u^{\prime}=0$. Having applied this argument
to our R1 and R2 simple rarefaction waves, we obtain the following PA
approximation to the required isentrope

$\displaystyle p(\bar{u}^{\prime})$ $\displaystyle=$
$\displaystyle\rho_{01}c_{01}^{2}\,\left(\bar{u}^{\prime}+b\bar{u}^{\prime
2}\right),$ (30) $\displaystyle z(\bar{u}^{\prime})$ $\displaystyle=$
$\displaystyle 1+2b\bar{u}^{\prime},$ (31)

where

$\bar{u}^{\prime}=\frac{u^{\prime}}{c_{01}},\quad
u^{\prime}=\left\\{\begin{array}[]{ll}u_{1}-u>0,&\mbox{R1 wave},\\\
u-u_{1}<0,&\mbox{R2 wave}.\end{array}\right.$ (32)

Below, we refer to the second terms on the right-hand sides of Eqs. (30) and
(31) as the _pu-nonlinearity corrections_. They characterize the nonlinearity
of the pertinent EOS and are proportional to the value of the fundamental
gasdynamic derivative Thompson (1971)

$\Gamma=\frac{c^{4}}{2v^{3}}\left(\frac{\partial^{2}v}{\partial
p^{2}}\right)_{s}=\frac{v^{3}}{2c^{2}}\left(\frac{\partial^{2}p}{\partial
v^{2}}\right)_{s}=2b$ (33)

in the state where $\bar{u}^{\prime}=0$. The accuracy of the PA approximation
to the GWEOS isentrope, passing through the normal state with parameters from
Table 1, is illustrated in Fig. 11, where the essentially positive
dimensionless ratio $p/(\rho_{0}c_{0}u^{\prime})$ is plotted versus
$\bar{u}^{\prime}$. The error in $p$ stays below 5% for
$\left|\bar{u}^{\prime}\right|\lesssim 0.15$, but exceeds 30% for
$\bar{u}^{\prime}\gtrsim 0.5$.

Figure 11: PA approximation to the release isentrope (dash-dotted), passing
through the normal state $\rho_{0}$, $p_{0}=0$, versus particle velocity
increment $u^{\prime}$. Shown also are the exact Hugoniot (solid gray) and
Poisson (solid blue) adiabats calculated from Eqs. (1)–(3)

If Eqs. (30), (31) are deemed to be exact along a considered isentrope, they
yield

$\displaystyle\rho(\bar{u}^{\prime})$ $\displaystyle=$
$\displaystyle\rho_{01}\left[1-(2b)^{-1}\ln(1+2b\bar{u}^{\prime})\right]^{-1}$
(34) $\displaystyle=$
$\displaystyle\rho_{01}\left[1+\bar{u}^{\prime}-(b-1)\bar{u}^{\prime
2}+\ldots\right],$ $\displaystyle c(\bar{u}^{\prime})$ $\displaystyle=$
$\displaystyle
c_{01}(1+2b\bar{u}^{\prime})\left[1-(2b)^{-1}\ln(1+2b\bar{u}^{\prime})\right]^{-1}$
(35) $\displaystyle=$ $\displaystyle
c_{01}\left[1+(2b-1)\bar{u}^{\prime}-b\bar{u}^{\prime 2}+\ldots\right].$

Because neither $\rho$ nor $z$ can be negative, the applicability of the above
formulae is in any case limited to the range

$-1/2b<\bar{u}^{\prime}<\left(e^{2b}-1\right)/2b,$ (36)

which is usually sufficient for application to spall experiments. In fact, if
separate values of $\rho(u^{\prime})$ or $c(u^{\prime})$ are needed, one can
justifiably restrict oneself to the zero- and first-order terms in powers of
$\bar{u}^{\prime}$.

###  IV.3 Evaluation of the strain rate

Introduction of nonlinear terms opens up a possibility to obtain an a priori
estimate of the strain rate at spallation without invoking the $u_{fs}(t)$
curve — the possibility not available in the acoustic approximation. Our
estimate is based on the assumption of time-independent linear velocity
profile

$\frac{\partial u(t,m)}{\partial m}\approx\frac{\partial u(t_{1m},m)}{\partial
m}\approx\frac{u_{1}}{m_{1m}}$ (37)

throughout the reflection zone R12 — the LU approximation; here the _mid-
reflection_ time $t_{1m}$ is defined as the moment when the characteristics
$C^{-}_{t1}$ and $C^{+}_{h2}$ in Fig. 10 intersect at point $M_{1}$, and the
zone O1 of constant flow behind R1 has its maximum extension $m_{0}-m_{1m}$
along the $m$ coordinate. Figure 12 confirms that this approximation remains
quite accurate even in the highly nonlinear case of $u_{p}=4.0$.

Figure 12: Normalized velocity profiles at the moment $t_{1m}$ of R1
reflection and at the moment $t_{sp}$ of spallation, calculated with the DEIRA
code for the two reference cases of $u_{p}=1.25$ (black) and $u_{p}=4.0$
(blue)

The LU approximation (37) can be supported by the observation that the piece-
wise linear profile of $u$, displayed in Fig. 9a, is in fact the exact
solution of the nonlinear Euler equation for $u(t,m)$, obtained by
substituting Eq. (30) into Eq. (22); in this solution the width $\Delta
m=(\rho_{1}c_{1}-\rho_{01}c_{01})t$ of the R1 wave grows linearly in time.
When the reflection of such a linear in $m$ rarefaction pulse from the $m=0$
boundary is treated in the acoustic limit of $c_{1}-c_{01}\ll c_{01}$, both
approximate equalities in (37) become exact. Also, the fluid state at mid-
reflection is then everywhere uniform with $p=0$ and $\rho=\rho_{01}$. Note
that the slope (37) remains finite as both $u_{1}$ and $m_{1m}$ become
infinitely small in the limit of $u_{p}\to 0$, $c_{1}\to c_{01}$.

To evaluate the mid-reflection coordinate $\bar{m}_{1m}=m_{1m}/m_{0}$, we use
the obvious fact that the integral $\oint d\bar{t}=0$ along any closed contour
in the $(\bar{m},\bar{t})$ plane. Then we combine this circulation rule with
Eq. (25) and apply it to the perimeter of the R1 simple-wave zone. Whereas
integration along straight characteristic segments, on which all flow
parameters remain constant, is trivial, for any curved segment — like
$\bar{t}_{1}$-$M_{1}$ in Fig. 10 — we assume the acoustic impedance to be
approximately constant and equal to its arithmetic mean $z=(z_{a}+z_{b})/2$
between the end values $z_{a}$ and $z_{b}$. Thus, the condition of zero $\oint
d\bar{t}=\oint(d\bar{m}/d\bar{t})^{-1}d\bar{m}$ along the triangular contour
1-$\bar{t}_{1}$-$M_{1}$-1 yields the equation

$z_{1}^{-1}+2(1+z_{1})^{-1}\bar{m}_{1m}-(1-\bar{m}_{1m})=0,$ (38)

from which we get

$\bar{m}_{1m}=\frac{z_{1}-z_{1}^{-1}}{z_{1}+3}=\frac{2\beta_{1}(1+\beta_{1})}{(2+\beta_{1})(1+2\beta_{1})},$
(39)

where, according to (31),

$z_{1}=1+2\beta_{1},\quad\beta_{1}=bu_{1}/c_{01}.$ (40)

Note that the acoustic limit corresponds to $\beta_{1}\ll 1$.

Now, with $\partial u/\partial m$ known from Eq. (37), the strain rate (10)
can be calculated from Eq. (21) as

$\dot{\varepsilon}=\rho\frac{\partial v}{\partial t}=\rho\frac{\partial
u}{\partial m}.$ (41)

Before finalizing our estimate of $\dot{\varepsilon}$, we recall that the
spatial derivatives of $u$ and $p$ are discontinuous on the $C^{-}_{t1}$
characteristic — as is manifested by the kinks at $\bar{m}=\bar{m}_{sp}$ on
the velocity profiles in Fig. 12. Hence, we have to choose between the two
different values $(\partial u/\partial m)_{\pm}$ obtained on the two sides
$m=m_{sp}\pm 0$ of the spall boundary. Because it is the material at
$m<m_{sp}$ which undergoes spall fracture, we make the logical choice of
$(\partial u/\partial m)_{-}=u_{1}/m_{1m}$. In contrast, the often used
acoustic formula (19) represents the value $(\partial u/\partial m)_{+}$ in
the unfractured material, which is exactly one half of $(\partial u/\partial
m)_{-}$ in the acoustic limit. Thus, having combined Eqs. (37)–(41), we obtain
the following PA estimates for the strain rates at mid-reflection,

$\dot{\varepsilon}_{1m}=\frac{\rho_{01}u_{1}}{m_{1m}}=\frac{c_{0}}{h_{0}z_{0}}\frac{(2+\beta_{1})(1+2\beta_{1})}{2b(1+\beta_{1})},$
(42)

and at spallation

$\dot{\varepsilon}_{sp}=\frac{\rho_{sp}}{\rho_{01}}\,\dot{\varepsilon}_{1m}=\left(1-\left|\bar{u}^{\prime}_{sp}\right|\right)\dot{\varepsilon}_{1m}.$
(43)

The small correction
$\left|\bar{u}^{\prime}_{sp}\right|\simeq\sigma_{th}/\rho_{0}c_{0}^{2}$,
accounting for the fluid expansion between the $t_{1m}$ and $t_{sp}$ moments,
is typically a few percent even for the highest theoretical values of
$\sigma_{th}$ Mayer and Mayer (2015); Cai, Wu, and Luo (2017); Mayer and Mayer
(2020) and can usually be neglected; the PA formula for
$|\bar{u}^{\prime}_{sp}|$ is given in the next section.

If time $t_{2}$ is known from the measurements of the free surface velocity
$u_{fs}(t)$, one readily obtains an alternative a posteriori LU estimate

$\displaystyle\bar{m}_{1m}$ $\displaystyle=$ $\displaystyle
1-\frac{1}{2}\frac{c_{0}t_{2}}{h_{0}z_{0}},$ (44)
$\displaystyle\dot{\varepsilon}_{1m}$ $\displaystyle=$
$\displaystyle\frac{c_{0}}{h_{0}}\frac{u_{1}/c_{01}}{z_{0}-c_{0}t_{2}/2h_{0}},$
(45)

which is a direct consequence of the LU approximation (37) and of the obvious
fact that $t_{1m}=t_{2}/2$. Note that, having the dimension of $t^{-1}$, the
strain rate $\dot{\varepsilon}_{1m}$ in Eqs. (42) and (45) scales in inverse
proportion to the initial target thickness $h_{0}$.

Regarding practical application of Eqs. (42) and (43), we will distinguish two
subcases of the PA approximation depending on the available thermodynamic
information for the material in question, namely, a cruder PA0 version and a
more accurate PA1 variant:

$\begin{array}[]{lll}\mbox{PA0:}&\rho_{0},c_{0},b_{0}\rightarrow&\mbox{known,}\\\
&\rho_{01},c_{01},b,u_{1}\rightarrow&\mbox{replaced by
}\rho_{0},c_{0},b_{0},u_{p};\\\\[3.0pt]
\mbox{PA1:}&\rho_{0},c_{0},\rho_{01},c_{01},b,u_{1}\rightarrow&\mbox{all
known.}\end{array}$

In the PA0 version, only the normal-state values of $\rho_{0}$, $c_{0}$,
$b_{0}$ are assumed to be known, while for the PA1 case one needs also the
parameters of the $p=0$ state $O_{1}$ on the release isentrope behind the
initial load shock. Analogously, the LU approximation is subdivided into a
cruder LU0 and a more accurate LU1 versions.

Table 3: Comparison of the strain rates, evaluated in different versions of the PA and LU approximations, with those obtained in numerical simulations with the DEIRA code | DEIRA | PA1 | PA0 | LU1 | LU0 | ac
---|---|---|---|---|---|---
$u_{p}=1.25$, $\Delta u_{fs}=1.111$
$\bar{m}_{1m}$ | 0.40 | 0.350 | 0.356 | 0.401 | 0.392 | –
$(h_{0}/c_{0})\dot{\varepsilon}_{1m}$ | 0.51 | 0.578 | 0.587 | 0.505 | 0.532 | –
$(h_{0}/c_{0})\dot{\varepsilon}_{sp}$ | 0.41 | 0.500 | 0.508 | 0.437 | 0.461 | 0.086
$u_{p}=4.0$, $\Delta u_{fs}=1.063$
$\bar{m}_{1m}$ | 0.72 | 0.603 | 0.612 | 0.717 | 0.633 | –
$(h_{0}/c_{0})\dot{\varepsilon}_{1m}$ | 0.89 | 0.963 | 1.090 | 0.811 | 1.054 | –
$(h_{0}/c_{0})\dot{\varepsilon}_{sp}$ | 0.66 | 0.812 | 0.947 | 0.683 | 0.915 | 0.125

Table 3 compares the values of $\bar{m}_{1m}$, $\dot{\varepsilon}_{1m}$, and
$\dot{\varepsilon}_{sp}$, evaluated in different versions of the PA and LU
approximations, with those obtained in the numerical DEIRA runs 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-
invariant dimensionless product $(c_{0}/h_{0})\dot{\varepsilon}$. All the
required EOS parameters are taken from Tables 1 and 2; the values of $\Delta
u_{fs}$, used to calculate $\dot{\varepsilon}_{sp,ac}$ and
$\bar{u}^{\prime}_{sp}$, as well as the times $t_{2}$ are inferred from the
DEIRA produced plots of $u_{fs}(t)$ as shown in Fig. 8 for the case of
$u_{p}=1.25$. The numbers in Table 3 demonstrate that the PA formulae (42),
(42) provide a robust estimate of the strain rate at spallation with an
accuracy of about 15–30%, which is not sensitive to the loading shock
strength. The conventional acoustic formula (19), on the contrary, tends to
underestimate $\dot{\varepsilon}_{sp}$ by a significant factor of about 3–5,
as was already noticed in MD simulations Luo  _et al._ (2009); Cai, Wu, and
Luo (2017). If one can reliably measure the arrival time $t_{2}$ of the R2
rarefaction wave, Eqs. (43)–(45) yield even more accurate values of
$\bar{m}_{1m}$ and $\dot{\varepsilon}_{sp}$, which confirms the high accuracy
of the LU approximation (37).

###  IV.4 Evaluation of the spall strength

Similar to the acoustic formula (17), our PA estimate of the spall strength is
based on the measured free-surface velocity pullback $\Delta u_{fs}$. We begin
by arguing that the outer boundary $\bar{m}_{sp}$ of the spall zone,
represented by the pb-shock starting point $M_{sp}$ in Fig. 10, should
generally lie on the R1 trailing characteristic $C^{-}_{t1}$, where $\partial
p/\partial m$ and $\partial u/\partial m$ have a jump. Firstly, we note that
everywhere above and to the right of the $M_{1}$-$\bar{t}_{12}$ segment of the
$C^{-}_{t1}$ curve, where the pressure is negative in the R2 simple-wave zone,
we have

$\mathcal{M}_{sp}=c/D_{pb}<1.$ (46)

This inequality, obtained by substituting Eq. (41) and the relation
$dp=\pm\rho c\,du$ into the definition (14) of the effective spall Mach number
$\mathcal{M}_{sp}$, implies that the pb shock must originate either from the
$C^{-}_{t1}$ characteristic or from the interior of the R12 reflection zone.
On the other hand, we know that in the acoustic limit $\beta_{1}\ll 1$ it
cannot start inside the R12 zone because in this limit $\partial p/\partial
m\to 0$ and $\mathcal{M}_{sp}\to\infty$ everywhere inside R12. This proves
that at least for not too large values of $u_{p}$ and $\beta_{1}$ the pb shock
can only start from the $C^{-}_{t1}$ characteristic. In strongly nonlinear
cases with $\beta_{1}\gtrsim 1$, this fact is confirmed by numerical
simulations (at least for the present EOS), as is shown in Fig. 6.

With $M_{sp}$ lying on the $C^{-}_{t1}$ characteristic above the $M_{1}$
point, we can use the analytic properties of the R2 simple wave, represented
by Eqs. (30)–(35), to relate $\sigma_{sp}$ to $\Delta u_{fs}$. The simple
acoustic formula (17) is based on the premise that the negative-pressure spall
pulse (see Fig. 9b), created at $M_{sp}$ with the initial velocity decrement

$u^{\prime}_{sp}\equiv
c_{01}\bar{u}^{\prime}_{sp}\,\stackrel{{\scriptstyle\mbox{\scriptsize
def}}}{{=}}\,u(t_{sp},m_{sp})-u_{1}<0,$ (47)

propagates without attenuation down to the free surface, where the well-known
(Zel’dovich and Raizer, 2012, § 11, Ch. XI) velocity doubling rule $\Delta
u_{fs}=-2u^{\prime}_{sp}$ applies; then (17) follows directly from Eq. (30) in
the limit of $b|\bar{u}^{\prime}|\ll 1$.

In the PA approximation, we take into account the attenuation of the spall
pulse, caused by the decreasing amplitude of the pb shock as it propagates up
the rising density profile of the R2 simple wave. For this, we trace the
curvilinear pb-shock trajectory up to its intersection with the right-
reflected characteristic $C^{-}_{h3}$ at point $M_{2}$ and apply the doubling
rule

$\Delta u_{fs}=-2u^{\prime}_{2}$ (48)

to the velocity decrement

$u^{\prime}_{2}\equiv
c_{01}\bar{u}^{\prime}_{2}\,\stackrel{{\scriptstyle\mbox{\scriptsize
def}}}{{=}}\,u(t_{2m},m_{2m})-u_{1}<0.$ (49)

at that point, which has coordinates $(\bar{m}_{2m},\bar{t}_{2m})$ in Fig. 10.
Having introduced the _pb-attenuation correction_

$\delta_{pb}\,\stackrel{{\scriptstyle\mbox{\scriptsize
def}}}{{=}}\,u^{\prime}_{sp}/u^{\prime}_{2}-1$ (50)

and invoking Eq. (30), we obtain the following PA formula for the spall
strength

$\displaystyle\sigma_{sp}$ $\displaystyle=$
$\displaystyle-p(\bar{u}^{\prime}_{sp})=\rho_{01}c_{01}^{2}|\bar{u}^{\prime}_{sp}|\left(1-b|\bar{u}^{\prime}_{sp}|\right)$
(51) $\displaystyle=$ $\displaystyle\frac{1}{2}\rho_{0}c_{0}\Delta
u_{fs}\,z_{0}^{-1}(1+\delta_{pb})\left[1-\beta_{2}(1+\delta_{pb})\right],$

where

$\beta_{2}=b|\bar{u}^{\prime}_{2}|=b\frac{\Delta
u_{fs}}{2c_{01}},\quad|\bar{u}^{\prime}_{sp}|=(1+\delta_{pb})\frac{\Delta
u_{fs}}{2c_{01}}.$ (52)

The attenuation correction

$\delta_{pb}=\frac{2(1-2\beta_{2})}{\omega+\left[\omega^{2}-2\beta_{2}(1-2\beta_{2})(3\mu+1)\right]^{1/2}}$
(53)

is found as the solution of the quadratic equation (71), where

$\mu=\frac{1+\beta_{1}}{1+\frac{3}{2}\beta_{1}}\,\frac{1-\frac{3}{2}\beta_{2}}{1-\frac{1}{2}\beta_{2}},\quad\omega=\frac{1}{2}+\mu\left(2-3\beta_{2}\right).$
(54)

Note that formula (51) may be considered as a generalization of the simple
extrapolation of the Hugoniot into the negative-pressure region, advocated in
Ref. Kanel, 2010, which takes a more consistent account of the first-order
nonlinear terms. For the fractional mass of the spall zone we obtain

$\bar{m}_{sp}=\bar{m}_{1m}\left(1-\frac{|u^{\prime}_{sp}|}{u_{1}}\right)=\bar{m}_{1m}\left[1-\frac{(1+\delta_{pb})\Delta
u_{fs}}{2u_{1}}\right],$ (55)

where $\bar{m}_{1m}$ is given by Eq. (39). The derivation of Eqs. (53)–(55) is
detailed in Appendix A.

Although the effect of the pb-shock attenuation is fairly well known and has
been discussed in literature Antoun  _et al._ (2003); Cai, Wu, and Luo (2017),
the explicit formula (53) for the corresponding correction $\delta_{pb}$ is an
important new result of this paper. For the particular considered problem of
symmetric plate collision, this correction is a function of two dimensionless
parameters $\beta_{1}$ and $\beta_{2}$ representing (in relative terms),
respectively, the intensity of the initial load and the amplitude of the spall
pulse. Figure 13, which displays the behavior of
$\delta_{pb}(\beta_{1},\beta_{2})$ in the allowed range of $0\leq\beta_{2}\leq
0.5$, reveals that for $\beta_{2}\lesssim 0.3$ the pb-attenuation correction
stays within a narrow range $0.4\leq\delta_{pb}<0.56$ for any value of
$\beta_{1}>0$, i.e. for arbitrary velocity of the impactor plate. This
observation is important for applications because in real experiments the
limit of $\beta_{2}\approx 0.3$ is hardly ever exceeded.

Figure 13: The pb-attenuation correction
$\delta_{pb}=\delta_{pb}(\beta_{1},\beta_{2})$ as a function of $\beta_{2}$
for three selected values of $\beta_{1}$

Of special interest is the seemingly paradoxical result that in the limit of
$\beta_{1},\beta_{2}\to 0$ the pb-attenuation correction does not vanish but
approaches the value

$\delta_{00}\equiv\delta_{pb}(0,0)=\frac{2}{5}$ (56)

that is not much lower than the upper limit of $\delta_{pb,max}\approx 0.56$.
The latter is explained by the fact that for $\beta_{1}\sim\beta_{2}\ll 1$ the
small amplitude of the pb shock has the same order of magnitude as the
difference between its propagation speed and the local sound velocity. In this
sense, attenuation of the pb shock is qualitatively similar to the attenuation
of weak saw-tooth-like pulses (Landau and Lifshitz, 1987, § 102) calculated in
the first-order nonlinear approximation and controlled by the value of the
fundamental gasdynamic derivative $\Gamma$ defined in Eq. (33). Remarkably, in
our case $\delta_{00}$ turns out to be a universal constant which depends
neither on $\Gamma$ nor on target parameters. Its value can be calculated
independently and more rigorously by direct integration of the equation (58)
for the pb-shock trajectory in the asymptotic limit of $\beta_{1},\beta_{2}\to
0$, which yields $\delta_{00}=\sqrt{2}-1$ that is practically
indistinguishable from $\delta_{00}=0.4$ obtained from Eq. (53).

In our PA formula (51) for the spall strength, the nonlinear corrections
appear as three different terms, namely, the modified (compared to the normal
state) acoustic impedance $\rho_{01}c_{01}=z_{0}^{-1}\rho_{0}c_{0}$, the pb-
attenuation correction $\delta_{pb}$, and the negative pu-nonlinearity
correction $-\beta_{2}(1+\delta_{pb})$. The relative importance of these
correction terms can be judged from the data in Tables 2 and 4. The negative
impedance correction $z_{0}^{-1}-1$ is small (a few percent) for moderate
loads near the spall threshold, but becomes significant for highly nonlinear
cases with large $u_{p}$. The positive pb-attenuation correction $\delta_{pb}$
is always significant and must never be ignored, at least not in the
experiments where the flyer and the target plates have comparable masses.

Table 4:  Comparison of the PA and acoustic (ac) approximations with the exact (EOS or DEIRA) results for the spall strength $\sigma_{th}$ and the spalled mass fraction $\bar{m}_{sp}$, calculated for the two reference values of $u_{p}$. The values of $\Delta u_{fs}$ are taken from the $u_{fs}(t)$ curves produced by the DEIRA runs | EOS/DEIRA | PA1 | PA0 | ac
---|---|---|---|---
$u_{p}=1.25$, $\Delta u_{fs}=1.111$
$\beta_{1}$ | – | 0.6442 | 0.6614 | –
$\beta_{2}$ | – | 0.2941 | 0.2939 | –
$\delta_{pb}$ | – | 0.4370 | 0.4378 | –
$\sigma_{sp}$ | 11.15 (EOS) | 10.65 | 10.80 | 13.01
$\bar{m}_{sp}$ | 0.131 (DEIRA) | 0.120 | 0.128 | -0.14
$u_{p}=4.0$, $\Delta u_{fs}=1.063$
$\beta_{1}$ | – | 2.031 | 2.117 | –
$\beta_{2}$ | – | 0.2879 | 0.2812 | –
$\delta_{pb}$ | – | 0.4740 | 0.4800 | –
$\sigma_{sp}$ | 8.219 (EOS) | 8.147 | 10.75 | 12.45
$\bar{m}_{sp}$ | 0.568 (DEIRA) | 0.477 | 0.492 | 0.277

The relative impact of the pu-nonlinearity term depends on the value of
$\beta_{2}$, which is insensitive to $u_{p}$ and reaches its maximum for the
maximum possible spall strength. But even then the value $\beta_{2}\approx
0.3$, obtained in our simulations and quoted in Table 4, appears as rather an
absolute upper limit; more typical for liquid metals at maximum possible
tension would be $\beta_{2}\lesssim 0.1$–0.15, like the value
$\beta_{2}\approx 0.12$ obtained in the MD simulations for liquid copper Cai,
Wu, and Luo (2017). The latter implies that the pu-nonlinearity should be
taken into account near the theoretical limit for $\sigma_{th}$, but can be
ignored in experiments where the spall failure occurs at tensions 3–10 times
below this limit.

A characteristic feature of the formula (51) is that the contributions of its
three PA correction terms tend to compensate for one another. For large
$u_{p}$, the decreasing impedance $\rho_{01}c_{01}$ of the zero-pressure
release state strongly suppresses the impact of the $1+\delta_{pb}$ factor;
for $\beta_{2}\approx 0.2$–0.25 the effect of the pb-attenuation is
practically canceled by that of the pu-nonlinearity term. The latter explains
why it is not easy to establish a clear tendency in the role of nonlinearity
corrections Cai, Wu, and Luo (2017), and why in practice the simple formula
(17) may frequently be more accurate than more complex expressions with only a
partial account for the nonlinearity effects. On the whole, the data in Table
4 confirm that Eqs. (51) and (55) provide a significant improvement over the
traditional acoustic formulae for interpretation of spall experiments with
symmetric initial loading.

Our final remark is on the applicability of the acoustic limit. Naively, one
might expect that in the limit of small amplitudes of both the loading shock
and the spall pulse, i.e. in the limit of $\beta_{1},\beta_{2}\to 0$,
$z_{0}\to 1$, the acoustic formula (17) should be asymptotically exact. Our
present analysis demonstrates that, whenever a pull-back shock is present,
this is not the case. Instead, for symmetrically loaded impact plates, the
correct small-amplitude limit is given by the formula

$\sigma_{sp}=\frac{1}{2}(1+\delta_{00})\rho_{0}c_{0}\Delta
u_{fs}=\frac{1}{\sqrt{2}}\,\rho_{0}c_{0}\Delta u_{fs},$ (57)

which exceeds the acoustic result by 41%, and which would be suitable for
interpretation of experiments where $\beta_{1}\lesssim 0.3$ and
$\beta_{2}\lesssim 0.05$.

##  V Conclusions

The paper presents the results of detailed analysis of the phase-flip (PF)
model of spallation in pure liquids, where spall failure is treated as an
infinitely fast irreversible liquid-gas phase transition. The spall strength
is set equal to its theoretical maximum determined by the thermodynamic
stability limit, i.e. by the negative pressure on the liquid branch of the
spinodal in a two-phase EOS. The main advantage of the proposed model is that
it can be easily and in a thermodynamically fully consistent way incorporated
into the framework of ideal hydrodynamics. Unlike in the approach proposed in
Ref. Grady, 1988, 1996, energy dissipation due to spall fracture is described
without invoking the concepts of surface tension and/or viscosity, simply as
the entropy jumps in the PF transitions and in shocks generated by the PF
pressure jumps.

The process of spall fracture in the proposed model develops as a sudden
appearance of one or more “fog” bubbles behind hypersonic spall fronts
(representing a new type of discontinuity in ideal hydrodynamics) that are
launched from one or more isolated inception points at local pressure minima.
In our context, the term “fog” refers to the equilibrium matter state, to
which a metastable liquid relaxes in the process of spinodal decomposition.
The volume and mass of each fog bubble continue to grow until the spall front,
precipitating the PF transition, is overtaken by the pull-back compression
shock.

To illustrate how the PF model could be used in practice, the results of
numerical modeling of a symmetric plate impact experiment are presented. In
most aspects, they look quite similar to those found in the MD simulation of
an analogous problem Cai, Wu, and Luo (2017) — but obtained at immeasurably
lower computational cost. The principal difference between the two approaches
manifests itself only in the microstructure of the fracture zone: while in the
MD case it is stipulated by the statistics and the inherent atomic scale, in
the fully deterministic PF approach the sizes and distribution of liquid
fragments are determined by small flow perturbations and inhomogeneities of
the initial state. When equipped with an adequate two-phase EOS, the PF model
has a potential to become a powerful tool for modeling experiments where
spallation and cavitation in liquids is an issue.

Using the PF numerical results as a reference solution, we derive new post-
acoustic formulae for evaluation of the strain rate, the fractured mass
fraction, and the spall strength in plate impact experiments where these
quantities are to be inferred from the measured free-surface velocity. Though
valid for the particular case of a symmetric plate impact, the new approximate
formula (51) for $\sigma_{th}$ reveals a non-trivial interplay of the three
nonlinear correction terms, which in reality often compensate for one another.
A notable new result is that the correction factor for the pb-shock
attenuation never becomes close to unity, and in the limit $u_{p}\ll c_{0}$ of
weak initial loads approaches a universal value $1+\delta_{00}=\sqrt{2}$ for
any equation of state.

The key assumption of instant phase relaxation, on which the proposed PF model
is based, is justified by explosive increase of the rate of homogeneous bubble
nucleation near the limit of thermodynamic stability. In reality, the
respective phase relaxation timescale $\tau_{pr}$ remains finite and limited
to $\tau_{pr}\gtrsim 1$–3 ps, as is indicated by pressure relaxation dynamics
in MD simulations. Hence, the hypersonic spall fronts, discussed in Sect.
II.3, have a finite width of $\sim u_{spf}\tau_{pr}$, where $u_{spf}$ is the
spall-front velocity in the comoving system. Consequently, the proposed model
becomes applicable on time scales $t\gg\tau_{pr}$ and length scales $l\gg
c_{0}\tau_{pr}\gtrsim 10$ nm, where $c_{0}$ is the sound velocity in the
relevant zero-pressure metastable state. Also, the length scale of emerging
fracture zones must exceed the typical size of critical vapor bubbles, which
again sets a limit of $l\gg 10$ nm. If the strain-rate dependence of the spall
strength is to be accounted for, the present assumption of the spinodal-
limited spall strength can be relaxed along the lines proposed in Refs. Dekel
_et al._ , 1998, Pisarev  _et al._ (2009), Faik  _et al._ (2012) and in other
similar works.

###### Acknowledgements.

This work was funded by the Keldysh Institute of Applied Mathematics of the
Russian Academy of Sciences under the state contract with the Ministry of
Science and Higher Education of the Russian Federation.

##  Appendix A Derivation of formula (53) for the pb-attenuation correction

We begin by ascertaining that the trajectory $\bar{m}_{pb}(\bar{t})$ of the pb
shock, as it propagates from point $M_{sp}$ to point $M_{2}$ in Fig. 10
against the background of the R2 simple wave where all flow characteristics
are functions of a single parameter $u^{\prime}<0$ [see Eqs. (30)–(35)], is
governed by the equation

$\frac{d\bar{m}_{pb}}{d\bar{t}}\equiv\zeta(\bar{u}^{\prime})=1-b|\bar{u}^{\prime}|.$
(58)

Indeed, the trajectory $m(t)$ of any shock front, traveling with a speed
$D_{if}$ relative to the fluid in an initial state $(v_{i},p_{i})$, obeys the
equation (Zel’dovich and Raizer, 2012, § 16, Ch. I)

$\frac{dm}{dt}=\rho_{i}D_{if}=\left(\frac{p_{f}-p_{i}}{v_{i}-v_{f}}\right)^{1/2},$
(59)

where $(v_{f},p_{f})$ is the final post-shock state. Important here is that
the shock impedance (59) remains invariant by interchanging the initial and
final states. In our context, where the compression Hugoniot for a final state
$p_{f}=p_{H}(u^{\prime})>0$ at $u^{\prime}>0$ and the initial state
$p_{i}=p_{H}(0)=0$ is given by Eqs. (27), (28), the normalized shock impedance
takes the form

$\frac{d\bar{m}}{d\bar{t}}=1+b\bar{u}^{\prime}.$ (60)

If we consider now the pb compression shock, whose initial pressure
$p_{i}=p_{H}(u^{\prime})<0$ is given by Eq. (27) [identical to Eq. (30)] for
some $u^{\prime}<0$ and the final pressure $p_{f}=p_{H}(0)=0$, we conclude
that, due to the $i\rightleftarrows f$ interchange symmetry, its impedance
must be equal to that of the _rarefaction shock_ (perhaps non-physical but
mathematically conceivable) represented by the extension of Eqs. (27), (28),
and (60) into the $u^{\prime}<0$ domain, i.e. by Eq. (58). Note that, unlike
its normalized impedance $\zeta$, the front speed of the pb shock

$D_{pb}=c_{01}-(b-1)|u^{\prime}|$ (61)

is not equal to the extension of Eq. (26) to negative $u^{\prime}$. The fact
that $\rho_{H}(u^{\prime})$ in Eq. (28) is not exactly equal to
$\rho(\bar{u}^{\prime})$ from Eq. (34) has no practical significance because,
when expanded in powers of $\bar{u}^{\prime}$, the two formulae begin to
diverge in $O(\bar{u}^{\prime 3})$ terms only, which can be safely ignored for
$|\bar{u}^{\prime}|\leq|\bar{u}^{\prime}_{sp}|\lesssim 0.1$–0.15 encountered
in reality. Of negligible effect is also the essentially positive pressure
behind the pb shock, which, being limited from above by $p_{EQ,sp}$ (see Sect.
II.3), can practically always be neglected in comparison with $\sigma_{sp}$,
as is seen from Table 2.

Next, we introduce an auxiliary characteristic $C^{+}_{o2}$ into the chart in
Fig. 10, which belongs to the $C^{+}$ fan of the R2 wave and crosses the
right-reflected characteristic $C^{-}_{h3}$ at the same point $M_{2}$ as the
pb shock. Because for $u^{\prime}<0$ the acoustic impedance (31) is smaller
than the pb-shock impedance $\zeta$, $C^{+}_{o2}$ is straight and lies below
the pb-shock trajectory everywhere between $M_{2}$ and its intersection with
the $C^{-}_{t1}$ characteristic at point $M_{o2}$ with coordinates
$(\bar{m}_{o2},\bar{t}_{o2})$, which lie in the range

$\bar{m}_{sp}<\bar{m}_{o2}<\bar{m}_{1m},\quad\bar{t}_{1m}<\bar{t}_{o2}<\bar{t}_{sp}.$
(62)

When applied to points $M_{o2}$ and $M_{sp}$, the LU approximation (37) takes
the form

$\frac{u_{1}}{m_{1m}}=\frac{u_{1}-|u^{\prime}_{2}|}{m_{o2}}=\frac{u_{1}-|u^{\prime}_{sp}|}{m_{sp}},$
(63)

which relates $\bar{m}_{o2}$ and $\bar{m}_{sp}$ to the known quantity
$\bar{m}_{1m}$, leading, in particular, to Eq. (55).

Now, the two unknown quantities $\delta_{pb}$ and $\bar{m}_{2m}$ can be
calculated from the system of two equations

$1-\bar{m}_{1m}+\frac{1-\bar{m}_{2m}}{(1+z_{2})/2}-\frac{\bar{m}_{2m}-\bar{m}_{o2}}{z_{2}}-\frac{\bar{m}_{1m}-\bar{m}_{o2}}{(1+z_{2})/2}=0,$
(64)
$\frac{\bar{m}_{2m}-\bar{m}_{o2}}{z_{2}}-\frac{\bar{m}_{2m}-\bar{m}_{sp}}{(\zeta_{2}+\zeta_{sp})/2}-\frac{\bar{m}_{o2}-\bar{m}_{sp}}{(z_{2}+z_{sp})/2}=0,$
(65)

obtained, respectively, by applying the $\oint d\bar{t}=\oint
z^{-1}\,d\bar{t}=0$ circulation rule to the two closed contours
$M_{1}$-$\bar{t}_{2}$-$M_{2}$-$M_{o2}$-$M_{1}$ and
$M_{o2}$-$M_{2}$-$M_{sp}$-$M_{o2}$ in the $(\bar{m},\bar{t})$ plane, and by
using the mean impedance values along the curvilinear characteristic segments
and the pb-shock trajectory. Here

$z_{2}=1-2\beta_{2},\quad\zeta_{2}=1-\beta_{2},$ (66)

are the normalized acoustic and pb-shock impedances at point $M_{2}$, and

$z_{sp}=1-2\beta_{2}(1+\delta_{pb}),\quad\zeta_{sp}=1-\beta_{2}(1+\delta_{pb})$
(67)

are the same impedances at point $M_{sp}$; $\beta_{2}$ is defined in Eq. (52).
Equations (64) and (65) are easily transformed to

$\displaystyle\bar{m}_{2m}-\bar{m}_{o2}$ $\displaystyle=$
$\displaystyle(1-\bar{m}_{1m})\frac{z_{2}(3+z_{2})}{1+3z_{2}},$ (68)
$\displaystyle\frac{\bar{m}_{2m}-\bar{m}_{o2}}{\bar{m}_{o2}-\bar{m}_{sp}}$
$\displaystyle=$
$\displaystyle\frac{2z_{2}(z_{2}+z_{sp}+\zeta_{2}+\zeta_{sp})}{(z_{2}+z_{sp})(\zeta_{2}+\zeta_{sp}-2z_{2})},$
(69)

while from Eq. (63) we get

$\bar{m}_{o2}-\bar{m}_{sp}=\bar{m}_{1m}(\beta_{2}/\beta_{1})\delta_{pb}.$ (70)

Having substituted Eqs. (66)–(68) together with Eqs. (70) and (39) into Eq.
(69), we obtain the quadratic equation

$(3\mu+1)\beta_{2}\delta_{pb}^{2}-2\omega\delta_{pb}+2(1-2\beta_{2})=0$ (71)

for the attenuation correction $\delta_{pb}$, where $\mu$ and $\omega$ are
defined in Eq. (54). The basic equations (64), (65) and the ensuing formulae
are physically meaningful only if the impedance $1\geq z_{2}\geq 0$, i.e.
under the condition

$0\leq\beta_{2}\leq\frac{1}{2}.$ (72)

It is not difficult to prove that the above condition guarantees positiveness
of the discriminant in Eq. (71).

## References

  * Antoun  _et al._ (2003) T. Antoun, L. Seaman, D. R. Curran, G. I. Kanel, S. V. Razorenov, and A. V. Utkin,  _Spall Fracture_ , Shock Wave and High Pressure Phenomena (Springer-Verlag New York, 2003).
  * Kanel, Fortov, and Razorenov (2007) G. I. Kanel, V. E. Fortov, and S. V. Razorenov, “Shock waves in condensed-state physics,” Physics-Uspekhi 50, 771–791 (2007).
  * Kanel  _et al._ (2015) G. I. Kanel, A. S. Savinykh, G. V. Garkushin, and S. V. Razorenov, “Dynamic strength of tin and lead melts,” JETP Letters 102, 548–551 (2015).
  * Zaretsky (2016) E. B. Zaretsky, “Experimental determination of the dynamic tensile strength of liquid Sn, Pb, and Zn,” Journal of Applied Physics 120, 025902 (2016).
  * de Rességuier  _et al._ (2007) T. de Rességuier, L. Signor, A. Dragon, M. Boustie, G. Roy, and F. Llorca, “Experimental investigation of liquid spall in laser shock-loaded tin,” Journal of Applied Physics 101, 013506 (2007).
  * Agranat  _et al._ (2010) 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 an aluminum melt at extremely high tension rates under the action of femtosecond laser pulses,” JETP Letters 91, 471–477 (2010).
  * Krivokorytov  _et al._ (2017) M. S. Krivokorytov, A. Y. Vinokhodov, Y. V. Sidelnikov, V. M. Krivtsun, V. O. Kompanets, A. A. Lash, K. N. Koshelev, and V. V. Medvedev, “Cavitation and spallation in liquid metal droplets produced by subpicosecond pulsed laser radiation,” Phys. Rev. E 95, 031101 (2017).
  * Stan  _et al._ (2016) 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 and spallation in water drops subjected to nanosecond shock waves,” J. Phys. Chem. Lett. 7, 2055–2062 (2016).
  * Basko  _et al._ (2017) 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, “Fragmentation dynamics of liquid-metal droplets under ultra-short laser pulses,” Laser Physics Letters 14, 036001 (2017).
  * Grigoryev  _et al._ (2018) 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. Sidelnikov, V. M. Krivtsun, K. N. Koshelev, and V. V. Medvedev, “Expansion and fragmentation of a liquid-metal droplet by a short laser pulse,” Phys. Rev. Applied 10, 064009 (2018).
  * Pisarev  _et al._ (2009) V. V. Pisarev, A. Y. Kuksin, G. E. Norman, V. V. Stegailov, and A. V. Yanilkin, “Microscopic theory and kinetic model of spall in liquids,” AIP Conference Proceedings 1195, 801–804 (2009).
  * Cai, Wu, and Luo (2017) Y. Cai, H. A. Wu, and S. N. Luo, “Spall strength of liquid copper and accuracy of the acoustic method,” Journal of Applied Physics 121, 105901 (2017).
  * Mayer and Mayer (2020) A. E. Mayer and P. N. Mayer, “Strain rate dependence of spall strength for solid and molten lead and tin,” International Journal of Fracture 222, 171–195 (2020).
  * Anisimov  _et al._ (1999) S. I. Anisimov, N. A. Inogamov, A. M. Oparin, B. Rethfeld, T. Yabe, M. Ogawa, and V. E. Fortov, “Pulsed laser evaporation: equation-of-state effects,” Applied Physics A 69, 617–620 (1999).
  * Colombier  _et al._ (2005) J. P. Colombier, P. Combis, F. Bonneau, R. L. Harzic, and E. Audouard, “Hydrodynamic simulations of metal ablation by femtosecond laser irradiation,” Phys. Rev. B 71, 165406 (2005).
  * Zhao  _et al._ (2011) N. Zhao, A. Mentrelli, T. Ruggeri, and M. Sugiyama, “Admissible shock waves and shock-induced phase transitions in a van der Waals fluid,” Phys. Fluids 23, 086101 (2011).
  * Basko (2018a) M. M. Basko, “Centered rarefaction wave with a liquid-gas phase transition in the approximation of “phase-flip” hydrodynamics,” Physics of Fluids 30, 123306 (2018a).
  * Boteler and Sutherland (2004) J. M. Boteler and G. T. Sutherland, “Tensile failure of water due to shock wave interactions,” Journal of Applied Physics 96, 6919–6924 (2004).
  * Grady (1988) D. E. Grady, “The spall strength of condensed matter,” Journal of the Mechanics and Physics of Solids 36, 353–384 (1988).
  * Grady (1996) D. E. Grady, “Spall and fragmentation 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.
  * Bakshi (2006) V. Bakshi,  _EUV Sources for Lithography_ , Press Monographs (SPIE Press, 2006).
  * Mizoguchi  _et al._ (2015) H. Mizoguchi, H. Nakarai, T. Abe, K. M. Nowak, Y. Kawasuji, H. Tanaka, Y. Watanabe, T. Hori, T. Kodama, Y. Shiraishi, T. Yanagida, G. Soumagne, T. Yamada, T. Yamazaki, S. Okazaki, and T. Saitou, “Performance of one hundred watt HVM LPP-EUV source,” Proc. SPIE 9422, 94220C–94220C–13 (2015).
  * Martynyuk (1991) M. M. Martynyuk, “Generalized van der Waals equation of state for liquids and gases,” Zh. Fiz. Khim. 65, 1716–1717 (1991).
  * Martynyuk (1993) M. M. Martynyuk, “Transition of liquid metals into vapor in the process of pulse heating by current,” International Journal of Thermophysics 14, 457–470 (1993).
  * Basko (2018b) M. M. Basko, “Generalized van der Waals equation of state for in-line use in hydrodynamic codes,” Keldysh Institute Preprints 112, 28 p. (2018b).
  * Haynes (2015) W. M. Haynes,  _CRC Handbook of Chemistry and Physics, 96th Edition_ , 100 Key Points (CRC Press, 2015).
  * Landau and Lifshitz (1996) L. D. Landau and E. M. Lifshitz,  _Statistical Physics_ , 3rd ed. (Butterworth Heinemann, 1996).
  * Landau and Lifshitz (1987) L. D. Landau and E. M. Lifshitz,  _Fluid Mechanics_ , Course of theoretical physics (Pergamon Press, 1987).
  * Kanel (2010) G. I. Kanel, “Spall fracture: methodological aspects, mechanisms and governing factors,” International Journal of Fracture 163, 173–191 (2010).
  * Imre  _et al._ (2008) 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 20, 244104 (2008).
  * Malyshev  _et al._ (2015) V. L. Malyshev, D. F. Marin, E. F. Moiseeva, N. A. Gumerov, and I. S. Akhatov, “Study of the tensile strength of a liquid by molecular dynamics methods,” High Temperature 53, 406–412 (2015).
  * Blander and Katz (1975) M. Blander and J. L. Katz, “Bubble nucleation in liquids,” AIChE Journal 21, 833–848 (1975).
  * Skripov and Skripov (1979) V. P. Skripov and A. V. Skripov, “Spinodal decomposition (phase transitions via unstable states),” Sov. Phys. Usp. 22, 389 (1979).
  * Skripov (1974) V. P. Skripov,  _Metastable Liquids_ (Wiley, New York, 1974).
  * Dekel  _et al._ (1998) 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 84, 4851–4858 (1998).
  * Faik  _et al._ (2012) S. Faik, M. M. Basko, A. Tauschwitz, I. Iosilevskiy, and J. A. Maruhn, “Dynamics of volumetrically heated matter passing through the liquid-vapor metastable states,” High Energy Density Phys. 8, 349–359 (2012).
  * Basko (2001) 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).
  * Basko (1990) M. M. Basko, “Spark and volume ignition of DT and D2 microspheres,” Nuclear Fusion 30, 2443–2452 (1990).
  * Fortov, Iakubov, and Khrapak (2006) V. Fortov, I. Iakubov, and A. Khrapak,  _Physics of Strongly Coupled Plasma_ , 1st ed., International series of monographs on physics (Clarendon Press, Oxford, 2006).
  * Thompson (1971) P. A. Thompson, “A fundamental derivative in gasdynamics,” Phys. Fluids (1958-1988) 14, 1843–1849 (1971).
  * Mayer and Mayer (2015) A. E. Mayer and P. N. Mayer, “Continuum model of tensile fracture of metal melts and its application to a problem of high-current electron irradiation of metals,” Journal of Applied Physics 118, 035903 (2015).
  * Luo  _et al._ (2009) 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 106, 013502 (2009).
  * Zel’dovich and Raizer (2012) Y. B. Zel’dovich and Y. P. Raizer,  _Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena_ , Dover Books on Physics (Dover Publications, 2012).