A multi-physics methodology for the simulation of the two-way interaction of reactive flow and elastoplastic structural response

We propose a numerical methodology for the numerical simulation of distinct, interacting physical processes described by a combination of compressible, inert and reactive forms of the Euler equations, multiphase equations and elastoplastic equations. These systems of equations are usually solved by coupling finite element and CFD models. Here we solve them simultaneously, by recasting all the equations in the same, hyperbolic form and solving them on the same grid with the same finite-volume numerical schemes. The proposed compressible, multiphase, hydrodynamic formulation can employ a hierarchy of five reactive and non-reactive flow models, which allows simple to more involved applications to be directly described by the appropriate selection. The communication between the hydrodynamic and elastoplastic systems is facilitated by means of mixed-material Riemann solvers at the boundaries of the systems, which represent physical material boundaries. To this end we derive approximate mixed Riemann solvers for each pair of the above models based on characteristic equations. The components for reactive flow and elastoplastic solid modelling are validated separately before presenting validation for the full, coupled systems. Multi-dimensional use cases demonstrate the suitability of the reactive flow-solid interaction methodology in the context of impact-driven initiation of reactive flow and structural response due to violent reaction in automotive (e.g. car crash) or defence (e.g. explosive reactive armour) applications. Several types of explosives (C4, Deetasheet, nitromethane, gaseous fuel) in gaseous, liquid and solid state are considered.


Introduction
The accurate numerical simulation of a wide range of industrial, automotive, aerospace and defence processes necessitates the consideration of gaseous or condensed-phase explosives initiated by the impact by and their interaction with fluid or elastoplastic solid materials. Examples include accidental fuel initiation in a car crash and fuel tank containment in the context of safety studies and explosive reactive armour (ERA). This article is concerned with the development of numerical methods for the simultaneous solution of multiphase, reactive, inert fluid and elastoplastic solid equations suitable for the numerical modelling of such problems. The methodology can also be used as part of the manufacturing process for optimising car (or other device) compartments in terms of shapes and materials.
An integrated numerical methodology for this kind of simulation has three elements; the formulation describing the elastoplastic solid (impactor and/or confiner), the formulation describing the gaseous or condensed phase fuel or explosive, including the explicit capturing of the reaction zone and the transition from reactants to products and the communication between the two systems through material coupling. By explicit capturing of the reaction zone we mean that the reactants and products are described by distinct equations of state and transition between them occurs in a numerically mixed zone leading to the generation of the reaction front and detonation wave. This is in contrast to programmed burn approaches where assumptions are made with regards to the location of the front and the energy deposed behind it. an inert component with the reactant-product mixture, through a diffused interface approach. This would allow, for example, for the inclusion of an air gap between the explosive stick and the metal confiner. Reduced versions of this formulation are also presented to model cases when the inert component is not present, or when the explicit modelling of the products of reaction is not required or even when the phases are all non-reacting and could form free-surfaces.
There are various approaches for the numerical solution of coupled solid-fluid problems. Traditionally, Lagrangian techniques are followed. These, however, present difficulties for large deformations of the materials as these are inherently translated to large deformations in the underlying mesh. Others include Smoothed Particle Hydrodynamics (SPH) and arbitrary Lagrangian-Eulerian (ALE). Also, traditionally, the modelling of solid-fluid problems is using one-way coupled finite element codes for modelling the solid part and CFD codes to model the fluid part. As a result, the two processes are solved in a 'co-simulation' environment, each on it's own grid with a distinct numerical method. This may lead to discretisation errors passing from one method to the other. Even though each class has its merits and shortcomings, in this work, we retain a regular, Cartesian mesh (and data) structure, as this is relatively easy to implement in our existing structured, hierarchically adaptive mesh refinement (AMR) framework.
Studies including solid-reactive coupling in the Eulerian frame include work by Miller and Colella [23] and Barton et al. [24]. The solid material in these studies is described by a full elastoplastic system, while the explosive system considered incorporates either a single-phase Euler formulation or program burn. Although these are suitable for gaseous combustion, they are not adequate for more complex or condensed phase explosives. Schoch et al. [25] couple the full elastoplastic system to a two-phase, five-equation model for condensed explosives, which is more complex and as a result more restrictive than the model used in this work. Examples of inert Euler-solid coupling can also be found, as for example in [26,27].
The coupling in multi-material simulations (including fluid-fluid, solid-solid or fluid-solid interaction) can follow inherently from the Lagrangian framework, as for example by Howell and Ball [28] who apply the coupling in solid-solid applications. In the Eulerian framework, interface fitting approaches can be followed, such as Volume of Fluid (VOF) and Ghost Fluid Methods (GFM) as in the work by Barton et al. [29] who use the modified GFM and Schoch et al. [25] who use the 'real GFM'. We define as multi-material the framework where material interfaces are tracked rather than captured. This includes different sets of equations being solved on either side of the material interface and special methods used to apply 'boundary conditions' across the material interface. In this work we consider the Riemann ghost fluid method as presented by Sambasivan and Udaykumar [30], based on the pioneering work by Fedkiw et al. [31] and we provide the coupling by deriving mixed-material Riemann solvers to be used at material interfaces. To this end, we derive characteristic equations that lead to the formulation of a linearised mixed Riemann solver, applied to the one side of the material interface for the full hydrodynamic formulation by Michael and Nikiforakis [22] and its reduced models. Depending on the application or the fuel/explosive in consideration, the full or the reduced versions of the model is employed.
In the remainder of this article, we first present the distinct mathematical formulations describing the explosive (including the reduced versions) and the elastoplastic solid. The coupling technique is then presented and derivations of the mixed Riemann solvers for fluid-fluid and fluid-solid coupling are presented for each explosives model considered. Then, we validate separately the solid and explosives components by invoking solid-only and explosive-only test problems with known solutions. The explosive-solid coupling is first tested in one dimension and then multi-material inert and reactive simulations are considered, illustrating the applicability of the coupling in impact-initiation of condensed-phase explosives in ERA examples and gaseous fuel in car-crash examples for fuel containment.

