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

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 mixture 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.


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 [48] or cut-cell approaches such as Gokhale et al. [27]. However these methods have their own associated problems. Recently Kemm et al. [36] 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 both the material interfaces and rigid body interfaces. The methods presented build on the work of Wallis et al. [57], 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. [61] to be used to keep interfaces sharp.
The work of Wallis et al. [57] used a reduced-equation Allaire-type [1] multi-material diffuse interface scheme, adapted by Barton [10] 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 2.1. Base System 2.1. 1

. Evolution Equations
The base system, devised by Barton [10], 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 ρ (l) , volume fraction φ (l) , symmetric left unimodular stretch tensor V e , velocity vector u, and internal energy 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, . . . , N materials: where E = E + |u| 2 /2 denotes the specific total energy, σ denotes the Cauchy stress tensor, β j = ∂V e k j /∂x k , and Φ represents the contribution from plastic effects. The plastic source term Φ limits the possible range of elastic deformation and initiates the growth of plastic deformation. This work follows the method of convex potentials, defining this source term as where χ is a plastic multiplier and ϕ is a convex potential, both of which are closure models. Further details of this source term can be found in Barton [10] and comprehensive details of numerous alternative closure models can be found in Asaro and Lubarda [2]. Some multi-physics closure models introduce a dependence on material history variables. For example, the equivalent plastic strain ε p(l) or the damage parameter D (l) in the case of fracture [57], or the reaction progress variable λ for reactive fluid mixtures [58]. For these variables, additional evolution equations are required: Here α (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 α represents a vector.

Thermodynamics
The internal energy E for each material is defined by an equation-of-state that conforms to the general form: where 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, E c (l) ρ (l) , α (l) , the contribution due to temperature deviations, E t (l) ρ (l) , T (l) , and the contribution due to shear strain E s (l) ρ (l) , dev (H e ) , α (l) . The cold compression energy will generally be provided by the specific closure model for each material. The thermal energy is given by where C V (l) is the heat capacity, T 0 (l) is a reference temperature, and θ D (l) (ρ (l) ) is the non-dimensional Debye temperature. The Debye temperature is related to the Grüneisen function, Γ(ρ (l) ), via 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 where G (l) (ρ (l) , α (l) ) is the shear modulus, and is the second invariant of shear strain. For each component, the Cauchy stress, σ, and pressure, p, are inferred from the second law of thermodynamics and classical arguments for irreversible elastic deformations: 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. It can be seen that, by presenting the internal energy in the form outlined, equation (15) can be written in the form: Here, E ref,(l) = E c (l) + E s (l) and p ref,(l) = ρ 2 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 in Section 2.1.3. Following the approach of Barton [10], rather than calculating the full acoustic tensor, the sound-speed c (l) for a material is approximated by: where where S is the entropy and is the transverse wave speed squared.

Closure Models
For the case of fluids like air, a simple ideal gas law will be used. In this case, p ref = 0, E ref = 0, and the material is entirely specified by the choice of the adiabatic index γ = Γ + 1 which is taken to be constant. This reduces the form of Equation 17 to: When elastoplastic solids are employed, this work follows previous studies by using the equation of state from Dorovskii et al. [21]. For this equation of state, Γ(ρ) = Γ 0 , where Γ 0 is a material dependent constant, and: 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 Γ(ρ) = Γ 0 and

Mixture Rules
Having defined the thermodynamics and closure models for each separate material, they may now be combined into a mixture. In order to do this, 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,10]: where c is the sound-speed and Y (l) = (φ (l) ρ (l) )/ρ is the mass fraction. Using these mixture rules and the form of the Mie-Grüneisen equation of state, the total hydrodynamic pressure and total stress can be obtained from the total energy. Initially, to calculate the total pressure, equation 29 is combined with the isobaric assumption that forms part of the mechanical equilibrium of the underlying Allaire et al. [1] model. First the internal energy is obtained: Then the Mie-Grüneisen form of the equation of state is used to rewrite equation 29, using the isobaric assumption: Finally, this equation can be rewritten in terms of the total pressure p: which may alternatively be written as a weighted sum of component pressures: To obtain the total stress, the total pressure is combined with the mechanical equilibrium assumption that all materials share the same deviatoric strain and the mixture rule for the shear modulus 30: The sound-speed mixture rule is used for the estimation of the global time-step and the calculation of the wave speeds required for the flux calculation, as detailed in Barton [10].

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 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. This gives the evolution equation:

Reduced System
Clearly not every one of the equations presented above will be relevant for a given application. Frequently, singlematerial 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: Here the stress tensor takes the simple form σ i j = −pδ i j , where p is the hydrodynamic pressure and δ 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.

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. [57]. 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 [63], which includes an implementation of the structured adaptive mesh refinement (SAMR) method of Berger and Colella [12] 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. ∆t l := ∆t l−1 /r l and ∆x l := ∆x l−1 /r l with r l ∈ N, r l ≥ 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 ∈ Z and each cell C l i jk has the dimensions ∆x l i jk . 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 and non-conservative terms from the rigid volume fraction update for each spatial dimension: Note that the non-conservative term will also include the non-conservative updates for the material volume fraction and stretch tensor if the full system is used. For the reduced case, this system is given by: The inhomogeneous system is integrated for time intervals [t n , t n+1 ], where the time-step ∆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 [18].

