# A unified diffuse interface method for the interaction of rigid bodies with
elastoplastic solids and multi-phase mixtures

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 multi-physics-compatible immersed rigid body method
for Eulerian finite-volume simulations. To achieve this, rigid bodies are
represented as a diffuse scalar field and an interface seeding method is
employed to mediate the interface boundary conditions. The method is based on
an existing multi-material diffuse interface method that is capable of
handling an arbitrary mix of fluids and elastoplastic solids. The underlying
method is general and can be extended to a range of different applications
including high-strain rate deformation in elastoplastic solids and reactive
fluid mixtures. As such, the new method presented here is thoroughly tested
through a variety of problems, including fluid-rigid body interaction,
elastoplastic-rigid body interaction, and detonation-structure interaction.
Comparison is drawn between both experimental work and previous numerical
results, with excellent agreement in both cases. The new method is
straightforward to implement, highly local, and parallelisable. This allows
the method to be employed in three dimensions with multiple levels of adaptive
mesh refinement using complex immersed geometries. The rigid body field can be
static or dynamic, with the THINC interface reconstruction method being used
to keep the interface sharp in the dynamic case.

© British Crown Owned Copyright 2021/AWE

###### keywords:

Rigid bodies , Diffuse interface , Multi-physics , Elastoplastic solids ,
Reactive fluids

##  1 Introduction

Immersed rigid boundaries feature in the simulation of a wide variety of
different physical applications. These range from fluid-solid interaction
problems to elastoplastic machining studies. The inclusion of a rigid body in
a simulation is most often to save computational expense; discretising inside
a solid body and solving a fully elastoplastic system of equations is
expensive for a number of reasons. In certain applications, the evolution
inside a solid body may not be relevant, and in this case it saves time and
expense to replace the solid with a rigid boundary.

This paper introduces a new diffuse interface method for the inclusion of
static and dynamic rigid bodies in Cartesian Eulerian finite volume
simulations. A variety of other methods exist to model rigid boundaries,
designed for a range of different problems including inviscid hypersonic
flows, viscous Navier-Stokes problems and turbulence modelling. Many fluid-
structure interaction studies consider body-fitted or unstructured meshes.
These methods, however, often feature some difficulty and expense in creating
suitable meshes around complex geometries in three dimensions. Cartesian
approaches also exist, offering a more straightforward approach to mesh
generation and the data structures that describe them. Such approaches include
level-set based Ghost Fluid methods such as Sambasivan and UdayKumar [26] or
cut-cell approaches such as Gokhale et al. [15]. However these methods have
their own associated problems. Recently Kemm et al. [22] also include rigid
bodies with the use of a diffuse interface model, incorporating the effect of
the rigid body as a forcing term in the momentum equation, using the ADER-WENO
approach.

This work chooses a diffuse interface representation for rigid body
interfaces. The methods presented build on the work of Wallis et al. [33],
where novel flux-modifiers and interface seeding routines were employed to
mediate void-opening, slip and fracture internal boundary conditions in a
diffuse interface context, without the use of level sets. The work at hand
follows a very similar approach by adding an additional scalar field $r$, the
rigid body volume fraction, into the system of equations to track rigid
bodies. From here, an analogous interface seeding method is derived to enable
rigid boundaries to be incorporated in the system. The rigid body volume
fraction is treated as any other volume fraction, allowing the THINC
reconstruction methods of Xiao et al. [36] to be used to keep interfaces
sharp.

The work of Wallis et al. [33] used a reduced-equation Allaire-type [1] multi-
material diffuse interface scheme, adapted by Barton [2] to include
elastoplastic solids. On top of this underlying multi-material method,
additional history parameters such as plastic strain, and scalar fields such
as the void volume fraction, were then used to include additional physics in a
straightforward manner.

As outlined by Wallis et al. [33], the major benefit of the new scheme is its
straightforward, practical approach. Very little additional numerical
machinery must be added to incorporate rigid bodies, and any desired physics.
For the hyperbolic update, the same set of partial differential equations is
solved across the entire domain, incorporating the effect of all materials in
the mixture. The numerical methods required to mediate the rigid boundary
conditions are highly local to the rigid body interface, avoiding the parallel
communication that can be required for other methods. This simplicity allows
the method to parallelise efficiently over several layers of adaptive mesh
refinement (AMR) and in three dimensions. Additionally, as rigid bodies are
represented as a continuous scalar field in a standard Cartesian mesh, no
additional effort must be made in creating an appropriate mesh, or re-meshing,
and the use of body fitted meshes is avoided. The new methods developed in
this work are entirely compatible with the previous works on which it is
based, as the same underlying multi-material system of equations employed by
Wallis et al. [33] is used here. This allows the new rigid body method to be
applied to a wide variety of applications, including single material fluid
flow, elastoplastic-rigid interaction, and reactive fluid mixtures (as
presented by Wallis et al. [34]). This work will therefore present a number of
validation cases covering this range, with both experimental and numerical
comparison.

##  2 Governing Theory

The system considered in this work will be laid out as follows. Firstly, the
full multi-material base system of equations and its thermodynamics will be
briefly outlined. The base system has been laid out in detail in previous
works [2, 33, 34] and is not affected by the later additions of a rigid body,
so shall not be repeated here in depth. Subsequently, the additional
components required for the rigid body will be introduced. From here, a
smaller subsystem will be derived, for use in single material fluid structure
interaction problems, which will be the primary focus of this work.

###  2.1 Base System

####  2.1.1 Evolution Equations

The base system, devised by Barton [2], is an Allaire-type [1] multi-material
diffuse interface system, capable of accounting for an arbitrary mixture of
both fluids and elastoplastic solids. The state of any material $l$ is
characterised by its 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 simplifying assumption is made that materials in a mixture
are in mechanical and thermal equilibrium, resulting in a reduced equation
system where materials share the same momentum, energy and deformation
equations. For $l=1,\ldots,N$ materials:

$\displaystyle\frac{\partial\phi_{(l)}}{\partial
t}+\frac{\partial\phi_{(l)}u_{k}}{\partial x_{k}}$
$\displaystyle=\phi_{(l)}\frac{\partial u_{k}}{\partial x_{k}}$ (1)
$\displaystyle\frac{\partial\phi_{(l)}\rho_{(l)}}{\partial
t}+\frac{\partial\phi_{(l)}\rho_{(l)}u_{k}}{\partial x_{k}}$ $\displaystyle=0$
(2) $\displaystyle\frac{\partial\rho u_{i}}{\partial t}+\frac{\partial(\rho
u_{i}u_{k}-\sigma_{ik})}{\partial x_{k}}$ $\displaystyle=0$ (3)
$\displaystyle\frac{\partial\rho E}{\partial t}+\frac{\partial(\rho
Eu_{k}-u_{i}\sigma_{ik})}{\partial x_{k}}$ $\displaystyle=0$ (4)
$\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=\frac{2}{3}\overline{\mathbf{V}}^{e}_{ij}\frac{\partial
u_{k}}{\partial x_{k}}-u_{i}\beta_{j}-\boldsymbol{\Phi}_{ij}\ ,$ (5)

