A Diffuse Interface Model of Reactive-fluids and Solid-dynamics
AA Di ff use Interface Model of Reactive-fluids and Solid-dynamics Tim Wallis a, ∗ , Philip T. Barton b , Nikolaos Nikiforakis a a Department of Physics, University of Cambridge, Cavendish Laboratory, JJ Thomson Avenue, CB3 0HE, UK b AWE Aldermaston, Reading, Berkshire, RG7 4PR, UK
Abstract
This article presents a multi-physics methodology for the numerical simulation of physical systems that involve the non-linearinteraction of multi-phase reactive fluids and elastoplastic solids, inducing high strain-rates and high deformations. Each state ofmatter is governed by a single system of non-linear, inhomogeneous partial di ff erential equations, which are solved simultaneouslyon the same computational grid, and do not require special treatment of immersed boundaries. To this end, the governing equationsfor solid and reactive multiphase fluid mechanics are written in the same mathematical form and are discretised on a regularCartesian mesh. All phase and material boundaries are treated as di ff use interfaces. An interface-steepening technique is employedat material boundaries to keep interfaces sharp whilst maintaining the conservation properties of the system. These algorithmsare implemented in a highly-parallelised hierarchical adaptive mesh refinement platform, and are verified and validated usingnumerical and experimental benchmarks. Results indicate very good agreement with experiment and an improvement of numericalperformance compared to certain existing Eulerian methods, without loss of conservation.c (cid:13) British Crown Owned Copyright 2020 / AWE
Keywords:
Multi-physics, Multi-phase, Reactive fluids, Elastoplastic Solids, Interface Sharpening, Di ff use interface
1. Introduction
The ability to model systems containing both solid andreactive-fluid materials is of interest in many industries and aca-demic disciplines. Examples include applications in explosivesafety, mining, explosive welding, and blast-structure interac-tion. It is desirable to model materials with di ff erent physicalproperties and equations of motion in a single framework, asthis ensures full physical interaction between materials and re-duces the complexity of simulations. Furthermore, these prob-lems often feature disparate length scales for which it is de-sirable to use adaptive mesh refinement (AMR) to better focuscomputational resources. Eulerian methods are well suited toproblems of this kind, where complex, high strain-rate deforma-tions can result in large topological changes that conventionalLagrangian methods cannot resolve. However, previous Eule-rian methods are not without drawbacks and have challengesassociated with their implementation in higher dimensions andwith AMR. Resolving material boundaries in Eulerian methodsis a non-trivial problem. This is further exacerbated by mate-rials such as condensed phase explosives that introduce addi-tional phase boundaries between their reactants and products.This paper outlines a new method that resolves multi-physicsand multi-materials simultaneously, and is straightforward tosolve on a structured AMR framework.Most existing Eulerian models capable of simulating fluidsand elastoplastic solids are based on tracking material interfacesand fall into one of three broad categories: ∗ Corresponding author
Email address: [email protected] (Tim Wallis ) • Methods based on homogenised mixed cells and volume-of-fluid reconstruction, along the lines of the approachdescribed in [6]. These methods underpin many well-established and legacy multi-material codes. • Ghost-fluid methods such as [27, 33, 23, 24]. Introducedmore recently, these methods capture internal boundaryconditions by extending the multi-fluid methods from[15, 31, 37]. The main drawback of these methods is thatthey are non-conservative. • Cut-cell methods [26, 4]. These methods resolve the ge-ometry of cells intersected by interfaces and apply a strictfinite volume discretisation.All of these approaches have the advantage that they maintainarbitrarily sharp interfaces, but to do so they involve complexinterface reconstructions and mixed cell algorithms. Further-more, applying load-balancing to these schemes when usingAMR is di ffi cult due to the non-uniform numerical methods.Robust implementations can therefore be challenging to con-struct. Some of these di ffi culties have led to the developmentof hybrid approaches such as Arbitrary Lagrangian Eulerian(ALE) methods and co-simulation methods, based on embed-ding finite element grids in fluid domains [such as 10]. Thesemethods are popular for fluid-structure interaction (FSI) prob-lems, but again introduce the natural complexities of mesh man-agement associated with Lagrangian methods.Di ff use interface methods are a practical alternative classof Eulerian techniques. These methods allow a finite-volumecomputational cell to contain a mixture of several di ff erentc (cid:13) British Crown Owned Copyright 2020 / AWE a r X i v : . [ phy s i c s . c o m p - ph ] S e p aterials, and are governed by single set of evolution equa-tions that encompasses the physics of all the components inthe mixture. These methods are well established for multi-fluidproblems [see for example 1, 32, and references therein], buthave emerged only recently for coupled solid-fluid dynamics[3, 14, 13, 16]. A particular benefit of di ff use interface modelsis that they provide conservation equations for mass, momen-tum and energy across interfaces. Previous work has shownthat these schemes have the considerable advantage over inter-face tracking methods in that they can support genuine fluidmixtures, allowing for the study of phenomena such as cavi-tation and chemical reaction which rely on physical mixtures[29]. In di ff use interface methods the complexities of inter-face interactions and multi-physics are built into the equationsthemselves, and only having to solve a single system of evolu-tion equations removes much of the di ffi culty of constructingnumerical methods. In the recent model from Barton [3] for in-stance, which this work builds upon, the numerical methods areonly marginally more complicated than a conventional shockcapturing method for inviscid gas dynamics.This paper extends the method of Barton [3] to include re-active fluids by following along the lines of Michael and Niki-forakis [22]. The resulting model is capable of resolving notonly material boundaries, but also boundaries between materialphases. A practical AMR-based numerical method is detailed,which includes a recently developed interface sharpening tech-nique to counter the di ff usion of material interfaces and bringthe interface dimensions on par with those of sharp interfacemethods. To demonstrate the potential of the method, resultsare provided for several challenging tests including a multi-dimensional simulation of an explosively formed elastoplasticjet.
2. Governing Theory
This system is based on the Allaire five-equation model [1],augmented with evolution equations for elastoplastic solids fol-lowing Barton [3], and reactive, physically mixing fluids fol-lowing Michael and Nikiforakis [22] (MiNi16). With the basemodel of conservation equations in place, the additional multi-physics components require closure models to convey specificmaterial behaviour (such as plastic strain rate or reaction rate),which this work will also outline.
In the interests of focusing on multi-physics, only a sum-mary of the core evolution equations is provided here. Full de-tails, including the derivation, are outlined by Barton [3]. Ma-terials are allowed to mix at their interfaces. A material’s con-tribution to a spatially averaged physical quantity is weightedby its volume fraction, φ , in that region. This mixing is referredto as numerical mixing, to distinguish it from the physical mix-tures produced by detonations. The state of any material l ischaracterised by the phasic density ρ ( l ) , volume fraction φ ( l ) ,symmetric left uni-modal stretch tensor V e , velocity vector u , and internal energy E . The model assumes mechanical equi-librium; materials in a mixture region share a single velocity,pressure and deviatoric strain. Mixture rules are provided for aconsistent definition of thermodynamically averaged quantitiesin these mixture regions. For l = , . . . , N materials: ∂φ ( l ) ∂ t + ∂φ ( l ) u k ∂ x k = φ ( l ) ∂ u k ∂ x k (1) ∂ρ ( l ) φ ( l ) ∂ t + ∂ρ ( l ) φ ( l ) u k ∂ x k = ∂ρ u i ∂ t + ∂ ( ρ u i u k − σ ik ) ∂ x k = ∂ρ E ∂ t + ∂ ( ρ Eu k − u i σ ik ) ∂ x k = ∂ V ei j ∂ t + ∂ (cid:16) V ei j u k − V ek j u i (cid:17) ∂ x k = V ei j ∂ u k ∂ x k − u i β j − Φ i j , (5)where E = E + | u | / β j = ∂ V ek j /∂ x k . Some multi-physics closure models introducea dependence on material history variables such as equivalentplastic strain, ε p ( l ) , and reactant mass fraction, λ ( l ) . For thesematerials, additional evolution equations are required: ∂ρ ( l ) φ ( l ) α ( l ) ∂ t + ∂ρ ( l ) φ ( l ) α ( l ) u k ∂ x k = ρ ( l ) φ ( l ) ˙ α ( l ) . (6)Here α ( l ) represents any such history parameter which is ad-vected and evolved with a material as time progresses. A mate-rial may have more than one history variable, in which case α represents a vector.The system presented here allows for the fully-coupledmulti-physics solution of problems involving the interaction ofelastoplastic solids with reactive fluid mixtures. The extensionto yet more multi-physics applications would follow straight-forwardly by the inclusion of additional history parameters andclosure relations. The internal energy E for each material is defined by anequation-of-state that conforms to the general form: E ( l ) (cid:0) ρ ( l ) , T ( l ) , dev ( H e ) (cid:1) = E c ( l ) (cid:0) ρ ( l ) (cid:1) + E t ( l ) (cid:0) ρ ( l ) , T ( l ) (cid:1) + E s ( l ) (cid:0) ρ ( l ) , dev ( H e ) (cid:1) , (7)wheredev ( H e ) = ln (cid:16) V e (cid:17) (8)is the deviatoric Hencky strain tensor and T is the temperature.The three terms on the right hand side are the contribution dueto cold compression or dilation, E c ( l ) (cid:0) ρ ( l ) (cid:1) , the contribution dueto temperature deviations, E t ( l ) (cid:0) ρ ( l ) , T ( l ) (cid:1) , and the contributiondue to shear strain E s ( l ) (cid:0) ρ ( l ) , dev ( H e ) (cid:1) . The cold compressionenergy will generally be provided by the specific closure modelfor each material, covered in Section 2.3. The thermal energyis given by E t ( l ) ( ρ ( l ) , T ) = C V( l ) (cid:16) T − T l ) θ D ( l ) (cid:0) ρ ( l ) (cid:1)(cid:17) , (9)c (cid:13) British Crown Owned Copyright 2020 / AWE C V( l ) is the heat capacity, T l ) is a reference temperature,and θ D( l ) ( ρ ( l ) ) is the non-dimensional Debye temperature. TheDebye temperature is related to the Gr¨uneisen function, Γ ( ρ ( l ) ),via Γ ( l ) ( ρ ( l ) ) = ∂ ln θ D( l ) ( ρ ( l ) ) ∂ ln (cid:0) /ρ ( l ) (cid:1) = ρ ( l ) θ D( l ) ( ρ ( l ) ) ∂θ D( l ) ( ρ ( l ) ) ∂ρ ( l ) . (10)The shear energy is given by E s ( l ) ( ρ ( l ) , dev ( H e )) = G ( l ) (cid:0) ρ ( l ) (cid:1) ρ ( l ) J (dev ( H e )) , (11)where G ( ρ ) is the shear modulus, and J (dev ( H e )) = tr (cid:16) dev ( H e ) · dev ( H e ) T (cid:17) (12)is the second invariant of shear strain. This form is chosen suchthat the resultant stresses are analogous to Hooke’s law. Bruhnset al. [8] find that this form provides good empirical agreementfor a range of materials and deformations.For each component, the Cauchy stress, σ , and pressure, p ,are inferred from the second law of thermodynamics and clas-sical arguments for irreversible elastic deformations: σ ( l ) = p ( l ) I + dev (cid:0) σ ( l ) (cid:1) (13) p ( l ) = ρ l ) ∂ E ( l ) ∂ρ ( l ) (14)dev (cid:0) σ ( l ) (cid:1) = G ( l ) · dev ( H e ) . (15)Although it might appear that the model describes solid materi-als, fluids can be considered a special case where the shear mod-ulus is zero, resulting in a spherical stress tensor and no shearenergy contribution. It is the equation-of-state for each materialthat ultimately distinguishes solids from fluids. This formula-tion lends itself well to di ff use interface modelling where di ff er-ent phases that share the same underlying model can combineconsistently in mixture regions.Mixture rules must be provided to represent the state ofregions containing multiple materials in a thermodynamicallyconsistent way. The following mixture rules are applied [22, 1,3]: 1 = N (cid:88) l = φ ( l ) (16) ρ = N (cid:88) l = φ ( l ) ρ ( l ) (17) ρ E = N (cid:88) l = φ ( l ) ρ ( l ) E ( l ) (18) G = (cid:80) Nl = (cid:0) φ ( l ) G ( l ) ( ρ ( l ) ) / Γ ( l ) (cid:1)(cid:80) Nl = (cid:0) φ ( l ) / Γ ( l ) (cid:1) (19) c = (cid:80) Nl = (cid:16) φ ( l ) Y ( l ) c l ) / Γ ( l ) (cid:17)(cid:80) Nl = (cid:0) φ ( l ) / Γ ( l ) (cid:1) (20)where c is the sound speed and Y ( l ) = φ ( l ) ρ ( l ) ρ is the mass fraction. For any N × N matrix M , dev ( M ) : = M − N tr( M ) I denotes the matrixdeviator, tr( M ) denotes the trace, and I denotes the identity matrix. It can be seen that, by writing the internal energy in the formoutlined, equation (14) can be written in the form: p ( l ) = p ref , ( l ) + ρ ( l ) Γ ( l ) (cid:0) E ( l ) − E ref , ( l ) (cid:1) . (21)Here, E ref , ( l ) = E c ( l ) + E s ( l ) and p ref , ( l ) = ρ l ) ∂ E ref , ( l ) ∂ρ ( l ) . This is the formof the standard Mie-Gr¨uneisen equation-of-state. This encom-passes a wide range of di ff erent materials, not limited to solids,depending on the choice of the reference curves E ref , ( l ) ( ρ ( l ) ), p ref , ( l ) ( ρ ( l ) ), and Γ ( l ) ( ρ ( l ) ). These additional freedoms incorpo-rate thermal, compaction and shear e ff ects. For example: • E ref , ( l ) = p ref , ( l ) = , Γ ( ρ ) = Γ = γ − • p ref , ( l ) = − γ p ∞ , E ref , ( l ) = e ∞ , Γ ( ρ ) = Γ = γ − ff ened gas equation of state, used to model denserfluids such as water. • Γ ( ρ ) = Γ = γ − p ref , ( l ) ( ρ ( l ) ) = A e −R ρ ρ ( l ) + B e −R ρ ρ ( l ) (22) E ref , ( l ) ( ρ ( l ) ) = AR ρ e −R ρ ρ ( l ) + BR ρ e −R ρ ρ ( l ) (23)gives the JWL equation-of-state, widely used for con-densed phase explosives or reaction products [22]. • Γ ( ρ ) = Γ = γ − E ref , ( l ) ( ρ ( l ) ) = K ρ ( l ) α (cid:32)(cid:32) ρ ( l ) ρ (cid:33) α − (cid:33) + E s ( l ) (24) G ( ρ ( l ) ) = G (cid:32) ρ ( l ) ρ (cid:33) β + (25) p ref , ( l ) ( ρ ( l ) ) = ρ l ) ∂ E ref , ( l ) ∂ρ ( l ) (26)gives the Romenskii equation-of-state, used for elasticsolids [12]. The reactive components considered in this work do notmodel the specific chemistry of any given reaction, ratherthe approach will be to model the e ff ects on the continuumscale. As such, following Michael and Nikiforakis [22], reac-tive species are modelled as a fluid mixture, with this mixturebeing treated as one material for the purposes of bookkeeping.In this mixture, ‘reactants’ will turn to ‘products’ in a simpleone-step exothermic reaction:R → P . (27)This reactive mixture is referred to as a physical mixture todistinguish it from the numerical mixtures inherent in di ff useinterface schemes. This distinction is drawn as both the ori-gin (physical reactive source terms) and the length scale of themixing are di ff erent. It is assumed for simplicity that physicalc (cid:13) British Crown Owned Copyright 2020 / AWE α and β , each ofwhich may be governed by di ff erent equations-of-state with dif-ferent parameters. The progress of the reaction will be trackedby following the reactant mass fraction, λ ; when only reactantsare present λ =
1, and when λ = ρ ( l ) = λ ( l ) ρ α + − λ ( l ) ρ β (28) E ( l ) = λ ( l ) E α + (1 − λ ( l ) ) E β (29)This exothermic reaction powers detonation waves. The en-ergy input will generally arise by including the reaction energyas a term in the energy of the products: E ref , ( l ) → E ref , ( l ) − Q , (30)where Q is the energy released by the reaction. The rate atwhich this reaction occurs is defined by the choice of reactionrate law. For example, this could be given by the Arrhenius ratelaw: ˙ λ = k λ = Ae − EaRT λ , (31)where k is the reaction rate in s − , A is the so-called ‘pre-exponential factor’, E a is the activation energy, R is the univer-sal gas constant and T is the temperature of the mixture. Thisequation can equally be written:˙ λ = A λ e − T ref T , (32)where T ref is an activation temperature. Other rate laws willalso be considered. This equation will be included to track themixture mass fraction as it evolves: ∂ρ ( l ) φ ( l ) λ ( l ) ∂ t + ∂ρ ( l ) φ ( l ) λ ( l ) u k ∂ x k = ρ ( l ) φ ( l ) λ ( l ) Ae − T ref T ( l ) . (33)The temperature is taken to be T ( l ) = p − p ref , ( l ) ρ ( l ) Γ ( l ) C V ( l ) , (34)where C V is the specific heat at constant volume for the mate-rial. The sound speed in a physical mixture is also calculatedusing the reactant mass fraction weighted form from Michaeland Nikiforakis [22]: c , ( l ) = p ρ l ) − (cid:16) ∂ E ( l ) ∂ρ ( l ) (cid:17) p (cid:16) ∂ E ( l ) ∂ p (cid:17) ρ ( l ) (35a) (cid:32) ∂ E ( l ) ∂ p (cid:33) ρ ( l ) = λ ( l ) (cid:32) ∂ E α ∂ p (cid:33) ρ ( l ) + (1 − λ ) (cid:32) ∂ E β ∂ p (cid:33) ρ ( l ) (35b) (cid:32) ∂ E ( l ) ∂ρ ( l ) (cid:33) p = λ ( l ) (cid:32) ∂ E α ∂ρ α (cid:33) p (cid:32) ∂ρ α ∂ρ ( l ) (cid:33) p + (1 − λ ) (cid:32) ∂ E β ∂ρ β (cid:33) p (cid:32) ∂ρ β ∂ρ ( l ) (cid:33) p . (35c) The introduction of plasticity through the source term Φ fol-lows the method of convex potentials (see for example [28]) andis therefore thermodynamically compatible. In this approach,the von Mises yield criterion forms the scalar potential whichleads to the plastic flow rule: Φ = χ (cid:114)
32 dev ( σ ) || dev ( σ ) || V e . (36)The plastic flow rate χ is a closure model and must be suitablefor arbitrary mixtures. For a mixture of N materials: χ = (cid:80) Nl = (cid:0) φ ( l ) χ ( l ) ( ρ ( l ) ) / Γ ( l ) (cid:1)(cid:80) Nl = (cid:0) φ ( l ) / Γ ( l ) (cid:1) . (37)If any material is does not obey a plasticity model, that materialcontributes χ ( l ) = χ relates the particular material flow model.This paper considers both ideal plasticity, where χ is a Heavi-side function such that the relaxation is non-zero only when thestress exceeds the yield surface σ Y : χ ( l ) = χ l ) H (cid:114) || dev( σ ) ( l ) || − σ Y , (38)and the rate sensitive isotropic work-hardening plasticity out-lined by Johnson and Cook [19]: χ ( l ) = χ l ) exp C (cid:113) || dev (cid:0) σ ( l ) (cid:1) || σ Y (cid:16) ε p , ( l ) (cid:17) − , (39)where χ > C controls the rate dependency. In the Johnson and Cook modelthe yield stress is given by: σ Y (cid:16) ε p , ( l ) (cid:17) = c + c ( ε p , ( l ) ) n , (40)where c , c and n are material dependent constants. This yieldsurface is a function of the accumulated plastic strain, ε p , mak-ing it necessary to add an additional evolution equation to thesystem to advect and evolve the plastic strain: ∂ρ ( l ) φ ( l ) ε ( l ) , p ∂ t + ∂ρ ( l ) φ ( l ) ε ( l ) , p u k ∂ x k = ρ ( l ) φ ( l ) χ ( l ) . (41)Using the Johnson-Cook model in this way results in a vis-coplastic flow rule, where plastic deformations can accumulatefrom the onset of loading. Note however that the parameter C is usually small such that χ ( l ) (cid:28) χ l ) for stresses much belowthe characteristic stress. Indeed, as C →
3. Numerical Approach
The model is solved on a Cartesian mesh with local res-olution adaptation in space and time. This is achieved usingc (cid:13)
British Crown Owned Copyright 2020 / AWE ff erentialequations (PDEs) of the form of equation (42). In this approach,cells of identical resolution are grouped into logically rectangu-lar sub-grids or ‘patches’. Refined grids are derived recursivelyfrom coarser ones, based upon a flagging criterion, to form ahierarchy of successively embedded levels. All mesh widths onlevel l are r l -times finer than on level l −
1, i.e. ∆ t l : = ∆ t l − / r l and ∆ x l : = ∆ x l − / r l with r l ∈ N , r l ≥ l > r = l by calling a single-grid update routine in a loop over all patches constituting thelevel. The discretisation of the constitutive models does not dif-fer between patches or levels, so for clarity the method shall bedescribed for a single sub-grid. Cell centres are denoted by theindices i , j , k ∈ Z and each cell C li jk has the dimensions ∆ x li jk .The system of equations can be written compactly in vectorform by separating it into various qualitatively di ff erent parts:a conservative hyperbolic part for each spatial dimension, non-conservative terms from the volume fraction and stretch tensorupdates, a source term due to plastic flow, a source term due toreactive species, and a source term to account for geometricale ff ects. This can be written as: ∂ q ∂ t + ∂ g k ∂ x k = s non-con. + s p + s r + s g . (42)Subject to the closure relations previously outlined, this is givenby: ∂∂ t φ ( l ) φ ( l ) ρ ( l ) φ ( l ) ρ ( l ) λ ( l ) φ ( l ) ρ ( l ) ε p , ( l ) ρ u i ρ E V ei j + ∂∂ x k φ ( l ) u k φ ( l ) ρ ( l ) u k φ ( l ) ρ ( l ) λ ( l ) u k φ ( l ) ρ ( l ) ε p , ( l ) u k ρ u i u k − σ ik ρ Eu k − u i σ ik V ei j u k − V ek j u i + · · · (43) = φ ( l ) ∂ u k ∂ x k V ei j ∂ u k ∂ x k − u i ∂ V ekj ∂ x k + φ ( l ) ρ ( l ) ˙ ε p , ( l ) Φ i j + φ ( l ) ρ ( l ) ˙ λ ( l ) + s g (44)When considering cylindrical symmetry, the conservative vari-ables and geometrical term are given by: q = φ ( l ) φ ( l ) ρ ( l ) φ ( l ) ρ ( l ) λ ( l ) φ ( l ) ρ ( l ) ε p , ( l ) ρ u r ρ u z ρ E V ei j , s g = − r φ ( l ) ρ ( l ) u r φ ( l ) ρ ( l ) λ ( l ) u r φ ( l ) ρ ( l ) ε p , ( l ) u r ρ u r − σ rr + σ θθ ρ u z u r − σ rz ρ Eu r − ( u r σ rr + u z σ zr ) V ei j u r − δ i θ V ei j u r . (45) The inhomogeneous system is integrated for time intervals[ t n , t n + ] where the time-step ∆ t = t n + − t n is chosen to be a frac-tion of the global maximum allowable time step required forstability of the hyperbolic update method. For high strain-rateor highly reactive applications, the plastic relaxation or reac-tions can occur over smaller time scales, which can lead to localsti ff ness. To address this issue, without resorting to forecastingsti ff zones and resolving the time scales of irreversible physicalprocesses, Godunov’s method of fractional steps is used, wherethe hyperbolic part is updated first, followed by serially addingin contributions from each source term: ∂ q ∂ t = − ∂ g k ∂ x k + s non-con. IC: q n ∆ t = ⇒ q (cid:73) (46) ∂ q ∂ t = s x IC: q (cid:73) ∆ t = ⇒ q (cid:70) , (47)where the result of the each step is used as the initial condition(IC) for the next. Here s x = s p , s r , s g . Once all source termshave been added, the last q (cid:70) becomes q n + . Employing the method of lines and replacing the spatialderivatives with a conservative approximation, the hyperbolicsystem can be writtendd t q i jk + D i jk ( q ) = , (48)where q i jk represents the vector of conservative variables storedat cell centres, and D i jk : = ∆ x i jk (cid:16)(cid:101) g i + / , jk − (cid:101) g i − / , jk (cid:17) + ∆ x i jk (cid:16)(cid:101) g i , j + / , k − (cid:101) g i , j − / , k (cid:17) + ∆ x i jk (cid:16)(cid:101) g i j , k + / − (cid:101) g i j , k − / (cid:17) − s non-con. , i jk , (49)where (cid:101) g lm ± / , for m = i , j , k , are the cell wall numerical fluxfunctions. The numerical fluxes are computed through suc-cessive sweeps of each spatial dimension and summed accord-ing to equation (49). Fluxes are computed using the HLLDsolver from Barton [3]. To achieve higher order spatial accu-racy the initial conditions for the Riemann solver are taken tobe MUSCL reconstruction of cell centred primitive variables. Itshould be noted that this procedure is generally carried out onthe conservative variables; this work instead follows Johnsenand Colonius [18], who showed that, for multi-material prob-lems, reconstruction of the primitive variables leads to feweroscillations around interfaces. An artificial interface recon-struction is applied to reduce numerical di ff usion around inter-faces. This is achieved using the Tangent of Hyperbola IN-terface Capturing (THINC) method: an algebraic interface re-construction technique that fits a hyperbolic tangent functionto variables inside a cell. In contrast to Barton [3] who usedthe original THINC method of Xiao et al. [38], the more re-cent MUSCL-BVD-THINC scheme of Deng et al. [11] is em-ployed. This more recent scheme provides an additional checkc (cid:13) British Crown Owned Copyright 2020 / AWE
5o minimise oscillations by comparing the reconstructed state’scell boundary variation with the previously calculated MUSCLreconstruction. THINC-reconstructed states are only acceptedwhen their total boundary variation is lower than that of theMUSCL scheme alone. This algorithm is provided in AppendixA for convenience. The non-conservative source terms that re-sult from the volume fraction and deformation tensor updatesare added in the hyperbolic step following a procedure similarto that used by several authors [22, 1, 3]. To achieve a highertemporal resolution in the update of the hyperbolic terms, a sec-ond order Runge-Kutta time integration is used. q (1) = q n − ∆ t D ( q n ) (50) q (2) = q (1) − ∆ t D (cid:16) q (1) (cid:17) (51) q (cid:73) = q n + q (2) . (52) The plastic update is performed after the hyperbolic stepas detailed above. The potentially sti ff ODEs governing theplasticity evolution are solved using an analytical technique asdetailed by Barton [3], which reduces the problem to a singleODE for each material in a given cell. The algorithm can besummarised as follows: (cid:16) V e (cid:17) (cid:70) = exp (cid:16) dev ( H e ) (cid:70) (cid:17) (53)dev ( H e ) (cid:70) = J (cid:70) J (cid:73) dev (cid:18) H e (cid:18) V e (cid:73) (cid:19)(cid:19) (54) J (cid:70) = (cid:80) Nl = (cid:0) φ ( l ) / Γ ( l ) (cid:1) J (cid:70) ( l ) (cid:80) Nl = φ ( l ) / Γ ( l ) (55) R (cid:16) J (cid:70) ( l ) (cid:17) = J (cid:70) ( l ) − J (cid:73) ( l ) + ∆ t (cid:114) χ ( l ) (cid:16) J (cid:70) ( l ) , ε (cid:70) p ( l ) (cid:16) J (cid:70) ( l ) (cid:17)(cid:17) (56)The last objective equation is solved using a simple bisectionalgorithm between the limits J (cid:70) ( l ) ∈ [0 : J (cid:73) ( l ) ]. As a resultnumerical error the stretch tensor to lose its symmetry and uni-modular properties. To mitigate this, at the end of every timestep, the symmetry must be reinstated by evaluating V e ← (cid:113) V e V e T (57)It is necessary to evaluate this algorithm, at least in part, every-where in a problem irrespective of whether a material adheresto a plasticity law or not, since the algorithm encompasses theenforcement of equation (57) to ensure that the stretch tensorremains symmetric, and the condition V e = I is true for fluids.The implementation of this algorithm by Barton [3] employssingular value decomposition (SVD) so that equation (53) andequation (54) are evaluated exactly. Since SVD is relatively ex-pensive, this part of the overall numerical scheme constitutes a primary overhead. To alleviate this computational burden, thefollowing modifications are made.Firstly, prior to proceeding with the algorithm the volumefractions are interrogated to check if a cell contains all or pre-dominantly fluid materials, in which case the condition V e = I is enforced directly. This simple additional step makes a pro-found di ff erence for problems where fluids occupy large regionsof a problem, but adds negligible cost to the original method.Secondly, the logarithmic strain tensor in equation (54) isevaluated using a variant of the approximation by Baˇzant [5]:dev ( H e ) ≈ H e (1) B (cid:18) V e V e T (cid:19) , (58)where H e ( m ) B (cid:16) V e (cid:17) = m (cid:18) V em − V e − m (cid:19) , (59)and the invariants evaluated from this in turn; the exponentialmatrix in equation (53) is evaluated using the first Pad´e approx-imant:( V e ) n + ≈ (cid:32) −
12 dev ( H e ) n + (cid:33) − (cid:32) +
12 dev ( H e ) n + (cid:33) . (60)The approximation of the strain was proposed in Barton [3] butwas used only where stresses are computed for the numericalflux functions. Use of the first Pad´e approximate is found tobe su ffi cient over second or higher variants since the secondnorm of H e is known a priori to not exceed su ffi ciently smallvalues for the materials of interest. These approximations avoidthe use of a SVD, and the algorithm reduces to simple matrixoperations. Further e ffi ciency can be achieved noting V e anddev ( H e ) are unimodular, meaning their inverses are equivalentto their cofactor, the evaluation of which avoids computing thedeterminant. The potentially sti ff reactive update consists of solving theequation:dd t (cid:0) ρ ( l ) φ ( l ) λ ( l ) (cid:1) = ρ ( l ) φ ( l ) ˙ λ ( l ) . (61)A second order Runge-Kutta integration was used for this task.In the course of computation it is necessary to be able to convertbetween the primitive variables and the conservative variables,of which the only non-trivial conversion is between p and E .This requires knowledge of the mixture phase densities ρ α and ρ β . As only the total phase density of the mixture is evolved,a non-linear root finding procedure is required to define thesecomponent densities before the conversion can be performed.The outline of this root finding procedure is included in Ap-pendix B. In this work, only a single physical mixture species isconsidered in a simulation, as including more than this requiresan unstable multi-dimensional root finding procedure. Oncethese component densities are defined, the conversion may be-gin. Starting with the single-phase equations of state: p ( l ) = p ref , ( l ) + ρ ( l ) Γ ( l ) (cid:0) E ( l ) − E ref , ( l ) (cid:1) , (62)c (cid:13) British Crown Owned Copyright 2020 / AWE igure 1: A detonation wave impacting an elastoplastic solid, transmitting ashock wave. This test was used by Schoch et al. [33] in a sharp interface con-text, and this paper’s results in a purely di ff use interface scheme agree wellthis author’s results. The images are taken at t = -10, 0, 0.4, 19.5 µ s after thecollision and show − σ xx . This test demonstrates the ability of the model tosimultaneously handle reactive fluid mixtures and elastoplastic solids. the mixing rule for internal energy, equation (29), gives: ρ E = N (cid:88) l = ρ ( l ) E ( l ) = (cid:88) Non-mixtures, x φ x ρ x E x + (cid:88) Mixtures, y φ y (cid:16) ρ y λ y E α, y + ρ y (1 − λ y ) E β, y (cid:17) . (63)Pressure equilibrium assumptions then allow these equations tobe combined, to allow ρ E to be calculated from p : ρ E = (cid:88) Non-mixture, x φ x (cid:32) ( p − p ref , x ) Γ x + E ref , x (cid:33) + (64) φ y (cid:32) ρ y λ y ( p − p ref ,α ) ρ α Γ α + ρ y λ y E ref ,α ρ α + (65) ρ y (1 − λ y )( p − p ref ,β ) ρ β Γ β + ρ y (1 − λ y ) E ref ,β ρ β (cid:33) . Alternatively, the equation can be rearranged to find p from ρ E : p = (cid:80) x φ x Γ x + φ y ρ y λ y ρ α Γ α + φ y ρ y (1 − λ y ) ρ β Γ β (cid:34) ρ E + · · · (cid:88) x φ x (cid:32) p ref , x Γ x − E ref , x (cid:33) + · · · φ y (cid:32) λ y ρ y p ref ,α ρ α Γ α − ρ y λ y E ref ,α ρ α + · · · (1 − λ y ) ρ y p ref ,β ρ β Γ β − ρ y (1 − λ y ) E ref ,β ρ β (cid:33) (cid:35) . (66) Figure 2: Comparison with experiment [34] for the LX-17 rate stick test. Goodagreement is observed with experimental values, and a marked di ff erence be-tween confined and unconfined detonation velocity is seen.
4. Validation and Evaluation
The interaction of reactive fluids with elastic solids was ex-amined to ensure this multiphase aspect can be handled by thisscheme. A one-dimensional test employed by Schoch et al. [33]was considered. This test involves a detonation wave travel-ling through a multiphase explosive and then hitting a slab ofpurely elastic copper. Here, x = [0 , .
4] m, CFL = N = x < .
15 m with explosive occupying the rest ofthe domain. The solid is governed by the Romenskii equation-of-state, with the parameters ρ = − , K = . G = . Γ = α = β =
3. The explosiveused is the idealised multiphase mixture used by Schoch et al.[33]. The reactants are governed by a sti ff ened gas equation-of-state with γ = . p ∞ = . − ,and the products are governed by an ideal gas equation-of-statewith γ = .
0. The reaction rate law used by Schoch et al. [33]is a simplified mass fraction based law. In terms of this model,this law can be expressed as:˙ λ = k √ λ , (67)where k = × s − . The reaction energy is Q = . × Jkg − . To ignite the detonation wave, the region 0 . < x < .
40 m has a pressure of 10 Pa, with the rest of the domainhaving a pressure of 10 Pa.Figure 1 shows this test. A correct, stable interaction be-tween the fluid detonation wave and the quiescent copper wasobserved, matching the sharp interface results of Schoch et al.[33].
Similar to the previous test is the confined rate stick test,commonly used to assess reactive models [2], [22]. This testc (cid:13)
British Crown Owned Copyright 2020 / AWE igure 3: Pressure ( top ), density and AMR patches ( bottom ) in the 20 mm LX-17 confined rate stick test. This image is taken at 1.55 × − s. This test demonstratesthe both the AMR working correctly and the dual multi-physics capabilities of plasticity and reactive fluids working simultaneously. features a cylinder of the condensed phase explosive LX-17 sur-rounded by an elastoplastic copper container. This scenario isthen compared to the unconfined case where the explosive issurrounded by air. LX-17 is a slightly non-ideal explosive, so adi ff erence in the confined and unconfined steady state detona-tion velocities is observed experimentally [34].The explosive’s products and reactants are governed by theJWL equation-of-state previously outlined. The reactants havethe parameters A = . × Pa , B = − . × Pa , R = . , R = . , Γ = . , C V = . − K − and theproducts have the parameters A = . × Pa , B = . × Pa , R = . , R = . , Γ = . , C V = . − K − [35, 36, 17]. The detonation energy is givenby Q = . × Jkg − .The ignition and growth rate law was used to describe theexplosive. More details of the use of this rate law for LX-17 areoutlined in Tarver [36]. This is a three-stage rate law, based onphenomenological experience of how detonations in condensedphase explosives evolve. The rate can be expressed as:˙ λ = I (1 − F ) b (cid:32) ρρ − − a (cid:33) x H ( F ig − F ) (68) + G (1 − F ) c F d p y H ( F G − F ) + G (1 − F ) e F g p z H ( F − F G )where F = − λ is the reacted fraction and H is theHeaviside function. Other parameters are material dependentconstants, given in this case by I = × s − , G = − s − , G = . × GPa − s − , a = . , b = . , c = . , d = , e = . , x = , y = , z = , F ig = . , F G = . , F G = . ρ = − .A cylindrically symmetric domain was used for these tests,with a reflective boundary condition on the axis. A CFL num-ber of 0.4 was used. The explosive cylinder had a radius vary- ing from r = r = [0 , r ], z = [0 , r ]. AMR is required for this test, as the resolution ofthe detonation front is crucial in measuring accurate detonationvelocities. Four refinement levels each with a refinement factorof 2 and base resolution of 24 ×
160 are used, giving an e ff ec-tive resolution of 384 × ff erence-based approach istaken, where when the fractional change in a refinement vari-able between a cell and its neighbour is larger than a specifiedamount, those cells are tagged for refinement in the procedureoutlined previously. Generally this flagging is performed on thevariables φ ( l ) , φ ( l ) ρ ( l ) , u i , and p .The solid was described by the Romenskii equation-of-statewith the parameters laid out in Section 4.1, with a yield stress of σ y = . γ = . × Pa and a thicknessof 0.5 r , with the other materials at atmospheric pressure. Theexplosive is initialised with a density of ρ = ρ = − .It is worth reiterating here that this test is performed with allmaterials having di ff use interfaces between them. There is nouse of level sets, only the THINC algorithm to keep the appro-priate interfaces sharp. This contrasts the mixed approach usedby both Schoch et al. [33] and Michael and Nikiforakis [23]Figure 2 shows the results of several confined and uncon-fined rate stick tests with comparison against experiment [34].The results are non-dimensionalised following the procedure inIoannou et al. [17]. The results show good agreement with theexperimental values, validating the reactive part of the model.Figure 3 shows the r =
20 mm confined test for reference. Thistest shows good qualitative agreement with previous numericalrate stick tests, such as the strong confinement case of Banksc (cid:13)
British Crown Owned Copyright 2020 / AWE
8t al. [2] and Michael and Nikiforakis [22]. The THINC recon-struction maintains sharp interfaces, even when using AMR.
To more rigorously test the coupling of the explosive andelastoplastic solids and the ability of the THINC algorithm tomaintain sharp interfaces over large deformations, the ‘explo-sive filled cooper vessel’ tests from Schoch et al. [33] were con-sidered. This test is commonly employed in detonation studies,and appears in a variety of forms [25], [22], [23].Both the booster- and flyerplate-ignited versions of thistest were considered. A cylindrically symmetric domain of r = [0 ,
20] cm, z = [0 ,
40] cm was used, with a CFL numberof 0.4 and a resolution of N = × − . For the booster ignited version, the partition is removedand the detonation is started by a 1 cm thick booster region witha pressure of 10 Pa.The explosive here is the same idealised explosive as wasused in the one-dimensional reactive test and the copper isagain governed by the Romenskii equation-of-state. Follow-ing Schoch et al. [33], this test uses ideal plasticity with a yieldstress of σ Y = .
07 GPa.Figure 4 shows how the pressure and density progressedthrough the flyerplate test. This method recovers the same be-haviour as that found by Schoch et al. [33]; the shock is able topass through the solid partition and reignite on the other side, ademonstration of the full multi-physics coupling of this model.Figure 6 shows the percentage mass change for the solidphase in the booster ignited test, compared to the results ob-tained by Schoch et al. [33]. The method presented here out-performed the Ghost Fluid method at all resolutions tested.This is due to the intrinsic conservation errors in sharp inter-face methods that this scheme avoids. On top of this, Figure 5shows how the THINC reconstruction helps to reduce numer-ical di ff usion; evidently this is required when the copper con-finer undergoes such large deformations. To strenuously test all parts of the model, the final test con-sidered is an explosively-initiated shock colliding with a hemi-spherical indentation in copper. This test has been performedexperimentally [20] and numerically [30, 9], making it a suit-able candidate for validation. Experimental research [20, 21]has found that the problem has the benefit of being charac-terised by the radius of the indentation, if material parametersare kept constant, allowing for straight forward comparison be-tween tests. When the shock hits the hemispherical indenta-tion, the resulting converging wave causes a very high speed ( ≈ Mach 10), high deformation, thin spike to jet out from the metalsurface. Experimentally, the shock is generated by the detona-tion of an explosive below the metal. Previous numerical sim-ulations have instead imposed shock conditions on the metal Test Jet Velocity / kms − Jet Radius / mm4mm RadiusExperimental [20] 2.5 0.8This work 2.58 ± ± ± ± ± ± Table 1: Results of the hemispherical indentation test from di ff erent sourcesfor two di ff erent radii: 4mm and 15mm. Excellent agreement with previoussimulation and experimental work is seen for the work at hand. without the use of an explosive. This work includes both thisand a multiphase explosive to drive the shock, with the multi-physics approach enabling the full experimental set-up to bemodelled.The thin jet is also a good test of the interface capturingmethod, and this work’s results can be compared to those ofSambasivan et al. [30] who performed this test with a GhostFluid method.The copper in these tests used the Romenskii equation ofstate previously outlined, but now following the Johnson andCook plasticity model (39) with the parameters C = . c = . c = .
177 GPa, n = . χ = . − . The air wasan ideal gas with γ = .
4, and the explosive was the idealisedcondensed phase explosive used in the confined rate stick test,now following the Arrhenius reaction rate to demonstrate thecapability of the model to handle more complex reaction rates.This rate law is given by:˙ λ = A λ e − T ref T , (69)with A = . × , T ref =
210 K, and a detonation energy Q = × Jkg − .Two radii were considered: 4mm and 15mm. The 15mmtest used the full multi-phase explosive ignition and is shown inFigure 9a. The 4mm test uses the simple shocked metal initialconditions and is shown in Figure 8.The initial conditions for these tests are outlined in Figure7. Both tests employed a cylindrically symmetrical domain. Inthe 4mm test, the shock is started by accelerating the bottom ofthe domain to 1000 ms − and in the 15mm test the explosiveis ignited by a 1mm booster region with a pressure of 1 GPaat its base. Adaptive mesh refinement is used in both of thesetests, with a base mesh of 64 ×
560 and one level of refinementfor the 4mm test, and a base mesh of 24 ×
352 and two levels ofrefinement for the 15mm test. A CFL number of 0.4 was usedfor these tests.Table 1 shows how these numerical tests compare with pre-vious results in terms of the radius and initial velocity of the jetproduced.c (cid:13)
British Crown Owned Copyright 2020 / AWE igure 4: The Schoch et al. [33] shock-induced copper can test. The images are taken at t =
5, 10, 19.5, 29, 33.5 and 54.5 µ s, chosen to correspond to the resultsof Schoch et al. [33]. In each frame, the top shows the pressure and the bottom shows the density. This test demonstrates a full coupling of the elastoplastic andreactive mixture components of the model, showing good agreement with the sharp interface results of Schoch et al. [33]Figure 5: A comparison of THINC reconstruction and no reconstruction for the solid volume fraction profile of the Schoch et al. [33] can test at 40 µ s. The blackpoints show the THINC profile and the blue points show the no reconstruction profile. c (cid:13) British Crown Owned Copyright 2020 / AWE igure 6: The percentage solid mass loss over the course of the Schoch explo-sive can test for di ff erent spatial resolutions. Data obtained by Schoch et al.[33] is plotted as points and this paper’s data is plotted as lines. Correspondingcolours represent the same resolution of test. In all cases, the mass loss incurredby this model is significantly lower than that of Schoch et al. [33]. ExplosiveBoosterSolidAir
SolidShockAir
Figure 7: ( left ) multi-physics and ( right ) simple shock initial conditions for thehemispherical indentation test. Figure 8: Experimental comparison with the work of Mali [20] for the 4mmindentation test. Images are taken at 4, 8, 16 and 20 mu s. Figure 8 shows the 4mm test alongside experimental imagesfrom Mali [20]. Good qualitative agreement is seen, backed upby the quantitative agreement in Table 1 for the jet radius andvelocity. Figure 9a shows 3D images of the 15mm test, show-ing the pressure in the solid and explosive, and a numericalschlieren of the shock waves in the surrounding air. Addition-ally, Figure 9b shows the e ff ect of the THINC interface sharp-ening. The jet is longer and wider, its interfaces are sharper andshock waves in the surrounding air are captured earlier, all forthe same resolution.
5. Conclusion
In this work, a multi-material, multi-phase, multi-physicsdi ff use interface scheme with crisp material interfaces has beencreated and numerically validated against experiment and previ-ous simulations. This was accomplished by combining di ff useinterface elastoplastic methods, reactive fluid mixture methods,and THINC interface reconstruction techniques. This com-bined system solves the same set of equations over the entiredomain, with all materials being defined globally, without theneed for a dividing level set. The model is general, capable ofhandling an arbitrary number of materials that obey the broadclass of Mie-Gr¨uneisen equations of state in three dimensions,some of which can be reactive fluid mixtures. This enables thec (cid:13) British Crown Owned Copyright 2020 / AWE a) The 15mm radius hemispherical indentation test. The images show the pressure in the explosive and solid,and a numerical schlieren of the density in the surrounding air. The images are taken at 0, 15, 30, 40, 60 and60 µ s. (b) A comparison of no reconstruc-tion (left) and THINC reconstruc-tion (right) for the 15mm radiushemispherical indentation test. Aclear di ff erence in the ability of themodel to capture the thin jet is ob-served. c (cid:13) British Crown Owned Copyright 2020 / AWE
ReferencesReferences [1] Gr´egoire Allaire, S´ebastien Clerc, and Samuel Kokh. A five-equationmodel for the numerical simulation of interfaces in two-phase flows.
Comptes Rendus de l’Academie des Sciences Series I Mathematics , 331(12):1017–1022, 2000. ISSN 0764-4442.[2] J. W Banks, W. D Henshaw, D. W Schwendeman, and A. K Kapila.A study of detonation propagation and di ff raction with compliant con-finement. Combustion Theory and Modelling , 12(4):769–808, 2008.ISSN 1364-7830. URL .[3] Philip T. Barton. An interface-capturing godunov method for the sim-ulation of compressible solid-fluid problems.
Journal of ComputationalPhysics , 2019. ISSN 0021-9991.[4] P.T. Barton, B. Obadia, and D. Drikakis. A conservative level-set basedmethod for compressible solid / fluid problems on fixed grids. Journal ofComputational Physics , 230(21):7867 – 7890, 2011. ISSN 0021-9991.doi: https: // doi.org / / j.jcp.2011.07.008.[5] Zdenˇek P. Baˇzant. Easy-to-compute tensors with symmetric inverse ap-proximating hencky finite strain and its rate. Journal of Engineering Ma-terials and Technology , 120(2):131–136, 04 1998. ISSN 0094-4289. doi:10.1115 / https://doi.org/10.1115/1.2807001 .[6] David J. Benson. Computational methods in lagrangian and eulerian hy-drocodes. Comput. Methods Appl. Mech. Eng. , 99(23):235394, Septem-ber 1992. ISSN 0045-7825. doi: https: // doi.org / / Journal of Computational Physics , 82(1):64–84, 1989.ISSN 0021-9991.[8] O. T. Bruhns, H. Xiao, and A. Meyers. Constitutive inequalities foran isotropic elastic strain-energy function based on hencky’s logarithmicstrain tensor.
Proceedings of the Royal Society A: Mathematical, Physi-cal and Engineering Sciences , 457(2013):2207–2226, 2001. ISSN 1364-5021.[9] S.R Cooper, D.J Benson, and V.F Nesterenko. A numerical exploration ofthe role of void geometry on void collapse and hot spot formation in duc-tile materials.
International Journal of Plasticity , 16(5):525–540, 2000.ISSN 0749-6419.[10] Ralf Deiterding, Raul Radovitzky, Sean Mauch, Ludovic Noels, Julian C.Cummings, and Daniel I. Meiron. A virtual test facility for the e ffi cientsimulation of solid material response under strong shock and detonationwave loading. Engineering with Computers , 22:325–347, 2006.[11] Xi Deng, Satoshi Inaba, Bin Xie, Keh-Ming Shyue, and Feng Xiao. Highfidelity discontinuity-resolving reconstruction for compressible multi-phase flows with moving interfaces.
Journal of Computational Physics ,371:945–966, 2018. ISSN 0021-9991.[12] V. Dorovskii, A. Iskol’dskii, and E. Romenskii. Dynamics of impulsivemetal heating by a current and electrical explosion of conductors.
Journalof Applied Mechanics and Technical Physics , 24(4):454–467, 1983. ISSN0021-8944.[13] N. Favrie and S.L. Gavrilyuk. Di ff use interface model for compressiblefluid compressible elasticplastic solid interaction. Journal of Computa-tional Physics , 231(7):2695–2723, 2012. ISSN 0021-9991.[14] N. Favrie, S.L. Gavrilyuk, and R. Saurel. Solidfluid di ff use interfacemodel in cases of extreme deformations. Journal of ComputationalPhysics , 228(16):6037–6077, 2009. ISSN 0021-9991.[15] Ronald P Fedkiw, Tariq Aslam, Barry Merriman, and Stanley Osher. Anon-oscillatory eulerian approach to interfaces in multimaterial flows (theghost fluid method).
Journal of Computational Physics , 152(2):457–492,1999. ISSN 0021-9991. [16] Sarah Hank, Sergey Gavrilyuk, Nicolas Favrie, and Jacques Massoni. Im-pact simulation by an eulerian model for interaction of multiple elastic-plastic solids and fluids.
International Journal of Impact Engineering ,109(C):104–111, 2017. ISSN 0734-743X.[17] Eleftherios Ioannou, Stefan Schoch, Nikolaos Nikiforakis, and LouisaMichael. Detonation propagation in annular arcs of condensed phase ex-plosives. 29(11), 2017. ISSN 10706631.[18] Eric Johnsen and Tim Colonius. Implementation of weno schemes incompressible multicomponent flow problems.
Journal of ComputationalPhysics , 219(2):715–732, 2006. ISSN 0021-9991.[19] Gordon R. Johnson and William H. Cook. Fracture characteristics ofthree metals subjected to various strains, strain rates, temperatures andpressures.
Engineering Fracture Mechanics , 21(1):31–48, 1985. ISSN0013-7944.[20] V. I. Mali. Flow of metals with a hemispherical indentation under theaction of shock waves.
Combustion, Explosion and Shock Waves , 9(2):241–245, Mar 1973. ISSN 1573-8345. doi: 10.1007 / BF00814821. URL https://doi.org/10.1007/BF00814821 .[21] V. I. Mali, V. V. Pai, and A. I. Skovpin. Investigation of the breakdownof flat jets.
Combustion, Explosion and Shock Waves , 10(5):676–682,Sep 1974. ISSN 1573-8345. doi: 10.1007 / BF01463984. URL https://doi.org/10.1007/BF01463984 .[22] L. Michael and N. Nikiforakis. A hybrid formulation for the numeri-cal simulation of condensed phase explosives.
Journal of ComputationalPhysics , 316:193–217, 2016. ISSN 0021-9991.[23] L. Michael and N. Nikiforakis. A multi-physics methodology for thesimulation of reactive flow and elastoplastic structural response.
Journalof Computational Physics , 367(C):1–27, 2018. ISSN 0021-9991.[24] Louisa Michael, Stephen T. Millmore, and Nikolaos Nikiforakis. A multi-physics methodology for four-states of matter. 2019.[25] G.H. Miller and P. Colella. A conservative three-dimensional eulerianmethod for coupled fluid-solid shock capturing.
Journal of Computa-tional Physics , 183(1), 2002. ISSN 0021-9991.[26] G.H. Miller and P. Colella. A conservative three-dimensional eulerianmethod for coupled solid-fluid shock capturing.
Journal of Computa-tional Physics , 183(1):26 – 82, 2002. ISSN 0021-9991. doi: https: // doi.org / / jcph.2002.7158.[27] A. L´opez Ortega, M. Lombardini, P.T. Barton, D.I. Pullin, and D.I. Me-iron. Richtmyermeshkov instability for elasticplastic solids in converginggeometries. Journal of the Mechanics and Physics of Solids , 76:291 –324, 2015. ISSN 0022-5096. doi: https: // doi.org / / j.jmps.2014.12.002.[28] Niels Saabye Ottosen. The mechanics of constitutive modeling / NielsSaabye Ottosen, Matti Ristinmaa.
Elsevier, Amsterdam ; London, 1st ed.edition, 2005. ISBN 008044606X.[29] Fabien Petitpas, Jacques Massoni, Richard Saurel, Emmanuel Lapebie,and Laurent Munier. Di ff use interface model for high speed cavitatingunderwater systems. International Journal of Multiphase Flow , 35(8):747–759, 2009. ISSN 0301-9322.[30] S Sambasivan, A Kapahi, and H.S Udaykumar. Simulation of high speedimpact, penetration and fragmentation problems on locally refined carte-sian grids.
Journal of Computational Physics , 235:334–370, 2013. ISSN0021-9991.[31] Shiv Kumar Sambasivan and H. S Udaykumar. Ghost fluid method forstrong shock interactions part 1: Fluid-fluid interfaces.
AIAA Journal , 47(12):2907–2922, 2009. ISSN 0001-1452.[32] R. Saurel and R. Abgrall. A simple method for compressible multifluidflow.
SIAM Journal on Scientific Computing , 21(3), 1999. ISSN 1064-8275.[33] Stefan Schoch, Kevin Nordin-Bates, and Nikolaos Nikiforakis. Aneulerian algorithm for coupled simulations of elastoplastic-solids andcondensed-phase explosives.
Journal of Computational Physics , 252:163–194, 2013. ISSN 0021-9991.[34] P. Clark Souers, Andrew Hernandez, Chris Cabacungan, Raul Garza, LisaLauderbach, SenBen Liao, and Peter Vitello. Air gaps, size e ff ect, andcornerturning in ambient lx17. Propellants, Explosives, Pyrotechnics , 34(1):32–40, 2009. ISSN 0721-3115.[35] Craig M Tarver. Corner turning and shock desensitization experimentsplus numerical modeling of detonation waves in the triaminotrinitroben-zene based explosive lx-17.
The journal of physical chemistry. A , 114(8):2727, 2010. ISSN 1089-5639. c (cid:13) British Crown Owned Copyright 2020 / AWE
36] CraigM. Tarver. Ignition and growth modeling of lx17 hockey puck ex-periments.
Propellants, Explosives, Pyrotechnics , 30(2):109–117, 2005.ISSN 0721-3115.[37] C. W. Wang, T. G. Liu, and B. C. Khoo. A real ghost fluid method for thesimulation of multimedium compressible flow.
SIAM Journal on Scien-tific Computing , 28(1):278–302, 2006. ISSN 1064-8275.[38] F. Xiao, Y. Honma, and T. Kono. A simple algebraic interface capturingscheme using hyperbolic tangent function.
International Journal for Nu-merical Methods in Fluids , 48(9):1023–1040, 7 2005. ISSN 1097-0363.doi: 10.1002 / fld.975. URL https://doi.org/10.1002/fld.975 .[39] David Youngs. An interface tracking method for a 3d eulerian hydrody-namics code. Technical report, AWRE, 1984.[40] Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, JohannesBlaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, DanielGraves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka,Michele Rosso, Samuel Williams, and Michael Zingale. AMReX: aframework for block-structured adaptive mesh refinement. Journal ofOpen Source Software , 4(37):1370, May 2019. doi: 10.21105 / joss.01370.URL https://doi.org/10.21105/joss.01370 . Appendix A. THINC Reconstruction Algorithm
The THINC algorithm models a thermodynamic quantity q in a cell i near an interface as obeying a tanh profile, given by: q i ( x ) THINC = q min + q max + θ tanh β x − x i − x i + − x i − − x . (A.1)Here q min = min( q i − , q i + ), q max = max( q i − , q i + ) − q min , θ = sgn( q i + − q i − ) and β is a parameter that controls the thick-ness of the interface. x is the unknown interface location asa fraction of the cell width, such that when x = x i − , and when x = x i + . This interfacelocation is found by solving q i = ∆ x (cid:90) x i + x i − q i ( x ) THINC d x . (A.2)Using this form also ensures conservation of the variable. Thecell-edge values are computed as q Li = q min + q max + θ A ) (A.3a) q Ri = q min + q max (cid:32) + θ tanh( β ) + A + A tanh( β ) (cid:33) , (A.3b)where A = B cosh( β ) − β ) , B = exp( θβ (2 C − , C = q i − q min + (cid:15) q max + (cid:15) , (A.4)are constants that arise in the integration of Equation A.2. Asmall positive quantity, (cid:15) , is introduced to prevent division byzero. A value of (cid:15) = − is used in this work.The decision of whether to apply this reconstruction isbased on criteria for whether the cell in question is a mixed cell ,defined as satisfying: δ < C < − δ, ( q i + − q i )( q i − q i − ) > , (A.5) where δ is a small positive value. This work uses δ = − .Following this initial reconstruction along the lines of Xiaoet al. [38], the BVD (Boundary Variation Diminishing) algo-rithm of Deng et al. [11] is then applied in conjunction with aMUSCL extrapolation scheme. The premise of this scheme in-volves minimising the variation between cells, hence reducingthe dissipation at interfaces, by comparing the THINC recon-struction with the MUSCL extrapolation and taking whichevergives a lower boundary variation. This is achieved by compar-ing the total boundary variation (TBV) of each scheme, givenby: TBV P = min( | q MUSCL i − , R − q P i , L | + | q P i , R − q MUSCL i + , L | , | q THINC i − , R − q P i , L | + | q P i , R − q THINC i + , L | , | q MUSCL i − , R − q P i , L | + | q P i , R − q THINC i + , L | , | q THINC i − , R − q P i , L | + | q P i , R − q MUSCL i + , L | ) , (A.6)where P stands for either MUSCL or THINC. Therefore thefinal criterion for whether the THINC reconstruction is appliedto a cell is given by: q BVD i = (cid:40) q THINC i if TBV THINC < TBV
MUSCL q MUSCL i otherwise . (A.7)Multidimensionality is accounted for following the proce-dure used by Xiao et al. [38], which takes β = β | n d | + . , (A.8)where β = . | n d | is the magnitude of thecomponent of the interface normal in a given direction. Thisnormal vector is calculated using Youngs’ method [39]. Appendix B. Root Finding procedure