Mathematical models
In this section, the distinct mathematical formulations describing the materials involved in solid-fluid impact and interaction applications are presented. To aid the description of the principal components and without loss of generality, we use as a reference configuration. As the accidental impact of a fuel tank in a car crash scenario involves a complex tank geometry, we chose that for explaining the mathematical formulations and numerical methods to use the simpler sandwiched-plate impact in an ERA configuration as illustrated in Fig. 1. An explosive (yellow) is residing between two steel plates (grey) and a steel projectile (grey) impacts this sandwiched configuration from below. The whole system could be surrounded by air or vacuum making this a solid-explosive-fluid or solid-explosive-vacuum configuration. The choice of the explosives model depends on the physical complexity of the material (e.g. porous, homogeneous etc) and the dominant timescales in the mechanical and thermodynamic properties of the scenario considered. The explosives model could therefore range from the simple single-EoS augmented Euler formulation (e.g. suitable for gaseous combustion) to the full Baer-Nunziato formulation (e.g. suitable for porous materials), with the possibility of using a model of intermediate complexity like the ones proposed by Banks et al. [10], Kapila et al. [18] and Saurel et al. [17]. In this work, the formulation proposed by Michael and Nikiforakis [22] (henceforth referred to as MiNi16) is used for the mathematical description of the explosive. This formulation additionally captures the air gap that can be included inbetween the explosive and the front plate (zoom in Figure 1). The steel plates and the steel projectile are described mathematically by an elastoplastic model in the Eulerian frame, as presented by Barton et al. [5]. The plates and projectile are treated numerically as separate materials in the multi-material framework. This would allow easily for the physical material of one or both plates as well as of the projectile to be changed. The surrounding material can be vacuum or air. If it is air, this is described by the compressible inert Euler equations.

The explosives model (Case 1)
In this section, we summarise the MiNi16 formulation [22] based on the Case 1 example illustrated in Figure 1 (zoom) of a sandwiched-plate impact scenario with an air gap. In this configuration, the explosive reactant is described as phase α, with density, velocity vector and pressure (ρ α , u α , p α ) and the products of reaction as phase β with (ρ β , u β , p β ) equivalently. These are assumed to form a homogeneous mixture, which we will hereafter call interchangeably the 'explosive mixture' or the 'reactant-products mixture' and denote this as phase 2 with mixture density, velocity vector and pressure (ρ 2 , u 2 , p 2 ). We can then denote by λ the mass fraction of the reactants, such that λ = 1 denotes fully unburnt material and λ = 0 denotes fully burnt material. The air gap (or any other inert material that could confine the explosive mixture) is denoted as phase 1 with density, velocity vector and pressure (ρ 1 , u 1 , p 1 ).
We denote by z a colour function, which can be considered to be the volume fraction of the air with respect to the volume of the total mixture of phases 1 and 2, with density ρ. Equivalently, the volume fraction of the reactant-product mixture with respect to the volume of the total mixture is given by 1 − z. For convenience, we denote z by z 1 and 1 − z by z 2 . Then, the closure condition z 1 + z 2 = 1 holds.
Velocity and pressure equilibrium applies between the all phases, such that u α = u β = u 1 = u 2 = u and p α = p β = p 1 = p 2 = p. Temperature equilibrium is only assumed between the phases of the explosive mixture, i.e. T α = T β , although other closure conditions can be found to be more suitable for other applications (see for example Stewart et al. [32]).
Then, the MiNi16 system is described as in [22] by: where u = (u, v, w) denotes the total vector velocity, i denotes space dimension, i = 1, 2, 3, ρ the total density of the system and E the specific total energy given by E = e + 1 2 i u 2 i , with e the total specific internal energy of the system.
K is a function giving the rate of conversion of reactants to products and for numerical purposes, it is considered a source term to the hyperbolic part of the system. In this work, depending on the reaction rate law form, term K usually depends on the temperature or pressure of the phases, as well as the total density, the density of phase 2 and λ. The model is not restrictive with the reaction rate law, therefore other types of reaction rates can also be used.
In this work, all fluid components described by the MiNi16 model are assumed to be governed by a Mie-Grüneisen equation of state, of the form: Material interfaces between the three phases are described by a diffused interface technique. Hence, mixture rules need to be defined for the diffusion zone, relating the thermodynamic properties of the mixture with those of the individual phases. The mixture rules for the specific internal energy, density and adiabatic index (γ) are: where e 1 , e α , and e β denote the specific internal energies of phases 1, α and β, ξ = 1/(γ − 1), and C vα and C v β denote the specific heat at constant volume for phases α and β.
The model in this form can solve multi-component problems involving two miscible phases (phase α and phase β) forming a mixture represented as phase 2 and one inert, immiscible component (phase 1 ). Only one of the phases α and β can be reactive.
The soundspeed also follows a mixture rule given as: where c i are the individual soundspeeds of phases 1 and 2 and c 2 depends on averaging procedures of energy and density derivatives of phases α and β. For more information on this as well as for the numerical evaluation of the total equation of state the reader is referred to [22].

Reduced MiNi16 reactive model (Case 2)
Suppose that the explicit modelling of the combustion products is not necessary and the fluid system comprises of the reacting explosive and an inert phase, as illustrated in Case 2 of Fig. 1. The full model then reduces to a two-phase, reacting, five-equation interface model as given by Michael and Nikiforakis [33]. In this case, phase β is dropped and the system comprises of phase 1 and phase 2 = phase α (effectively quantities ( ) 2 = ( ) α ).

Reduced MiNi16 inert model (Case 3)
Suppose that phase 2 does not comprise a homogeneous mixture of reactant and products, but instead it represents a single, inert constituent (i.e. λ = 0 or λ = 1 and K = 0 everywhere), as illustrated in Case 3 of Fig. 1. Then, the full model reduces to the inert five-equation interface model, as given by Allaire et al. [34] (effectively quantities ( ) 2 = ( ) α = ( ) β and equation (6) becomes inactive): The mixture rules of the full model also reduce to the ones given by Allaire et al. [34].

Augmented reactive Euler model (Case 4)
Suppose that the reactant-product mixture is not confined by an inert material (in the limit of z 1 → 0) and it is considered its own entity, as illustrated in Case 4 of Fig. 1. The mixture rules then reduce to the mixture rules found in [10] and the system reduces to the two-component fluid-mixture model (or two-phase, in the context of this work), as described in the same work:

The elastoplastic model
In this work, we use the elastic solid model described by Barton et al. [5] based on the formulation by Godunov and Romenskii [35]. The plasticity is included following the work of Miller and Colella [23].
Consider the steel plate of Fig. 1 in isolation. In an Eulerian frame, which we employ here, there is no mesh distortion that can be used to describe the solid material deformation. Thus the material distortion needs to be accounted for in a different way. Here, this is done by defining the elastic deformation gradient as: which maps the coordinate X in the initial configuration to the coordinate x in the deformed configuration.
The state of the solid is characterised by the elastic deformation gradient, velocity u i and entropy S. Following the work by Barton et al. [5], the complete three-dimensional system forms a hyperbolic system of conservation laws for momentum, strain and energy: with the vector components · i and tensor components · ij . The first two equations along with the densitydeformation gradient relation: where ρ 0 is the density of the initial unstressed medium, essentially evolve the solid material hydrodynamically. Here, σ is the stress, E the total energy such that E = 1 2 |u| 2 + e, with e the specific internal energy and κ the scalar material history that tracks the work hardening of the material through plastic deformation. We denote the source terms associated with the plastic update as P ij .
The system is closed by an analytic constitutive model relating the specific internal energy to the deformation gradient, entropy and material history parameter (if applicable): The stress tensor is given by: For F e to represent a physical deformations, the equations for deformation gradient satisfy three compatibility constraints: ∂ρF e kj ∂x k = 0, j = 1, 2, 3, which hold true for t > 0 if true for initial data. This is based on the fact that F e is defined as a gradient. The deformation is purely elastic until the physical state is evolved beyond the yield surface (f > 0), which in this work is taken to be: where σ Y is the yield stress and the matrix norm ||.|| the Shur norm (||σ|| 2 = tr(σ T σ)).
As this identifies the maximum yield allowed to be reached by an elastic-only step, a predictor-corrector method is followed to re-map the solid state onto the yield surface. Assuming that the simulation timestep is small, this is taken to be a straight line, using the associative flow rate (˙ p = η ∂F ∂σ ), satisfying the maximum plastic dissipation principle (i.e. the steepest path). In general, this is re-mapping procedure is governed by the dissipation law: where Σ = GσF and : is the double contraction of tensors (e.g. σ : σ = tr(σ T σ)). The initial prediction is F = F e and F p = I, where F is the specific total deformation tensor and F p the plastic deformation tensor that contains the contribution from plastic deformation. This is then relaxed to the yield surface according to the procedure of Miller and Colella [23]. The explosive and solid mathematical formulations described in this section are solved numerically using high-resolution, shock-capturing, Riemann-problem based methods and structured, hierarchical adaptive mesh refinement, as described in previous work [25,21,22,33,36,37].

The multi-material approach
Ghost fluid methods, in combination with level set methods, provide a robust and efficient technique for tracking interfaces and boundaries, such as the material interfaces between solid and fluid materials. In this work, we use level set methods to track the solid-explosive 1 interface. Each component, e.g. the solid or the explosive is called a material in this framework. Such methods only give the location of the interface; they do not affect the evolution of the material components. The behaviour of the material components at the interface is modelled by the implementation of dynamical boundary conditions with the aid of the Riemann ghost fluid method and the devise of mixed-material Riemann solvers to solve the interfacial Riemann problems between materials.
A signed distance function φ(x, y), called the level set function is used, with the zero contour given by Φ = (x, y)|φ(x, y) = 0. The sign of φ i , where i is the cell index, determines which material is present in that cell. The evolution of φ(x, y), assuming no mass transfer and continuous velocity across the interface, is given by the advection equation: The Riemann ghost fluid method by Sambasivan and Udaykumar [30,38] is utilised to model the behaviour of the material component at the interface. This method, in contrast to the original ghost fluid method [31] and the 'modified' ghost fluid method [39] uses a Riemann solver to predict ghost-cell states adjacent to the interface. For every cell i adjacent to the interface the following procedure is used: 1. locate the interface within the cell at the point P = i + φ∇φ 2. project two probes into the adjacent materials, reaching the points P 1 = P + n · ∆x and P 2 = P − n · ∆x 3. interpolate states at each point using information from the surrounding cells 4. solve a mixed Riemann problem (as described in Sec. 3.1) between the two states to extract the state of the real-material cells, adjacent to the interface (left star state W * L , in Fig. 2) 5. replace the state in cell i by the computed star state.
After the above procedure is followed for each material, a fast-marching method is used to fill in the ghost cells for each material.

Mixed Riemann solvers
In this section we describe how the mixed-material Riemann problem at material interfaces is solved (step 4 in the procedure described above). As we are considering five hydrodynamic models and one elastoplastic model, there are a lot of ways of combining them. In these section, we derive mixed Riemann solvers for each pair that can be encountered, although some extensions of the different combinations are straightforward. The Riemann solver at the material interface takes as input two states from the two different materials that are modelled by different mathematical models, providing one-sided estimations of the interface-adjacent (star) states. These estimations are based on the characteristic equations deduced from the mathematical system describing the left material and by invoking appropriate 'boundary conditions' between the two materials at the interface. That is why for each model described in Sec. 2 to be coupled with another fluid or solid model, a new mixed Riemann solver has to be derived. In this section, we first describe how the full MiNi16 model and the reduced versions of it are coupled with a simple Euler system and then with a full elastoplastic solid system (including how the elastoplastic system is coupled with the Euler equations). The remainder combinations should be directly deductible from these.

Reactive MiNi16 model coupled with the Euler model
Consider a cell which contains a material boundary. Without loss of generality, we assume that on the left side of the interface lies the material governed by the MiNi16 model of Sec. 2.1 (hereby called the MiNi16 material) and on the right side of the interface lies a material governed by the Euler equations. Hereforth we assume that we are currently solving for the MiNi16 material (in GFM terminology, the 'real' material). At the material boundary a Riemann problem is solved between the left MiNi16 and the right Euler system to provide the star state for the real material. We use a Riemann solver that takes into account the two different materials and all the wave patterns in the MiNi16 system, as described in this section.
We write the MiNi16 model described in Sec. 2.1 by equations (1)-(6) in primitive form as: where is the mass fraction 2 of material k, k = 1, 2, with respect to the mixture of phases 1 and 2 (i.e. the total fluid).
As we are solving for the MiNi16 model, we only look for the left star state, W * L (see Fig. 2). To obtain star values on the left of the interface, we use characteristic equations.
Characteristics define directions dx dt = µ j , in which So, along dx dt = µ 1 = u: β 1 (0, 0, 0, 0, 0, 1) · (dρ, dY, du, dw, dp, dz) giving, dλ = 0 (39) Figure 2: The Riemann problem at the material interface between any combination of fluid hydrodynamic and solid elastoplastic materials. The shaded region represents the initial ghost region and the white region the initial 'real' material region in GFM terminology. We are solving for the material in the white, left region. If this is a fluid, the one-sided solution consists one or more degenerate and hence overlapping waves (left) and if it is a solid, the solution can consist with more than one non-overlapping waves (right).
Applying the same along dx dt = µ j = u for j = 2, 3, 4, 5 we obtain, respectively: dp − c 2 dρ = 0, Finally, along dx dt = µ 6 = u − c and dx dt = µ 7 = u + c, we obtain: dp − ρc du = 0 (44) dp + ρc du = 0 For a single phase described by the ideal gas or stiffened gas equation of state, the characteristic equations can be integrated directly, as the expressions for the sound speed are simple. However, for more complex equations of state or MiNi16 formulations where the sound speed is the 'mixture' of the individual phase sound speeds (see equation (11)) this might not be possible. In such case, one can obtain an approximate mixed Riemann solver by replacing the differentials with the difference of the initial and the final state without integrating (i.e. across the characteristics), as presented below.
Using dp + ρc du = 0, we connect the states W * L and W L to obtain: Using dp − ρc du = 0, we connect the states W * R and W R to obtain: Pressure and velocity don't change across the material interface, hence p * L = p * R = p * and u * L = u * R = u * . Applying these conditions to (46) and (47) we obtain two expressions for p * and u * : Solving the above two simultaneously, we obtain an expression for the pressure in the star region: where To calculate the left MiNi16 fluid state, connect states W * L and W L , using equation (46) and dp−c 2 dρ = 0 to obtain: Using the remaining characteristic equations we obtain: Using the above and the definition for and so: Equations (49) If the real fluid is described by the reduced two-phase reactive system given in Sec. 2.1.1 then the same procedure as here is used, as the same governing equations (and hence characteristic equations) hold.