where $E=\mathscr{E}+|\mathbf{u}|/2$ denotes the specific total energy,
$\boldsymbol{\sigma}$ denotes the Cauchy stress tensor,
$\beta_{j}=\partial\overline{\mathbf{V}}^{e}_{kj}/\partial x_{k}$, and
$\boldsymbol{\Phi}$ represents the contribution from plastic effects (details
of which can be found in Barton [2]).

Some multi-physics closure models introduce a dependence on material history
variables. For example, the equivalent plastic strain $\varepsilon_{p(l)}$ or
the damage parameter $D_{(l)}$ in the case of fracture [33], or the reaction
progress variable $\lambda$ for reactive fluid mixtures [34]. For these
variables, 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)}\ .$ (6)

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.

####  2.1.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)\
,$ (7)

where

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

is the deviatoric 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)\
,$ (9)

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)}}\
.$ (10)

The specific form of the Grüneisen function for each material must be
provided, but it is generally taken to be constant. The elastic 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)\
,$ (11)

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)$
(12)

is the second invariant of shear strain.

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

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

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, 2]:

$\displaystyle 1$ $\displaystyle=$ $\displaystyle\sum_{l=1}^{N}\phi_{(l)}$
(16) $\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\phi_{(l)}\rho_{(l)}$ (17)
$\displaystyle\rho\mathscr{E}$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N}\phi_{(l)}\rho_{(l)}\mathscr{E}_{(l)}$ (18)
$\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)}$
(19) $\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)}\
,$ (20)

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

####  2.1.3 Closure Models

It can be seen that, by presenting the internal energy in the form outlined,
equation (14) 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)\
.$ (21)

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, and it is
therefore the choice of the reference curves
$p_{\text{\scriptsize{ref}}},\mathscr{E}_{\text{\scriptsize{ref}}}$ and
$\Gamma$ which will convey specific material properties. This form is not
limited to solely fluids or solids, with the choice of reference curves
capable of incorporating cold-compression, thermal, and shear effects. A
variety of different materials will be considered in the validation tests,
which are detailed here.

For the case of fluids like air, a simple ideal gas law will be used. In this
case, $p_{\text{\scriptsize{ref}}}=0,\mathscr{E}_{\text{\scriptsize{ref}}}=0$,
and the material is entirely specified by the choice of the adiabatic index
$\gamma=\Gamma+1$ which is taken to be constant. This reduces the form of
Equation 21 to:

$\displaystyle p=(\gamma-1)\rho\mathscr{E}\ .$ (22)

When elastoplastic solids are employed, this work follows previous studies by
using the equation of state from Dorovskii et al. [13]. For this equation of
state, $\Gamma(\rho)=\Gamma_{0}$, where $\Gamma_{0}$ is a material dependent
constant, and:

$\displaystyle\mathscr{E}_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle=\mathscr{E}_{(l)}^{c}+\mathscr{E}^{s}_{(l)}=\frac{K_{0}}{2\rho_{(l)}\alpha^{2}}\left(\left(\frac{\rho_{(l)}}{\rho_{0}}\right)^{\alpha}-1\right)+\frac{G_{(l)}\left(\rho_{(l)},\alpha_{(l)}\right)}{\rho_{(l)}}\mathcal{J}^{2}\left(\mbox{dev}\left({\mathbf{H}^{e}}\right)\right)$
(23) $\displaystyle G(\rho_{(l)})$
$\displaystyle=G_{0}\left(\frac{\rho_{(l)}}{\rho_{0}}\right)^{\beta+1}$ (24)
$\displaystyle p_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle=\rho^{2}_{(l)}\frac{\partial\mathscr{E}_{\text{\scriptsize{ref}},(l)}}{\partial\rho_{(l)}}\
.$ (25)

Finally, when reactive fluids are considered, the Jones-Wilkins-Lee (JWL)
equation of state is employed for both the reactants and products of the
explosive. For this equation of state $\Gamma(\rho)=\Gamma_{0}$ and

$\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)}}}$ (26)
$\displaystyle p_{\text{\scriptsize{ref}},(l)}(\rho_{(l)})$
$\displaystyle=\rho^{2}_{(l)}\frac{\partial\mathscr{E}_{\text{\scriptsize{ref}},(l)}}{\partial\rho_{(l)}}={\cal
A}e^{-{\cal R}_{1}\frac{\rho_{0}}{\rho_{(l)}}}+{\cal B}e^{-{\cal
R}_{2}\frac{\rho_{0}}{\rho_{(l)}}}$ (27)

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 18 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}|\ .$ (28)

Then the Mie-Grüneisen form of the equation of state is used to rewrite
equation 18, 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)\
.$ (29)

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)}}}\
.$ (30)

###  2.2 Rigid bodies

Rigid bodies are represented as scalar fields. This field, $r$, varies from 0
to 1, with the $r=0.5$ contour marking the rigid body interface. When the
rigid body is dynamic it evolves under its own velocity field
$\mathbf{v}_{r}$, which can be any required function of space and time, but
the velocity does not mutually interact with the flow at large; it is simply
imposed on the rigid body, giving the evolution equation:

$\displaystyle\frac{\partial r}{\partial t}+(\mathbf{v}_{r})_{k}\frac{\partial
r}{\partial x_{k}}$ $\displaystyle=0\ .$ (31)

###  2.3 Reduced System

Clearly not every one of the equations presented above will be relevant for a
given application. Frequently, single-material fluid tests are used to examine
the performance of a rigid body method through fluid-solid interaction tests.
In this case, a reduced system of equations can be employed, simply consisting
of the Euler equations, plus the rigid body evolution law:

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial\rho
u_{k}}{\partial x_{k}}$ $\displaystyle=0$ (32)
$\displaystyle\frac{\partial\rho u_{i}}{\partial t}+\frac{\partial(\rho
u_{i}u_{k}-\sigma_{ik})}{\partial x_{k}}$ $\displaystyle=0$ (33)
$\displaystyle\frac{\partial\rho E}{\partial t}+\frac{\partial(\rho
Eu_{k}-u_{i}\sigma_{ik})}{\partial x_{k}}$ $\displaystyle=0$ (34)
$\displaystyle\frac{\partial r}{\partial t}+(\mathbf{v}_{r})_{k}\frac{\partial
r}{\partial x_{k}}$ $\displaystyle=0\ .$ (35)

Here the stress tensor takes the simple form $\sigma_{ij}=-p\delta_{ij}$,
where $p$ is the hydrodynamic pressure and $\delta$ is the Dirac delta
function.

This system is intended to clarify the method in simple fluid structure
interaction problems. However, the full system of equations is still employed
when considering more complex situations. For example, the deformation tensor
and plastic history variables must also be included alongside this system in
elastoplastic-rigid interaction studies.

##  3 Numerical Approach

