# A Flux-enriched Godunov Method for Multi-material Problems with Interface
Slide and Void Opening

Tim Wallis<EMAIL_ADDRESS>Philip T. Barton  Nikolaos Nikiforakis  Department
of Physics, University of Cambridge, Cavendish Laboratory, JJ Thomson Avenue,
CB3 0HE, UK  AWE Aldermaston, Reading, Berkshire, RG7 4PR, UK

###### Abstract

This work outlines a new three-dimensional diffuse interface finite volume
method for the simulation of multiple solid and fluid components featuring
large deformations, sliding and void opening. This is achieved by extending an
existing reduced-equation diffuse interface method by means of a number of
novel flux-modifiers and interface seeding routines that enable the
application of different material boundary conditions. The method allows for
slip boundary conditions across solid interfaces, material-void interaction,
and interface separation. The method is designed to be straightforward to
implement, inexpensive and highly parallelisable. This makes it suitable for
use in large, multi-dimensional simulations that feature many complex
materials and physical processes interacting over multiple levels of adaptive
mesh refinement. Furthermore, the method allows for the generation of new
interfaces in a conservative fashion and therefore naturally facilitates the
simulation of high-strain rate fracture. Hence, the model is augmented to
include ductile damage to allow for validation of the method against demanding
physical experiments. The method is shown to give excellent agreement with
both experiment and existing Eulerian interface tracking algorithms that
employ sharp interface methods.

© British Crown Owned Copyright 2020/AWE

###### keywords:

Multi-physics , Elastoplastic Solids , Diffuse interface , Fracture , Slide ,
Void opening

##  1 Introduction

Multi-material, multi-physics simulation techniques are vital for the accurate
modelling of many real-world physical systems. These problems are dynamic and
often highly non-linear. Such simulations often feature multiple complex
materials with different equations of state, large deformations and strain
rates, and additional physical effects such as material fracture or
reactivity. The combination of all of these factors makes a multi-physics
approach crucial; modelling the system in a single framework ensures full
physical interaction between all the components. To this end, realistic
physical models, suitable mathematical solution techniques, and efficient
computer algorithms must be brought together in a single framework to enable a
complete system to be modelled.

Importantly, real-life problems of interest often involve large, three
dimensional systems evolving over long time-scales. This makes the development
of a practical numerical scheme crucial, and provides the motivation for the
work at hand. This work will be focused on providing a straightforward
algorithm that can form a base for large, 3D, multi-material, multi-physics
simulations. As such, particular care will be taken to ensure the resulting
algorithm is direct, prioritising the ease with which additional physics can
be incorporated.

A variety of mature methods for solving multi-material shock-capturing
compressible flow problems already exist. These are divided into a number of
approaches, the broadest of which being the division between Lagrangian and
Eulerian methods. In essence, this choice comes down to choosing the frame of
reference in which to solve the equations of motion. Lagrangian methods
(typically finite element methods) consider a frame that is fixed with the
material which deforms as the simulation progresses. Eulerian methods
(typically finite volume methods) consider a frame fixed in space, through
which material flows. These methods both have their advantages and
disadvantages, but in the high-deformation, high-strain-rate problems
considered in this work, large topological changes can cause Lagrangian
meshing methods to fail due to the severe mesh distortion. As Eulerian meshes
do not deform with the material, they can smoothly handle large material
deformations and complex topological changes, and so are well suited to the
task at hand.

However, Eulerian methods are not without drawbacks and have challenges
associated with their implementation in higher dimensions and with AMR. Most
existing Eulerian models capable of simulating fluids and elastoplastic solids
are based on tracking material interfaces and fall into one of three broad
categories. The first is methods based on homogenised mixed cells and volume-
of-fluid reconstruction, such as Benson [12]. These methods underpin many
well-established and legacy multi-material codes. The next is Ghost Fluid
methods, such as Schoch et al. [53], Michael and Nikiforakis [39], Michael et
al. [40], Barton and Drikakis [8] and Barton [6]. These are a fairly recent
development, based on the methods of Fedkiw et al. [26], Sambasivan and
Udaykumar [51] and Wang et al. [60]. Ghost Fluid methods capture material
boundaries with the use of a level set and apply dynamic boundary conditions
to allow for material interaction. These methods have the downside of being
non-conservative at material interfaces. Lastly, Cut-cell methods such as
Gokhale et al. [27], Miller and Colella [41] and Barton et al. [10] resolve
the geometry of cells intersected by interfaces and apply a strict finite
volume discretisation to capture the flux across material boundaries. All of
these approaches have the advantage that they maintain arbitrarily sharp
interfaces, but to do so they involve complex interface reconstructions and
mixed-cell algorithms. Furthermore, applying load-balancing to many of these
schemes when using AMR is difficult due to the non-uniform numerical methods.
Robust implementations can therefore be challenging to construct. Some of
these difficulties have led to the development of hybrid approaches such as
Arbitrary Lagrangian Eulerian methods such as Barlow et al. [4] and co-
simulation methods, based on embedding finite element grids in fluid domains,
see for example Deiterding et al. [21]. These methods are popular for fluid-
structure interaction problems, but introduce the natural complexities of mesh
management associated with Lagrangian methods.

Diffuse interface methods represent a different way to handle material
interfaces in Eulerian methods. These methods track materials through the
volume fraction that they occupy in a given computational cell. This volume
fraction is continuous and multiple materials can be present in a single cell.
This numerical mixture is then evolved using a single system of partial
differential equations that encompasses the physics of all the components in
the mixture. A particular benefit of diffuse interface models is that they
provide conservation equations for mass, momentum and energy across
interfaces. Previous work has shown that these schemes have the considerable
advantage over interface tracking methods in that they can support genuine
fluid mixtures, allowing for the study of phenomena such as cavitation and
chemical reaction which rely on physical mixtures [47, 59].

Within diffuse interface methods, there are a range of different approaches.
These different approaches are often due to equilibrium assumptions made about
the physical model in question. The most general approach is to assume full
non-equilibrium between mixture components. This results in each material in
the mixture having its own density, velocity, pressure and temperature, with
relaxation terms augmenting these equations to allow for materials to interact
and exchange conserved quantities. An example of such a model is the work of
Baer and Nunziato [2]. However, experiments reveal that the scales over which
both pressure and velocity equilibrate are very small, and evolving a velocity
and pressure for every material is an encumbrance. It is therefore common [34,
47] to perform an asymptotic reduction of non-equilibrium models, applying
equilibrium conditions to variables such as velocity, to obtain a reduced-
equation model. This process leaves scope for how many equilibrium assumptions
are to be applied, with authors such as Allaire et al. [1] assuming both
pressure and velocity equilibrium, whereas authors such as Favrie et al. [25]
retain pressure non-equilibrium. These approaches lead to components of a
numerical mixture sharing a single equation for quantities such as momentum or
energy. This in turn results in having to solve fewer equations overall, and
is considered beneficial in terms of both efficiency and stability. However,
these equilibrium assumptions often come at a price. In some scenarios, such
as when considering porous materials like Baer and Nunziato [2], it is wholly
inappropriate to assume velocity equilibrium between phases. At an even more
basic level, velocity equilibrium prevents the use of slip boundary conditions
at material interfaces and prevents material separation. This can be seen
straightforwardly, as both of these problems require a velocity discontinuity
at the interface: tangentially for slip and normally for separation. Such
discontinuities are not permitted by a velocity equilibrium model, which will
average the two velocities at the interface, effectively applying a stick
boundary condition.

Diffuse interface methods are well established for multi-fluid problems, see
for example Allaire et al. [1], Saurel and Abgrall [52], Massoni et al. [37].
By contrast, relatively few authors have applied diffuse interface methods to
coupled solid-fluid dynamics, including Barton [7], Favrie et al. [25], Favrie
and Gavrilyuk [24], Hank et al. [28] and Tavelli et al. [54]. This is likely
due to a number of factors. Firstly, Eulerian diffuse interface methods
require qualitatively different numerical methods compared to those
traditionally applied to solid dynamics, especially in Lagrangian schemes.
Secondly, the diffuse nature of the material interfaces has, in the past,
proved to be a limiting factor on the accuracy of such methods. This issue
has, however, recently been addressed by the development of various interface
sharpening methods, such as the THINC method designed by Xiao et al. [61],
implemented by Barton [7]. Lastly, diffuse interface methods have not been
generally applied to solid dynamics due the limiting nature of the enforced
stick boundary condition when using reduced-equation models. Unlike with
fluids, it is natural for materials with strength to be able to slide along
each other and separate from each other. This problem is addressed by the work
at hand.

This paper details a new diffuse interface method for the simulation of
coupled elastoplastic solids and fluids. Namely, this work introduces methods
to allow for slip boundary conditions, material-void interaction, and material
fracture in a diffuse interface context. The work extends the method of Barton
[7], who developed a model for coupled solid-fluid dynamics based on the
Allaire five-equation multi-fluid model [1]. The new method introduces a novel
set of flux-modifiers and interface seeding routines that will enable the
desired boundary conditions to be applied on top of the base model. Alongside
this, a damage model is incorporated to enable material damage and fracture to
be studied.

As has been alluded to, many other methods exist that can model multi-material
systems, and these will be mentioned by way of contrast and comparison in the
validation of the new method. However, the new method presents a number of
specific advantages – it is designed to be straightforward to implement,
inexpensive and highly parallelisable, making it suitable for use in large,
multi-dimensional simulations that feature multiple complex materials and
physical processes over several levels of adaptive mesh refinement (AMR). As
such, the method can be considered an ‘engineering solution’ to the problems
presented by multi-material interactions.

The highly-parallelisable nature of the method stems from the locality of the
algorithm. There is no need to use the expensive level set reinitialisation
and extrapolation routines used by Ghost Fluid methods, which require multi-
patch and multi-level communication. The method is also highly numerically
homogeneous. One single set of partial differential equations that encompasses
the physics of all materials present is solved across the entire domain and
material boundary conditions are facilitated by a number of simple flux-
modifiers. These flux-modifiers only change how the flux between cells is
calculated, rather than requiring any drastically different numerical
machinery.

How well a given algorithm can scale with the number of materials present is
also important as the number of materials in a simulation gets increasingly
large. As has been mentioned, non-equilibrium equation systems allow for
material shear by providing every material with its own complete set of
variables. By way of example, Ndanou et al. [45] facilitate a form of fracture
by giving each material its own deformation and stress tensor. However, this
method gets increasingly expensive as the number of materials grows. In a
scenario where compaction dynamics and velocity relaxation times are not of
specific interest, the new method has the advantage of allowing material shear
while only using a single velocity, momentum, deformation and stress. This
greatly reduces how the algorithmic complexity scales with the number of
materials. It should be stressed, however, that the void and fracture
components of the method could still be used in a non-equilibrium type
equation system, if so desired.

The problems of interest often occur over disparate scales, with both regions
of fine detail with strong interaction and regions with smoothly varying
features. It is therefore desirable to use AMR to better focus computational
resources where they are needed. With computational efficiency in mind, the
algorithm outlined in this work is designed to work straightforwardly with
hierarchical AMR algorithms, saving computational resources by allowing high
resolution around regions of interest and lower resolution in smooth regions.

Finally, the method facilitates a more conservative fracture algorithm than
can be achieved with many other sharp interface methods, due to the
continuous, diffuse representation of interfaces. It is commonplace in many
existing methods for simulating fracture to introduce cracks by eroding
material from critically damaged cells [e.g. 18, 16, 6]. This approach is
highly non-conservative and leads to issues if the conserved state of eroded
material is simply redistributed to its neighbours. The new method avoids this
problem, instead allowing a material interface to open continuously.

##  2 Governing Theory

###  2.1 Evolution Equations

In the interests of focusing on multi-physics, only a summary of the core
evolution equations is provided here. Full details, including the derivation,
are outlined by Barton [7].
Materials are allowed to mix at their interfaces. A material’s contribution to
a spatially averaged physical quantity is weighted by its volume fraction,
$\phi$, in that region. This mixing is referred to as numerical mixing, to
distinguish it from the physical mixtures such as those produced by reactive
fluids. The state of any material $l$ is characterised by the phasic density
$\rho_{(l)}$, volume fraction $\phi_{(l)}$, symmetric left unimodular stretch
tensor $\overline{\mathbf{V}}^{e}$, velocity vector $\mathbf{u}$, and internal
energy $\mathscr{E}$.

The left stretch tensor $\mathbf{V}^{e}$ is related to the deformation tensor
$\mathbf{F}$ by the polar decomposition:

$\mathbf{F}=\mathbf{V}^{e}\mathbf{R},$ (1)

after which it is normalised to obtain $\overline{\mathbf{V}}^{e}$:

$\overline{\mathbf{V}}^{e}=\det\left(\mathbf{V}^{e}\right)^{-1/3}\mathbf{V}^{e}.$
(2)

The base model of Barton [7] assumes mechanical equilibrium; materials in a
mixture region share a single velocity, pressure and deviatoric strain. This
assumption will be maintained, in that only a single equation for each of
these quantities will be evolved, but flux-modifiers in the numerical method
will allow for the correct material boundary conditions to be imposed on top.

For $l=1,\ldots,N$ materials:

$\displaystyle\frac{\partial\phi_{(l)}}{\partial
t}+\frac{\partial\phi_{(l)}u_{k}}{\partial x_{k}}$ $\displaystyle=$
$\displaystyle\phi_{(l)}\frac{\partial u_{k}}{\partial x_{k}}$ (3)
$\displaystyle\frac{\partial\rho_{(l)}\phi_{(l)}}{\partial
t}+\frac{\partial\rho_{(l)}\phi_{(l)}u_{k}}{\partial x_{k}}$ $\displaystyle=$
$\displaystyle 0$ (4) $\displaystyle\frac{\partial\rho u_{i}}{\partial
t}+\frac{\partial(\rho u_{i}u_{k}-\sigma_{ik})}{\partial x_{k}}$
$\displaystyle=$ $\displaystyle 0$ (5) $\displaystyle\frac{\partial\rho
E}{\partial t}+\frac{\partial(\rho Eu_{k}-u_{i}\sigma_{ik})}{\partial x_{k}}$
$\displaystyle=$ $\displaystyle 0$ (6)
$\displaystyle\frac{\partial\overline{\mathbf{V}}^{e}_{ij}}{\partial
t}+\frac{\partial\left(\overline{\mathbf{V}}^{e}_{ij}u_{k}-\overline{\mathbf{V}}^{e}_{kj}u_{i}\right)}{\partial
x_{k}}$ $\displaystyle=$
$\displaystyle\frac{2}{3}\overline{\mathbf{V}}^{e}_{ij}\frac{\partial
u_{k}}{\partial x_{k}}-u_{i}\beta_{j}-\bm{\Phi}_{ij}\ ,$ (7)

where $E=\mathscr{E}+|\mathbf{u}|/2$ denotes the specific total energy,
$\bm{\sigma}$ denotes the Cauchy stress tensor,
$\beta_{j}=\partial\overline{\mathbf{V}}^{e}_{kj}/\partial x_{k}$, and
$\bm{\Phi}$ represents the contribution from plastic effects. Some multi-
physics closure models introduce a dependence on material history variables
such as the equivalent plastic strain, $\varepsilon_{p(l)}$, or the damage
parameter $D_{(l)}$. For these materials, additional evolution equations are
required:

$\displaystyle\frac{\partial\rho_{(l)}\phi_{(l)}\alpha_{(l)}}{\partial
t}+\frac{\partial\rho_{(l)}\phi_{(l)}\alpha_{(l)}u_{k}}{\partial x_{k}}$
$\displaystyle=$ $\displaystyle\rho_{(l)}\phi_{(l)}\dot{\alpha}_{(l)}\ .$ (8)

Here $\alpha_{(l)}$ represents any such history parameter which is advected
and evolved with a material as time progresses. A material may have more than
one history variable, in which case $\alpha$ represents a vector.

Although it might appear that the model describes solid materials, it will be
shown in the next section that fluids can be considered as a special case. It
is the equation-of-state for each material that ultimately distinguishes
solids from fluids. The system allows for the fully-coupled multi-physics
solution of problems, which in this work includes the interaction of
damageable elastoplastic solids, fluid mixtures and vacuum. The extension to
yet more multi-physics applications follows straightforwardly by the inclusion
of additional history parameters and closure relations, for example reactive
fluid mixtures as shown in Wallis et al. [59].

###  2.2 Thermodynamics

The internal energy $\mathscr{E}$ for each material is defined by an equation-
of-state that conforms to the general form:

$\mathscr{E}_{(l)}\left(\rho_{(l)},T_{(l)},\mbox{dev}\left({\mathbf{H}^{e}}\right),\alpha_{(l)}\right)=\mathscr{E}_{(l)}^{c}\left(\rho_{(l)},\alpha_{(l)}\right)+\mathscr{E}_{(l)}^{t}\left(\rho_{(l)},T_{(l)}\right)+\mathscr{E}_{(l)}^{s}\left(\rho_{(l)},\mbox{dev}\left({\mathbf{H}^{e}}\right),\alpha_{(l)}\right)\
,$ (9)

where

$\mbox{dev}\left({\mathbf{H}^{e}}\right)=\ln\left(\overline{\mathbf{V}}^{e}\right)\
$ (10)

is the deviatoric111 For any $N\times N$ matrix $\mathbf{M}$,
$\mbox{dev}\left({\mathbf{M}}\right):=\mathbf{M}-\frac{1}{N}\tr(\mathbf{M})\mathbf{I}\
$ denotes the matrix deviator. Here, $\tr(\mathbf{M})$ denotes the trace, and
$\mathbf{I}$ denotes the identity matrix.  Hencky strain tensor and $T$ is the
temperature. The three terms on the right hand side are the contribution due
to cold compression or dilation,
$\mathscr{E}_{(l)}^{c}\left(\rho_{(l)},\alpha_{(l)}\right)$, the contribution
due to temperature deviations,
$\mathscr{E}_{(l)}^{t}\left(\rho_{(l)},T_{(l)}\right)$, and the contribution
due to shear strain
$\mathscr{E}_{(l)}^{s}\left(\rho_{(l)},\mbox{dev}\left({\mathbf{H}^{e}}\right),\alpha_{(l)}\right)$.
The cold compression energy will generally be provided by the specific closure
model for each material. The thermal energy is given by

$\displaystyle\mathscr{E}_{(l)}^{t}(\rho_{(l)},T)$ $\displaystyle=$
$\displaystyle
C_{(l)}^{\text{V}}\left(T-T_{(l)}^{0}\theta_{(l)}^{D}\left(\rho_{(l)}\right)\right)\
,$ (11)

where $C_{(l)}^{\text{V}}$ is the heat capacity, $T_{(l)}^{0}$ is a reference
temperature, and $\theta_{(l)}^{\text{D}}(\rho_{(l)})$ is the non-dimensional
Debye temperature. The Debye temperature is related to the Grüneisen function,
$\Gamma(\rho_{(l)})$, via

$\Gamma_{(l)}(\rho_{(l)})=\frac{\partial\ln\theta_{(l)}^{\text{D}}(\rho_{(l)})}{\partial\ln(1/\rho_{(l)})}=\frac{\rho_{(l)}}{\theta_{(l)}^{\text{D}}(\rho_{(l)})}\frac{\partial\theta_{(l)}^{\text{D}}(\rho_{(l)})}{\partial\rho_{(l)}}\
.$ (12)

The specific form of the Grüneisen function for each material must be
provided, but it is generally taken to be constant. The shear energy is given
by

$\mathscr{E}_{(l)}^{s}(\rho_{(l)},\mbox{dev}\left({\mathbf{H}^{e}}\right),\alpha_{(l)})=\frac{G_{(l)}\left(\rho_{(l)},\alpha_{(l)}\right)}{\rho_{(l)}}\mathcal{J}^{2}\left(\mbox{dev}\left({\mathbf{H}^{e}}\right)\right)\
,$ (13)

where $G(\rho_{(l)},\alpha_{(l)})$ is the shear modulus, and

$\mathcal{J}^{2}(\mbox{dev}\left({\mathbf{H}^{e}}\right))=\tr\left(\mbox{dev}\left({\mathbf{H}^{e}}\right)\cdot\mbox{dev}\left({\mathbf{H}^{e}}\right)^{\text{T}}\right)$
(14)

is the second invariant of shear strain. This form is chosen such that the
resultant stresses are a non-linear analogy of Hooke’s law. Bruhns et al. [17]
find that this form provides good empirical agreement for a range of materials
and deformations. It is remarked that other extensions of infinitesimal strain
theory to finite deformations have taken the density divisor to be the
reference value $\rho^{0}_{(l)}$. This formulation can be replicated here
without further modification to the model or numerical method by setting
$G=\left(\frac{\rho_{(l)}}{\rho^{0}_{(l)}}\right)\overline{G}$, and redefining
$\overline{G}$ as the actual shear modulus.

For each component, the Cauchy stress, $\bm{\sigma}$, and pressure, $p$, are
inferred from the second law of thermodynamics and classical arguments for
irreversible elastic deformations:

$\displaystyle\bm{\sigma}_{(l)}$ $\displaystyle=$ $\displaystyle
p_{(l)}\mathbf{I}+\mbox{dev}\left({\bm{\sigma}_{(l)}}\right)$ (15)
$\displaystyle p_{(l)}$ $\displaystyle=$
$\displaystyle\rho^{2}_{(l)}\frac{\partial\mathscr{E}_{(l)}}{\partial\rho_{(l)}}$
(16) $\displaystyle\mbox{dev}\left({\bm{\sigma}_{(l)}}\right)$
$\displaystyle=$ $\displaystyle
2G_{(l)}\cdot\mbox{dev}\left({\mathbf{H}^{e}}\right)\ .$ (17)

It is reiterated that, although it might appear that the model describes solid
materials, fluids can be considered a special case where the shear modulus is
zero, resulting in a spherical stress tensor and no shear energy contribution.
This formulation lends itself well to diffuse interface modelling where
different phases that share the same underlying model can combine consistently
in mixture regions.

Mixture rules must be provided to represent the state of regions containing
multiple materials in a thermodynamically consistent way. The following
mixture rules are applied [1, 7]:

$\displaystyle 1$ $\displaystyle=$ $\displaystyle\sum_{l=1}^{N}\phi_{(l)}$
(18) $\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\phi_{(l)}\rho_{(l)}$ (19)
$\displaystyle\rho\mathscr{E}$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\phi_{(l)}\rho_{(l)}\mathscr{E}_{(l)}$ (20)
$\displaystyle G$ $\displaystyle=$
$\displaystyle\frac{\sum_{l=1}^{N}\left(\phi_{(l)}G_{(l)}(\rho_{(l)},\alpha_{(l)})/\Gamma_{(l)}\right)}{\sum_{l=1}^{N}\left(\phi_{(l)}/\Gamma_{(l)}\right)}$
(21) $\displaystyle c^{2}$ $\displaystyle=$
$\displaystyle\frac{\sum_{l=1}^{N}\left(\phi_{(l)}Y_{(l)}c_{(l)}^{2}/\Gamma_{(l)}\right)}{\sum_{l=1}^{N}\left(\phi_{(l)}/\Gamma_{(l)}\right)}\
,$ (22)

where $c$ is the sound speed and $Y_{(l)}=({\phi_{(l)}\rho_{(l)}})/{\rho}$ is
the mass fraction.

###  2.3 Closure Models

It can be seen that, by presenting the internal energy in the form outlined,
equation (16) can be written in the form:

$\displaystyle
p_{(l)}=p_{\text{\scriptsize{ref}},(l)}+\rho_{(l)}\Gamma_{(l)}\left(\mathscr{E}_{(l)}-\mathscr{E}_{\text{\scriptsize{ref}},(l)}\right)\
.$ (23)

Here,
$\mathscr{E}_{\text{\scriptsize{ref}},(l)}=\mathscr{E}^{c}_{(l)}+\mathscr{E}^{s}_{(l)}$
and
$p_{\text{\scriptsize{ref}},(l)}=\rho^{2}_{(l)}\frac{\partial\mathscr{E}_{\text{\scriptsize{ref}},(l)}}{\partial\rho_{(l)}}$.
This is the form of the standard Mie-Grüneisen equation-of-state. This
encompasses a wide range of different materials, not limited to solids,
depending on the choice of the reference curves
$\mathscr{E}_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$,
$p_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$, and $\Gamma_{(l)}(\rho_{(l)})$.
These additional freedoms incorporate thermal, compaction and shear effects.
For example:

  * •

$\mathscr{E}_{\text{\scriptsize{ref}},(l)}=p_{\text{\scriptsize{ref}},(l)}=0,\
\Gamma(\rho)=\Gamma_{0}=\gamma-1$, where $\gamma$ is the adiabatic index,
gives the ideal gas law used for relatively simple gases such as air.

  * •

$p_{\text{\scriptsize{ref}},(l)}=-\gamma p_{\infty},\
\mathscr{E}_{\text{\scriptsize{ref}},(l)}=e_{\infty},\
\Gamma(\rho)=\Gamma_{0}=\gamma-1$ gives the stiffened gas equation of state,
used to model denser fluids such as water.

  * •

$\Gamma(\rho)=\Gamma_{0}=\gamma-1$ and

$\displaystyle p_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle={\cal A}e^{-{\cal R}_{1}\frac{\rho_{0}}{\rho_{(l)}}}+{\cal
B}e^{-{\cal R}_{2}\frac{\rho_{0}}{\rho_{(l)}}}$ (24)
$\displaystyle\mathscr{E}_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle=\frac{{\cal A}}{{\cal R}_{1}\rho_{0}}e^{-{\cal
R}_{1}\frac{\rho_{0}}{\rho_{(l)}}}+\frac{{\cal B}}{{\cal
R}_{2}\rho_{0}}e^{-{\cal R}_{2}\frac{\rho_{0}}{\rho_{(l)}}}$ (25)

gives the JWL equation-of-state, widely used for condensed phase explosives or
reaction products [38].

  * •

$\Gamma(\rho)=\Gamma_{0}=\gamma-1$ and

$\displaystyle\mathscr{E}_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle=\frac{K_{0}}{2\rho_{(l)}\alpha^{2}}\left(\left(\frac{\rho_{(l)}}{\rho_{0}}\right)^{\alpha}-1\right)+\mathscr{E}^{s}_{(l)}$
(26) $\displaystyle G(\rho_{(l)})$
$\displaystyle=G_{0}\left(\frac{\rho_{(l)}}{\rho_{0}}\right)^{\beta+1}$ (27)
$\displaystyle p_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle=\rho^{2}_{(l)}\frac{\partial\mathscr{E}_{\text{\scriptsize{ref}},(l)}}{\partial\rho_{(l)}}$
(28)

gives the Romenskii equation-of-state, used for elastic solids [23].

The form of the Mie-Grüneisen equation of state now allows the total pressure
to be obtained from the total energy. To achieve this, equation 20 is combined
with the isobaric assumption that forms part of the mechanical equilibrium of
the underlying Allaire et al. [1] model. Firstly the internal energy is
obtained:

$\rho\mathscr{E}=\rho E-\frac{1}{2}\rho|\mathbf{u}|\ .$ (29)

Then the Mie-Grüneisen form of the equation of state is used to rewrite
equation 20, using the isobaric assumption:

$\rho\mathscr{E}=\sum_{l=1}^{N}\phi_{(l)}\rho_{(l)}\left(\mathscr{E}_{\text{\scriptsize{ref}},(l)}+\frac{(p-p_{\text{\scriptsize{ref}},(l)})}{\rho_{(l)}\Gamma_{(l)}}\right)\
.$ (30)

Finally, this equation can be rewritten in terms of $p$:

$p=\frac{\rho\mathscr{E}-\sum_{l}\phi_{(l)}\left(\rho_{(l)}\mathscr{E}_{\text{\scriptsize{ref}},(l)}-\frac{p_{\text{\scriptsize{ref}},(l)}}{\Gamma_{(l)}}\right)}{\sum_{l}\frac{\phi_{(l)}}{\Gamma_{(l)}}}\
.$ (31)

####  2.3.1 Plasticity

The introduction of plasticity through the source term $\bm{\Phi}$ follows the
method of convex potentials (see for example Ottosen [46]) and is therefore
thermodynamically compatible. In this approach, the von Mises yield criterion
forms the scalar potential which leads to the plastic flow rule:

$\displaystyle\bm{\Phi}$
$\displaystyle=\chi\sqrt{\frac{3}{2}}\frac{\mbox{dev}\left({\bm{\sigma}}\right)}{||{\mbox{dev}\left({\bm{\sigma}}\right)}||}\cdot\overline{\mathbf{V}}^{e}\
.$ (32)

The plastic flow rate $\chi$ is a closure model and must be suitable for
arbitrary mixtures. For a mixture of $N$ materials:

$\displaystyle\chi=\frac{\sum_{l=1}^{N}\left(\phi_{(l)}\chi_{(l)}(\rho_{(l)})/\Gamma_{(l)}\right)}{\sum_{l=1}^{N}\left(\phi_{(l)}/\Gamma_{(l)}\right)}\
.$ (33)

If any material does not obey a plasticity model, such as fluids, then that
material contributes $\chi_{(l)}=0$.

The form of $\chi$ relates the particular material flow model. This work
considers both ideal plasticity, where $\chi$ is a Heaviside function such
that the relaxation is non-zero only when the stress exceeds the yield surface
$\sigma_{Y}$:

$\displaystyle\chi_{(l)}=\chi_{(l)}^{0}H\left[\sqrt{\frac{3}{2}}||\text{dev}({\bm{\sigma}})_{(l)}||-\sigma_{Y}\right]\
,$ (34)

and the rate sensitive isotropic work-hardening model outlined by Johnson and
Cook [30]:

$\chi_{(l)}=\chi_{(l)}^{0}\exp\left[\frac{1}{c_{3}}\left(\frac{\sqrt{\frac{3}{2}}||\mbox{dev}\left({{\bm{\sigma}}_{(l)}}\right)||}{\sigma_{Y}\left(\varepsilon_{p,(l)}\right)}-1\right)\right],$
(35)

where ${\chi}_{0}>0$ is the reference plastic strain-rate and the constant
$c_{3}$ controls the rate dependency. In the Johnson and Cook model the yield
stress is given by:

$\displaystyle\sigma_{Y}\left(\varepsilon_{p,(l)}\right)=\left(c_{1}+c_{2}(\varepsilon_{p,(l)})^{n}\right)\left(1-\left(\frac{T-T_{0}}{T_{\text{melt}}-T_{0}}\right)^{m}\right)\
,$ (36)

where $c_{1}$ is the yield stress, $c_{2}$ is the strain hardening factor, $n$
is the strain hardening exponent, $T_{\text{melt}}$ is the melting temperature
of the material, $T_{0}=298$ K is a reference temperature, and $m$ is the
thermal softening exponent. This yield surface at the reference strain rate
$\chi_{(l)}^{0}$ is a function of the accumulated plastic strain,
$\varepsilon_{p}$, making it necessary to add an additional evolution equation
to the system to advect and evolve the plastic strain:

$\displaystyle\frac{\partial\rho_{(l)}\phi_{(l)}\varepsilon_{(l),p}}{\partial
t}+\frac{\partial\rho_{(l)}\phi_{(l)}\varepsilon_{(l),p}u_{k}}{\partial
x_{k}}=\rho_{(l)}\phi_{(l)}\chi_{(l)}\ .$ (37)

Using the Johnson-Cook model in this way results in a viscoplastic flow rule,
where plastic deformations can accumulate from the onset of loading. Note
however that the parameter $c_{3}$ is usually small such that
$\chi_{(l)}\ll\chi_{(l)}^{0}$ for stresses much below the characteristic
stress. Indeed, as $c_{3}\rightarrow 0$ the plastic flow becomes rate
independent and the stress becomes bounded by a yield surface.

###  2.4 Damage

Modelling damage is important if realistic material behaviour is to be
captured for solid materials undergoing finite deformation. This work focuses
on ductile damage, a significant mode of damage in most metals, where
substantial plastic deformation occurs before the onset of fracture, as
opposed to brittle fracture, where fracture directly proceeds from elastic
loading.

Physically, damage reduces the load carrying ability of a solid material
until, at a critical level, the material loses all strength and fractures.
Various branches of physics deal with damage on different scales, depending on
the application and type of damage at hand. The increase of damage in a
material can be explained on the micro-scale by the nucleation and coalescence
of micro-voids in the structure of the material and is intimately linked to
plastic effects. At a macro-scale, this presents itself as the formation of
visible fractures within a material.

A continuum damage mechanics (CDM) model is used to simulate ductile fracture.
CDM is a thermodynamically consistent way of modelling the effect of damage at
the continuum scale, where thermodynamic forces are derived using the second
law in order to advect and evolve damage in a given material [44]. Damage is
assumed to be isotropic, and can thus can be quantified and tracked by the
addition of a single scalar field for each damageable material, $D_{(l)}$:

$\displaystyle\frac{\partial\rho_{(l)}\phi_{(l)}D_{(l)}}{\partial
t}+\frac{\partial\rho_{(l)}\phi_{(l)}D_{(l)}u_{k}}{\partial
x_{k}}=\rho_{(l)}\phi_{(l)}\dot{D}_{(l)}\ .$ (38)