Iterative extension
To improve accuracy of all of the explosive-solid mixed Riemann solvers, in the expressions for p * , u * and ρ * , the expressions for c 2 L and ρ L c L (denoted now as c L 2 and ρ L c L ) can approximated not just by the ( ) L state but by an average between the original and predicted star states, i.e.
. This can be considered as an average between the material interior state and the interface state. For other options for constructing the ( ) state see Schoch et al. [25].
This, however generates an implicit problem; the unknown interface state now depends on itself, i.e. W = W(W L , W * L ). A predictor-corrector method is used to iterate through and repeatedly update the star states until convergence (based on p * ). The initial guesses for the iteration process are the primitive states of the two materials.

Reduced MiNi16 inert model coupled with the Euler model
Suppose now that the material on the left side of the interface is governed by the MiNi16 inert model of Sec. 2.1.2 (the 'real material), using two phases only for demonstration purposes, while the material on the right side is governed by the Euler equations.

Augmented Euler reactive model coupled with the Euler model
Suppose now that the material on the left side of the interface is governed by the augmented Euler reactive model of Sec. 2.1.3, while the material on the right side is governed by the Euler equations.
We write the reactive fluid model in primitive form as given by equation (33) with The Jacobian matrix A(W) has eigenvalues The right eigenvectors are and the left eigenvectors are l 1 = −c 2 , 0, 0, 1, 0 , l 2 = (0, 0, 1, 0, 0) , l 3 = (0, 0, 0, 0, 1) , Thus, the new system is directly reduced from the full MiNi16 system and the same procedure can be used to derive characteristic equations and connect states across characteristics to obtain left star values. As the reduced system does not contain a z−equation, the characteristic equations are given by equations (39)-(41) and (44)

Solving for the Euler system
When we are solving for the Euler system (i.e. the Euler fluid is the 'real' fluid), irrespective of which system is on the right side of the interface, we repeat the above process and only compute the quantities (p * , u * , w * and ρ * L ).
3.1.6. MiNi16 model coupled with the elastoplastic solid model Consider a cell which contains a material boundary. Without loss of generality, we assume that on the left side of the interface lies the material governed by the MiNi16 equations (the muti-phase material) and on the right side of the interface lies a material governed by the elastoplastic solid equations (the solid material). We develop a Riemann solver that takes into account the two different materials to determine the star state in the fluid material.
We follow a similar procedure to that in Sec. 3.1.1. Referring to Fig. 2, W L corresponds to the original MiNi16 state, W R to the original elastoplastic state and W * L to the MiNi16 star state that we are looking to compute in this Riemann solver. Since we are solving for the fluid as the real material, the Riemann problem still has three types of waves (two non-linear and four overlapping linear). The same characteristic relations (39)- (45) are defined as before and we use the approach of representing the differentials with the state difference. Connecting fluid states W * L and W L using dp + ρc du = 0 and solid states W * R and W R using dp − ρc du = 0 we obtain a mixed-material expression for p * : where Q is an orthogonal matrix and D is the diagonal matrix of positive eigenvalues for the solid system, such that the acoustic tensor is defined by: Considering (49). Then, using dp + c 2 dρ = 0 to connect states W * L and W L we obtain equations (50)-(51) and values for u * L and ρ * L , using the remaining characteristic equations we obtain equations (52)-(55) and values for w * , λ * , Y k * and z * and using the definition of Y k we obtain ρ * k via equations (56)-(57). At the interface, we apply conditions To improve accuracy, an iterative approach as described in Sec. 3.1.2 is used. If the real fluid is described by the reduced two-phase reactive system given in Sec. 2.1.1 then the same procedure as here is used, as the same governing equations (and hence characteristic equations) hold.

Reduced MiNi16 inert model coupled with the elastoplastic solid model
Suppose now that on the left side of the interface lies the material governed by the MiNi16 inert equations of Sec. 2.1.2 (here we use two phases only here for convenience) and on the right side of the interface still lies a material governed by the elastoplastic solid equations. We follow a similar procedure to that in Sec. 3.1.6 to determine the star state (W * L ) in the fluid material. The characteristic relations for this system are defined by equations (40)- (45). By connecting the appropriate states we obtain the mixed material expression for p * given by equation (64) and the remaining star states by equations (50)-(52) and (54)-(57). The difference now is that the ( ) L quantities come from the two-phase inert model rather than the full MiNi16 model. At the interface, we apply the conditions given by (66) as before.
As with the MiNi16-elastoplastic coupling, the expressions for c F 2 and ρ F c F can be taken to be averages of the original and predicted star states and hence a predictor-corrector method is required to update the star states until convergence (based on p * ).

Augmented reactive Euler model coupled with the elastoplastic solid model
If now that on the left side of the interface lies the material governed by the augmented reactive Euler equations of Sec. 2.1.3 and on the right side of the interface still lies a material governed by the elastoplastic solid equations. We follow a similar procedure to that in Sec. 3.1.6 to determine the star state (W * L ) in the fluid material.
The characteristic relations for this system are defined by equations (39)- (41) and (44)- (45). By connecting the appropriate states we obtain the mixed material expression for p * given by equation (64) and the remaining star states by equations (50)-(53) and (56)-(57). The difference now is that the ( ) L quantities come from the reactive augmented Euler model rather than the full MiNi16 model. At the interface, we apply the conditions given by (66) as before.
As with the MiNi16-elastoplastic coupling, the expressions for c F 2 and ρ F c F can be taken to be averages of the original and predicted star states and hence a predictor-corrector method is required to update the star states until convergence (based on p * ).

Solving for the elastoplastic solid model
Suppose we again have a cell interface which is a material boundary but on the left side of the interface lies an elastoplastic solid material. On the right side of the interface lies a fluid material which can be a MiNi16 material or any of its reduced versions (including the unaugmented Euler equations). We follow a similar procedure as with the explosive-fluid mixed Riemann solvers. We give a brief outline on the derivation of the solid approach and refer the reader to Barton et al. [24] for more details.
The aim here is to determine the star state for the solid cell adjacent to the solid-fluid interface. The elastoplastic solid model can be written in primitive, quasilinear form for primitive variables W = (u, F e T , S): where where E ij represents the unit dyads E ij = e i × e T j , I is the identity matrix and The eigenvalues and corresponding right and left eigenvectors are computed and used to obtain characteristic relations, allowing for the computation of the primitive elastic star state (the plastic deformation gradient F p is not computed at this point): wherer are the right eigenvectors of the elastic system. At the material interface, slip conditions are used which specify the continuity of the normal components of velocity and traction and zero tangential stresses in the solid: A more accurate star state can be determined if an iterative procedure is used. Starting with initial guesses for the left and right star states the primitive solid and the primitive fluid states, the star state computation described above is iterated until convergence of p * .

Validation of the elastoplastic component
In this section, the implementation of the solid-solid formulation is validated in isolation, without the influence of the fluid formulation or solid-fluid mixed Riemann solvers.

Aluminium plate impacting an aluminium target
In this test, a projectile block of aluminium moving at 400m/s impacts a stationary target block of aluminium. This was originally presented by Howell and Ball [28].  [8], consisting of a two-term hydrodynamic component and a one-term shear deformation component. The full equation is given by: where I 1 , I 2 , I 3 are the invariants of the Finger tensor, K 0 = c 2 0 − 4 3 b 2 0 is the squared bulk speed of sound, B 0 is the reference shear wave speed, c v is the specific heat capacity at constant volume, T 0 is the reference temperature, α, γ are exponents determining the the non-linear dependence of the sound speed and temperature on density respectively and β an exponent determining the non-linear dependence of this shear wave speed on density.
The constitutive model parameters for the aluminium considered here are given in Table 1. Perfect plasticity is assumed, with a yield stress of 0.4 GPa. Between the two aluminium plates, a slip condition is assumed. The simulation is performed at an effective resolution of ∆x = ∆y = 50 µm up to a final time of 5 µs.   Fig. 3 illustrates the computed wave structure in the two plates at times t = 0.5, 1, 3 and 5 µs. Upon impact, shock waves are generated that travel upwards into the projectile and downwards into the target, as seen in Fig. 3(a). These waves are of the same strength since the impact is symmetric (projectile and target plate are made of the same material). In Fig. 3(b), the shock wave travelling in the target plate is seen to have split into an elastic precursor and a trailing plastic wave. The shock wave travelling into the projectile reaches the rear end of the plate, where it interacts with the solid/vacuum interface. This interaction generates a release wave travelling backwards into the projectile. It then crosses into the target plate generating a region of high tensile stress in the x-direction. This weakens the plastic wave traversing the target (Fig. 3(c)). By t = 4 µs, the elastic wave reaches the rear of the target plate, generating a downwards-moving release wave. Subsequently, this release wave interacts with the still rightward-travelling plastic wave, producing a new set of waves that travel upwards and downwards and continue reflecting on the rear end of the target plate. The deformation of the two plates is also apparent in the sequence of Fig. 3.

Hyperelastic and
An excellent match is seen between our results and the computation by Howell and Ball [28]. To compare further with their results we consider gauges embedded in the originally stationary target block to detect strain, velocity, pressure, and density. These gauges are allowed to move with the flow as the target block deforms. Five equally spaced gauges are placed along the centreline of the target plate. The first one is placed at 1.8125 mm from the original impact position and a distance of 3.625 mm is allowed between consecutive gauges. The time histories for each gauge for the x-wise velocity, pressure, density and x-wise stress are seen in Fig. 4. The arrival times of the waves and their amplitude compare well with the results of Howell and Ball [28], for all gauges. The split of the shock wave generated at impact in traversing the target is clearly seen and the two waves appear to be steeper in our results. We choose to run the simulation longer than Howell and Ball, to capture the reflection of the elastic wave at the end of the target plate and the interaction of the reflected wave with the travelling elastic wave. These phenomena are also seen in the time-histories.

Validation of the reactive, hydrodynamic component
In this section, the implementation of the MiNi16 formulation is validated in isolation, without the influence of the elastoplastic formulation or elastoplastic-hydrodynamic (or simply also referred to as solidfluid) mixed Riemann solvers. The tests are based on C4, as this is the explosive used in the later solid-fluid -0.09 e 0.667 Table 2: Scaled JWL parameters for C4 and its detonation products (left) and scaled Ignition and Growth parameters for C4 (right) [42]. configurations.

Equation of state and reaction rate
In this work, C4 and its products are modelled by the JWL equation of state, with parameters as given in Table 2. The JWL equation of state is a Mie-Grüneisen equation of state, with general form given by with reference pressure and energy curves given by: and Grüneisen coefficient given as: The JWL equation of state is usually used to model reaction products but it has been extensively used to model the unreacted phase of explosives as well [40,41]. To describe the conversion of reactants to products we use the the Ignition and Growth model developed by Lee and Tarver [44] and given by: where φ = 1 − λ is the mass fraction of the products, H is the Heaviside function and I, G 1 , G 2 , a, b, c, d, e, g, x, y and z are constants chosen for a particular explosive and a specific regime of the detonation process (i.e. initiation or propagation of established detonation). The constants φ IGmax , φ G1max and φ G2max determine for how long each of the three terms is dominant. The form of this reaction rate was constructed so that each term represents one of the three stages of reaction observed during the shock initiation and detonation of pressed solid explosives. The three terms can be interpreted differently depending which regime of the detonation process is studied. For shock initiation modelling, the first term represents the amount of reaction due to the formation and ignition of hot-spots which are generated by several mechanisms in heterogeneous explosives but mainly by void compression. The second term then describes the reaction due to the growth of the hot-spots and the third term the completion of reaction and transition to detonation. For detonation modelling, the first term still describes the amount of reaction at the initiation stage, due to the generation of hot spots. The second term describes the fast growth of the reaction as the reactant is converted to products. The third term describes the relatively slow diffusion-limited process of carbon formation [40]. The reaction rate law parameters for C4 are given in Table  2 which have been rendered non-dimensional using the CJ state for the explosive.

C4 ZND
Consider a one-dimensional slab of C4 at ambient conditions, initiated by a booster with pressure of 30 GPa. We model this computationally in a domain that spans [0,6.432]cm and the booster resides in the region [0,0.0402]cm. An effective resolution of ∆x = 33.5 µm is used. The initial data in the ambient region are: (ρ reactant , ρ product , u, p, λ) = (1590kg m −3 , 1590kg m −3 , 0m s −1 , 1 × 10 5 Pa, 1) and in the booster region are: The explosive transits very quickly to detonation, which in-turn settles down to steady state. The steady detonation structure is shown in Fig. 6a.
The computed CJ and von Neumann pressure values agree with the values predicted in the previous section and fall within the range of values found in the literature.

C4 Pop-plot
In this section we validate the the explosives model against Pop-plot data of run distance to detonation versus input pressure. This is equivalent to overtake time versus input pressure that is usually used as, for Pop-plot purposes, detonation is defined to be the overtake of precursor shock wave by the generated reactive wave.
In Fig. 6b we compare our numerical results for run distance to detonation to experimental data by Urtiew et al. [42] for different input pressures, where a very good match is observed.

Validation and evaluation of the hydrodynamic-elastoplastic coupling
The final step towards simulating the non-linear, two-way interaction of explosives and elastoplastic materials is the validation of the solid-explosive coupling for one-dimensional and cylidrically symmetric test cases. Thereafter, the coupled system is used to demonstrate example applications.

Stressed copper impacting quiescent PBX 9404
This test considers a stressed copper component for x ∈ [0, 0.005]m, impacting an initially quiescent, reacted PBX9404 gas for x ∈ [0.005, 1.0]m, as described in [24]. The PBX9404 is modelled by the ideal gas equation of state, with γ = 2.83 and the copper by the elastic Romenskii equation of state (72) with parameters as given in Table 1 The numerical and the exact solutions for density and tangential stress components σ xy and σ xz at t = 0.9 µs are given in Fig. 7, where good agreement for all quantities is demonstrated. A small density error at the interface can be seen, which is, however, visibly smaller than the error in [24].

Sandwich-plate impact: inert Detasheet confined by steel, impacted by steel
The next validation is in cylindrical symmetry for an elastoplastic-elastoplastic-hydrodynamic configuration and we compare our results against existing numerical solutions. The test considers a cylindrical flyer plate impact and the subsequent response of the explosive residing behind the target plate. Specifically, a steel projectile of 18 mm diameter and 50 mm length impacts a steel target plate of thickness 3.18 mm and diameter 100 mm. Behind the target plate sits a block of Detasheet explosive with the same diameter as the target plate. Two explosive thicknesses are considered; 6.35 mm and 3.18 mm. Another steel plate is considered to sit behind the explosive, with the same diameter and thickness as the front plate. The simulation is done in axial symmetry. Initially the projectile and target plate are separated by 1 mm. A schematic of the setup is shown in Fig. 8, residing in vacuum.
The Detasheet is modelled by the linear Hugoniot equation of state [45], which is of Mie-Grüneisen form (73), with reference functions given by: Figure 7: Comparison of our numerical solutions for total density and tangential xy and xz stress against the exact solution for the solid-fluid test problem described in Sec. 6.1, with an initially moving and stressed copper (green) and quiescent PBX9404 at ambient pressure (red). The material interface between the solid and fluid component is captured as a discontinuity, low numerical diffusion is generally observed and the solution is free from any spurious oscillations.
and Grüneisen coefficient given as: The hydrodynamic, elastic and plastic behaviour of the steel is described by a linear Hugoniot equation as well, in combination with a constant shear model with G = 95.4 × 10 9 Pa and the Johnson Cook plasticity model [46]. The Hugoniot parameters for steel are: The Johnson-Cook plasticity model assumes yield stress given by: with A the base yield stress, B the hardening constant multiplier, C the strain rate dependence constant multiplier, n the hardening exponent, m the temperature dependence exponent,ε 0 the reference strain rate, T 0 the reference temperature and T m the reference melt temperature. For this test, the Johnson-Cook parameters for steel are: The wave pattern in this scenario is complex so a wave diagram is included as Fig. 9 to keep track of the waves generated. The solid black lines denoted by 'S' represent the shock waves, the pairs of red lines denoted by 'R' the rarefaction fans and the black dotted lines the impulsively accelerated material interfaces. Upon impact, two shocks are generated; one travelling in the target plate (S 1 ) and one travelling back into the projectile (S 2 ). As a result of the impulse, the face of the plate is accelerated. The shock S 1 reaches the front plate/explosive interface, a high impedance/low impedance interface (HI-LI), where it is transmitted in the explosive as a shock (S 3 ) and reflected back into the front plate as a release wave (R 1 ). Again, the material interface is accelerated. The shock S 3 reaches the end of the explosive and encounters a low impedance-high impedance interface (LI-HI). Hence, it is transmitted into the rear plate as shock S 4 and reflected into the explosive as shock S 5 . It should be noted that if there was no rear plate, the wave diagram would stop after the generation of S 3 and no other waves would have affected the explosive. The shock S 5 reaches the rear end of the explosive, now a LI-HI interface, where two shocks are generated; S 6 , which is reflected back into the explosive and S 7 , which is transmitted into the front target plate. In the meantime, shock S 4 reaches the end of the rear plate and, since there is vacuum on the other side of the plate, the shock is reflected as a rarefaction wave (R 2 ) only, into the plate. The release wave R 2 reaches the rear plate/explosive interface which is now a HI-LI interface where it is reflected as a shock (S 8 ) back into the plate and transmitted as  a rarefaction wave (R 3 ) into the explosive. The release wave (R 3 ) interacts with the shock (S 6 ) within the explosive, leading to its weakening and making the wave pattern much more complex from here on. It is worth noting that if the rear plate was infinitely thick there would be no release waves coming into the explosive and a series of shock reflections would occur (like for S 5 and S 6 ), continuously increasing the pressure in the explosive. This is the general one-dimensional behaviour along the centreline of the experiment. Additional effects are generated due to the deformation of the solid materials and the weld assumption between the solid components. Pressure gauges are placed at the point of impact (station 1), at the front plate/explosive interface (station 2), at equally spaced points in the body of the explosive (stations [3][4][5] and at the explosive/rear plate interface (station 6). The pressure at each of these gauges over time is seen in Fig. 10 for both explosive thicknesses and we compare our results to the results by Lynch [45]. The station 1 curve (station 4 in [45]) shows the shock (S 1 ) as it is generated at the front plate/explosive interface, increasing the pressure at 40 GPa. The wave S 3 is seen traversing the explosive in the stations 2-5 (stations 6-10 in [45]), followed by release waves. In the meantime, the release wave R 1 in the front plate lowers considerably the pressure in the steel and the weld condition allows the pressure (or rather, stress) to go as low as −10 GPa. The wave S 5 is seen in stations 5 and 6. The release wave R 3 travels behind the shock S 5 within the explosive leading to its weakening.
A similar behaviour is seen for the explosive of thickness 3.18 mm, though the several features are seen to happen a lot faster (of course the gauges have also been moved). Moreover, the strength of the shock S 5 is greater than before and also the minimum stress achieved in the front plate is higher than before.
For both thicknesses our results agree well with the results by Lynch [45]. We note that in the results by Lynch [45] some oscillations are seen that are attributed to the numerical scheme used therein.

Sandwich-plate impact: reactive C4 confined by steel, impacted by steel
The previous test is repeated with a different explosive, namely C4, which is now chemically active. The explosive has been tested and validated individually in Sec. 5. We use the JWL equation of state as given by equations (73)-(76) and the ignition and growth reaction rate model (77) to represent the explosive. The parameters for these are found in Table 2.
The initial data for steel are as in the previous section, with the exception of the impact velocity which here is taken to be u z = 700 m s −1 . The initial data for the explosive are: C4: ρ = 1590 kg m −3 , u r = u z = 0 m s −1 , p = 10 5 Pa, λ = 1.0. (85) In Fig. 11 (top), the mass fraction of the reactants is seen on the left and the pressure distribution on the right, at times t = 2.4 and 4.9 µs. At t = 1.4 µs, minimal reaction, of the order of λ = 0.91, is observed due to the shock wave generated at impact (S 3 ). The effect of the shock S 5 on the reaction is larger, leading to λ = 0.65 at t = 2.4 µs at the ignition site near the explosive/rear plate interface and to λ = 0.55 at t = 4.9 µs in the same site.  The same test is repeated with the rear plate removed and the explosive extended to have thickness of 10.82 mm to demonstrate the effect of the rear plate on the ignition of the explosive. In Fig. 11 (bottom), the mass fraction of the reactants is seen on the left and the pressure distribution on the right, at times t = 2.4 and 4.9 µs. At t = 1.4 µs, minimal reaction, of the order of λ = 0.91, is observed due to the shock wave generated at impact (S 3 ). This is the same amount of reaction observed in the previous case as well, as up to this point, no waves have reached the rear plate. As mentioned earlier, if the rear plate is not present, only the wave S 3 affects the ignition process. At time t = 2.4 µs, reaction of the order λ = 0.85 is seen and at t = 4.9 µs reaction of the order λ = 0.83. The combination of the results of these two cases demonstrates the effect of the rear plate on accelerating the reaction.

Sandwich-plate impact with front air gap
The same test as in Sec. 6.3 is repeated with the inclusion of an air gap initially at atmospheric conditions (modelled as ideal gas with γ = 1.4) of width 3.18 mm between the front plate and the explosive. This demonstrates the full use of the MiNi16 model (air is phase 1, the explosive reactant is phase α and the explosive products is phase β; the last two form phase 2) coupled with the elastoplastic formulation. This problem is difficult to handle numerically, due to the big differences in the properties of the materials (reactants, products, air, solid) and the strong conditions involved. In this scenario we demonstrate how the presence of the air gap affects the ignition of the explosive so a direct comparison to the results of Sec. 6.3 is carried out. In both scenarios the impact generates a wave that travels in the front plate. In the air gap scenario, this reaches the air gap where a shock of the order of p = 0.5MPa is transmitted in the air. After a short propagation in the air, at t = 1.5µs this wave reaches the explosive-air interface where it is transmitted in the explosive as a shock of the order of p = 0.5MPa. The strength of this shock is considerably lower than the strength of the wave generated upon impact on the no-air gap case (p = 2.1GPa). As a result, this wave does not generate any reaction at all in the explosive. At t = 2.3µs the deformed front plate has squeezed away the air gap and reaches the lower explosive face. Upon impact a new shock wave is generated in the explosive (and an opposite one in the front target plate) of the order of p = 2.5GPa. The strength of this wave is comparable to the strength of the wave generated upon impact on the no-air gap case. Comparing the air-gap and non-air-gap cases, at 1µs after the first wave impacts the explosive, a reaction of the order of λ = 0.91 is observed in the non-air-gap case, whereas no reaction is seen in the air-gap scenario. At 2µs after the first wave impacts the explosive, a reaction of the order of λ = 0.65 is observed in the non-air-gap case, and of the order of λ = 0.71 in the air-gap scenario, demonstrating how the air gap hinders the reaction process. Figure 12: The mass fraction of the reactants (left) and the pressure distribution for the solid plates, the explosive and air (right) are presented at times t = 2.6 and 4.6 µs, for the steel-confined explosive with an air gap, impacted by steel problem of Sec. 6.4. Note that the figures have been stretched in the y-direction by a factor of 2, to allow for a clear view of the air gap.

Rod impact leading to detonation in a copper vessel
In this section, we present another showcase as a high-speed impact example involving multiphase reactive flow (explosive), elastoplastic structural (solid) response and inert air response. Consider a copper can of outer diameter of 14 cm, inner diameter 10 cm (resulting in a 2 cm wall thickness) and height of 49.2 cm. The can is filled with a nitromethane charge of diameter 10 cm and height of 45.2 cm (i.e. no air gaps between the charge and the confiner). The can is impacted by a copper projectile travelling at 2000 m s −1 with the same diameter as the explosive charge and a semi-infinite length, starting with length 6 cm but extending also out of the domain. The configuration resides in air, of outer diameter 18 cm and height 59.2 cm. All materials are initially at atmospheric conditions. This test is similar to the confined explosion tests by Miller and Colella [23], Barton et al. [24] and Schoch et al. [25]. However, these tests describe the initiation of the explosive with a booster whereas here we consider directly the impact resulting to the initiation and detonation of the explosive. Moreover, in the aforementioned studies the explosive was either simpler and described with simpler EoS and/or the entire configuration resided in vacuum.
In this work, we consider the explosive to be nitromethane, described by the Cochran-Chan EoS, which is of a Mie-Grüneisen form given by Eq. 73 with reference curves: reference energy given by Grüneisen coefficient Γ(ρ) = Γ 0 and parameters Γ 0 = 1.19, A = 0.819 GPa, B = 1.51 GPa, E 1 = 4.53, E 2 = 1.42, ρ 0 = 1134 kg m −3 , c v = 1714 J kg −1 K −1 and Q = 4.48 MJ kg −1 . The chemical reaction follows a singlestep, temperature dependent reaction rate law: with C = 6.9 × 10 10 s −1 and activation temperature T A = 11 350 K [47] and nitromethane temperature recovered as . We neglect the explicit description of products so the explosives model reduces to augmented Euler. The copper is described by the elastic Romenskii equation of state with parameters as given in Table 1 and perfect plasticity with yield of 70 MPa. The entire configuration is residing in air, described as an ideal gas with γ = 1.4, so we are able to visualise the transmission of the waves from the explosive, to the solid and finally to the gas. The simulation is performed with a base spatial resolution of ∆x = ∆y = 1 mm and two levels of refinement each with refinement factor 2, resulting in an effective resolution of ∆x = ∆y = 0.25 mm. Figure 13 shows the pressure distribution in the explosive, the solid and the surround gas, as well as the AMR grids, at times t = 0, 15, 35 and 60 µs. At the start of the simulation, the explosive is at ambient conditions, with pressure at 10 5 Pa. The impact sends a rightward-moving shock wave into the explosive and a leftward-moving wave in the projectile. Reflections of the waves in all materials also occur due to the confiner and the projectile not having the same width as the casing. The explosive is initiated and the reaction wave transits to a steady detonation. The detonation wave (red region in explosive) induces a shock wave in the solid (dark grey region in solid) which, upon reaching the copper-air boundary sends a shock wave in the air (dark blue region in air) and a release wave back into the solid. The repeated reflections at the material interfaces generate alternate regions of compression and tension in the solid and pressure waves in the explosive (often referred to as 'ringing'). The effect of the detonation on the deformation of the can is seen by comparing the shape of the can in the early and late stages of the event.

Rod impact leading to fuel ignition in a car fuel tank
In this section, we present a final showcase as a low-speed impact example involving reactive flow (fuel), elastoplastic structural (solid) response and inert air response. Consider a fuel tank as presented in Fig. 14a. The steel tank is filled with gaseous explosive and is impacted by a steel projectile travelling at 45 m s −1 , which can be considered to be the speed of a head-on car collision. The steel rod has diameter 3 cm and a semi-infinite length, starting with length 5.5 cm but extending also out of the domain. The configuration resides in air, in a domain with dimensions 23 cm × 23 cm with reflective right boundary condition and transmissive conditions at all other boundaries. The simulation is performed with a base spatial resolution of ∆x = ∆y = 1 cm and one level of refinement with refinement factor 4, resulting in an effective resolution of ∆x = ∆y = 0.25 cm.
This use case aims to demonstrate the ability of the algorithm to handle low-speed impact scenarios involving multiple materials and combustion. The fuel is modelled by an ideal gas EoS with γ = 1.4, c v = 718 J kg −1 K −1 and Q = 0.497 MJ kg −1 and its combustion by a single-step Arrhenius law as per Eq. 88 with C = 9.1 × 10 11 s −1 and activation temperature T A = 7974.68 K. The fuel tank casing and the projectile are both taken to be mild steel, described by the Hugoniot elastic equation of state with parameters: and perfect plasticity with yield of 137 MPa. The entire configuration is residing in air, described as an ideal gas with γ = 1.4. All materials are initially at atmospheric conditions. The left halves of the images in Fig. 14 show the pressure distribution in the fuel, in the casing and the impacting rod, while the right halves show the distribution of the AMR grids. Similarly, left halves of the images in Fig. 14 show the evolution of the reaction progress variable (λ) at late stages of the impact. After several reflections of the pressure waves at the casing walls, initiation of the gas is observed in the blue region of Fig. 15a. The reaction region grows as depicted in Figs. 15b-c. In Fig. 15c a secondary initiation zone is seen to be generated at the leftmost (and symmetrically rightmost) bottom part of the tank. In Fig. 15f the distinct reaction fronts coalesce and move to deplete the remaining fuel in the tank. It should be noted that the elastoplastic equations are solved both in the projectile and the fuel tank shell. In Fig. 15e we zoom in on the fuel tank shell to show waves that are travelling in this material. Being able to solve for the evolution of the fuel tank shell signifies that one can use this methodology to optimise the shape and material of the fuel tank shell for accident prevention. Similarly, as the impactor can be thought of as the car panel, our techniques can be used for optimising this car part as well.

Conclusions
In this work we present a methodology for the numerical simulation of the two-way interaction of reactive flow and elastoplastic structural response for automotive and defence applications. We present the coupling of a MiNi16 formulation and reduced models suitable for modelling gaseous fuels and condensed-phase explosives, with fluid (compressible Euler) and elastoplastic solid formulations. The coupling utilises level set  and Riemann ghost fluid methods to achieve communication between the materials described by different systems of equations. The fuel/explosive, solid and fluid models are all formulated as hyperbolic conservation laws and thus they can be solved by the same or similar numerical methods, facilitating the communication at material interfaces. The communication is achieved by solving mixed-material Riemann problems at the interfaces. To this end, we derive mixed Riemann solvers for each formulation pair considered in this work. These are based on the characteristic equations derived for each system pair, allowing for the computation of the star state in the 'real' material, by taking into account the two different systems on either side of the interface and applying appropriate interface boundary conditions. The elastoplastic solid model and the explosives model are in the first instance validated separately. This assesses the implementation of each model without the influence of the other formulation and the ghost fluid boundaries. The validation test chosen for the solid is at low enough speed to observe splitting of the elastic and plastic components of waves. The explosives model is validated for C4, as this is the explosive of interest in later applications. Pop-plots, as well as ZND and CJ conditions are demonstrated to match well with experiments.
The solid-explosive coupling is validated against several tests from the literature. Firstly, a one-dimensional mixed Riemann problem is utilised. We compare our approximate solutions for this test against exact solutions and other numerical solutions found in the literature. Then, we validate the coupling in cylindrical symmetry, in the context of sandwich-plate impact tests. The tests consider a non-reacting Detasheet material or a reacting C4 explosive residing between two steel plates. The sandwiched configuration is impacted by a steel projectile. In the inert scenario, we present and analyse the generated wave pattern and wave interaction. Pressure gauges are placed in Detasheet and pressure histories are compared with gauge results found in the literature, observing a very good match between the two. In the reacting scenario (ERA), we study how the generated wave pattern leads to the ignition of C4. We compare ignition times between the sandwich-plate configuration and a configuration without a rear steel plate, illustrating the effect of the plate in accelerating the ignition process.
We demonstrate the capability of the methodology to simulate high-and low-speed impact leading to the initiation of fuel/explosive in three use cases. The first use case, extends the ERA scenario by including an air gap between the front plate and the explosive and we demonstrate how the air gap hinders the initiation. The second use case is the impact of a copper can filled with reactive nitromethane. The impact initiates the nitromethane which transits to detonation. The detonation wave interacts with the confiner and shows wave-ringing in the confiner as reported in the literature. The deformation of the copper can is also highlighted, as well as the transmission of waves from the projectile to the explosive, to the confiner and finally to the surrounding air. The last use case is concerned with the accidental initiation of the fuel in the car fuel tank in crash scenarios. The speed of impact in this case is much lower than in the ERA scenario and it is demonstrated that the methodology and the code can handle this often difficult to handle scenario. The impact generates waves that reflect multiple times on the fuel tank boundaries until the fuel cannot sustain the excess energy and ignites. Multiple ignition sites are seen. The multi-dimensional use cases illustrate the accurate implementation of mixed Riemann solvers for the different systems considered in this work, describing explicitly explosive, fluid and elastoplastic materials. The methodology can also be used as part of the manufacturing process for optimising car (or other device) compartments in terms of shapes and materials.