Apart from the new algorithms required for rigid boundaries, the numerical
method for the full multi-material system is the same as presented in Wallis
et al. [33]. Therefore, the bulk of the method shall not be repeated, but the
core algorithms are summarised here. For clarity, the methods are described in
the context of the reduced equation system.

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 [38], which includes an implementation of the
structured adaptive mesh refinement (SAMR) method of Berger and Colella [4]
for solving hyperbolic systems of partial differential equations. 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 reduced
system of equations in compact vector form by separating it into various
qualitatively different parts: a conservative hyperbolic part for each spatial
dimension and non-conservative terms for the rigid body volume fraction:

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

This is given by:

$\displaystyle\partialderivative{t}\matrixquantity(\rho\\\ \rho u_{i}\\\ \rho
E\\\ r\\\ )+\partialderivative{x_{k}}\matrixquantity(\rho u_{k}\\\ \rho
u_{i}u_{k}-\sigma_{ik}\\\ \rho Eu_{k}-u_{i}\sigma_{ik}\\\ 0\\\
)=\matrixquantity(0\\\ 0\\\ 0\\\ (\mathbf{v}_{r})_{k}\frac{\partial
r}{\partial x_{k}}\\\ )\ .$ (37)

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 [10]. The hyperbolic part is updated first, followed by updating
the rigid body field if required.

####  3.0.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,$
(38)

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)$
(39) $\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),$

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 (39), resulting in a spatially unsplit method.

Fluxes are computed using the HLLC solver developed by Toro et al. [32] (note
that the full multiphase system uses the HLLD solver from Barton [2]). To
achieve higher order spatial accuracy, the initial conditions for the Riemann
solver are taken to be the piecewise linear (MUSCL) reconstruction of the cell
centred primitive variables. Primitive reconstruction is employed here so as
to be consistent with the full multiphase system, but conservative
reconstruction is also possible for the reduced system. An interface
reconstruction method is also 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. Specifically,
this work uses the MUSCL-BVD-THINC scheme of Deng et al. [11]. This 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. This
procedure is applied to all variables, including the rigid volume fraction
field. The rigid volume fraction field is finally updated using a simple
upwind method, based on the sign of $\mathbf{v}_{r}$ and the reconstructed
cell-wall values. If greater accuracy is required for the advection of the
rigid body field, this upwind method can be replaced by any desired scheme,
such as volume of fluid methods. Thermodynamic variables do not need to be
updated in rigid ($r>0.5$) cells, saving computational expense.

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)$
(40) $\displaystyle\mathbf{U}^{(2)}$ $\displaystyle=$
$\displaystyle\mathbf{U}^{(1)}-\Delta
t\mathcal{D}\left({\mathbf{U}}^{(1)}\right)$ (41)
$\displaystyle\mathbf{U}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbf{U}^{n}+\frac{1}{4}\mathbf{U}^{(2)}$ (42)
$\displaystyle\mathbf{U}^{(4)}$ $\displaystyle=$
$\displaystyle\mathbf{U}^{(3)}-\Delta
t\mathcal{D}\left({\mathbf{U}}^{(3)}\right)$ (43)
$\displaystyle\mathbf{U}^{n+1}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\mathbf{U}^{n}+\frac{2}{3}\mathbf{U}^{(4)}\ .$ (44)

###  3.1 Rigid body method

Rigid bodies introduce undeforming, perfectly reflective boundaries within the
domain, also known as immersed boundaries, which are governed by the evolution
equation (31). The boundary condition at rigid body interfaces has two
components: no-penetration and either slip or no-slip. At the interface, the
no-penetration condition requires:

$(\mathbf{u}-\mathbf{v}_{r})\cdot\mathbf{\hat{n}}_{r}=0\ ,$ (45)

where $\mathbf{\hat{n}}_{r}$ is the rigid body normal vector. This is enforced
by reflecting normal-aligned vector- and tensor-components inside the rigid
body. For the slip boundary condition, scalar and tangential vector quantities
remain unchanged, and are copied into the rigid body from outside, resulting
in a zero-tangential stress condition. The no-slip condition is not employed
here, but is enforced by also reflecting tangential vector quantities inside
the body, resulting in the condition $\mathbf{u}=\mathbf{v}_{r}$ at the
interface.

The method presented here shall be applicable to both the reduced system and
the full multi-material system. Analogously to the method in Wallis et al.
[33], the rigid body boundary conditions are mediated by flux-modifiers and
interface seeding routines. However, by comparison, no new flux modifiers are
required: the extant void-opening flux modifier is applied to rigid-solid
interfaces to allow solid bodies to collide with and separate from rigid
bodies, just as they would with any other solid body in the domain. The flux
modifier is applied when the rigid volume fraction ‘changes sign’ between the
two cells being used to calculate the flux:

$(0.5-r_{L})(0.5-r_{R})<0\ .$ (46)

Although no additional flux-modifiers are required, rigid bodies must be
seeded with suitable states to enable the correct boundary interaction. This
is done using an interface seeding routine, similar to those presented in
Wallis et al. [33]. This procedure is carried out at the start of every time-
step, before the calculation of the hyperbolic update.

In the following, it is necessary to distinguish between variables that are
reflected across a rigid body interface and those that are not, denoted by
$\mathbf{U}^{\mbox{\scriptsize{ref}}}$ and
$\mathbf{U}^{\mbox{\scriptsize{const}}}$ respectively. All scalar quantities
are constant. For vector and tensor quantities, it is easiest to define these
quantities in the reference frame where $\mathbf{\hat{n}}_{r}$ points in the
$x$-direction $\mathbf{\hat{x}}$ (although analogous results can be obtained
for any direction). States that have been rotated into this reference frame
will be denoted with a tilde, $\mathbf{U}\rightarrow\tilde{\mathbf{U}}$, and
are found in the normal way, where for a vector $\mathbf{v}$ and matrix M:

$\displaystyle\mathbf{\tilde{v}}$ $\displaystyle=\mathbf{R}\mathbf{v}$ (47)
$\displaystyle\tilde{\mbox{M}}$
$\displaystyle=\mathbf{R}\mbox{M}\mathbf{R}^{\text{T}}\ .$ (48)

The rotation matrix $\mathbf{R}$ is detailed in Wallis et al. [33].

Having rotated a state into this new frame, the reflected and constant parts
of the velocity and stretch tensor are defined as:

The interface seeding routine proceeds as follows:

The rigid seeding routine: 1. In all rigid cells (with $r>0.5$), the outward
pointing interface normal, $\mathbf{\hat{n}}_{r}$, is calculated. 2. A probe
is sent out along the normal direction. The probe can either simply be taken
to be a constant length, such as 1.5$\differential x$, or sample points can
chosen at various distances stepping out in the normal direction (from
$\differential x$ to 4$\differential x$), and when the $2^{\text{DIM}}$
enclosing cell centres around the sample point are valid ($r<0.5$), this point
is chosen. Once the probe point is chosen, material values are interpolated at
that point to give a new state $\mathbf{U}_{\text{Interp}}$. If a valid probe
position cannot be found, no further action is taken. 3. The interpolated cell
is rotated so the normal lies in the $\mathbf{\hat{x}}$ direction, giving
$\tilde{\mathbf{U}}_{\text{Interp}}$. 4. A new state,
$\tilde{\mathbf{U}}^{\mbox{\scriptsize{new}}}$, is then created by reflecting
the suitable components of the interpolated cell:
$\tilde{\mathbf{U}}^{\mbox{\scriptsize{new}}}=\tilde{\mathbf{U}}_{\text{Interp}}^{\mbox{\scriptsize{const}}}-\tilde{\mathbf{U}}_{\text{Interp}}^{\mbox{\scriptsize{ref}}}\
.$ If the rigid body is dynamic, then the normal component of the rigid body
velocity must also be included. This is done adding
$2\mathbf{\tilde{v}}_{r}^{\mbox{\scriptsize{ref}}}=2(\mathbf{v}_{r}\cdot\mathbf{\hat{n}}_{r})$
to the reflected velocity of $\tilde{\mathbf{U}}^{\mbox{\scriptsize{new}}}$.
Therefore, in total:
$\mathbf{\tilde{u}}^{\mbox{\scriptsize{new}},\mbox{\scriptsize{ref}}}=2\mathbf{\tilde{v}}_{r}^{\mbox{\scriptsize{ref}}}-\mathbf{\tilde{u}}_{\text{Interp}}^{\mbox{\scriptsize{ref}}}\
.$ 5. Here, the method then diverges with two options: the original or Riemann
problem based version. (a) For the original method, this new state is then
rotated back into the original frame of reference and replaces the original
state. $\mathbf{U}\leftarrow\mathbf{U}^{\mbox{\scriptsize{new}}}\ .$ (b) For
the Riemann problem based method, the HLLC(D) star states are calculated using
the same method as for the flux calculation, with the left and right states
being $\tilde{\mathbf{U}}^{\mbox{\scriptsize{new}}}$ and
$\tilde{\mathbf{U}}_{\text{Interp}}$ respectively. This gives the intermediate
states $\tilde{\mathbf{U}}_{L}^{*}$ and $\tilde{\mathbf{U}}_{R}^{*}$. The
$\tilde{\mathbf{U}}_{R}^{*}$ state is then rotated back into the original
frame of reference and replaces the original state.
$\mathbf{U}\leftarrow\mathbf{U}_{R}^{*}\ .$

####  3.1.1 Normal Estimation

Two different methods are considered for the calculation of the rigid body
normal from the volume fraction field $r$. The first is from Youngs [37].
Youngs’ method is a first order approximation to the normal, attained using a
simple difference formula for the local $3^{\text{\scriptsize{DIM}}}$ cells.
The method is cheap, and is easily extended into three dimensions. Following
the notation of Youngs [37], the normal of a field $f$ is estimated as:

$\displaystyle\hat{n}=-\frac{\gradient f}{|\gradient f|}\ ,$ (49)

with

$\displaystyle\partialderivative{f}{x}$
$\displaystyle=\frac{f_{E}-f_{W}}{2\differential x}$ (50)
$\displaystyle\partialderivative{f}{y}$
$\displaystyle=\frac{f_{N}-f_{S}}{2\differential y}$ (51)
$\displaystyle\partialderivative{f}{z}$
$\displaystyle=\frac{f_{T}-f_{B}}{2\differential z}\ ,$ (52)

where symmetric difference formulae must be used for each direction. These
formulae are outlined by Youngs [37] and Pilliod and Puckett [24].

Although Youngs’ method is efficient, it is only first order, and it is not
guaranteed to reproduce linear interfaces exactly.

The second order ELVIRA method is also considered. This method is outlined by
Pilliod and Puckett [24], and produces a second order accurate interface
normal. The method attempts to find the best normal from a set of candidate
normals that are again obtained using simple difference formulae. The measure
of a candidate normal’s quality found by calculating the volume fractions the
candidate normal would produce in neighbouring cells, and comparing these
volume fractions to the true cell volume fractions. Whichever candidate normal
has the closest match is then selected. This results in a more accurate
calculation of the normal, with linear interfaces being reconstructed exactly
in cells where the interface passes through the centre cell [24]. Naturally,
the ELVIRA method is more expensive as it requires the calculation of multiple
candidate normals and calculation of geometric informations from the
interface. As such, ELVIRA is only employed in two dimensional validation
tests, with Youngs’ method being used in three dimensions.

##  4 Validation and Verification

Throughout all the tests, the rigid body is plotted in solid grey for all
rigid ($r>0.5$) cells.

###  4.1 Shock Reflection From a Wedge

Initially, two shock-wedge interaction tests are considered. These tests
feature a shock wave in an ideal gas colliding with an angled plane. The first
test features a Mach 1.7 shock colliding with a 25∘ wedge, producing a single
reflection. The second test features a Mach 10 shock colliding with a 30∘
wedge, producing a more complex double reflection. These tests examine the
ability of the model to give accurate solutions across different shock
strengths. The Mach 1.7 test can be compared to the similar test in Kemm et
al. [22], and the Mach 10 test can be compared to the similar tests in
Woodward and Colella [35] and Forrer and Jeltsch [14].

Both of these tests are non-dimensionalised, with a domain $x=[0:3]$,
$y=[0:2]$, a base resolution of 300 $\times$ 200, and two levels of AMR, both
refinement factor 2. The gas is a $\gamma=1.4$ ideal gas, with a quiescent
state of $\rho=1.4,p=1.0,\mathbf{u}=\mathbf{0}$. The rigid body wedge starts
at $x=0.5$ and the shock is initially located at $x=0.25$ in both cases. A CFL
of 0.6 is used, with the first test being run to a time of 1.2 and the second
test being run to a time of 0.2. These tests are shown in Figure 1. The method
handles the test well, replicating the results of previous numerical
simulations. In this case, the second order ELVIRA normal estimation method is
preferred over the first order Young’s method. The ELVIRA method produces a
much more accurate representation of the normal, which is constant from cell
to cell, whereas the first order method varies slightly, causing oscillations
near the interface.

Figure 1: Different shock reflections problems. The images show (left) the
density and (right) a numerical schlieren with AMR levels. (top) The Mach 1.7
single shock reflection test at $t=1.2$ and (bottom) the Mach 10 double shock
reflection test at $t=0.2$. The rigid body is depicted in grey.

###  4.2 Shardin Wedge

Similar to the previous test is the Shardin wedge problem, for which there is
a range of experimental and numerical comparison. The original study by
Schardin [27] took high-frequency schlieren photographs of a shock hitting a
triangular wedge. This work examines the near identical case from Chang and
Chang [8], in which a Mach 1.34 shock hits an equilateral triangle. The
quantitative experimental comparison for the test is obtained by tracking the
positions of different flow features over the course of the simulation. The
positions are non-dimensionalised with respect to the triangle length and
width to enable comparison to a broader range of tests. Similar experiments
can be found in the literature: [6, 9, 22, 26, 29].