It is remarked that other formulations of damage can be derived, including
fully anisotropic models such as Barton [5] that quantify damage using a
tensor quantity. However, the efficacy of the damage model itself is outside
the scope of this paper. Rather, the intention is to provide a general
framework that demonstrates the capability of the new numerical method,
regardless of what specific constitutive material models are employed.

The contributions to the internal energy from cold compression and shear
energy for damageable materials are affected by the damage parameter, and
taken to be

$\displaystyle\mathscr{E}_{(l)}^{c}\left(\rho_{(l)},D_{(l)}\right)$
$\displaystyle=$
$\displaystyle(1-D_{(l)})\mathscr{E}_{(l)}^{c,\text{u}}\left(\rho_{(l)}\right)$
(39)
$\displaystyle\mathscr{E}_{(l)}^{s}\left(\rho_{(l)},\mbox{dev}\left({\mathbf{H}^{e}}\right),D_{(l)}\right)$
$\displaystyle=$
$\displaystyle(1-D_{(l)})\mathscr{E}_{(l)}^{s,\text{u}}\left(\rho_{(l)},\mbox{dev}\left({\mathbf{H}^{e}}\right)\right)\
.$ (40)

where $\mathscr{E}_{(l)}^{c,\text{u}}$ and $\mathscr{E}_{(l)}^{s,\text{u}}$
are the functions corresponding to the undamaged state, and can be any of the
forms highlighted in section 2.3. From this proposition, the Cauchy stress
tensor derived from the Clausius-Duhem inequality becomes:

$\displaystyle\bm{\sigma}_{(l)}$ $\displaystyle=$ $\displaystyle
p_{(l)}\mathbf{I}+\mbox{dev}\left({\bm{\sigma}_{(l)}}\right)$ (41)
$\displaystyle p_{(l)}$ $\displaystyle=$
$\displaystyle(1-D_{(l)})\rho^{2}_{(l)}\frac{\partial\mathscr{E}_{(l)}}{\partial\rho_{(l)}}$
(42) $\displaystyle\mbox{dev}\left({\bm{\sigma}_{(l)}}\right)$
$\displaystyle=$
$\displaystyle(1-D_{(l)})2G_{(l)}\cdot\mbox{dev}\left({\mathbf{H}^{e}}\right)\
.$ (43)

The presence of damage causes a reduction in load carrying area, with
$(1-D_{(l)})$ being the reduction factor, leading to the effective-stress
principal first proposed by Kachanov [32]. In this approach, damage $D_{(l)}$
is quantified by the fractional area of micro-voids on a plane intersecting a
representative volume element.

In a similar way to plasticity, the Clausius-Dunhem inequality gives the form
of the evolution equation for damage in terms of its conjugate thermodynamic
force, and the method of convex potentials is then used to derive the
evolution of the damage parameter:

$\displaystyle\dot{D}_{(l)}=-\dot{\lambda}_{(l)}\partialderivative{F^{D}_{(l)}}{Y_{(l)}}\
,$ (44)

where $\lambda_{(l)}$ is the Lagrange multiplier,
$Y_{(l)}={\partial\mathscr{E}_{(l)}}/{\partial D_{(l)}}$ is the elastic energy
release rate due to damage (the thermodynamic force conjugate to damage), and
$F^{D}_{(l)}$ is the damage dissipation potential. The associative flow rule
for damaged materials is

$\dot{\lambda}_{(l)}=\dot{\varepsilon}_{p,(l)}(1-D_{(l)})\ .$ (45)

The damage dissipation potential is taken from Bonora [14], who propose the
following potential to account for micro-mechanical processes:

$\displaystyle
F^{D}=\left[\frac{1}{2}\left(-\frac{Y}{S}\right)\frac{S}{1-D_{(l)}}\right]\frac{(D_{\text{crit},(l)}-D_{(l)})^{\frac{\alpha-1}{\alpha}}}{\varepsilon_{p}^{\frac{2+n}{n}}}\
.$ (46)

Under certain assumptions on the equation of state (detailed in the appendix)
this gives:

$\displaystyle\dot{D}_{(l)}=\alpha\frac{(D_{\text{crit},(l)}-D_{0,(l)})^{\frac{1}{\alpha}}}{\ln(\frac{\varepsilon_{\text{crit}}}{\varepsilon_{\text{thresh}}})}R_{t}\left(\frac{p}{\sigma_{\text{eq}}}\right)(D_{\text{crit},(l)}-D_{(l)})^{\frac{\alpha-1}{\alpha}}\left(\frac{\dot{\varepsilon}_{p}}{\varepsilon_{p}}\right)\
,$ (47)

where
$D_{\text{crit}},D_{0},\alpha,\varepsilon_{\text{crit}},\varepsilon_{\text{thresh}}$
are material parameters and $R_{t}$ is the stress triaxiality function:

$R_{t}\left(\frac{p}{\sigma_{\text{eq}}}\right)=\frac{2E_{(l)}(1-D_{(l)})^{2}\mathscr{E}_{(l)}^{c,\text{u}}\rho_{(l)}}{\sigma^{2}_{\text{eq}}}+\frac{E_{(l)}}{3G_{(l)}^{\text{u}}}$
(48)

where $G_{(l)}^{\text{u}},E_{(l)}^{\text{u}}$ are respectively the shear and
Young’s modulus in the undamaged state. A full derivation of the damage model
is provided in the appendix.

Pirondi and Bonora [48] extended this model to deal with cyclic loading,
proposing the following addition:

$\displaystyle\varepsilon_{p}^{-}$
$\displaystyle=\left\\{\matrixquantity{\varepsilon_{p}&&\text{if: }\ p<0\\\
0&&\text{else}\\\ }\right.$ (49) $\displaystyle\dot{D}$
$\displaystyle=\left\\{\matrixquantity{\dot{D}&&\text{if: }\ p<0\ \text{and}\
\varepsilon_{p}^{-}>\varepsilon_{p,\text{thresh}}\\\ 0&&\text{else}\\\
}\right.$ (50)

which conveys the underlying assumption that damage can only accrue in states
of tension, not under compression. Finally, threshold plastic strain
$\varepsilon_{p,\text{thresh}}$ is taken to depend on the stress triaxiality
using the model from Bonora et al. [15]:

$\displaystyle\varepsilon_{p,\text{thresh}}=\frac{\varepsilon_{\text{thresh}}^{1/R_{t}}}{R_{t}^{1/n}}\left(\frac{\varepsilon_{\text{crit}}^{2n}-\varepsilon_{\text{thresh}}^{2n}}{\varepsilon_{\text{crit}}^{2n/R_{t}}-\varepsilon_{\text{thresh}}^{2n/R_{t}}}\right)^{\frac{1}{2n}}\
,$ (51)

where $n$ is the strain hardening exponent.

##  3 Numerical Approach

The model is solved on a Cartesian mesh with local resolution adaptation in
space and time. This is achieved using the AMReX software from Lawrence
Berkeley National Laboratory [63], which includes an implementation of the
structured adaptive mesh refinement (SAMR) method of Berger and Colella [13]
for solving hyperbolic systems of partial differential equations of the form
of equation (52). In this approach, cells of identical resolution are grouped
into logically rectangular sub-grids or ‘patches’. Refined grids are derived
recursively from coarser ones, based upon a flagging criterion, to form a
hierarchy of successively embedded levels. All mesh widths on level $l$ are
$r_{l}$-times finer than on level $l-1$, i.e. $\Delta t_{l}:=\Delta
t_{l-1}/r_{l}$ and $\Delta\mathbf{x}_{l}:=\Delta\mathbf{x}_{l-1}/r_{l}$ with
$r_{l}\in\mathbb{N},r_{l}\geq 2$ for $l>0$ and $r_{0}=1$. The numerical scheme
is applied on level $l$ by calling a single-grid update routine in a loop over
all patches constituting the level. The discretisation of the constitutive
models does not differ between patches or levels, so for clarity the method
shall be described for a single sub-grid. Cell centres are denoted by the
indices $i,j,k\in\mathbb{Z}$ and each cell $C^{l}_{ijk}$ has the dimensions
$\Delta\mathbf{x}^{l}_{ijk}$.

It is convenient when describing the numerical method to express the system of
equations in compact vector form by separating it into various qualitatively
different parts: a conservative hyperbolic part for each spatial dimension,
non-conservative terms from the volume fraction and stretch tensor updates, a
source term due to plastic flow and associated damage accumulation, and a
source term to account for geometrical effects. This can be written as:

$\frac{\partial\mathbf{U}}{\partial t}+\frac{\partial\mathbf{F}_{k}}{\partial
x_{k}}=\mathbf{s}_{\text{\scriptsize{non-
con.}}}+\mathbf{s}_{p}+\mathbf{s}_{g}\ .$ (52)

Subject to the closure relations previously outlined, this is given by:

$\displaystyle\partialderivative{t}\matrixquantity(\phi_{(l)}\\\
\phi_{(l)}\rho_{(l)}\\\ \phi_{(l)}\rho_{(l)}D_{(l)}\\\
\phi_{(l)}\rho_{(l)}\varepsilon_{p,(l)}\\\ \rho u_{i}\\\ \rho E\\\
\overline{\mathbf{V}}^{e}_{ij}\\\
)+\partialderivative{x_{k}}\matrixquantity(\phi_{(l)}u_{k}\\\
\phi_{(l)}\rho_{(l)}u_{k}\\\ \phi_{(l)}\rho_{(l)}D_{(l)}u_{k}\\\
\phi_{(l)}\rho_{(l)}\varepsilon_{p,(l)}u_{k}\\\ \rho u_{i}u_{k}-\sigma_{ik}\\\
\rho Eu_{k}-u_{i}\sigma_{ik}\\\
\overline{\mathbf{V}}^{e}_{ij}u_{k}-\overline{\mathbf{V}}^{e}_{kj}u_{i}\\\
)=\matrixquantity(\phi_{(l)}\partialderivative{u_{k}}{x_{k}}\\\ 0\\\ 0\\\ 0\\\
0\\\ 0\\\
\frac{2}{3}\overline{\mathbf{V}}^{e}_{ij}\partialderivative{u_{k}}{x_{k}}-u_{i}\partialderivative{\overline{\mathbf{V}}^{e}_{kj}}{x_{k}}\\\
)+\matrixquantity(0\\\ 0\\\ \phi_{(l)}\rho_{(l)}\dot{D}_{(l)}\\\
\phi_{(l)}\rho_{(l)}\dot{\varepsilon}_{p,(l)}\\\ 0\\\ 0\\\ \bm{\Phi}_{ij}\\\
)+\mathbf{s}_{g}\ .$ (53)

When considering cylindrical symmetry, the conservative variables and
geometrical term are given by:

$\displaystyle\mathbf{U}=\matrixquantity(\phi_{(l)}\\\ \phi_{(l)}\rho_{(l)}\\\
\phi_{(l)}\rho_{(l)}D_{(l)}\\\ \phi_{(l)}\rho_{(l)}\varepsilon_{p,(l)}\\\ \rho
u_{r}\\\ \rho u_{z}\\\ \rho E\\\ \overline{\mathbf{V}}^{e}_{ij}\\\ ),\ \
\mathbf{s}_{g}=-\frac{1}{r}\matrixquantity(0\\\ \phi_{(l)}\rho_{(l)}u_{r}\\\
\phi_{(l)}\rho_{(l)}D_{(l)}u_{r}\\\
\phi_{(l)}\rho_{(l)}\varepsilon_{p,(l)}u_{r}\\\ \rho
u_{r}^{2}-\sigma_{rr}+\sigma_{\theta\theta}\\\ \rho u_{z}u_{r}-\sigma_{rz}\\\
\rho Eu_{r}-(u_{r}\sigma_{rr}+u_{z}\sigma_{zr})\\\
\frac{1}{3}\overline{\mathbf{V}}^{e}_{ij}u_{r}-\delta_{i\theta}\overline{\mathbf{V}}^{e}_{ij}u_{r})\
.$ (54)

In what follows, the underpinning method of Barton [7] is summarised first
before detailing the new modifications that facilitate the resolution of
interface slide, separation, and void.

###  3.1 Underpinning Method

Following Barton [7], the inhomogeneous system is integrated for time
intervals $[t^{n},t^{n+1}]$, where the time-step $\Delta t=t^{n+1}-t^{n}$ is
chosen to be a fraction of the global maximum allowable time step required for
stability of the hyperbolic update method [19]. For high strain-rate
applications, the plastic relaxation can occur over smaller time scales, which
can lead to local stiffness. To address this issue, without resorting to
forecasting stiff zones and resolving the time scales of irreversible physical
processes, Godunov’s method of fractional steps is used, where the hyperbolic
part is updated first, followed by serially adding in contributions from each
source term:

$\displaystyle\frac{\partial\mathbf{U}}{\partial t}$
$\displaystyle=-\frac{\partial\mathbf{F_{k}}}{\partial
x_{k}}+\mathbf{s}_{\text{\scriptsize{non-con.}}}$
$\displaystyle\text{IC:}\quad\mathbf{U}^{n}\xRightarrow{\Delta
t}{\mathbf{U}}^{\text{\char 73}}$ (55)
$\displaystyle\frac{\partial\mathbf{U}}{\partial t}$
$\displaystyle=\mathbf{s}_{x}$
$\displaystyle\text{IC:}\quad\mathbf{U}^{\text{\char 73}}\xRightarrow{\Delta
t}{\mathbf{U}}^{\bigstar}\ ,$ (56)

where the result of the each step is used as the initial condition (IC) for
the next. Here $\mathbf{s}_{x}=\mathbf{s}_{p},\mathbf{s}_{g}$. Once all source
terms have been added, the last $\mathbf{U}^{\bigstar}$ becomes
$\mathbf{U}^{n+1}$.

####  3.1.1 Hyperbolic Update

Replacing the spatial derivatives with a discretised conservative
approximation, the hyperbolic system can be written

$\frac{\text{d}}{\text{d}t}{\mathbf{U}}_{ijk}+\mathcal{D}_{ijk}\left({\mathbf{U}}\right)=0,$
(57)

where ${\mathbf{U}}_{ijk}$ represents the vector of conservative variables
stored at cell centres, and

$\displaystyle\mathcal{D}_{ijk}:=$ $\displaystyle\frac{1}{\Delta
x^{x}_{ijk}}\left(\widetilde{\mathbf{F}}^{x}_{i+1/2,jk}-\widetilde{\mathbf{F}}^{x}_{i-1/2,jk}\right)$
(58) $\displaystyle+$ $\displaystyle\frac{1}{\Delta
x^{y}_{ijk}}\left(\widetilde{\mathbf{F}}^{y}_{i,j+1/2,k}-\widetilde{\mathbf{F}}^{y}_{i,j-1/2,k}\right)$
$\displaystyle+$ $\displaystyle\frac{1}{\Delta
x^{z}_{ijk}}\left(\widetilde{\mathbf{F}}^{z}_{ij,k+1/2}-\widetilde{\mathbf{F}}^{z}_{ij,k-1/2}\right)-\mathbf{s}_{\text{\scriptsize{non-
con.}},ijk},$

where $\widetilde{\mathbf{F}}^{d}_{m\pm 1/2}$, for $m=i,j,k$, are the cell
wall numerical flux functions in the direction $d=x,y,z$. The numerical fluxes
are computed through successive sweeps of each spatial dimension and summed
according to equation (58).

Fluxes are computed using the HLLD solver from Barton [7] (but with the
additional flux-modifiers that shall be set out next). To achieve higher order
spatial accuracy, the initial conditions for the Riemann solver are taken to
be the MUSCL reconstruction of the cell centred primitive variables. It should
be noted that this procedure is generally carried out on the conservative
variables; the method instead follows Johnsen and Colonius [29], who showed
that, for multi-material problems, reconstruction of the primitive variables
leads to fewer oscillations around interfaces. An artificial interface
reconstruction is applied to reduce numerical diffusion around interfaces.
This is achieved using the Tangent of Hyperbola INterface Capturing (THINC)
method: an algebraic interface reconstruction technique that fits a hyperbolic
tangent function to variables inside a cell. In contrast to Barton [7], who
used the original THINC method of Xiao et al. [61], the more recent MUSCL-BVD-
THINC scheme of Deng et al. [22] is employed. This more recent scheme provides
an additional check to minimise oscillations by comparing the reconstructed
state’s cell boundary variation with the previously calculated MUSCL
reconstruction. THINC-reconstructed states are only accepted when their total
boundary variation is lower than that of the MUSCL scheme alone.

The non-conservative source terms that result from the volume fraction and
deformation tensor updates are added in the hyperbolic step following a mid-
point rule based procedure similar to that used by several authors [38, 1, 7].

To achieve a higher temporal resolution in the update of the hyperbolic terms,
the third order Runge-Kutta time integration scheme is used:

$\displaystyle\mathbf{U}^{(1)}$ $\displaystyle=$
$\displaystyle\mathbf{U}^{n}-\Delta t\mathcal{D}\left({\mathbf{U}}^{n}\right)$
(59) $\displaystyle\mathbf{U}^{(2)}$ $\displaystyle=$
$\displaystyle\mathbf{U}^{(1)}-\Delta
t\mathcal{D}\left({\mathbf{U}}^{(1)}\right)$ (60)
$\displaystyle\mathbf{U}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbf{U}^{n}+\frac{1}{4}\mathbf{U}^{(2)}$ (61)
$\displaystyle\mathbf{U}^{(4)}$ $\displaystyle=$
$\displaystyle\mathbf{U}^{(3)}-\Delta
t\mathcal{D}\left({\mathbf{U}}^{(3)}\right)$ (62)
$\displaystyle\mathbf{U}^{\text{\char 73}}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\mathbf{U}^{n}+\frac{2}{3}\mathbf{U}^{(4)}\ .$ (63)

It is remarked that the second order Runge-Kutta method could also be used,
but in practice it is found that the third-order method offers additional
robustness in some cases.

####  3.1.2 Plastic and Damage Update

The plastic update is performed after the hyperbolic step as detailed above
and consists of solving the equations:

$\displaystyle\dot{\varepsilon_{p}}$ $\displaystyle=\chi$ (64)
$\displaystyle\dot{\overline{\mathbf{V}}^{e}}$
$\displaystyle=-\chi\sqrt{\frac{3}{2}}\frac{\mbox{dev}\left({\sigma}\right)}{||{\mbox{dev}\left({\sigma}\right)}||}\overline{\mathbf{V}}^{e}\
.$ (65)

The potentially stiff ODEs governing the plasticity evolution are solved using
an analytical technique as detailed by Barton [7], which reduces the problem
to a single ODE for each material in a given cell. The algorithm can be
summarised as follows:

$\left(\overline{\mathbf{V}}^{e}\right)^{\bigstar}=\exp\left(\mbox{dev}\left({\mathbf{H}^{e}}\right)^{\bigstar}\right)$
(66)
$\mbox{dev}\left({\mathbf{H}^{e}}\right)^{\bigstar}=\frac{\mathcal{J}^{\bigstar}}{\mathcal{J}^{\text{\char
73}}}\mbox{dev}\left({\mathbf{H}^{e}\left({\overline{\mathbf{V}}^{e}}^{\text{\char
73}}\right)}\right)$ (67)
$\mathcal{J}^{\bigstar}=\frac{\sum_{l=1}^{N}\left(\phi_{(l)}/\Gamma_{(l)}\right)\mathcal{J}^{\bigstar}_{(l)}}{\sum_{l=1}^{N}\phi_{(l)}/\Gamma_{(l)}}$
(68)
$\mathscr{R}\left(\mathcal{J}_{(l)}^{\bigstar}\right)=\mathcal{J}_{(l)}^{\bigstar}-\mathcal{J}_{(l)}^{\text{\char
73}}+\Delta
t\sqrt{\frac{3}{2}}\chi_{(l)}\left(\mathcal{J}_{(l)}^{\bigstar},\varepsilon_{p_{(l)}}^{\bigstar}\left(\mathcal{J}_{(l)}^{\bigstar}\right)\right)=0\
.$ (69)

The last objective equation is solved using a simple bisection algorithm
between the limits
$\mathcal{J}_{(l)}^{\bigstar}\in[0:\mathcal{J}_{(l)}^{\text{\char 73}}]$.

This algorithm can be evaluated, at least in part, everywhere in a problem
irrespective of whether a material adheres to a plasticity law or not. This is
because the stretch tensor may lose its symmetry and unimodular properties as
a result of numerical error, and the algorithm encompasses the enforcement of

$\overline{\mathbf{V}}^{e}\leftarrow\sqrt{\overline{\mathbf{V}}^{e}\overline{\mathbf{V}}^{e^{\text{T}}}}\
,$ (70)

to ensure that the stretch tensor remains symmetric, and the condition
$\overline{\mathbf{V}}^{e}=\mathbf{I}$ is true for fluids. However, for tests
that do not feature materials that obey a plasticity law, performing this
algorithm is not necessarily required, as not enforcing the symmetry condition
was not found to significantly affect the results.

The implementation of this algorithm by Barton [7] employs singular value
decomposition (SVD) so that equation (66) and equation (67) are evaluated
exactly. Since SVD is relatively expensive, this part of the overall numerical
scheme constitutes a primary overhead. To alleviate this computational burden,
the following modifications are made. The Hencky strain tensor in equation
(67) is evaluated using a variant of the approximation by Bažant [11]:

$\mbox{dev}\left({\mathbf{H}^{e}}\right)\approx\frac{1}{2}\mathbf{H}_{B}^{e^{(1)}}\left(\overline{\mathbf{V}}^{e}{\overline{\mathbf{V}}^{e}}^{\text{T}}\right)\
,$ (71)

where

$\mathbf{H}_{B}^{e^{(m)}}\left(\overline{\mathbf{V}}^{e}\right)=\frac{1}{2m}\left({\overline{\mathbf{V}}^{e}}^{m}-{\overline{\mathbf{V}}^{e}}^{-m}\right)\
,$ (72)

and the invariants are evaluated from this in turn. The exponential matrix in
equation (66) is evaluated using the first diagonal Padé approximant:

$(\overline{\mathbf{V}}^{e})^{n+1}\approx\left(1-\frac{1}{2}\mbox{dev}\left({\mathbf{H}^{e}}\right)^{n+1}\right)^{-1}\left(1+\frac{1}{2}\mbox{dev}\left({\mathbf{H}^{e}}\right)^{n+1}\right)\
.$ (73)

This approximation of the strain was proposed in Barton [7] but was used only
where stresses are computed for the numerical flux functions. Use of the first
Padé approximate is found to be sufficient compared to second or higher
variants since the second norm of $\mathbf{H}^{e}$ is known from the outset
not to exceed sufficiently small values for the materials of interest. These
approximations avoid the use of a SVD, and the algorithm reduces to simple
matrix operations. The use of these approximations is further discussed in B.

Alongside this algorithm, the plastic strain and damage parameters are
updated:

$\displaystyle\dot{\varepsilon}_{p,(l)}$
$\displaystyle=\chi=\sqrt{\frac{2}{3}}\left(\mathcal{J}_{(l)}^{\text{\char
73}}-\mathcal{J}_{(l)}^{\bigstar}\right)/\Delta t$ (74)
$\displaystyle\dot{D}_{(l)}$
$\displaystyle=\dot{D}_{(l)}(\dot{\varepsilon}_{p,(l)})\ .$ (75)

###  3.2 Flux-modifiers

Since the underpinning model does not account for certain interface boundary
conditions, the fluxes in the hyperbolic update are modified in selected cells
in such a way that the desired boundary conditions are captured. This is the
crux of the new method, and the primary distinction compared to sharp
interface Ghost Fluid methods. In Ghost Fluid methods, the same numerical flux
is used across the entire domain, meaning ghost-regions of the domain must be
filled with states that provide suitable material boundary interaction. By
contrast, the new method uses a modified flux at the material interface,
encoding the effect of the material boundary condition into the flux method
itself, thereby not requiring a ghost-region to hold any additional states.
The following sections will therefore explain how appropriate material
boundary conditions are included in the modified fluxes.

In some cases when a flux-modifier is present, two different fluxes are
required across a cell interface. This two-flux approach has an approximately
equivalent effect to introducing a Ghost Fluid on either side, where each
side’s flux is based on the cell on one side and its own Ghost Fluid cell on
the other side. The presence of two fluxes modifies the conservative update
formula: what was

$\frac{\text{d}}{\text{d}t}{\mathbf{U}}_{ijk}+\frac{1}{\Delta
x^{x}_{ijk}}\left(\widetilde{\mathbf{F}}^{x}_{i+1/2,jk}-\widetilde{\mathbf{F}}^{x}_{i-1/2,jk}\right)=0\
,$ (76)

becomes

$\frac{\text{d}}{\text{d}t}{\mathbf{U}}_{ijk}+\frac{1}{\Delta
x^{x}_{ijk}}\left(\widetilde{\mathbf{F}}^{x}_{L,i+1/2,jk}-\widetilde{\mathbf{F}}^{x}_{R,i-1/2,jk}\right)=0\
.$ (77)

The method is illustrated in Figure 1.

The method sacrifices some degree of strict conservation at interfaces at the
benefit of allowing shear and void opening. However, it is emphasised that the
fluxes are modified only in a compact region around material interfaces – when
the flux-modifiers are not present the two fluxes at each cell interface are
the same, returning the method to the standard conservative update formula.
Moreover, in the case of shear, this flux imbalance applies predominantly to
the tangential vector quantities, retaining a greater degree of conservation
for scalar quantities such as mass. In the case of void opening, the full
effect of this imbalance is not felt as the flux into a void region is of
comparatively less significance as it will be over-written by the void seeding
routine. A two-flux approach such as this is not without precedent; the volume
fraction parts of the original Allaire et al. [1] model and the stretch tensor
parts of Barton [7] employ non-conservative update formulae. The extent of
these conservation loss effects will be quantified in Section 4.

The fluxes in the following sections are found by considering the Riemann
problem between two states. Various different flux methods are outlined, so
the notation $\mathbf{F}^{\bullet}(\mathbf{U}_{1},\mathbf{U}_{2})$ will denote
the flux between states $\mathbf{U}_{1}$ and $\mathbf{U}_{2}$ using the
$\bullet$ flux method, where $\bullet$ can be substituted for $n$
(normal/standard), $\nu$ (void) or $s$ (shear). The higher order methods and
interface reconstruction of the underpinning method can be included in the
usual way in the definition of the initial states $\mathbf{U}_{L}$,
$\mathbf{U}_{R}$ between which the flux is to be calculated.

Figure 1: The two-flux approach. The fluxes that the cell in the middle ‘sees’
for the purposes of the conservative update formula are highlighted in blue.

####  3.2.1 Shear

Physically, slip boundary conditions are produced at material interfaces where
the lack of adhesion between materials causes them to be able to slide along
each other in directions tangential to the interface normal. This is
represented mathematically as the tangential components of the interface
stress tensor being zero. Vector components normal to the interface and scalar
quantities are not affected by this condition, but tangential vector
components such as the transverse momentum/velocity are. A similar approach
has been used Favrie et al. [25] and Favrie and Gavrilyuk [24] for allowing
slip between solid-gas interfaces, but this was not extended to the case of
solid-solid interaction. Another successful approach has been led by Vitali et
al. [56], called the XEM method.

The underpinning algorithm proceeds until it is required to find the flux
across a material interface. Thus the two states between which the flux is to
be found, $\mathbf{U}_{L},\mathbf{U}_{R}$ have already been obtained.

Material interfaces are defined as where the volume fraction for any material
‘changes sign’ between the two cells being considered. The sign of the volume
fraction is defined as:

$\displaystyle\text{sign}(0.5-\phi)\ .$ (78)

The material for which this is the case is denoted the mixture material. When
three or more materials are present, there is a chance that several material
interfaces can coincide between the two cells being considered. In this case,
this method approximates the solution by selecting only one mixture material,
taken to be whichever of the volume fractions that change sign has the largest
difference between the two cells. This approximation was not found to affect
the results, but could be further examined in future.

In the method, a distinction is drawn between the ‘normal’ and ‘tangential’
state variables. The normal state variables consist of all scalar variables
(such as density) and the parts of tensorial quantities that are aligned with
the interface, and the tangential parts are the remainder. The normal parts of
a state $\mathbf{U}$ are denoted $\mathbf{U}^{\parallel}$ and the tangential
parts are denoted $\mathbf{U}^{\perp}$. For vector quantities such as momentum
and velocity, these different parts are relatively easy to define:

$\displaystyle\mathbf{u}^{\parallel}$
$\displaystyle=(\mathbf{u}\cdot\mathbf{\hat{n}})\mathbf{\hat{n}}$ (79)
$\displaystyle\mathbf{u}^{\perp}$
$\displaystyle=\mathbf{u}-\mathbf{u}^{\parallel}\ ,$ (80)

where $\mathbf{\hat{n}}$ is the outward pointing unit normal vector to the
mixture interface. For the stretch tensor, however, the process is slightly
more involved. It proves easier to rotate the tensor into the frame where the
interface normal lies along the $\mathbf{\hat{x}}$ direction (or any other
coordinate basis aligned frame). In what follows, states that have been
rotated into this reference frame will be denoted with a tilde,
$\mathbf{U}\rightarrow\tilde{\mathbf{U}}$. In this frame, the components of
the $x-$row and $x-$column, except for the component on the diagonal, are the
tangential parts, and the rest are the normal parts:

(81)

The normal and tangential parts in the original reference frame are then found
by reversing the rotation:

$\displaystyle\mathbf{\overline{\mathbf{V}}^{e}}^{\parallel}$
$\displaystyle=\mathbf{R}^{\text{T}}\tilde{\mathbf{V}}^{\parallel}\mathbf{R}$
(82) $\displaystyle\mathbf{\overline{\mathbf{V}}^{e}}^{\perp}$
$\displaystyle=\mathbf{R}^{\text{T}}\tilde{\mathbf{V}}^{\perp}\mathbf{R}\ .$
(83)

Here $\mathbf{R}$ is the rotation matrix which maps the normal vector
$\mathbf{\hat{n}}$ to the $x$-axis $\mathbf{\hat{x}}$. In two dimensions,
$\mathbf{R}$ is a simple rotation around the $z$-axis, given in terms of the
components of the normal vector by:  The shear flux-modifier: 1. First, the
outward pointing interface normal, $\mathbf{\hat{n}}$, is estimated using
Young’s method [62] applied to the mixture material volume fraction, and
averaged between the two cells. 2. $\mathbf{U}_{L}$ and $\mathbf{U}_{R}$ are
rotated into the reference frame where the normal $\mathbf{\hat{n}}$ points
along $\mathbf{\hat{x}}$. 3. Two new states are constructed:
$\tilde{\mathbf{U}}_{L}^{s},\tilde{\mathbf{U}}_{R}^{s}$. These contain the
normal parts of one side, combined with the tangential parts of the other
side: $\displaystyle\tilde{\mathbf{U}}_{L}^{s}$ $\displaystyle=$
$\displaystyle\tilde{\mathbf{U}}_{L}^{\perp}+\tilde{\mathbf{U}}_{R}^{\parallel}$
$\displaystyle\tilde{\mathbf{U}}_{R}^{s}$ $\displaystyle=$
$\displaystyle\tilde{\mathbf{U}}_{R}^{\perp}+\tilde{\mathbf{U}}_{L}^{\parallel}\
.$ 4. The appropriate stress tensor is then calculated in these new states,
where the tangential stress components are zero. This can be done in two ways,
see Section 3.2.2. 5. All the aforementioned cells are then rotated back into
the original frame of reference:
$\displaystyle\mathbf{U}=\mathbf{R}^{\text{T}}\tilde{\mathbf{U}}\mathbf{R}\ .$
6. The standard flux is then calculated between each pair of states:
$\displaystyle\mathbf{F}^{s}_{L}$ $\displaystyle=$
$\displaystyle\mathbf{F}^{n}(\mathbf{U}_{L},\mathbf{U}_{L}^{s}),$
$\displaystyle\mathbf{F}^{s}_{R}$ $\displaystyle=$
$\displaystyle\mathbf{F}^{n}(\mathbf{U}_{R}^{s},\mathbf{U}_{R}).$ This
procedure predominantly affects the tangential components of any state.
However, in two and three dimensions, the interface will not generally be
coordinate basis aligned, so the procedure should modify components that are
normal to a given hyperbolic sweep direction, as they will have some component
pointing along the interface normal.

As it stands, the shear flux-modifier alone would not be sufficient. If only
this flux-modifier were applied, the tangential components would not move with
the interface normal velocity. Before every time step, a seeding routine is
first performed. This routine propagates the correct tangential velocity
components into mixed cell regions. Mixed cells are defined as those either
next to a material interface or having a volume fraction for any material that
is in the range $0.5<\phi<1.0-\delta$ where delta is some small amount, taken
in this work to be $1\times 10^{-4}$. The volume fraction for which this
inequality applies is then the mixture material. The following procedure is
then employed:  The shear seeding routine: 1. The outward pointing interface
normal, $\mathbf{\hat{n}}$, is calculated using Young’s method [62]. 2. A
probe is sent out along the normal direction using the method described in
Section 3.3, and the material values are interpolated at that point to give a
new state $\mathbf{U}_{\text{Interp}}$. 3. If the interpolated mixture
material volume fraction is higher than in the original cell, both the
original and the interpolated cell are rotated so the normal lies in the
$\mathbf{\hat{x}}$ direction, giving $\tilde{\mathbf{U}}$ and
$\tilde{\mathbf{U}}_{\text{Interp}}$. 4. A new state is then created
containing the tangential components of the interpolated cell and the normal
components of the original cell, in the same way as the flux-modifier:
$\tilde{\mathbf{U}}^{s}=\tilde{\mathbf{U}}^{\parallel}+\tilde{\mathbf{U}}_{\text{Interp}}^{\perp}\
.$ 5. This new state is then rotated back into the original frame of reference
and replaces the original state. $\mathbf{U}\leftarrow\mathbf{U}^{s}\ .$ The
procedure is found to adequately seed the mixture regions with the suitable
tangential components. If correctly applied, this procedure should only affect
mixture regions, and with the application of interface reconstruction, these
regions are generally small and the process remains localised.

####  3.2.2 Strain-free and Stress-free Tangential Conditions

To calculate the correct interface stress tensor in the shear flux-modifier
states $\tilde{\mathbf{U}}_{L}^{s},\tilde{\mathbf{U}}_{R}^{s}$ there are two
options. The first is to leave the stretch tensor unchanged and simply set the
tangential components of the stress tensor to zero. This corresponds to the
stress-free condition. In this intuitive approximation, the strain-free nature
of the interface will emerge naturally as the simulation progresses. The
second option is to calculate the strain-free state that will result in a
stress-free interface. This approach is somewhat more involved, but more
physically grounded. Here, a new stretch tensor is calculated that results in
the tangential stress being zero when it is calculated in the ordinary way.
This strain-free state is calculated in a similar way to how the stretch
tensor is reconstructed during the plastic update. The deviatoric Hencky
strain is obtained, its tangential components are set to zero, and then it is
reconstructed into the stretch tensor, following the procedure laid out in
Section 3.1.2. In the course of this work, it was found that both techniques
produce satisfactory results, but the second, more physical method is
preferred.

####  3.2.3 Void

To enable interaction with voids such as those produced when two solid
materials separate, an additional scalar variable, the void volume fraction,
$\nu$, is introduced. The void volume fraction is advected in exactly the same
way as material volume fractions in the standard flux, by solving

$\displaystyle\frac{\partial\nu}{\partial t}+\frac{\partial\nu u_{k}}{\partial
x_{k}}=\nu\frac{\partial u_{k}}{\partial x_{k}}\ .$ (86)

The void flux is then modified in two scenarios:

  1. 1.

Either cell $\mathbf{U}_{L}$, $\mathbf{U}_{R}$ already contains some void –
the ‘pre-existing void’ case.

  2. 2.

The cells are at a material interface or are critically damaged and are coming
apart. The ‘void generation’ case.

In the case where the void already exists and doesn’t need to be generated,
the stress tensor in each cell is weighted by the void volume fraction:

$\bm{\sigma}\leftarrow(1-\nu)\bm{\sigma}$ (87)

and the standard flux can be used. This in itself is insufficient to model
realistic void interaction and a void seeding routine is employed in a similar
way to the shear flux method. This will be outlined below.

In addition, this work presents two options for the void generation flux-
modifier. The first is a ‘simple’ method, akin to the original Ghost Fluid
method of Fedkiw et al. [26], where quantities in the ghost cells are simply
copied across from the real material. The second is a Riemann-solver-based
method akin to the later Ghost Fluid Methods of Sambasivan and Udaykumar [51]
and Wang et al. [60]. Both methods are found to work satisfactorily, with the
simple method being cheaper and generally more robust, but the Riemann-
problem-based method retaining greater conservation for strong shocks.

Other methods with which this work can be compared and contrasted include
Klima et al. [35] and Barlow et al. [3]

The void generation flux-modifier is applied in all interfacial cells that are
coming apart. Material interfaces are flagged in the same way as the shear
flux-modifier, namely when a volume fraction ‘changes sign’ between the two
cells. Damage interfaces are defined as where either cell between which the
flux is to be calculated is critically damaged (this fracture criterion is
explained in detail in the next section). Interfaces are deemed to be moving
apart when they are under tension:

$\displaystyle\mathbf{\hat{n}}^{\text{T}}\cdot\left(\bm{\sigma}_{R}+\bm{\sigma}_{L}\right)\cdot\mathbf{\hat{n}}>0\
.$ (88)

In this case, the following procedure is carried out:  The void flux-modifier:
1. The interface normal $\mathbf{\hat{n}}$ is calculated using Young’s method
[62]. 2. Two new states are constructed:
$\mathbf{U}_{L}^{\nu},\mathbf{U}_{R}^{\nu}$. • For the ‘simple’ method, these
states are directly copied from their counterpart:
$\displaystyle\mathbf{U}_{L}^{\nu}$ $\displaystyle=$
$\displaystyle\mathbf{U}_{R}$ $\displaystyle\mathbf{U}_{R}^{\nu}$
$\displaystyle=$ $\displaystyle\mathbf{U}_{L}\ .$ • For the Riemann-problem-
based method, these contain the star-states of the interaction
$\mathbf{U}_{L}$ and $\mathbf{U}_{R}$ with void, found with the HLLD solver
using a void boundary condition: $\displaystyle\mathbf{U}_{L}^{\nu}$
$\displaystyle=$ $\displaystyle\mathbf{U}_{L}^{**}$
$\displaystyle\mathbf{U}_{R}^{\nu}$ $\displaystyle=$
$\displaystyle\mathbf{U}_{R}^{**}\ .$ 3. The void volume fraction in these new
states is then set to 1, and the appropriate weighted stress tensor is then
calculated. 4. The standard flux is then calculated between each pair of
states: $\displaystyle\mathbf{F}^{\nu}_{L}$ $\displaystyle=$
$\displaystyle\mathbf{F}^{n}(\mathbf{U}_{L},\mathbf{U}_{L}^{\nu}),$
$\displaystyle\mathbf{F}^{\nu}_{R}$ $\displaystyle=$
$\displaystyle\mathbf{F}^{n}(\mathbf{U}_{R}^{\nu},\mathbf{U}_{R}).$

To enable the correct interaction with void states, the thermodynamic
variables in areas with high void volume fraction present are seeded with
values that mimic the interaction of the flow with void. This is somewhat akin
to the Ghost Fluid method’s filling of the ghost cells with the material
boundary condition. The following procedure is used:  The void seeding
routine: 1. In cells with a void volume fraction above a specified threshold
(generally this work takes $\nu=0.9$, but a range of values is found to work),
the void normal $\mathbf{\hat{n}}_{\nu}$ is calculated using Young’s method
[62]. 2. A probe is sent out along the normal direction using the method
described in Section 3.3 and the material values are interpolated at that
point to give a new state $\mathbf{U}_{\text{Interp}}$. 3. If the interpolated
void volume fraction is lower than the void volume fraction of the original
cell, a new state $\mathbf{U}^{\nu}$ is calculated: • For the ‘simple’ method,
this new state is directly copied from the interpolated cell:
$\mathbf{U}^{\nu}=\mathbf{U}_{\text{Interp}}.$ • For the Riemann-problem-based
method, this new state is the star-state of the interaction of
$\mathbf{U}_{\text{Interp}}$ with void. This is found as before by rotating
$\mathbf{U}_{\text{Interp}}$ into the direction where the normal
$\mathbf{\hat{n}}_{\nu}$ points along $\mathbf{\hat{x}}$, and then using the
HLLD solver with a void boundary condition. This star-state then becomes the
new state: $\tilde{\mathbf{U}}^{\nu}=\tilde{\mathbf{U}}^{**}\ .$ This state
$\tilde{\mathbf{U}}^{\nu}$ is then rotated back to the original reference
frame. 4. The new state is then linearly combined with the old state, but
keeping history variables such as damage, plastic strain and volume fractions
constant: $\mathbf{U}\leftarrow(1-\nu)\mathbf{U}+\nu\mathbf{U}^{\nu}\ .$ Again
this seeding routine only affects void mixture regions, and with interface
reconstruction this procedure is kept localised.

###  3.3 Probe Calculation

The seeding routines require probes to be sent out in different directions in
order to interpolate material values at different points. The calculation of
these probes is explained in this section. In order to start, the current cell
centre position $\mathbf{c}$, the interface normal (probe direction)
$\mathbf{\hat{n}}$, and vector of cell dimensions $\mathbf{\differential x}$,
are assumed to be known. In two and three dimensions, if the normal does not
point along a coordinate basis direction, the distance a probe has to travel
to reach another cell is greater, as the probe travels along the diagonal of
the cell. Therefore, in order to ensure that probe directions along the
diagonal can still reach suitable material values outside the current cell,
the probe vector is scaled depending on the normal direction. This is
performed in a dimension-independent way by stretching the unit sphere defined
by all normal directions to fill a cube with side length 2 (in two dimensions
this corresponds to stretching the unit circle to fill a square with side
length 2). This stretching is denoted
$\mathbf{\hat{n}}\rightarrow\mathbf{\hat{n}}^{\prime}$. Figure 2 demonstrates
this mapping in two dimensions.

Figure 2: Mapping the unit circle to a square with side length 2. This stretch
ensures that the probe can reach exterior cell centres in diagonal directions.

To perform the mapping, firstly the component of $\mathbf{\hat{n}}$ with the
largest absolute value (either positive or negative) is denoted
$n_{\text{max}}$. The components of $\mathbf{\hat{n}}^{\prime}$ are then given
by:

$\displaystyle n_{i}^{\prime}=\frac{n_{i}}{|n_{\text{max}}|}\ .$ (89)

Once this stretch has been performed, the final probe position is given by:

$\mathbf{p}=\mathbf{c}+l\cdot\left(\mathbf{\hat{n}}^{\prime}\odot\mathbf{\differential
x}\right)\ ,$ (90)

where $\odot$ is element-wise multiplication. The factor $l$ modifies the
probe length, and values between 1 and 1.5 were found to work well for the
problems presented.

###  3.4 Fracture

One advantage of the new method is that fracture is handled naturally by the
same routines that mediate interface separation. Critically damaged cells,
defined as any cell where:

$D_{(l)}\geq D_{\text{crit},(l)}\ ,$ (91)

are simply treated as any other material interface where the void generation
flux-modifier can be applied. No additional algorithm needs to be applied to
initiate fracture. In this case, the interface normal can be taken as the
current sweep direction. This approach also removes the need for the highly
non-conservative approach of deleting critically damaged cells to form
fracture, as is commonplace in finite element and level set based codes [18],
[16], [6]. This also avoids the issue of potentially having to redistribute
the conserved quantities such a mass and momentum to the neighbouring cells
around the deleted cells.

However, this method does come with its own associated challenges; the
diffusive nature of all Eulerian codes means that damaged cells may undergo
artificial healing. As an initially critically damaged material is advected,
diffusion can ‘heal’ the damage by causing it to fall under the critical
damage threshold. This causes a number of issues:

  * •

Regions that should be able to come apart and form a crack are prevented from
separating.

  * •

The transition to critically damaged and back is discontinuous; quantities
such as the shear modulus and yield strength go suddenly to zero.

  * •

When a history-dependent plasticity model such as strain-hardening is
employed, the plastic strain remains unchanged through healing. This means a
material can be damaged, heal, and then continue to accrue plastic strain,
leading to unphysical increased hardening.

  * •

It prevents the use of a stochastic initial damage profile to mimic the effect
of defects in the structure of real materials, as high-frequency modes in the
initial distribution are quickly smeared out.

Previous works, such as Vitali and Benson [57] and Barton [6], have employed
semi-Lagrangian updates to avoid these issues. This technique could also be
applied to the current method, but is beyond the scope of the current study.
In this work, THINC reconstruction is applied to the damage field and this is
found to be sufficient to allow the large strain rate ductile fracture in
problems at hand. Additionally, this work does not apply a stochastic initial
condition to the damage profile to mimic material defects and ageing. Rather,
the aim is to demonstrate the validity of the underlying fracture mechanics,
while acknowledging that the method can be extended with techniques similar to
Vitali and Benson [57] and Barton [6] if so desired.

##  4 Validation and Verification

This section provides both numerical and experimental comparison for the new
method. In two- and three-dimensional tests, materials are delineated from
void by the $\nu=0.5$ volume fraction contour, with void areas shown in white
for clarity.

Material |  $\rho_{0}$ kgm-3 |  $K_{0}$ GPa |  $G_{0}$ GPa | $\alpha$ | $\beta$ | $\Gamma_{0}$
---|---|---|---|---|---|---
Aluminium | 2703.0 | 76.3 | 26.36 | 0.627 | 2.288 | 1.484
Al 5083 H32 | 2670.0 | 72.2 | 25.8 | 0.63 | 2.29 | 1.48
Copper | 8930.0 | 136.45 | 39.38 | 1.0 | 3.0 | 2.0
CuBe | 8370.0 | 131.3 | 53.6 | 1.0 | 3.0 | 2.0
Steel | 8030.0 | 156.2 | 77.2 | 0.569 | 2.437 | 1.563
Polycarbonate | 1193.0 | 5.52 | 1.18 | 1.49 |  * | 0.61
AerMet | 7940.0 | 156.2 | 77.2 | 0.569 | 2.437 | 1.84
Table 1: Romenskii equation of state material parameters. *A constant shear
modulus variant was used for this material.

###  4.1 Solid-Solid Riemann Problem

Firstly, the shear algorithms are verified in one dimension. This test is a
Riemann problem taken from Barton [7] involving the collision of a stressed
piece of aluminium with a stressed piece of copper. The test is challenging,
representing a highly unrealistic initial stress configuration. Moreover, the
problem features 7 different waves: 2 longitudinal, 4 transverse and one
contact. The problem has an exact solution with which this work can be
compared to find an estimate for the convergence; further details on the exact
solution are provided in C. The initial conditions are given by:

where $F$ is the deformation tensor. Note that the stretch tensor,
$\overline{\mathbf{V}}^{e}$, can be found from $F$ using a singular value
decomposition. A full, exact singular value decomposition can be used for this
task as this is only required once at the start of a simulation. The test is
run at a resolution of 500 cells on a domain $x=[0,0.01]$ m with an interface
at $x=0.005$ m, for a time of 0.5 $\mu$s, using a CFL of 0.9. For this initial
test it is assumed that the materials behave purely elastically and both
materials are governed by the Romenskii equation of state used by Barton [7]
(see section 2.3), using the material parameters laid out in Table 1.

The results are shown in Figure 3, where comparison is drawn between both slip
and stick boundary conditions. It can be seen that the method performs well in
this tough test, with a clean tangential velocity discontinuity, smooth zero
tangential stress, and unaffected longitudinal waves. The slip boundary
condition also removes the small THINC-induced oscillations in the $xz$-stress
noticed by Barton [7].

Table 2 shows the convergence for both stick and slip cases for the density
and velocities. In both cases, the density and $x$-velocity show a roughly
first order convergence, as expected for such a shock-dominated problem. This
also verifies that the flux-modifier does not affect the convergence in the
normal direction. The convergence for the two transverse velocities are
slightly lower, but importantly it is better in the slip case where the flux-
modifier is applied, even despite the additional contact wave that must be
resolved. This suggests that flux-modifiers do not degrade the convergence of
the underlying solution.

Figure 3: The solid-solid Riemann problem. This test is used to provide one-dimensional shear verification, comparing the case with the shear flux-modifier and the case without to their exact solutions. Excellent agreement is observed in both cases for this strenuous test. Method | Resolution | Density |  $x$-velocity |  $y$-velocity |  $z$-velocity
---|---|---|---|---|---
$L_{1}$ Error | Order |  $L_{1}$ Error | Order |  $L_{1}$ Error | Order |  $L_{1}$ Error | Order
Stick |  |  |  |  |  |  |  |  |
100 | 0.05653 |  | 0.01548 |  | 0.00143 |  | 0.00116 |
200 | 0.02428 | 1.21917 | 0.00793 | 0.96479 | 0.00086 | 0.74022 | 0.00084 | 0.46520
400 | 0.01286 | 0.91657 | 0.00427 | 0.89176 | 0.00050 | 0.78438 | 0.00054 | 0.64268
800 | 0.00710 | 0.85778 | 0.00198 | 1.11011 | 0.00031 | 0.69452 | 0.00033 | 0.70210
1600 | 0.00343 | 1.05106 | 0.00120 | 0.71821 | 0.00020 | 0.60264 | 0.00022 | 0.61283
3200 | 0.00205 | 0.74423 | 0.00058 | 1.06062 | 0.00014 | 0.54894 | 0.00014 | 0.66939
Slip |  |  |  |  |  |  |  |  |
100 | 0.05656 |  | 0.01551 |  | 0.00240 |  | 0.00304 |
200 | 0.02427 | 1.22096 | 0.00787 | 0.97837 | 0.00135 | 0.82735 | 0.00163 | 0.90084
400 | 0.01282 | 0.92021 | 0.00418 | 0.91192 | 0.00078 | 0.78909 | 0.00086 | 0.92802
800 | 0.00722 | 0.82829 | 0.00204 | 1.03288 | 0.00051 | 0.61001 | 0.00047 | 0.86087
1600 | 0.00339 | 1.09041 | 0.00111 | 0.87972 | 0.00031 | 0.71911 | 0.00027 | 0.77833
3200 | 0.00205 | 0.72749 | 0.00057 | 0.95303 | 0.00019 | 0.74642 | 0.00017 | 0.66865
Table 2: Convergence for the solid-solid Riemann problem. The proposed flux-
modifier does not degrade the convergence of the underlying method.

Additionally, this test can be run in two dimensions with the interface normal
rotated by and angle of $45^{\circ}$ to verify that the algorithm can equally
cope with non-grid-aligned interfaces. The initial conditions for the test are
the same, except the velocity and deformation tensors are rotated to account
for the rotated interface. The test is then run using a CFL of 0.6 and
resolution of $250\times 250$ cells, with one layer of AMR of refinement
factor 2. Figure 4 shows this rotated test, where a cut has been taken through
the domain along the line $(x,y)=(1,1)$ and the spatial axis has been scaled
by a factor of $\sqrt{2}$ to make the results comparable with the exact
solution.

Figure 4: The rotated solid-solid Riemann problem. This test is run in 2D
using an interface normal of $(x,y)=(1,1)/\sqrt{2}$ to test the shear flux-
modifier in the non-grid-aligned case. A slice is then taken through the
domain to compare to the 1D case and the exact solution (Figure 3), and the
results match well.

###  4.2 Solid-Vacuum Riemann Problem

The interaction of solids with vacuum in one dimension is now examined. This
is done using a test from Barton [7] that features a high-speed stressed block
of aluminium in contact with vacuum. The initial conditions are given by:

The problem features a rarefaction wave, 2 transverse shocks in the solid, and
a material-vacuum contact wave. The test is run for 0.6 $\mu$s on a domain
$x=[0:10]$ mm using a CFL 0.9 with the initial discontinuity at $x=5$ mm.

Figure 5: The solid-vacuum Riemann problem. This problem tests the interaction
of a stressed solid with an area of pre-existing void. The work at hand is
compared to the exact solution, and matches very well, keeping a sharp
material-void interface while resolving the multiple waves within the
material.

The results are shown in Figure 5. Again, the results are compared with the
exact solution. The method handles the strong velocity discontinuities across
the material interface well, for both normal and tangential components. This
test also matches the results in Barton [7] well, where the void was instead
replaced with air. Due to the large material differences between air and
aluminium it is expected that the wave structure inside the solid is very
similar in both the air- and vacuum-case, which is observed in the results of
this test. This is an example of where this method can be particularly useful
- the solution inside the air has little effect on the wave structure in the
solid, and replacing the air with void saves computational resources. The
effect of resolution is also examined, as shown in Figure 6 and Table 3.
Figure 6 compares multiple resolutions to the exact solution, derived using
the method laid out in C. For practical resolutions, it can be seen that the
method matches the exact solution very well. Additionally, Table 3 shows the
void algorithms also maintain the convergence to the exact solution, even in
this strenuous test.

Figure 6: The solid-vacuum Riemann problem shown for various different resolutions. For practical resolutions, the method matches the exact solution extremely well. Resolution | Density | $\sigma_{xx}$ |  $x$-velocity |  $y$-velocity
---|---|---|---|---
$L_{1}$ Error | Order |  $L_{1}$ Error | Order |  $L_{1}$ Error | Order |  $L_{1}$ Error | Order
100 | 0.02609 |  | 0.24831 |  | 0.03114 |  | 0.00112 |
200 | 0.01024 | 1.34911 | 0.12300 | 1.01344 | 0.01312 | 1.24658 | 0.00053 | 1.06342
400 | 0.00554 | 0.88558 | 0.05958 | 1.04568 | 0.00698 | 0.91040 | 0.00032 | 0.74531
800 | 0.00329 | 0.75150 | 0.02677 | 1.15449 | 0.00396 | 0.81883 | 0.00019 | 0.72190
1600 | 0.00221 | 0.57811 | 0.01210 | 1.14490 | 0.00253 | 0.64459 | 0.00013 | 0.60614
3200 | 0.00110 | 0.99971 | 0.00593 | 1.02980 | 0.00134 | 0.91885 | 0.00009 | 0.48300
Table 3: Convergence for the solid-vacuum Riemann problem.

###  4.3 Gas-Vacuum Riemann Problem

The void components are validated in one dimension using an ideal gas vacuum
Riemann problem, with a $\gamma=1.4$ ideal gas on the left and vacuum on the
right. This is a notoriously strenuous test of the material-void interaction
which many approximate Riemann solvers struggle to capture, especially in
regards to the interface velocity and internal energy of the gas. The test
features only two waves: a strong rarefaction back into the material and the
material-void contact discontinuity. This test is difficult because the
density of the gas at the material interface is incredibly small, while its
velocity is large. For an ideal gas, the test has an exact solution with which
this method can be compared to find an estimate for the convergence rate. The
test is non-dimensional and is run until a time of $t=0.1$ over a domain of
$x=[0:1.5]$, with varying resolution using a CFL of 0.9. The initial
conditions for the problem are:

$\displaystyle\mbox{Left: }\ (x<x_{0})\quad\rho_{L}=1,\quad
p_{L}=2.5,\quad\mathbf{u}_{L}=\mathbf{0},\quad\nu_{L}=0$ (125)
$\displaystyle\mbox{Right: }\ (x>x_{0})\quad\rho_{R}=0,\quad
p_{R}=0,\quad\mathbf{u}_{R}=\mathbf{0},\quad\nu_{R}=1\ ,$ (126)

where the interface $x_{0}$ is taken to be 0.3. The exact solution is given
by:

$\displaystyle\matrixquantity(\rho&&\mathbf{u}&&p)$
$\displaystyle=\left\\{\matrixquantity{\matrixquantity(\rho_{L}&&\mathbf{u}_{L}&&p_{L})&&\text{if}\quad\frac{x-x_{0}}{t}<-c_{0}\\\
\matrixquantity(\rho^{*}&&\mathbf{u}^{*}&&p^{*})&&\text{if}\quad\frac{x-x_{0}}{t}<\frac{2c_{0}}{\gamma-1}\\\
\matrixquantity(\rho_{R}&&\mathbf{u}_{R}&&p_{R})&&\text{if}\quad\frac{x-x_{0}}{t}>\frac{2c_{0}}{\gamma-1}\\\
}\right.\ ,$ (127) where: $\displaystyle\mathbf{u}^{*}$
$\displaystyle=2\frac{\frac{x-x_{0}}{t}+c_{0}}{\gamma-1}$ (128)
$\displaystyle\rho^{*}$
$\displaystyle=\left(\left(\mathbf{u}^{*}-\frac{x-x_{0}}{t}\right)^{2}\left(\frac{\rho_{L}^{\gamma}}{\gamma
p_{L}}\right)\right)^{\frac{1}{\gamma-1}}$ (129) $\displaystyle p^{*}$
$\displaystyle=p_{L}\left(\frac{\rho^{*}}{\rho_{L}}\right)^{\gamma}\ .$ (130)

Figure 7: The gas-vacuum Riemann problem. This difficult test examines the interaction of a fluid with a pre-defined area of void, comparing this method to the analytical solution, and good agreement is seen in the density and pressure profiles. As with all the tests presented here, THINC reconstruction is applied to the diffuse gas-vacuum boundary to keep it sharp. Figure 8: The density convergence rate for the gas-vacuum Riemann problem, comparing the work at hand to a Ghost Fluid method. Both methods show first order convergence. $N$ is the number of cells used in the test. Method | Resolution | Density | Pressure | Internal Energy
---|---|---|---|---
$L_{1}$ Error | Order |  $L_{1}$ Error | Order |  $L_{1}$ Error | Order
This work |  |  |  |  |  |  |
100 | 0.02567 |  | 0.04737 |  | 0.56130 |
200 | 0.01344 | 0.93393 | 0.02226 | 1.08937 | 0.37705 | 0.57400
400 | 0.00656 | 1.03531 | 0.01034 | 1.10613 | 0.20617 | 0.87096
800 | 0.00329 | 0.99276 | 0.00521 | 0.99035 | 0.09981 | 1.04659
1600 | 0.00165 | 0.99761 | 0.00262 | 0.99298 | 0.04236 | 1.23633
3200 | 0.00083 | 0.99907 | 0.00131 | 0.99588 | 0.02160 | 0.97203
6400 | 0.00041 | 0.99931 | 0.00066 | 0.99649 | 0.01725 | 0.32432
Ghost Fluid |  |  |  |  |  |  |
100 | 0.02578 |  | 0.04528 |  | 0.55784 |
200 | 0.01666 | 0.62938 | 0.02666 | 0.76402 | 0.44393 | 0.32953
400 | 0.00944 | 0.82010 | 0.01390 | 0.93990 | 0.33020 | 0.42699
800 | 0.00485 | 0.96168 | 0.00668 | 1.05690 | 0.21812 | 0.59821
1600 | 0.00243 | 0.99756 | 0.00329 | 1.02356 | 0.13510 | 0.69105
3200 | 0.00122 | 0.99393 | 0.00164 | 1.00159 | 0.08101 | 0.73792
6400 | 0.00061 | 0.99681 | 0.00082 | 0.99978 | 0.04695 | 0.78708
Table 4: The gas-vacuum Riemann problem test convergence. Here, the method at
hand is compared to a level set Ghost Fluid method, and matches the
convergence well.

This test is shown in Figure 7. This method is also compared with a Ghost
Fluid method, and Figure 8 and Table 4 demonstrate both of these methods show
first order convergence to the exact solution in the density and pressure, as
would be expected from this kind of interface problem. However, the internal
energy converges slower for both methods, and the velocity will eventually
overshoot the exact solution in both cases. It should be stressed that
detailed fluid-void interaction studies are not the intended purpose of this
method; other exact methods may be more appropriate in these scenarios.
Despite this, this test shows that the method performs very capably compared
to existing sharp interface methods in a strenuous test, maintaining its
convergence.

###  4.4 Sliding Blocks Test

For straightforward two dimensional shear verification, a simple sliding
blocks test was considered. In this test, two copper blocks slide past each
other, surrounded by vacuum. The velocities of the two blocks are taken along
the diagonal of the domain to verify there is no grid-dependence in the shear
algorithm. The blocks are both 1.18 mm $\times$ 2.36 mm with a chamfer of 0.17
mm and are given initial velocities:

$\displaystyle\mathbf{u}=\matrixquantity(\pm 200\\\ \pm 200\\\ 0)\
\text{ms}^{-1}\ .$ (149)

The test is run using a base resolution of 100 $\times$ 100 cells with 3
layers of AMR, using a CFL of 0.6, for a time of 5 $\mu$s. Results are shown
in Figure 9. The blocks can be seen to slide smoothly along each other,
despite the diagonal interface, and separate cleanly. Figure 10(a) depicts the
conservation error over the course of this test. The error remains low, and is
predominantly generated by the quasi-non-conservative void volume fraction
update, rather than the flux-modifier. Figure 10(b) shows the $x$-velocity
along a cut normal to the interface, where the THINC reconstruction and flux-
modifiers help to keep the interfaces sharp.

Figure 9: The sliding blocks test. Times are shown in $\mu$s. The blocks are
surrounded by void, and the shear flux-modifier allows them to slide past each
other, separating cleanly at the end.

(a) The conservation error for the sliding blocks test, showing mass, momentum
and energy loss over the course of the simulation. The method at hand performs
very well in terms of conservation, even at this relatively low resolution.

(b) A cut of the $x$-velocity along a line normal to the interface in the
sliding blocks test, taken at a time of 1.25 $\mu$s. The plot shows the THINC
reconstruction and shear flux-modifier keeping the velocity discontinuities
sharp, despite the diffuse interface representation.

Figure 10: Further detail on the sliding blocks test.

###  4.5 Bouncing Blocks and Cylinders

For a simple two dimensional verification of the void-generation parts of the
algorithm, two elastic solid impact tests were considered; a bouncing blocks
test and a bouncing cylinders test. These straightforward tests features two
purely elastic aluminium materials governed by the Romenskii equation of state
colliding and separating in a vacuum. The velocities of the two bodies in
either case are taken along the diagonal of the domain to verify there is no
grid-dependence in the void-generation algorithm, and the two different
geometries test both the planar and non-planar interface separation
capabilities. The blocks are both 1 mm $\times$ 2.35 mm, with an initial
separation of 0.35 mm. The cylinders both have radii of 0.83 mm, with their
centres having an initial separation of 1.53 mm. In both tests, the two bodies
are given initial velocities:

$\displaystyle\mathbf{u}=\matrixquantity(\pm 500\\\ \pm 500\\\ 0)\
\text{ms}^{-1}\ .$ (150)

The tests are run using a base resolution of 100 $\times$ 100 cells with 3
layers of AMR, using a CFL of 0.6, for a time of 1.6 $\mu$s. As can be seen in
Figures 11 and 12, the two bodies separate cleanly in both cases. It is in
this kind of test that the void components of this algorithm give large
computational savings. This test would otherwise have to be performed in air,
with the simulation having to capture all of the complex but insignificant
waves patterns in the air. Moreover, there would likely be large pressure
spikes in the air that would be trapped between the blocks, leading to time
step restrictions. Both of these problems are avoided with the work at hand.

Figures 13(a) and 13(b) depict the conservation errors of this scheme over the
course of the two simulations. In both cases, the mass and momentum
conservation is very good, with neither test going above or below a 0.2%
change. The total energy conservation is more erratic, but well within the
expected bounds for this kind of Eulerian scheme. The total energy
conservation error is larger in the case of the bouncing blocks both because
there is a larger contact surface, giving increased use of the flux-modifiers,
and the geometry for the blocks results in a stronger shock wave hitting the
back of the blocks compared to the circular case.

Figure 11: The bouncing blocks test, depicting the diagonal velocity. The
times shown are in $\mu$s. This test examines the void-generation flux-
modifier; as the blocks collide and bounce in vacuum, the interface between
closes and opens cleanly. Figure 12: The bouncing cylinders test, depicting
the diagonal velocity. The times shown are in $\mu$s. This test examines the
void-generation flux-modifier; as the cylinders collide and bounce in vacuum,
the interface between closes and opens cleanly. This test additionally
demonstrates the ability of the method to handle non-planar interface
configurations.

(a) Bouncing blocks.

(b) Bouncing cylinders.

Figure 13: The conservation error for both the bouncing blocks and bouncing cylinders tests, showing mass, momentum and energy loss over the course of the simulation. The conservation error for these tests remains low even despite the low resolution, with the bouncing blocks case showing larger error due to the larger interface area, and correspondingly increased use of the flux-modifier. Material |  $c_{1}$ GPa |  $c_{2}$ GPa | $c_{3}$ | $n$ | $m$ |  $T_{\text{melt}}$ K
---|---|---|---|---|---|---
Copper | 0.09 | 0.292 | 0.025 | 0.31 | 1.09 | 1356.0
CuBe | 1.041 | 0.0 | 0.025 | 0.31 | 1.09 | 1286.0
Table 5: Johnson Cook plasticity material parameters. Material | $D_{\text{crit}}$ | $\epsilon_{p,\text{thresh}}$ | $\epsilon_{p,\text{crit}}$ | $\alpha$ | $\nu$
---|---|---|---|---|---
Copper | 0.85 | 0.34 | 1.04 | 0.631 | 0.355
CuBe | 0.85 | 0.08 | 0.16 | 0.631 | 0.324
AerMet | 0.2 | 0.129 | 0.64 | 0.22 | 0.302
Table 6: Damage model parameters.

###  4.6 Tribological Pair Test

For more complex two dimensional shear validation, the tribological pair test
from Juanicotena [31] was considered. The tribological pair test is commonly
used to assess friction between different solid materials and allows for a
rigorous examination of the ability of the model to support multi-dimensional,
oblique sliding. As outlined by Barton and Drikakis [8], this test
demonstrates a marked difference between different material boundary
conditions. The test consists of a flyerplate colliding with a cylinder
comprised of two different materials, as shown in Figure 14(a). Owing to the
different material parameters of the inner and outer parts of the cone, the
shock wave produced by the collision of the flyerplate travels down the cone
at different speeds in the different materials. This produces a tangential
velocity discontinuity at the angled interface in the case of the slip
boundary conditions. When the wave reaches the bottom of the cone, the slip
conditions allow for the centre of the cone to protrude from the bottom, where
the initial shock wave arrives first, whereas the stick conditions do not.
These differences are captured in the free-surface velocity profile of the
test in Figure 14(b).

In this test, the work at hand is compared to the Ghost Fluid method of Barton
and Drikakis [8]. To match the test performed by Barton and Drikakis [8], the
materials are assumed to behave completely elastically. The initial conditions
for the tribological pair test are shown in Figure 14(a). The test is modelled
in 2D using cylindrical symmetry with a reflective boundary condition on the
$z$-axis. All materials are modelled using the Romenskii equation of state,
with parameters outlined in Table 1. The solid materials are surrounded by
void. The copper flyerplate has a downward velocity of 202 ms-1. The
computational domain spans $r=[0:3.2]\text{ cm},\ z=[0:4.28]\text{ cm}$. The
test is run at a resolution of 80 $\times$ 100 cells with one level of AMR
with a refinement factor 2 in order to match the simulation of Barton and
Drikakis [8]. The test is run at a CFL of 0.6, for 7.5 $\mu$s.

The results of the test are shown in figure 15, where both stick and slip
conditions are shown for comparison. The results are qualitatively examined
using the free surface velocity profile, taken on the axis at the bottom of
the cylinder. Figure 14(b) shows this velocity profile, and it can be seen
that the results of this method match those of Barton and Drikakis [8] very
well. This demonstrates the method’s ability to handle multi-material, multi-
dimension oblique shear well, matching an advanced elastoplastic Ghost Fluid
method. This test also has the added benefit of validating the material void
interaction and interface separation. The collision of the strong shock with
the void interface at the end of the cone is handled stably and the material
is allowed to protrude in the case of slip boundary conditions. Moreover, the
test is repeated with an additional layer of AMR, and run for a time of 13
$\mu$s to show how the different materials separate from each other at late
times, as shown in Figure 16.

(a) Initial conditions for the tribological pair test.

(b) A comparison of the free surface velocity profile in the tribological pair
test between this method and Barton and Drikakis [8] for both slip and stick
boundary conditions. Excellent agreement is observed.

Figure 14: Details of the tribological pair test. Figure 15: Comparison of the
stick and slip boundary conditions for the tribological pair test. In each
frame, the left hand side depicts the y-velocity and the right hand side shows
the pressure contours. The left column is with stick boundary conditions and
the right column is with slip boundary conditions. Times shown are in $\mu$s.
The method works well in both cases, with the slip boundary condition
correctly predicting the protruding cone and separating interfaces - features
that cannot be demonstrated without the suitable boundary conditions. Figure
16: The tribological pair test with slip boundary conditions at a higher
resolution at a time of $13\ \mu$s, showing the interface separation between
materials. This test was run without plasticity, meaning the interface
separation is reduced compared to experiment. Nonetheless, it validates the
correct working of the void-generation flux-modifier.

###  4.7 Fracture Tests with Experimental Comparison

####  4.7.1 One Dimensional Spallation

The fracture components of the method are validated in one dimension using an
experimental spallation test from Millett [42]. In this test, a block of
aluminium alloy 5083 H32 collides with a block of copper alloy CuBe TF00,
causing the copper to undergo spallation fracture. The test is run for 5
$\mu$s on a domain $x=[0:15]$ mm using a CFL 0.9 with a varying number of
cells. The initial conditions for the test are shown in Figure 17. Both
materials have plasticity, with the copper using the Johnson Cook model with
parameters laid out in Table 5, and the aluminium using ideal plasticity with
a yield stress of 0.275 GPa. The copper is damageable, with parameters laid
out in Table 6.

Figure 17: The 1D spallation test initial conditions.

The results are compared to experiment using the back surface velocity profile
of the copper block, shown in Figure 18. It can be seen that this method
matches the experiment well. However, it should be noted that the test is very
sensitive to the material damage parameters chosen.

Figure 18: The free-surface velocity profile for the 1D spallation test,
comparing to experiment. This test demonstrates combination of damage,
fracture and void generation working well in one dimension. A range of
resolutions is shown here, demonstrating the method’s convergence towards the
experimental data.

####  4.7.2 Plate Impact Spallation

To simultaneously validate all the components of the method, a test containing
material shear, separation and fracture is required. Such a test is the plate
impact spallation test from Thomas et al. [55]. This test features a copper
flyer-plate impacting another block of copper, the target, with sufficient
force that a spallation crack appears in the target. Spall fracture is induced
by the interaction of two strong tensile waves inside the target. As has been
mentioned, spall damage is important for any fracture model, as it requires
void-generation completely inside one material, which cannot physically be
modelled by the artificial injection of a low density gas. The initial set-up
for this test is shown in Figure 19.

Figure 19: The initial conditions for the plate impact spallation test.

The flyer-plate and target are treated as separate materials, both obeying the
Romenskii equation of state for copper. Both materials are damageable, with
the damage parameters given in Table 6.

The computational domain spans $r=[0:20\text{mm}],\ z=[0:7\text{mm}]$, and the
test is run for 3 $\mu$s at a CFL of 0.6. Figure 20 shows this test, where a
resolution of 500 $\times$ 180 cells with 2 levels of AMR is employed, giving
an effective resolution of 10 $\mu$m. The results show the method can handle
this large spallation fracture, combined with the slide and separation of the
two materials.

Similar to the 1D spallation case, the experimental validation for this test
comes from the free surface velocity profile, this time of the scab that is
formed after spallation. To match experiment, the profile is taken at $r=0.5$
mm. Comparison is also made with the Ghost Fluid method of Brauer et al. [16],
who performed the same test with a variety of damage models, using a level set
based approach. Four different resolutions are considered: 80, 40, 20 and 10
(Figure 20) $\mu$m. These velocity profiles are shown in Figure 21. It can be
seen that this work matches the experimental oscillation frequency well, and
improves with resolution. Moreover, this work matches the qualitative shape of
the material separation and the damaged region found by Brauer et al. [16]
with their level set based approach.

For other similar spallation tests on different materials, see for example
Kanel et al. [33] and Czarnota et al. [20].

Figure 20: The 2D plate impact spallation test. The plots show the downward
velocity in $m\cdot s^{-1}$, with the times given in $\mu$s. The test
demonstrates the shear, void-generation and fracture algorithms all working
together to enable this tough test to be modelled.  Figure 21: The plate
impact spallation test free surface velocity profile. Four different
resolutions are compared to experiment and previous numerical results and are
found to match well, demonstrating good convergence to the experimental data.

####  4.7.3 Tube Expansion Test

Figure 22: The 3D tube expansion test. Times shown are in $\mu$s. The steel
tube is depicted in grey and the two polycarbonate blocks are depicted in red
and blue. The test examines the three-dimensional capabilities of the method.
The work at hand proves to be a practical method to facilitate this kind of
large, anisotropic fracture test. Figure 23: A slice of the 3D tube expansion
test. Each frame shows the (left) $z-$velocity, (right) damage evolution and
(both) the AMR mesh. Times shown are in $\mu$s. This slice demonstrates the
flux-modifiers working well to allow the blocks to slide and separate inside
the tube, as well as depicting the formation of damage within the steel.
Figure 24: Experimental comparison for the tube expansion test. The radial
displacement profiles of the tube wall are compared in 2D and 3D versions of
this test. The 3D profiles are closer to the experimental results, reflecting
the anisotropic nature of the fracture in this problem and demonstrating the
necessity of a three dimensional algorithm for fracture studies.

To validate the method in three dimensions, the tube expansion test from
Vogler et al. [58] was considered. The test provides validation of all the
components of the algorithm, utilising shear, void and fracture. The test
consists of a tube containing a target plug and a projectile. As the
projectile hits the target, the tube around them expands and fractures. Both
the projectile and target consist of a cylinder of polycarbonate material,
with the tube being made of AerMet steel. The tube has a length of 50.8 mm, an
inner radius of 6.35 mm and a wall thickness of 3 mm. The two cylinders are
identical, both having a radius of 6.35 mm and a length of 25.4 mm. The upper
cylinder is given an initial downward velocity of 1.92 kms-1. The material
parameters for the test are laid out in Tables 1 and 6. Both the polycarbonate
and steel are assumed to obey ideal plasticity, with yield stresses of 0.08
and 2.56 GPa respectively. Exploiting the symmetry of the problem, a quarter
of the tube is modelled, with reflective boundaries on the $x-$, $y-$ and
$z-$axes. The test is run on a domain of $x=[0:20]$ mm, $y=[0:20]$ mm,
$z=[0:52]$ mm, using a CFL 0.3 for a time of 30 $\mu$s. A resolution of
$32\times 32\times 80$ cells with 2 levels of AMR, each of refinement factor
2, is used.

The full test is shown in Figure 22, where each material is delineated by its
$\phi=0.5$ volume fraction surface, with the steel tube shown in grey, and the
polycarbonate blocks in red and blue. Figure 23 shows the $z-$velocity, damage
evolution and AMR mesh on a vertical slice of the domain. Taking a slice both
allows the damage accumulation to be observed, and demonstrates the slide and
void opening required for the polycarbonate blocks.

The experimental comparison for this test is the radial displacement profile
measured by Vogler et al. [58], who give results up to 20 $\mu$s at three
different positions along the length of the tube using both VISAR probes and
photography. This work compares these measured values with those obtained here
in Figure 24, and good agreement is observed. The 3D results are obtained by
averaging the radial displacement over 20 different positions around the tube.
Moreover, these results demonstrate the ability of the method to facilitate
three dimensional simulations. This is especially important as many tests
often exhibit significant differences when modelled using a fully three
dimensional approach, compared to when modelled in two dimensions with radial
source terms. This effect is shown in Figure 24, where it can be seen that the
3D profiles match the experimental results better than running the same test
using the 2D model with radial symmetry. The ease with which the method can be
extended to three dimensions enables this difference to be captured. At an
even more fundamental level, a 2D approach can never capture the strain- and
fracture-localisation that a 3D model can allow.

##  5 Conclusions

This work has outlined a new method to enable interface slide, void-opening,
and material fracture in reduced-equation diffuse interface models. The model
has been validated in one, two and three dimensions against experiment and
previous simulations. The new model facilitates strenuous multi-material,
multi-physics simulations in a pragmatic, scalable fashion. This is achieved
by a set of flux-modifiers and interface seeding routines that enable material
boundary conditions to be applied across diffuse interfaces. The model is
general, capable of handling an arbitrary number of materials, in three
dimensions, with a broad range of flux methods. The method is broadly
applicable to any diffuse interface conservative update scheme, and it is not
limited to the system of equations outlined here – the method could be
suitably extended to any multi-phase system. Future work will focus on the
inclusion of different boundary conditions such as static and dynamic rigid
bodies and frictional boundaries, and the inclusion of additional multi-
physics effects such as reactive fluids.

##  Appendix A Damage Model Derivation

From the Clausius-Duhem inequality, irreversible damage requires

$-Y_{(l)}\dot{D}_{(l)}\geq 0\ .$ (158)

Using the definitions for the internal energies, the elastic energy release
rate is defined as

$-Y_{(l)}=\mathscr{E}_{(l)}^{c,\text{u}}+\frac{1}{\rho_{(l)}}G_{(l)}^{\text{u}}\mathcal{J}^{2}\
.$ (159)

The equivalent stress can be expressed as

$\sigma_{eq}^{2}=\frac{3}{2}||\text{dev}\bm{\sigma}||^{2}=6G^{2}_{(l)}\mathcal{J}^{2}=\left(1-D_{(l)}\right)^{2}{\sigma}_{eq}^{\text{u,2}}\
,$ (160)

where ${\sigma}_{eq}^{\text{u}}$ is the effective (undamaged) equivalent
stress. Substituting this into the above:

$-Y_{(l)}=\frac{\sigma_{eq}^{2}}{\left(1-D_{(l)}\right)^{2}}\frac{1}{2E_{(l)}^{\text{u}}}\left(\frac{2E_{(l)}^{\text{u}}\left(1-D_{(l)}\right)^{2}{\mathscr{E}}_{(l)}^{c,\text{u}}}{\sigma_{eq}^{2}}+\frac{E_{(l)}^{\text{u}}}{3\rho_{(l)}G_{(l)}^{\text{u}}}\right)\
,$ (161)

where the Youngs modulus $E_{(l)}^{\text{u}}$ has been introduced without loss
of generality.

It is not immediately apparent that this matches the well-recognised
expressions for $Y$ from the original works of Bonora [14], Lemaitre and
Desmorat [36] and others, which assume the St. Venant-Kirchhoff model for the
strain-energy density function, so here proof is provided.

Assume the undamaged cold energy is taken to be the logarithmic equation of
state from [49], which closely resembles the St. Venant-Kirchhoff model:

${\mathscr{E}}_{(l)}^{c,\text{u}}=\frac{K_{(l)}^{\text{u}}}{2\rho^{(l)}_{0}}\left[\ln\left(\frac{\rho_{(l)}}{\rho_{0}^{(l)}}\right)\right]^{2},\qquad{p}_{(l)}^{c,\text{u}}=K_{(l)}^{\text{u}}\left(\frac{\rho_{(l)}}{\rho_{0}^{(l)}}\right)\ln\left(\frac{\rho_{(l)}}{\rho_{0}^{(l)}}\right)\
,$ (162)

so that

${\mathscr{E}}_{c}^{(l)}=\frac{1}{2K}\frac{\rho_{0}^{(l)}}{\rho_{(l)}^{2}}{{p}_{(l)}^{c,\text{u}}}^{2}\
.$ (163)

Finally, using the relationships between elastic moduli

$\frac{E_{(l)}^{\text{u}}}{K_{(l)}^{\text{u}}}=3\left(1-2\nu_{(l)}^{\text{u}}\right),\qquad\frac{E_{(l)}^{\text{u}}}{2G_{(l)}^{\text{u}}}=1+\nu_{(l)}^{\text{u}}\
,$ (164)

and substituting all of the above, the elastic damage energy can be written

$-Y_{(l)}=\frac{1}{\rho}\frac{\sigma_{eq}^{2}}{\left(1-D_{(l)}\right)^{2}}\frac{1}{2E_{(l)}^{\text{u}}}R_{t}\left(\frac{p}{\sigma_{eq}}\right),\qquad
R_{t}\left(\frac{p}{\sigma_{e}}\right)=\left(3\left(1-2\nu_{(l)}^{\text{u}}\right)\frac{\rho^{(l)}_{0}}{\rho_{(l)}}\left(\frac{p}{\sigma_{eq}}\right)^{2}+\frac{2}{3}\left(1+\nu_{(l)}^{\text{u}}\right)\right)$
(165)

as required.

The variable $R_{t}$ is commonly referred to as the stress-triaxiality
function. Note that $R_{t}=1$ for uniaxial loading where $p/\sigma_{e}=\pm
1/3$. This result proves that the formulation conforms to the strain
equivalence hypothesis [36], which states that the constitutive model for
strain in the damaged state is equivalent to that for undamaged material
except the stress is replaced with the effective stress, and implies the
linear degradation of shear and bulk moduli:

${G}_{(l)}=(1-D_{(l)}){G}_{(l)}^{\text{u}},\qquad{K}_{(l)}=(1-D_{(l)}){K}_{(l)}^{\text{u}}\
.$ (166)

##  Appendix B Approximations for the Hencky Strain

Throughout the course of the algorithm presented here, such as in Section
3.1.2, it is necessary to calculate the strain from the deformation, and
reconstruct the deformation from the strain. In this work the Hencky strain
measure is employed. Although this measure admits a lot of benefits (as
described by Bažant [11] and Rezaee-Hajidehi et al. [50]), its primary
downside is its numerical complexity. This derives from the requirement to
solve matrix exponential/matrix logarithm problems.

The author of Barton [7] employs the approximation of Bažant [11] to calculate
the matrix logarithm required to evaluate the strain:

$\mbox{dev}\left({\mathbf{H}^{e}}\right)=\ln\left(\overline{\mathbf{V}}^{e}\right)\approx\frac{1}{2}\mathbf{H}_{B}^{e^{(1)}}\left(\overline{\mathbf{V}}^{e}{\overline{\mathbf{V}}^{e}}^{\text{T}}\right)\
,$ (167)

where

$\mathbf{H}_{B}^{e^{(m)}}\left(\overline{\mathbf{V}}^{e}\right)=\frac{1}{2m}\left({\overline{\mathbf{V}}^{e}}^{m}-{\overline{\mathbf{V}}^{e}}^{-m}\right)\
.$ (168)

To calculate the plastic update Barton [7] employs a singular value
decomposition (SVD), thus avoiding the need to calculate the matrix
exponential. This gives the numerically exact solution, but is computationally
expensive.

This work also uses the approximation of Bažant [11] to calculate the strain,
but employs a $(1,1)$ Padé approximation to the matrix exponential to perform
the plastic update:

$\overline{\mathbf{V}}^{e}=\exp\left(\mbox{dev}\left({\mathbf{H}^{e}}\right)\right)\approx\left(\mathbf{I}-\frac{1}{2}\mbox{dev}\left({\mathbf{H}^{e}}\right)\right)^{-1}\left(\mathbf{I}+\frac{1}{2}\mbox{dev}\left({\mathbf{H}^{e}}\right)\right)\
.$ (169)

For simulations involving plasticity, this is found to be sufficient since the
norm $||\mbox{dev}\left({\mathbf{H}^{e}}\right)||$ is known not to exceed
sufficiently small values for the materials of interest in physically
realistic scenarios [43].

If it is required to perform simulations without the physical limitations of
plasticity, thus allowing the norm
$||\mbox{dev}\left({\mathbf{H}^{e}}\right)||$ to grow unchecked, the
combination of approximations outlined here may be insufficient, and will lead
to the accumulation of errors. In this case, either the SVD can be used in
conjunction with the approximation of Bažant [11], or if the computational
burden of a SVD is still to be avoided, then this work suggests the use of the
$(2,2)$ diagonal Padé approximation for both the calculation of the strain and
the reconstruction of the deformation. To reconstruct the deformation, this
work follows Moler and Loan [43] in taking:

$\displaystyle\overline{\mathbf{V}}^{e}\approx\left(\mathbf{I}-\frac{1}{2}\mbox{dev}\left({\mathbf{H}^{e}}\right)+\frac{1}{12}\mbox{dev}\left({\mathbf{H}^{e}}\right)^{2}\right)^{-1}\left(\mathbf{I}+\frac{1}{2}\mbox{dev}\left({\mathbf{H}^{e}}\right)+\frac{1}{12}\mbox{dev}\left({\mathbf{H}^{e}}\right)^{2}\right)\
.$ (170)

To calculate the strain, this work follows Rezaee-Hajidehi et al. [50] in
taking:

$\displaystyle\mbox{dev}\left({\mathbf{H}^{e}}\right)\approx\frac{3}{2}\left(\mathbf{B}^{2}-\mathbf{I}\right)\left(\mathbf{B}^{2}+4\mathbf{B}+\mathbf{I}\right)^{-1}\
,$ (171)

where $\mathbf{B}=\overline{\mathbf{V}}^{e}{\overline{\mathbf{V}}^{e}}^{T}$ is
the left Cauchy-Green tensor. These higher order approximations are suitable
for all of the tests presented here, even those which do not employ
plasticity.

Bažant [11] also presents another alternative method for calculating higher
order approximations to the matrix logarithm, by taking linear combinations of
different approximations:

$\displaystyle\mbox{dev}\left({\mathbf{H}^{e}}\right)$
$\displaystyle\approx\frac{1}{2}\mathbf{H}_{B}^{e^{(1,2)}}\left(\overline{\mathbf{V}}^{e}{\overline{\mathbf{V}}^{e}}^{\text{T}}\right)=\frac{1}{2}\left(1.307\mathbf{H}_{B}^{e^{(1)}}-0.307\mathbf{H}_{B}^{e^{(2)}}\right)$
(172) $\displaystyle\mbox{dev}\left({\mathbf{H}^{e}}\right)$
$\displaystyle\approx\frac{1}{2}\mathbf{H}_{B}^{e^{(1/2,1)}}\left(\overline{\mathbf{V}}^{e}{\overline{\mathbf{V}}^{e}}^{\text{T}}\right)=\frac{1}{2}\left(1.326\mathbf{H}_{B}^{e^{(1/2)}}-0.326\mathbf{H}_{B}^{e^{(1)}}\right)\
.$ (173)

The latter of these was found to give results similar to the $(2,2)$ Padé
approximation, but the Padé approximations are preferred so as to avoid the
use of a matrix square root.

##  Appendix C Exact Solutions for Elastic Initial Value Problems

Exact solutions to the initial value problems in Sections 4.1 and 4.2 are
found using the method proposed in Barton et al. [9] and Barton and Drikakis
[8].

From the initial left and right states, an initial estimate of the states
between the waves (a.k.a star-states) is determined using numerical solutions.
The initial solution is used to determine the wave types (distinguishing
shocks from rarefactions by checking for entropy jumps and converging or
diverging characteristic wave speeds), and then the six wave speeds $S_{j}$
for $1\leq j\leq 6$ (from the Rankine-Hugoniot conditions for shocks, or the
eigenvalue evaluated at the tail of the wave for a rarefaction). New star-
states either side of the contact can then be determined exactly by solving
across the three waves on the left, from left to right, then the three waves
on the right, from right to left. In the case of shocks, the solution is found
from the Rankine-Hugoniot conditions:

$\mathbf{F}(\mathbf{U}^{+})-\mathbf{F}(\mathbf{U}^{-})-{S}\left(\mathbf{U}^{+}-\mathbf{U}^{-}\right)=\mathbf{0}.$
(174)

where $\pm$ indicate upstream/downstream states. For rarefaction waves the
solution is found by evaluating

$\frac{\partial\mathbf{W}}{\partial\xi}=\frac{\mathbf{r}_{j}(\mathbf{W})}{\mathbf{r}_{j}(\mathbf{W})\cdot\nabla_{\mathbf{W}}\lambda_{j}(\mathbf{W})},$
(175)

where $\lambda_{j}$ denotes the $j$-th eigenvalue, the speed $\xi$ is defined
within the range
$\lambda_{j}(\mathbf{W}^{-})\leq\xi=x/t\leq\lambda_{j}(\mathbf{W}^{+})$,
$\nabla_{\mathbf{W}}$ denotes the gradient operator with respect to components
of the vector of primitive variables, $\mathbf{W}$, and $\mathbf{r}_{j}$
denotes the $j$-th eigenvector. Full details of how to solve the wave
equations can be found in [9]

From the star-states next to the contact wave, the residual errors,
$\mathscr{R}$, in the required boundary conditions across the contact wave
give a measure of the error in the wave speeds (since the number of boundary
conditions across the contact equals the number of waves in the Riemann fan).
The exact solution method follows an iterative procedure which, given the
initial left and right states, seeks the exact wave speeds so as to minimise
$\mathscr{R}$. In the case of a solid/solid problem with stick conditions the
required conditions are equal velocity, normal stress and traction. In the
case of slip it is equal normal velocity and stress, and zero traction. In the
case of the solid/vacuum problem there are only three waves, and the boundary
conditions are zero normal stress and traction. The exact solution is
determined by solving the non-linear system
$\mathscr{R}(S_{1},S_{2},\ldots,S_{6})=\mathbf{0}$ for the wave speeds using
the Newton-Raphson method. For each new guess, the above procedure is repeated
to determine a new value for the residuals until they fall below the
threshold, which was here set to $10^{-9}$.

## Acknowledgements

This work was funded by AWE PLC. Tim Wallis acknowledges additional support by
a grant from the UK Engineering and Physical Sciences Research Council (EPSRC)
EP/L015552/1 for the Centre for Doctoral Training (CDT) in Computational
Methods for Materials Science. The authors would like to acknowledge the
valuable help of Dr Philip Blakely throughout the development of the method.

## References

## References

  * Allaire et al. [2000] Grégoire Allaire, Sébastien Clerc, and Samuel Kokh.  A five-equation model for the numerical simulation of interfaces in two-phase flows.  _Comptes Rendus de l’Academie des Sciences Series I Mathematics_ , 331(12):1017–1022, 2000.  ISSN 0764-4442.
  * Baer and Nunziato [1986] M.R. Baer and J.W. Nunziato.  A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials.  _International Journal of Multiphase Flow_ , 12(6):861–889, 1986.  ISSN 0301-9322.
  * Barlow et al. [2018] Andrew Barlow, Matej Klima, and Mikhail Shashkov.  Constrained optimization framework for interface-aware sub-scale dynamics models for voids closure in lagrangian hydrodynamics.  _Journal of computational physics_ , 371:914–944, 2018\.  ISSN 0021-9991.
  * Barlow et al. [2016] Andrew J Barlow, Pierre-Henri Maire, William J Rider, Robert N Rieben, and Mikhail J Shashkov.  Arbitrary lagrangian–eulerian methods for modeling high-speed compressible multimaterial flows.  _Journal of computational physics_ , 322(C):603–665, 2016.  ISSN 0021-9991.
  * Barton [2016] Philip T. Barton.  An Eulerian method for finite deformation anisotropic damage with application to high strain-rate problems.  _International Journal of Plasticity_ , 83:225–251, 2016\.  ISSN 0749-6419.
  * Barton [2018] Philip T Barton.  A level-set based Eulerian method for simulating problems involving high strain-rate fracture and fragmentation.  _International Journal of Impact Engineering_ , 117:75–84, 2018.  ISSN 0734-743X.
  * Barton [2019] Philip T. Barton.  An interface-capturing Godunov method for the simulation of compressible solid-fluid problems.  _Journal of Computational Physics_ , 2019.  ISSN 0021-9991.
  * Barton and Drikakis [2010] Philip T Barton and Dimitris Drikakis.  An Eulerian method for multi-component problems in non-linear elasticity with sliding interfaces.  _Journal of Computational Physics_ , 229(15):5518–5540, 2010.  ISSN 0021-9991.
  * Barton et al. [2009] P.T. Barton, D. Drikakis, E. Romenski, and V.A. Titarev.  Exact and approximate solutions of Riemann problems in non-linear elasticity.  _Journal of Computational Physics_ , 228(18):7046–7068, 2009.  ISSN 0021-9991.
  * Barton et al. [2011] P.T. Barton, B. Obadia, and D. Drikakis.  A conservative level-set based method for compressible solid/fluid problems on fixed grids.  _Journal of Computational Physics_ , 230(21):7867 – 7890, 2011.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2011.07.008.  URL http://www.sciencedirect.com/science/article/pii/S002199911100413X.
  * Bažant [1998] Zdeněk P. Bažant.  Easy-to-compute tensors with symmetric inverse approximating Hencky finite strain and its rate.  _Journal of Engineering Materials and Technology_ , 120(2):131–136, 04 1998.  ISSN 0094-4289.  doi: 10.1115/1.2807001.  URL https://doi.org/10.1115/1.2807001.
  * Benson [1992] David J. Benson.  Computational methods in Lagrangian and Eulerian hydrocodes.  _Comput. Methods Appl. Mech. Eng._ , 99(2–3):235–394, September 1992.  ISSN 0045-7825.  doi: 10.1016/0045-7825(92)90042-I.  URL https://doi.org/10.1016/0045-7825(92)90042-I.
  * Berger and Colella [1989] M.J Berger and P Colella.  Local adaptive mesh refinement for shock hydrodynamics.  _Journal of Computational Physics_ , 82(1):64–84, 1989.  ISSN 0021-9991.
  * Bonora [1997] N. Bonora.  A nonlinear CDM model for ductile failure.  _Engineering Fracture Mechanics_ , 58(1):11 – 28, 1997.  ISSN 0013-7944.  doi: https://doi.org/10.1016/S0013-7944(97)00074-X.  URL http://www.sciencedirect.com/science/article/pii/S001379449700074X.
  * Bonora et al. [2009] N. Bonora, A. Ruggiero, L. Esposito, and G. Iannitti.  Damage development in high purity copper under varying dynamic conditions and microstructural states using continuum damage mechanics.  _AIP Conference Proceedings_ , 1195(1):107–110, 2009.  doi: 10.1063/1.3294988.  URL https://aip.scitation.org/doi/abs/10.1063/1.3294988.
  * Brauer et al. [2018] A. Brauer, N. K. Rai, M. E. Nixon, and H. S. Udaykumar.  Modeling impact‐induced damage and debonding using level sets in a sharp interface Eulerian framework.  _International Journal for Numerical Methods in Engineering_ , 115(9):1108–1137, 2018.  ISSN 0029-5981.
  * Bruhns et al. [2001] O. T. Bruhns, H. Xiao, and A. Meyers.  Constitutive inequalities for an isotropic elastic strain-energy function based on Hencky’s logarithmic strain tensor.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 457(2013):2207–2226, 2001.  ISSN 1364-5021.
  * Chen et al. [2010] Qianyi Chen, Jingtao Wang, and Kaixin Liu.  Improved CE/SE scheme with particle level set method for numerical simulation of spall fracture due to high-velocity impact.  _Journal of Computational Physics_ , 229(19):7503–7519, 2010.  ISSN 0021-9991.
  * Courant et al. [1928] R. Courant, K. Friedrichs, and H. Lewy.  Über die partiellen differenzengleichungen der mathematischen physik.  _Mathematische Annalen_ , 100(1):32–74, 1928\.  ISSN 0025-5831.
  * Czarnota et al. [2008] C Czarnota, N Jacques, S Mercier, and A Molinari.  Modelling of dynamic ductile fracture and application to the simulation of plate impact tests on tantalum.  _Journal of the mechanics and physics of solids_ , 56(4):1624–1650, 2008.  ISSN 0022-5096.
  * Deiterding et al. [2006] Ralf Deiterding, Raul Radovitzky, Sean Mauch, Ludovic Noels, Julian C. Cummings, and Daniel I. Meiron.  A virtual test facility for the efficient simulation of solid material response under strong shock and detonation wave loading.  _Engineering with Computers_ , 22:325–347, 2006.
  * Deng et al. [2018] Xi Deng, Satoshi Inaba, Bin Xie, Keh-Ming Shyue, and Feng Xiao.  High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces.  _Journal of Computational Physics_ , 371:945–966, 2018\.  ISSN 0021-9991.
  * Dorovskii et al. [1983] V. Dorovskii, A. Iskol’dskii, and E. Romenskii.  Dynamics of impulsive metal heating by a current and electrical explosion of conductors.  _Journal of Applied Mechanics and Technical Physics_ , 24(4):454–467, 1983.  ISSN 0021-8944.
  * Favrie and Gavrilyuk [2012] N. Favrie and S.L. Gavrilyuk.  Diffuse interface model for compressible fluid – compressible elastic–plastic solid interaction.  _Journal of Computational Physics_ , 231(7):2695–2723, 2012.  ISSN 0021-9991.
  * Favrie et al. [2009] N. Favrie, S.L. Gavrilyuk, and R. Saurel.  Solid–fluid diffuse interface model in cases of extreme deformations.  _Journal of Computational Physics_ , 228(16):6037–6077, 2009.  ISSN 0021-9991.
  * Fedkiw et al. [1999] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher.  A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method).  _Journal of Computational Physics_ , 152(2):457–492, 1999.  ISSN 0021-9991.
  * Gokhale et al. [2018] Nandan Gokhale, Nikos Nikiforakis, and Rupert Klein.  A dimensionally split cartesian cut cell method for hyperbolic conservation laws.  _Journal of Computational Physics_ , 364:186–208, 2018\.  ISSN 0021-9991.
  * Hank et al. [2017] Sarah Hank, Sergey Gavrilyuk, Nicolas Favrie, and Jacques Massoni.  Impact simulation by an Eulerian model for interaction of multiple elastic-plastic solids and fluids.  _International Journal of Impact Engineering_ , 109(C):104–111, 2017.  ISSN 0734-743X.
  * Johnsen and Colonius [2006] Eric Johnsen and Tim Colonius.  Implementation of WENO schemes in compressible multicomponent flow problems.  _Journal of Computational Physics_ , 219(2):715–732, 2006.  ISSN 0021-9991.
  * Johnson and Cook [1985] Gordon R. Johnson and William H. Cook.  Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures.  _Engineering Fracture Mechanics_ , 21(1):31–48, 1985.  ISSN 0013-7944.
  * Juanicotena [2006] A. Juanicotena.  Experimental investigation of dynamic friction at high contact pressure applied to an aluminium/stainless steel tribo pair.  _Journal de Physique IV_ , 134(1):559–564, August 2006.  doi: 10.1051/jp4:2006134086.
  * Kachanov [1958] L. Kachanov.  On the creep fracture time.  _Izv Akad, Nauk USSR Otd Tech._ , 8:26–31, 1958.
  * Kanel et al. [1997] G.I Kanel, S.V Razorenov, A Bogatch, A.V Utkin, and Dennis E Grady.  Simulation of spall fracture of aluminum and magnesium over a wide range of load duration and temperature.  _International journal of impact engineering_ , 20(6):467–478, 1997.  ISSN 0734-743X.
  * Kapila et al. [2001] A. K. Kapila, R. Menikoff, J. B. Bdzil, S. F. Son, and D. S. Stewart.  Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations.  _Physics of Fluids_ , 13(10):3002–3024, 2001\.  ISSN 1070-6631.
  * Klima et al. [2020] Matej Klima, Andrew Barlow, Milan Kucharik, and Mikhail Shashkov.  An interface-aware sub-scale dynamics multi-material cell model for solids with void closure and opening at all speeds.  _Computers & fluids_, 208(C):104578, 2020.  ISSN 0045-7930.
  * Lemaitre and Desmorat [2005] J. Lemaitre and R. Desmorat.  _Engineering damage mechanics_.  Springer, 2005.