Hyperbolic Update
Replacing the spatial derivatives with a discretised conservative approximation, the hyperbolic system can be written where U i jk represents the vector of conservative variables stored at cell centres, and where F d m±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 (45), resulting in a spatially unsplit method.
The fluxes for the thermodynamic variables are computed using the HLLC solver developed by Toro et al. [56] (note that the full multiphase system uses the HLLD solver from Barton [10]). 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. [19]. 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. © British Crown Owned Copyright 2021/AWE 6 The fluxes for the rigid body volume fraction are calculated using a simple upwind method, based on the sign of v r and the cell-wall reconstructed values of r: r L and r R . For example, in the x-direction, let v r be the x-component of the rigid body velocity at the cell-wall, (v r ) i+1/2, jk , then the flux for the rigid volume fraction is given by: If greater accuracy is required for the advection of the rigid body field, this method can be replaced by any desired scheme, such as volume of fluid methods. The rigid body volume fraction is not a thermodynamic variable, so it may be updated completely separately to the other variables. Thermodynamic variables do not need to be updated in rigid (r > 0.5) cells, saving computational expense, as they will be replaced by the interface seeding method in Section 3.2.
Finally, the non-conservative term for the rigid body volume fraction is calculated. This term is added using a simple difference formula, similar to previous works [10] for consistency. This term can also be written as a sum of 'fluxes'. In what follows, variables with subscript L or R are evaluated at the cell wall on either side of a cell in a given direction. For example, in the x direction L → (i − 1/2, jk) and R → (i + 1/2, jk). All other variables are evaluated at the cell centre (i jk). The non-conservative term can then be written: where with similar definitions for the other directions. When the full multiphase system is used, the same form for the non-conservative update is used for the material volume fraction and stretch tensor update. For the material volume fraction, in direction d: where u * is the velocity at the cell-wall, found from the HLLD Riemann solver in the d direction. For the stretch tensor: where V e * is the stretch tensor at the cell-wall, found from the HLLD Riemann solver in the d direction. To achieve a higher temporal resolution in the update of the hyperbolic terms, the third order Runge-Kutta time integration scheme is used:

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 (37). 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: wheren r is the rigid body normal vector. In other words, the velocity aligned with the rigid interface normal component must equal the rigid body velocity in that direction. For the slip boundary condition, all other tangential velocity components remain unchanged across the interface. The no-slip condition is not employed here, but is enforced by the entire velocity vector to equal the rigid body velocity at the interface, resulting in the condition u = v r . Therefore, a method is required that will fill the rigid cells (r > 0.5) with a suitable thermodynamic state such as to reproduce these boundary conditions. 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. [57], the rigid body boundary conditions are mediated by flux-modifiers and interface seeding routines. However, by comparison, no new flux modifiers are required: only an interface seeding method is required for rigid body interaction. This is because the rigid body boundary only involves one real material, with the flux inside the rigid body not being relevant, therefore meaning that a two-flux method such as those presented by Wallis et al. [57] is not required.
Although no additional flux-modifiers are required, rigid bodies must be seeded with suitable states to enable the correct boundary interaction, as previously mentioned. This is done using an interface seeding routine, similar to those presented in Wallis et al. [57]. 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 will be reflected across a rigid body interface and those that will not, denoted by U ref and U const respectively. All scalar quantities are constant across the rigid interface. For vector and tensor quantities, it is easiest to define these quantities in the reference frame wheren r points in the x-direction 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, U →Ũ, and are found in the normal way, where for a vector v and matrix M: The rotation matrix R is detailed in Wallis et al. [57].
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,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.5dx, or sample points can chosen at various distances stepping out in the normal direction (from dx to 4dx), and when the 2 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 U 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 thex direction, givingŨ Interp . 4. A new state,Ũ new , is then created by reflecting the suitable components of the interpolated cell: If the rigid body is dynamic, then the normal component of the rigid body velocity must also be included. This is done adding 2ṽ ref r = 2(v r ·n r ) to the reflected velocity ofŨ new . Therefore, in total:

5.
Here, the method then diverges with two options: the original or Riemann problem based version.
• For the original method, this new state is then rotated back into the original frame of reference and replaces the original state.
U ← U new .
• 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Ũ new andŨ Interp respectively. This gives the intermediate statesŨ * L andŨ * R . TheŨ * R state is then rotated back into the original frame of reference and replaces the original state.
This procedure will then fill a few layers of cells inside the rigid body interface with a state that will produce the desired rigid body boundary condition, and the hyperbolic update can proceed as normal.
Although it has been mentioned that no new flux-modifiers are required for the method, the existing void-opening flux modifier can be 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. This then allows for void opening in solid-rigid collisions. The flux modifier is applied when the rigid volume fraction 'changes sign' between the two cells being used to calculate the flux:

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 [62]. Youngs' method is a first order approximation to the normal, attained using a simple difference formula for the local 3 DIM cells. The method is cheap, and is easily extended into three dimensions. Following the notation of Youngs [62], the normal of a field f is estimated as: where symmetric difference formulae must be used for each direction. These formulae are outlined by Youngs [62] and Pilliod and Puckett [44]. 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 [44], 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 [44]. 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.

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

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. [36], and the Mach 10 test can be compared to the similar tests in Woodward and Colella [60] and Forrer and Jeltsch [25].
Both of these tests are non-dimensionalised, with a domain x = [0 : 3], y = [0 : 2], a base resolution of 300 × 200, and two levels of AMR, both refinement factor 2. The gas is a γ = 1.4 ideal gas, with a quiescent state of ρ = 1.4, p = 1.0, u = 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.4 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.

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 [49] took high-frequency schlieren photographs of a shock hitting a triangular wedge. This work examines the near identical case from Chang and Chang [16], 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: [14,17,36,48,52].
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.4. A base resolution of 150 × 100 was used, along with 4 layers of AMR, each of refinement factor 2. The gas is again a γ = 1.4 ideal gas, with the quiescent state given by: ρ = 1 kgm −3 , p = 0.05 MPa, u = 0.
The results of the test and the experimental comparison are shown in Figure 3 and Figure 4. Excellent agreement with the results of Chang and Chang [16] is observed.

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 [15,17,35,48,64], with both numerical and experimental comparison, making it a suitable validation case.
The test is non-dimensionalised, with a domain spanning x = This test suffers from the carbuncle instability if the standard HLLC approximate Riemann solver is employed. To avoid this issue, the small HLLC-LM modification from Fleischmann et al. [24] is used in this test.
The test is shown in Figure 5, 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.

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 γ = 1.4 ideal gas with ρ = 1.4, p = 1, u = 0. The domain for this test is x = [0 : 5], y = [−2 : 2], with a resolution of 400 × 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 v r = 3x. The test is run until t = 1 with a CFL of 0.4.
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, and the test demonstrates the ability of the method to handle supersonic dynamic rigid body motion.

Solid Extrusion with a Piston
This test, proposed in Hill et al. [31], 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 presence of the elastoplastic solid in this case requires the use of the full multi-phase system of equations. 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 (22) with parameters given in Table 1. Following Hill et al. [31], the aluminium has an ideal plasticity law with the yield stress given by σ 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 × 216 is used, with 3 levels of AMR, refinement factor 2. The test is run for 20 µs with a CFL of 0.4. The piston is initially located at x = −20 mm and given a constant velocity of v r = 100x 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.

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. [28]. 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. [28] and other works in the literature [13,33,34,42,47]. 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   [17] and Zółtak and Drikakis [64], and the experimental studies of Bryson and Gross [15] and Kaca [35]. Each axis is non-dimensionalised with respect to the diameter of the cylinder, D.  et al. [57]. 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 9a. The copper is given an initial downward velocity of 400 ms −1 and is governed by the equation of state equation (22), 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 × 696 is used, with two layers of AMR, each of refinement factor 2. The test is run until t = 60 µs using a CFL of 0.4, and is shown in Figure 10.
The results agree well with the experimental results [28], shown in Figure 9b, 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. [47] and the range 6-9 quoted by Gray et al. [28]. 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.      [57]. The times are shown in µs, and the images depict the y-velocity and plastic strain. The results agree well with the experimental studies of Gray et al. [28].

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. [53], Deoliveira et al. [20] and Tarver [54] 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 multiphysics aspects to this test. Firstly, the explosive LX-17 is modelled using the multi-phase formulation detailed by Wallis et al. [58], 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 λ is added to the full multi-material system of equations in order to track the reaction: The explosive's products and reactants are governed by the JWL equation-of-state, with parameters given in Table 3. This approach further allows for the use of the ignition-and-growth model, as employed by Deoliveira et al. [20], 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: where F = 1−λ is the reacted fraction and H is the Heaviside function. Other parameters are material dependent constants, given in Table 4. Full details of the method are given by Wallis et al. [58]. Secondly, an explosive desensitisation model is also implemented. The model was used by Deoliveira et al. [20] 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 ϕ which tracks the desensitisation of the explosive: where A = 8.36 × 10 −3 is a material dependent constant, p is the pressure and = 10 −3 is a small value designed to increase ϕ 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: 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 ,min = λ c ϕ, with λ c = 0.01 being another calibration parameter. This means the reaction rate is now: 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 (22) with the parameters given in Table 1, obeying a ideal plasticity law with a yield stress of σ 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 = 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 11 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. [20].

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 [29], 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 [51] 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. [8]. The volume fraction is then set to be the rigid body volume fraction.
The test is non-dimensionalised, with a Mach 1. 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 with complex three-dimensional geometries.

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.