The initial set-up for this test is shown in Figure 2. The test is run for
0.35 ms, with a CFL of 0.6. A base resolution of $150\times 100$ was used,
along with 4 layers of AMR, each of refinement factor 2. The gas is again a
$\gamma=1.4$ ideal gas, with the quiescent state given by: $\rho=1\
\text{kgm}^{-3},p=0.05\ \text{MPa},\mathbf{u}=\mathbf{0}$.

Figure 2: The initial conditions for the Schardin wedge problem.

The results of the test and the experimental comparison are shown in Figure 4
and Figure 4. Excellent agreement with the results of Chang and Chang [8] is
observed.

Figure 3: The Schardin wedge test. The images show a numerical schlieren of
the solution, with times given in $1\times 10^{-4}$s. The last image shows the
AMR levels used. Experimental comparison is drawn by tracking the position of
different flow features: the upper shock triple point (TP1), the lower shock
triple point (TP2) and the vortex core (V).

Figure 4: Comparison of the present work with the previous experiments and
simulations of Chang and Chang [8]. The locations of the upper triple point
(TP1), lower triple point (TP2) and the vortex core are tracked over the
course of the simulation. Excellent agreement is observed. Each axis is non-
dimensionalised with respect to the dimensions of the triangle, $a$ and $b$,
outlined in Figure 2.

###  4.3 Cylinder Test

Another shock-body interaction problem is that of a cylinder being hit by a
shock wave. This test has been studied by many other authors [7, 9, 21, 26,
39], with both numerical and experimental comparison, making it a suitable
validation case.

The test is non-dimensionalised, with a domain spanning $x=[0:2]$, $y=[0:1]$
and a reflective boundary condition along both the upper and lower $x$-axis.
The cylinder has a radius of 0.105 and is centred at $(0.5,0)$. The test is
run until $t=0.4$ with a CFL of 0.6. A base resolution of 200 $\times$ 100 is
used, with 4 levels of AMR, all refinement factor 2. A Mach 2.81 shock is
initialised at $x=0.385$, with the quiescent state being a $\gamma=1.4$ ideal
gas with $\rho=1,p=1,\mathbf{u}=\mathbf{0}$.

The test is shown in Figure 6, where the images depict a numerical schlieren
of the density profile at different times. Figure 6 shows the numerical and
experimental comparison with previous works, and excellent agreement is
observed.

Figure 5: The cylinder test. The images show a numerical schlieren of the
solution, with the last image showing the AMR levels used. Experimental and
numerical comparison is drawn by tracking the position of different flow
features: the upper shock triple point (TP1) and the lower shock triple point
(TP2).

Figure 6: Comparison of cylinder test with the previous experiments and
simulations. The locations of the upper triple point (TP1) and lower triple
point (TP2) are tracked over the course of the simulation. Excellent agreement
is observed with both the previous simulations of Chaudhuri et al. [9] and
Zółtak and Drikakis [39], and the experimental studies of Bryson and Gross [7]
and Kaca [21]. Each axis is non-dimensionalised with respect to the diameter
of the cylinder, $D$.

###  4.4 Moving Cylinder Test

To test the ability of the model to handle dynamic rigid bodies, a moving
cylinder is considered. The test is again non-dimensionalised, with the
surrounding gas being a $\gamma=1.4$ ideal gas with
$\rho=1.4,p=1,\mathbf{u}=\mathbf{0}$. The domain for this test is $x=[0:5],\
y=[-2:2]$, with a resolution of 400 $\times$ 320\. The rigid cylinder has a
radius of 0.5 and is initially located at $(x,y)=(0.6,0)$. The velocity of the
rigid body is constant at $\mathbf{v}_{r}=3\ \hat{x}$. The test is run until
$t=1$ with a CFL of 0.6. For this test, the Riemann problem version of the
rigid seeding method is employed.

The test is shown in Figure 7, depicting density, pressure, $x$-velocity,
$y$-velocity, a numerical schlieren and the rigid volume fraction. The rigid
body volume fraction is kept sharp thanks to the THINC reconstruction
performed in its update.

Figure 7: The Riemann problem based seeding moving cylinder test. The images show (from left to right and top to bottom) density, pressure, $x$-velocity, $y$-velocity, a numerical schlieren and the rigid volume fraction. 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
Copper | 8960.0 | 130.1 | 43.33 | 1.0 | 3.0 | 2.0
Steel | 8030.0 | 156.2 | 77.2 | 0.569 | 2.437 | 1.563
Table 1: Equation of state material parameters for equation (23). Material |  $c_{1}$ GPa |  $c_{2}$ GPa | $c_{3}$ | $n$ | $m$ |  $T_{\text{melt}}$ K |  $C_{V}$ J kg-1K-1
---|---|---|---|---|---|---|---
Copper | 0.09 | 0.292 | 0.025 | 0.31 | 1.09 | 1358.0 | 385.0
Table 2: Johnson Cook plasticity parameters.

###  4.5 Solid Extrusion with a Piston

This test, proposed in Hill et al. [18], examines the interaction of both
static and dynamic rigid bodies with an elastoplastic solid. The test consists
of a rigid piston pushing a block of aluminium through a rigid wedge,
representing an extrusion problem. The test is more complex than the fluid-
wedge interaction problems, both because of the driving piston and the
separation of the elastic precursor from the plastic wave. The aluminium in
this test is given by the equation of state equation (23) with parameters
given in Table 1. Following Hill et al. [18], the aluminium has an ideal
plasticity law with the yield stress given by $\sigma_{Y}=0.2976$ GPa.

The domain for this test is $x=[-21:120]$ mm, $y=[0:120]$ mm, with a
reflective upper boundary. A base resolution of $256\times 216$ is used, with
3 levels of AMR, refinement factor 2. The test is run for 20 $\mu$s with a CFL
of 0.6. The piston is initially located at $x=-20$ mm and given a constant
velocity of $\mathbf{v}_{r}=100\hat{x}$ ms-1. The wedge has an angle of 15∘
with respect to the $x$-axis and starts at $x=0$. The method copes well with
this strenuous problem, as can be seen from the results in Figure 8.

Figure 8: The solid extrusion with a piston test at t = 20 $\mu$s. The images
show (from top to bottom, left to right) the density, AMR grids and density
contours, $x$-velocity and $y$-velocity. The test demonstrates the stability
of the method when handling a combination of elastoplastic solids and both
static and dynamic rigid boundaries.

###  4.6 Solid Impact Extrusion

Another example of the use of a rigid body and elastoplastic solid formulation
together is the impact extrusion test from Gray et al. [16]. In this test, a
copper sphere impacts an extrusion cone, resulting in a jet of metal being
ejected. The jet has a much higher velocity than the sphere, and undergoes
necking, breaking into fragments. The test can be compared with both the
results of Gray et al. [16] and other works in the literature [5, 19, 20, 23,
25].

This test is an example of the particular usefulness of the formulation at
hand. In experiment, the hardened steel extrusion cone is not expected to
deform significantly, meaning it can be replaced with a rigid body in
numerical studies, greatly reducing the cost of the test. Additionally, the
experiment is run surrounded by void using the model of Wallis et al. [33].
This removes the need to track the high speed waves in any surrounding fluid
which ultimately have little effect on the motion of the elastoplastic solid.

The set-up for the test is shown in Figure 9(a). The copper is given an
initial downward velocity of 400 ms-1 and is governed by the equation of state
equation (23), with parameters given in Table 1. The test is run using the
cylindrically symmetric form, with a reflective boundary condition along the
$r=0$ axis. The domain spans $r=[0:5]$ mm, $z=[-40:22.5]$ mm. A base
resolution of 56 $\times$ 696 is used, with two layers of AMR, each of
refinement factor 2. The test is run until $t=60\ \mu$s using a CFL of 0.6,
and is shown in Figure 10.

The results agree well with the experimental results [16], shown in Figure
9(b), despite the lack of a material damage model. Additionally, the maximum
plastic strain of 9.327 agrees well with the value of 9.3 quoted by Sambasivan
et al. [25] and the range 6-9 quoted by Gray et al. [16]. This test also
serves to demonstrate the application of more complex plasticity models, such
as the Johnson Cook model employed here, correctly producing necking effects
in the solid jet. However, the choice of plasticity model and parameters has a
significant effect on the relative fragment sizes, and a more in-depth
comparison of different models would be needed to be able to draw further
experimental comparison.

(a) The elastoplastic impact extrusion test initial conditions. Lengths shown
are in mm. In the domain, the base of the rigid body is positioned at $z=0$.

(b) The experimental results of elastoplastic impact extrusion test performed
by Gray et al. [16], presenting three different copper grain sizes. Reproduced
from Gray et al., Influence of Shock Prestraining and Grain Size on the
Dynamic‐Tensile‐Extrusion Response of Copper: Experiments and Simulation, AIP
Conference Proceedings, 2006 [16], with the permission of AIP Publishing.

Figure 9: The impact extrusion test. Figure 10: The elastoplastic impact
extrusion test. An elastoplastic sphere of copper collides with a rigid
extrusion cone (shown in grey), producing a jet of ejected metal. The metal is
surrounded by vacuum, using the void model outlined by Wallis et al. [33]. The
times are shown in $\mu$s, and the images depict the $y$-velocity and plastic
strain. The results agree well with the experimental studies of Gray et al.
[16].

I /s-1 | G1 / ($10^{11}$ Pa)-y s-1 | G2 / ($10^{11}$ Pa)-z s-1
---|---|---
$4\times 10^{12}$ | 4500 $\times 10^{6}$ | 30 $\times 10^{6}$
Table 4: The ignition and growth reaction rate parameters for LX-17.

Figure 11: The rigid LX-17 corner turning test, depicting the pressure in the
explosive. Times are shown in $\mu$s. The dead-zones are tracked by overlaying
a plot of the reaction progress variable with varying opacity; fully reacted
material is transparent and unreacted material is slightly opaque.

Figure 12: The elastoplastic solid LX-17 corner turning test, depicting the
pressure in both the explosive and the solid. Times are shown in $\mu$s. The
dead-zones are tracked by overlaying a plot of the reaction progress variable
with varying opacity; fully reacted material is transparent and unreacted
material is slightly opaque. Figure 13: The 3D STL geometry test. The test
features a Mach 1.5 flow over a complex three-dimensional geometry, defined
purely as volume fraction field. The images show the density and AMR levels at
$t=2,4$. The method handles this difficult test well, and is straightforward
to parallelise, even in three dimensions with AMR.

###  4.7 LX-17 Corner Turning

To demonstrate the multi-physics capabilities of the current model, the
desensitised LX-17 corner turning tests studied by Souers et al. [30],
Deoliveira et al. [12] and Tarver [31] are considered. In this test, the
domain is filled with the explosive LX-17, with a ‘hockey-puck’ of solid or
rigid material in the centre. A booster region directly above the puck ignites
the explosive, and the subsequent detonation wave travels around the corner of
the puck. There are several multi-physics aspects to this test. Firstly, the
explosive LX-17 is modelled using the multi-phase formulation detailed by
Wallis et al. [34], where both the reactants and products are explicitly
modelled as a mixture with distinct equations of state. This mixture is then
treated as one material for the purposes of bookkeeping, and the reaction
progress variable $\lambda$ is added to the full multi-material system of
equations in order to track the reaction:

$\displaystyle\frac{\partial\phi_{(l)}\rho_{(l)}\lambda_{(l)}}{\partial
t}+\frac{\partial\phi_{(l)}\rho_{(l)}\lambda_{(l)}u_{k}}{\partial x_{k}}$
$\displaystyle=\phi_{(l)}\rho_{(l)}\dot{\lambda}_{(l)}\ .$ (66)

The explosive’s products and reactants are governed by the JWL equation-of-
state, with parameters given in Table LABEL:tab:JWL. This approach further
allows for the use of the ignition-and-growth model, as employed by Deoliveira
et al. [12], to simulate the rate of reaction of the explosive. This is a
three-stage rate law, based on phenomenological experience of how detonations
evolve in condensed phase explosives. The rate can be expressed as:

$\displaystyle\dot{\lambda}=$ $\displaystyle
I(1-F)^{b}\left(\frac{\rho}{\rho_{0}}-1-a\right)^{x}H(F_{ig}-F)$ (67)
$\displaystyle+G_{1}(1-F)^{c}F^{d}p^{y}H(F_{G_{1}}-F)$
$\displaystyle+G_{2}(1-F)^{e}F^{g}p^{z}H(F-F_{G_{2}})\ ,$

where $F=1-\lambda$ is the reacted fraction and $H$ is the Heaviside function.
Other parameters are material dependent constants, given in Table 4.6. Full
details of the method are given by Wallis et al. [34].

Secondly, an explosive desensitisation model is also implemented. The model
was used by Deoliveira et al. [12] to study the formation of ‘dead-zones’
around the corner of the puck. These are regions where the explosive is not
ignited as the detonation diffracts around the corner, due to the shock-
compression of the material. This model is implemented by including an
additional history variable $\varphi$ which tracks the desensitisation of the
explosive:

$\displaystyle\frac{\partial\phi_{(l)}\rho_{(l)}\varphi_{(l)}}{\partial
t}+\frac{\partial\phi_{(l)}\rho_{(l)}\varphi_{(l)}u_{k}}{\partial
x_{k}}=\phi_{(l)}\rho_{(l)}Ap(1-\varphi_{(l)})(\varphi_{(l)}+\epsilon)\ ,$
(68)

where $A=8.36\times 10^{-3}$ is a material dependent constant, $p$ is the
pressure and $\epsilon=10^{-3}$ is a small value designed to increase
$\varphi$ from zero. When the explosive is shocked, the desensitisation goes
from 0 (no desensitisation) to 1 (fully desensitised). This parameter is then
fed into the ignition-and-growth reaction model to slow the reaction in areas
that have been desensitised. To this end, parameter $a$ in the ignition-and-
growth model is modified to:

$\displaystyle a(\varphi)=a_{0}(1-\varphi)+a_{1}\varphi\ ,$ (69)

where $a_{0}$ was the original parameter, and $a_{1}=0.5$ is a parameter used
to calibrate the desensitisation. Additionally, the first growth term,
$G_{1}$, is now only activated when the reaction progress variable is greater
than a threshold, $F_{G_{1},\text{min}}=\lambda_{c}\varphi$, with
$\lambda_{c}=0.01$ being another calibration parameter. This means the
reaction rate is now:

$\displaystyle\dot{\lambda}=$ $\displaystyle
I(1-F)^{b}\left(\frac{\rho}{\rho_{0}}-1-(a_{0}(1-\varphi)+a_{1}\varphi)\right)^{x}H(F_{ig}-F)$
(70)
$\displaystyle+G_{1}(1-F)^{c}F^{d}p^{y}H(F_{G_{1}}-F)H(F-\lambda_{c}\varphi)$
$\displaystyle+G_{2}(1-F)^{e}F^{g}p^{z}H(F-F_{G_{2}})\ .$

Thirdly, this test can either be performed with a rigid body or with an
elastoplastic solid, demonstrating how the methodology at hand is capable of
transitioning between either approach with no drastic change of numerical
scheme. The solid in this case is steel, governed by the equation of state
equation (23) with the parameters given in Table 1, obeying a ideal plasticity
law with a yield stress of $\sigma_{Y}=1.37$ GPa.

Both of tests have the same set-up. The tests are run in the cylindrically
symmetric formulation, adding radial source terms to model the effect of the
geometry and including a reflective boundary condition at the $r=0$ axis. The
domain spans $r=[0:44]$ mm, $z=[-15:40]$ mm. A base resolution of $272\times
344$ is employed, with 3 levels of AMR each of refinement factor 2. The tests
are run to a final time of 5 $\mu$s with a CFL of 0.6.

The puck has a height of 15 mm, a radius of 19.05 mm and a rounded corner with
a radius of 0.5 mm. The detonation is ignited by a hemispherical booster
region on top of the puck, centred at (0,0) with a radius of 7.68 mm.

Both the rigid and elastoplastic tests are shown in Figure 12 and 12. The
tests depict the pressure profile over the course of the test, showing how the
detonation diffracts around the corner. The dead-zones are tracked by
overlaying a plot of the reaction progress variable with varying opacity;
fully reacted material is transparent and unreacted material is slightly
opaque. Both cases handle the test well and demonstrate the marked difference
between discretising inside the solid, compared to assuming a rigid body. The
rigid tests agree well with the results of Deoliveira et al. [12].

###  4.8 Stereolithography (STL) Geometry

The method is designed to be easily extendible to large 3D simulations with
multiple levels of AMR. This example is designed to test the ability of the
model to handle arbitrary three-dimensional geometries, such as would be used
in research applications. In the test, the geometry in question is a scale
model of the NASA space shuttle [17], manually edited to form a watertight
geometry. The volume fraction field of the body is extracted from the
triangulated surface mesh by first using the VTK [28] object
vtkImplicitPolyDataDistance to obtain a signed distance function. From here,
the volume of the body in a cell is computed as the volume of the polyhedron
formed by the intersection of the signed distance function zero contour with
the cell, using the formula from Barton et al. [3]. The volume fraction is
then set to be the rigid body volume fraction.

The test is non-dimensionalised, with a Mach 1.5 flow over the body
initialised as $\rho=1.4,p=1,\mathbf{u}=1.5\hat{x}$. The domain spans
$x=[-10.0:10.0]$, $y=[-7.5:7.5]$, $z=[-5.0:5.0]$, with a base resolution of
$176\times 128\times 88$ and 2 levels of AMR, each of refinement factor 2. The
test is run until $t=4$, using a CFL of 0.3.

The test is shown in Figure 13, depicting the density and AMR levels at
$t=2,4$. The test performs well, with the method working well even in with
complex three-dimensional geometries.

##  5 Conclusions

This work has demonstrated a straightforward approach to the inclusion of
immersed rigid boundaries in a variety of multi-physics applications. The
example test cases have included single fluid tests with both static and
dynamic rigid bodies, the interaction of elastoplastic solids with rigid
bodies, multi-phase reactive fluids, and complex three-dimensional geometries.
All of these test cases are modelled using the same diffuse interface model.
The rigid body is included as a scalar field, and all material interaction is
incorporated into the underlying system of equations which can handle an
arbitrary number of solid and fluid components. The model is shown to be
extendible to three dimensions and adaptive mesh refinement. This is thanks to
the local, parallelisable nature of both the hyperbolic update and the
interface seeding method which mediates the rigid body boundary conditions.

## Acknowledgements

This work was funded by AWE PLC. Additionally, Tim Wallis is supported 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.

## Author declaration

The authors have no conflicts of interest to disclose.

## Data availability statement

The data that support the findings of this study are available within the
article.

## 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.
  * 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 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.
  * 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 et al. [2015] N Bonora, G Testa, A Ruggiero, G Iannitti, N Mortazavi, and M Hörnqvist.  Numerical simulation of dynamic tensile extrusion test of ofhc copper.  _Journal of dynamic behavior of materials_ , 1(2):136–152, 2015.  ISSN 2199-7446.
  * Boukharfane et al. [2018] Radouan Boukharfane, Fábio Henrique Eugênio Ribeiro, Zakaria Bouali, and Arnaud Mura.  A combined ghost-point-forcing / direct-forcing immersed boundary method (ibm) for compressible flow simulations.  _Computers & fluids_, 162:91–112, 2018.  ISSN 0045-7930.
  * Bryson and Gross [1961] A. E Bryson and R. W. F Gross.  Diffraction of strong shocks by cones, cylinders, and spheres.  _Journal of fluid mechanics_ , 10(1):1–16, 1961\.  ISSN 0022-1120.
  * Chang and Chang [2000] Se-Myong Chang and Keun-Shik Chang.  On the shock–vortex interaction in schardin’s problem.  _Shock waves_ , 10(5):333–343, 2000.  ISSN 0938-1287.
  * Chaudhuri et al. [2011] Arnab Chaudhuri, Abdellah Hadjadj, and Ashwin Chinnayya.  On the use of immersed boundary methods for shock/obstacle interactions.  _Journal of computational physics_ , 230(5):1731–1748, 2011.  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.
  * 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.
  * Deoliveira et al. [2006] G. Deoliveira, A. K. Kapila, and D. W. Schwendeman.  Detonation diffraction, dead zones and the ignition-and-growth model.  In _The Thirteenth Symposium (International) on Detonation._ , volume 16, 2006.
  * 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.
  * Forrer and Jeltsch [1998] Hans Forrer and Rolf Jeltsch.  A higher-order boundary treatment for cartesian-grid methods.  _Journal of computational physics_ , 140(2):259–277, 1998.  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.
  * Gray et al. [2006] G. T. Gray, E. Cerreta, C. A. Yablinsky, L. B. Addessio, B. L. Henrie, B. H. Sencer, M. Burkett, P. J. Maudlin, S. A. Maloy, C. P. Trujillo, and M. F. Lopez.  Influence of shock prestraining and grain size on the dynamic‐tensile‐extrusion response of copper: Experiments and simulation.  _AIP Conference Proceedings_ , 845(1):725–728, 2006.  doi: 10.1063/1.2263424.
  * Handbook [2009] NASA Johnson Space Center Space Educators’ Handbook.  Shuttle model 1, 2009.  URL https://nasa3d.arc.nasa.gov/detail/sts.
  * Hill et al. [2010] D.J Hill, D Pullin, M Ortiz, and D Meiron.  An eulerian hybrid weno centered-difference solver for elastic–plastic solids.  _Journal of computational physics_ , 229(24):9053–9072, 2010.  ISSN 0021-9991.
  * Hörnqvist et al. [2015] M. Hörnqvist, N. Mortazavi, M. Halvarsson, A. Ruggiero, G. Iannitti, and N. Bonora.  Deformation and texture evolution of ofhc copper during dynamic tensile extrusion.  _Acta Materialia_ , 89:163–180, 2015.  ISSN 1359-6454.  doi: https://doi.org/10.1016/j.actamat.2015.01.053.  URL https://www.sciencedirect.com/science/article/pii/S135964541500066X.
  * Iannitti et al. [2012] Gianluca Iannitti, Nicola Bonora, Andrew Ruggiero, and Simone Dichiaro.  Modeling ductile metals under large strain, pressure and high strain rate incorporating damage and microstructure evolution.  _AIP Conference Proceedings_ , 1426(1):1027–1030, 2012.  doi: 10.1063/1.3686453.
  * Kaca [1988] Jedrzej Kaca.  An interferometric investigation of the diffraction of a planar shock wave over a semicircular cylinder.  Technical report, UTIAS Technical Note, No. 269, 1988.
  * Kemm et al. [2020] Friedemann Kemm, Elena Gaburro, Ferdinand Thein, and Michael Dumbser.  A simple diffuse interface approach for compressible flows around moving solids of arbitrary shape based on a reduced baer–nunziato model.  _Computers & fluids_, 204:104536, 2020.  ISSN 0045-7930.
  * Park et al. [2013] Kyung-Tae Park, LeeJu Park, Hak Jun Kim, Seok Bong Kim, and Chong Soo Lee.  Dynamic tensile extrusion behavior of coarse grained and ultrafine grained ofhc cu.  _Materials Science and Engineering: A_ , 569:61–70, 2013\.  ISSN 0921-5093.  doi: https://doi.org/10.1016/j.msea.2013.01.015.  URL https://www.sciencedirect.com/science/article/pii/S0921509313000567.
  * Pilliod and Puckett [2004] James Edward Pilliod and Elbridge Gerry Puckett.  Second-order accurate volume-of-fluid algorithms for tracking material interfaces.  _Journal of computational physics_ , 199(2):465–502, 2004.  ISSN 0021-9991.
  * Sambasivan et al. [2013] S Sambasivan, A Kapahi, and H.S Udaykumar.  Simulation of high speed impact, penetration and fragmentation problems on locally refined Cartesian grids.  _Journal of Computational Physics_ , 235:334–370, 2013\.  ISSN 0021-9991.
  * Sambasivan and UdayKumar [2009] Shiv Kumar Sambasivan and H. S UdayKumar.  Ghost fluid method for strong shock interactions part 2: Immersed solid boundaries.  _AIAA journal_ , 47(12):2923–2937, 2009.  ISSN 0001-1452.
  * Schardin [1957] H. Schardin.  High frequency cinematography in the shock tube.  _The Journal of Photographic Science_ , 5(2):17–19, 1957.  doi: 10.1080/00223638.1957.11736588.
  * Schroeder [2006] Will Schroeder.  _The visualization toolkit : an object-oriented approach to 3D graphics / Will Schroeder, Ken Martin, Bill Lorensen, with special contributors Lisa Sobierajski Avila, Rick Avila, C. Charles Law._ Kitware, New York], 4th ed. edition, 2006.  ISBN 9781930934191.
  * Sivier et al. [1992] S Sivier, E Loth, J Baum, and R Löhner.  Vorticity produced by shock wave diffraction.  _Shock waves_ , 2(1):31–41, 1992.  ISSN 0938-1287.
  * Souers et al. [2004] P. Clark Souers, Henry G Andreski, Charles F Cook, Raul Garza, Ron Pastrone, Dan Phillips, Frank Roeske, Peter Vitello, and John D Molitoris.  Lx‐17 corner‐turning.  _Propellants, explosives, pyrotechnics_ , 29(6):359–367, 2004.  ISSN 0721-3115.
  * Tarver [2010] Craig M Tarver.  Corner turning and shock desensitization experiments plus numerical modeling of detonation waves in the triaminotrinitrobenzene based explosive LX-17.  _The journal of physical chemistry. A_ , 114(8):2727, 2010.  ISSN 1089-5639.
  * Toro et al. [1994] E. Toro, M. Spruce, and W. Speares.  Restoration of the contact surface in the HLL-Riemann solver.  _Shock Waves_ , 4(1):25–34, 1994.  ISSN 0938-1287.
  * Wallis et al. [2021a] Tim Wallis, Philip T. Barton, and Nikolaos Nikiforakis.  A flux-enriched godunov method for multi-material problems with interface slide and void opening.  _Journal of Computational Physics_ , 2021a.  doi: https://doi.org/10.1016/j.jcp.2021.110499.
  * Wallis et al. [2021b] Tim Wallis, Philip T. Barton, and Nikolaos Nikiforakis.  A diffuse interface model of reactive-fluids and solid-dynamics.  _Computers & Structures_, 2021b.  doi: https://doi.org/10.1016/j.compstruc.2021.106578.
  * Woodward and Colella [1984] Paul Woodward and Phillip Colella.  The numerical simulation of two-dimensional fluid flow with strong shocks.  _Journal of Computational Physics_ , 54(1):115–173, 1984.  ISSN 0021-9991.
  * Xiao et al. [2005] F. Xiao, Y. Honma, and T. Kono.  A simple algebraic interface capturing scheme using hyperbolic tangent function.  _International Journal for Numerical Methods in Fluids_ , 48(9):1023–1040, 7 2005.  ISSN 1097-0363.  doi: 10.1002/fld.975.  URL https://doi.org/10.1002/fld.975.
  * Youngs [1984] David Youngs.  An interface tracking method for a 3d Eulerian hydrodynamics code.  Technical report, AWRE, 1984.
  * Zhang et al. [2019] Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale.  AMReX: a framework for block-structured adaptive mesh refinement.  _Journal of Open Source Software_ , 4(37):1370, May 2019.  doi: 10.21105/joss.01370.  URL https://doi.org/10.21105/joss.01370.
  * Zółtak and Drikakis [1998] J Zółtak and D Drikakis.  Hybrid upwind methods for the simulation of unsteady shock-wave diffraction over a cylinder.  _Computer methods in applied mechanics and engineering_ , 162(1):165–185, 1998.  ISSN 0045-7825.