A unified Eulerian framework for multimaterial continuum mechanics

TL;DR

This paper introduces a unified Eulerian framework for simulating interactions among multiple continua, using a common set of equations and a novel interface handling method, enabling accurate and versatile multimaterial simulations.

Contribution

It presents a new multimaterial simulation framework based on the GPR model and a Riemann Ghost Fluid Method that simplifies and unifies the treatment of different material interfaces.

Findings

Accurately reproduces analytical Riemann problem solutions.

Successfully models heat conduction across interfaces.

Demonstrates impact-induced deformation and detonation scenarios.

Abstract

A framework for simulating the interactions between multiple different continua is presented. Each constituent material is governed by the same set of equations, differing only in terms of their equations of state and strain dissipation functions. The interfaces between any combination of fluids, solids, and vacuum are handled by a new Riemann Ghost Fluid Method, which is agnostic to the type of material on either side (depending only on the desired boundary conditions). The Godunov-Peshkov-Romenski (GPR) model is used for modeling the continua (having recently been used to solve a range of problems involving Newtonian and non-Newtonian fluids, and elastic and elastoplastic solids), and this study represents a novel approach for handling multimaterial problems under this model. The resulting framework is simple, yet capable of accurately reproducing a wide range of different physical scenarios. It is demonstrated here to accurately reproduce analytical results for known Riemann problems, and to produce expected results in other cases, including some featuring heat conduction across interfaces, and impact-induced deformation and detonation of combustible materials. The framework thus has the potential to streamline development of simulation software for scenarios involving multiple materials and phases of matter, by reducing the number of different systems of equations that require solvers, and cutting down on the amount of theoretical work required to deal with the interfaces between materials.

Figures

Tables

Table 1. Table 1: EOS parameters for different fluids (using SI units)

	${\rho }_0$	$c_v$	$\gamma$	$p_{\infty }$	$c_s$	$c_t$	$\mu$	$P_r$
Air	1.18	718	1.4	-	55	50	$1.85\times {10}^{-5}$	0.714
Helium	0.163	3127	5/3	-	55	50	$1.99\times {10}^{-5}$	0.688
Water	997	950	4.4	$6\times {10}^8$	1	1	${10}^{-3}$	7

Table 2. Table 2: Initial conditions for the helium bubble test

	$\rho$	$p$	$𝒗$	$A$	$𝑱$
	$1.3333$	$1.5\times {10}^5$		${\left(\frac{1.3333}{1.18}\right)}^{\frac{1}{3}}{I}_3$	$\mathrm{𝟎}$
	$1$	${10}^5$	$\mathrm{𝟎}$	${\left(\frac{1}{1.18}\right)}^{\frac{1}{3}}{I}_3$	$\mathrm{𝟎}$
	$0.1379$	${10}^5$	$\mathrm{𝟎}$	${\left(\frac{0.1379}{0.163}\right)}^{\frac{1}{3}}{I}_3$	$\mathrm{𝟎}$
$0.6\le x\le 1$	$1$	${10}^5$	$\mathrm{𝟎}$	${\left(\frac{1}{1.18}\right)}^{\frac{1}{3}}{I}_3$	$\mathrm{𝟎}$

Table 3. Table 3: Initial conditions for the water-air shock tube test

	$\rho$	$p$	$𝒗$	$A$	$𝑱$
	$1000$	${10}^9$	$\mathrm{𝟎}$	$I_3$	$\mathrm{𝟎}$
$0.7\le x\le 1$	$50$	${10}^5$	$\mathrm{𝟎}$	$\sqrt[3]{50}\cdot I_3$	$\mathrm{𝟎}$

Table 4. Table 4: Initial conditions for the heat conduction test

	$\rho$	$p$	$𝒗$	$A$	$𝑱$
	$2$	$1$	$\mathrm{𝟎}$	$\sqrt[3]{2}\cdot I_3$	$\mathrm{𝟎}$
$x\ge 0$	$0.5$	$1$	$\mathrm{𝟎}$	$\frac{1}{\sqrt[3]{2}}\cdot I_3$	$\mathrm{𝟎}$

Table 5. Table 5: Convergence Rates for the Water-Air Test, PBX-Copper Test, and Aluminium-Vacuum Test

	Water-Air Test		PBX-Copper Test		Aluminium-Vacuum Test
# Cells	Error	Rate	Error	Rate	Error	Rate
50	17.67419797		20.27348796		0.04232714
100	7.98596866	1.146	10.14165497	0.999	0.02180830	0.957
150	5.20404021	1.056	6.57517319	1.069	0.01365700	1.154
200	4.19411427	0.750	4.99100631	0.958	0.01005346	1.065
250	3.32431075	1.042	3.73732754	1.296	0.00793378	1.061
300	2.77425968	0.992	3.10851219	1.010	0.00571356	1.800

Equations

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v_k\right)}{\partial x_k}$$

$$\frac{\partial \left(\rho v_i\right)}{\partial t}+\frac{\partial (\rho v_i{v}_k+p{\delta }_{ik}-{\sigma }_{ik})}{\partial x_k}$$

$$\frac{\partial A_{ij}}{\partial t}+\frac{\partial \left(A_{ik}{v}_k\right)}{\partial x_j}+v_k\left(\frac{\partial A_{ij}}{\partial x_k}-\frac{\partial A_{ik}}{\partial x_j}\right)$$

$$\frac{\partial \left(\rho J_i\right)}{\partial t}+\frac{\partial \left(\rho J_i{v}_k+T{\delta }_{ik}\right)}{\partial x_k}$$

$$\frac{\partial \left(\rho E\right)}{\partial t}+\frac{\partial \left(\rho E{v}_k+\left(p{\delta }_{ik}-{\sigma }_{ik}\right)v_i+q_k\right)}{\partial x_k}$$

$$\{\begin{array}{ll}p={\rho }^2{\frac{\partial E}{\partial \rho }|}_{s,A}& \sigma =-\rho A^T{\frac{\partial E}{\partial A}|}_{\rho ,s}\\ T={\frac{\partial E}{\partial s}|}_{\rho ,A}& \mathbf{q}=T\frac{\partial E}{\partial \mathbf{J}}\end{array}$$

$${\tau }_1=\{\begin{array}{ll}\infty & elasticsolids\\ 0& inviscidfluids\end{array}$$

$$E=E_1\left(\rho ,s\right)+E_2\left(\rho ,s,A,\mathbf{J}\right)+E_3\left(\mathbf{v}\right)$$

$$E_3=\frac{1}{2}{\Vert \mathbf{v}\Vert }^2$$

$$E_2=\frac{c_s{\left(\rho ,s\right)}^2}{4}{\Vert \mathrm{dev}\left(G\right)\Vert }_F^2+\frac{c_t{\left(\rho ,s\right)}^2}{2}{\Vert \mathbf{J}\Vert }^2$$

$$\mathrm{dev}\left(G\right)=G-\frac{1}{3}\mathrm{tr}\left(G\right)I$$

$$c_h=\frac{c_t}{\rho }\sqrt{\frac{T}{c_v}}$$

$$E_1=\frac{p+\gamma p_{\infty }}{\rho \left(\gamma -1\right)}$$

$$E_1=\frac{p_{ref}}{2}\left(\frac{1}{{\rho }_0}-\frac{1}{\rho }\right)+\frac{p-p_{ref}}{{\Gamma }_0{\rho }_0}$$

$$p_{ref}=\frac{c_0^2\left(\frac{1}{{\rho }_0}-\frac{1}{\rho }\right)}{{\left(\frac{1}{{\rho }_0}-s\left(\frac{1}{{\rho }_0}-\frac{1}{\rho }\right)\right)}^2}$$

$$E_1=\frac{p}{\rho \gamma }+\frac{c_0^2}{2{\alpha }^2}\left({\left(\frac{\rho }{{\rho }_0}\right)}^{\alpha }-1\right)\left(\left(1-\frac{2\alpha }{\gamma }\right){\left(\frac{\rho }{{\rho }_0}\right)}^{\alpha }-1\right)$$

$$c_s=b_0{\left(\frac{\rho }{{\rho }_0}\right)}^{\beta }$$

$$E_1=\left(\frac{A}{{\rho }_0{R}_1}{e}^{-\frac{R_1{\rho }_0}{\rho }}+\frac{B}{{\rho }_0{R}_2}{e}^{-\frac{R_2{\rho }_0}{\rho }}\right)+\frac{p-\left(A{e}^{-\frac{R_1{\rho }_0}{\rho }}+B{e}^{-\frac{R_2{\rho }_0}{\rho }}\right)}{\rho {\Gamma }_0}$$

$${\theta }_1$$

$${\theta }_2$$

$$\sigma$$

$$\mathbf{q}$$

$$-\frac{\psi }{{\theta }_1({\tau }_1)}$$

$$-\frac{\rho \mathbf{H}}{{\theta }_2\left({\tau }_2\right)}$$

$$\det \left(A\right)=\frac{\rho }{{\rho }_0}$$

$${\Gamma }_c=\left\{x:f\left(x\right)=c\right\}$$

$$\frac{\partial f}{\partial t}=v\left|\nabla f\right|$$

$$\boldsymbol{P_{j}^{\left(1\right)}}=\begin{cases}\boldsymbol{P_{j}^{n}}&j\leq i\\ \left(\begin{array}[]{ccc}\rho\left(s_{i}^{n},p_{j}^{n},\gamma_{i}^{n}\right)&v_{j}^{n}&p_{j}^{n}\end{array}\right)&j>i\end{cases}$$

$$\boldsymbol{P_{j}^{\left(2\right)}}=\begin{cases}\boldsymbol{P_{j}^{n}}&j\geq i+1\\ \left(\begin{array}[]{ccc}\rho\left(s_{i+1}^{n},p_{j}^{n},\gamma_{i+1}^{n}\right)&v_{j}^{n}&p_{j}^{n}\end{array}\right)&j<i+1\end{cases}$$

$$\rho \left(s,p,\gamma \right)={\left(\frac{p}{s}\right)}^{\frac{1}{\gamma }}$$

$${\mathbf{P}}_{\mathbf{j}}^{\mathbf{n}+1}=\{\begin{array}{ll}{\mathbf{P}}_{\mathbf{j}}^{\left(1\right)}& j\le k\\ {\mathbf{P}}_{\mathbf{j}}^{\left(2\right)}& j>k\end{array}$$

$$\mathbf{P}\left(x,0\right)=\{\begin{array}{ll}{\mathbf{P}}_{\mathbf{L}}& x<0\\ {\mathbf{P}}_{\mathbf{R}}& x>0\end{array}$$

$$\frac{\partial \mathbf{Q}}{\partial t}+\nabla \cdot \mathbf{F}\left(\mathbf{Q}\right)+\mathbf{B}\left(\mathbf{Q}\right)\cdot \nabla \mathbf{Q}=\mathbf{S}\left(\mathbf{Q}\right)$$

$$\frac{\partial \mathbf{Q}}{\partial t}+\nabla \cdot \mathbf{F}\left(\mathbf{Q}\right)+\mathbf{B}\left(\mathbf{Q}\right)\cdot \nabla \mathbf{Q}$$

Full text

The Godunov-Peshkov-Romenski (GPR) model is used for modeling the continua (having recently been used to solve a range of problems involving Newtonian and non-Newtonian fluids, and elastic and elastoplastic solids), and this study represents a novel approach for handling multimaterial problems under this model.

The resulting framework is simple, yet capable of accurately reproducing a wide range of different physical scenarios. It is demonstrated here to accurately reproduce analytical results for known Riemann problems, and to produce expected results in other cases, including some featuring heat conduction across interfaces, and impact-induced deformation and detonation of combustible materials. The framework thus has the potential to streamline development of simulation software for scenarios involving multiple materials and phases of matter, by reducing the number of different systems of equations that require solvers, and cutting down on the amount of theoretical work required to deal with the interfaces between materials.

††journal: Journal of Computational Physics

1 Background

1.1 Multimaterial Models

This study concerns problems featuring immiscible materials. There are many approaches available to solve these problems, broadly including (but not limited to): Lagrangian and Arbitrary-Lagrangian-Eulerian methods [62, 14], volume of fluid methods [52, 45], diffuse interface methods [61], and level-set methods (including the ghost fluid approach [47, 23]).

Solids models tend to come in Lagrangian form, and often these are combined with ALE forms for the fluid phases, so that the fluid meshes may deform to match the deformation of the solid (see, for example, Pin et al. [51]). These schemes tend to be very accurate, but like all Lagrangian schemes, they fail if the meshes become highly contorted. Thus, adaptive remeshing is often necessary. Some authors have coupled a Lagrangian solid scheme with an Eulerian fluid scheme, but extra care must be taken when applying the boundary conditions to the interface, which corresponds to the intersection of the Eulerian and Lagrangian meshes (see Legay et al. [35] for an implementation using level sets, or Fedkiw [24] for a GFM coupling). Some authors, such as Ryzhakov et al. [58], have found success in using the common Lagrangian formulations for the solid, and a reformulated Lagrangian model for the fluid, implementing the necessary adaptive remeshing. Yet another option is to model both the fluid and the solid in an Eulerian framework, although this now necessitates a level set method or volume of fluid method [29] to track the interfaces. Also, these methods are more prone to losing small-scale geometric features of the media, unless methods such as AMR are employed to combat this [30].

Gavrilyuk, Favrie, et al have presented thermodynamically consistent schemes where solid-fluid interfaces are modeled with the diffuse interface method, with transverse velocities found using a ghost fluid method [21]. The fluid is governed by the compressible Euler equations, and the solid by a conservative hyperelastic model. This was later extended to encompass solids conforming to the visco-plastic model of Maxwell type materials [20], and later still to an arbitrary number of interacting hyperelastic solids and fluids governed by the compressible Euler equations [44].

In a recently submitted paper, Michael and Nikiforakis [41] (building on the work of Schoch et al. [63]) couple various Eulerian models of reactive and inert fluids and solids by use of a Riemann Ghost Fluid Method, with the ghost states calculated using specialised mixed-material Riemann solvers for each interaction (see Section 1.3). Whilst these techniques do not suffer from the mesh contortion issues inherent in Lagrangian formulations of continuum mechanics, a fair amount of theoretical work needs to be done to derive analytical relations describing the interactions between every pair of models used.

If it were possible to describe all phases with the same Eulerian model, this method could be used, with only one type of Riemann solver needed to cope with any multiphase problem posed. The GPR model represents such an opportunity. As will be seen, the model also includes terms for heat conduction, which do not appear in the basic formulations of many of the common models used in multiphase systems (e.g. the Euler equations, or the Godunov-Romenski equations of solid mechanics). Heat conduction is often ignored in multiphase modeling, but such a framework based on the GPR model would almost unavoidably include it.

At present, there is no way of dealing with material interfaces in the GPR model, however. In this study, a modification of Barton’s [2] application of Sambasivan and Udaykumar’s Riemann Ghost Fluid Method [59, 60] is devised for the GPR model, enabling the simulation of material interfaces. This new method is tested on a variety of interface problems.

It is interesting to note that de Brauer et al [13, 10] have presented a method for multimaterial modelling of a similar system (including the distortion tensor of the GPR model, but excluding the heat conduction terms). This method is based on level sets, similar to the method presented in this paper. It should be noted, however, that de Brauer et al do not apply their method to the modelling of viscous flows.

The following two subsections outline the theory behind the GPR model and ghost fluid methods. In Section 2 we explore the eigenstructure of the GPR model, and use it to derive a Riemann solver for the Riemann problem at the interfaces between different materials, which is able to incorporate the boundary conditions that we wish to use. Section 3 presents results of fluid-fluid, solid-solid, fluid-solid, and solid-vacuum problems, including some multidimensional cases, and some incorporating heat conduction across the interfaces. Conclusions are drawn in section 4, along with discussion of potential limitations to the method presented here, and ideas for further avenues of enquiry.

1.2 The Model of Godunov, Peshkov and Romenski

The GPR model takes the following form (see [49, 16, 9, 48]):

[TABLE]

[TABLE]

where $\rho$ is density, $\mathbf{v}$ is velocity, $\delta$ is the Kronecker delta, $p$ is pressure, $\sigma$ is the sheer stress tensor, $A$ is the distortion tensor, $\mathbf{J}$ is the thermal impulse vector, $T$ is temperature, $E$ is total energy, and $\mathbf{q}$ is heat flux. $\psi=\frac{\partial E}{\partial A}$ and $\boldsymbol{H}=\frac{\partial E}{\partial\boldsymbol{J}}$, and $\theta_{1}$ and $\theta_{2}$ are positive functions (given below for the problems at hand). Additionally, we have:

[TABLE]

where $s$ is the entropy of the system.

The GPR model represents the same set of equations as the model of elastoplastic deformation originally proposed by Godunov and Romenski. Peshkov and Romenski first subsequently proposed that these are the equations of motion for an arbitrary continuum - not just a solid. In doing so, they were able to apply the model to fluids too. Note that material elements have not only finite size, but also internal structure (encoded in the distortion), unlike in previous continuum models.

The idea of $\tau_{1}$ - the strain dissipation time - has its roots in Frenkel’s “particle settled life time” (see [25]). $\tau_{1}$ represents a continuous analogue of Frenkel’s object. It can be thought of as the characteristic time taken for a particle to move by a distance roughly the same as the particle’s size. Thus, the typical time taken for a material element to rearrange with its neighbors is characterized by $\tau_{1}$ . As long as a continuum description is appropriate for the material at hand, it is thus that the GPR model seeks to describe all three major phases of matter. For example, we have the following relations:

[TABLE]

The equation governing $\boldsymbol{J}$ - and its contribution to the system’s total energy - are derived from Romenski’s model of hyperbolic heat transfer, (see [40, 57]). These concepts were later implemented in [56, 53]. The entropy flux is the derivative of the specific internal energy with respect to $\boldsymbol{J}$, and it is in this way that $\mathbf{J}$ is defined (as the variable conjugate to the entropy flux). As remarked by Romenski, it is more convenient to evolve $\boldsymbol{J}$ and $E$ than $\mathbf{q}$ or the entropy flux, and thus the equations take the form given here. Similarly to $\tau_{1}$, $\tau_{2}$ is a relaxation time, characterizing the average speed of relaxation of thermal impulse due to heat exchange between neighboring particles.

$E$ must be specified to close the system. The energy contains contributions from the micro-, meso-, and macro-scale:

[TABLE]

$E_{3}$ is the usual specific kinetic energy per unit mass:

[TABLE]

$E_{2}$ takes the quadratic form:

[TABLE]

where

[TABLE]

$\left\|\cdot\right\|_{F}$ is the Frobenius norm111The Frobenius norm is defined by: $\left\|X\right\|_{F}=\sqrt{\sum_{i,j}\left|X_{ij}\right|^{2}}$ and $G=A^{T}A$ is the Gramian matrix of the distortion222$G$ is known as the Finger tensor in the solid mechanics community, and $\operatorname{dev}\left(G\right)$ is the deviator (trace-free part) of $G$. $c_{s}$ is the characteristic velocity of transverse perturbations. $c_{t}$ is related to the characteristic velocity of propagation of heat waves333Note that [16] denotes this variable by $\alpha$, which is avoided here due to a clash with a parameter of one of the equations of state used.:

[TABLE]

In previous studies (e.g. [16, 9]), $c_{t}$ has been taken to be constant, as it will be in this study.

In this study, $E_{1}$ is taken to be one of the following forms:

1. The stiffened gas EOS:

[TABLE]

2. The shock Mie-Gruneisen EOS:

[TABLE]

where

[TABLE]

3. The Godunov-Romenski EOS (see [4]):

[TABLE]

with

[TABLE]

4. The JWL EOS:

[TABLE]

Tabulated equations of state are common-place in the field (see [11, 36, 50] for a range of use cases). There is no a priori reason why they cannot be used for $E_{1}$ under the framework presented in this study, in the same manner as other hydrodynamic systems (e.g. see [8, 28]). This is out of scope of this paper, however, and presents an avenue of future research. Note that $E_{1},c_{s},c_{t}$ are permitted to depend upon $\rho,p$ instead of $\rho,s$ (as is the case in this study), or indeed $\rho,T$ if the material requires it (such as materials whose shear modulus depend on temperature).

The following forms are taken:

[TABLE]

[TABLE]

The justification of these choices is that classical Navier–Stokes–Fourier theory is recovered in the stiff limit $\tau_{1},\tau_{2}\rightarrow 0$ (see [16]). The rules for power-law fluids and elastoplastic solids are based on material from [33] and [6], respectively.

Finally, it is straightforward to verify that as a consequence of (2) we have the following relations:

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

Finally, the following constraint holds (see [49]):

[TABLE]

See [33] for an interpretation of the physical meaning of the relaxation times $\tau_{1},\tau_{2}$ and thermal impulse vector $\boldsymbol{J}$.

1.3 Ghost Fluid Methods

Ghost fluid methods, combined with level set methods, are used to model the evolution of interfaces between different materials. They are detailed here, as it is with such a method that this study proposes to model the interfaces between different materials described by the GPR model.

1.3.1 Level Set Methods

Given a scalar function $f$ on $\mathbb{R}^{n}$, the level set of $f$ at level $c$ is defined as:

[TABLE]

Given normal direction speed $v$, $f$ is advected according to the level set equation [46]:

[TABLE]

The advection of a point in a fluid with velocity $v$ can be modeled by taking $f=\left|\mathbf{x-x_{0}}\right|$ where $\mathbf{x_{0}}$ is the position of the point at time $t=0$, and tracking $\Gamma_{0}$. (20) is solved by an appropriate numerical method. The numerical methods used in this study are described in Chapter 2. $f$ will usually have to be renormalized at every time step, to avoid unwanted distortions such as becoming a multivalued function.

1.3.2 The Original Ghost Fluid Method

The Original Ghost Fluid Method of Fedkiw et al. [22] (an adaptation of the work of Glimm et al. [26]) is a numerical method for the Euler equations for simulating interfaces between multiple materials. The primitive variables for the Euler equations in 1D are given by \boldsymbol{P}=\left(\begin{array}[]{ccc}\rho&v&p\end{array}\right)^{T}.

Suppose the interface between two fluids is modeled on spatial domain $\left[0,1\right]$, divided into $N$ cells with width $\Delta x=\frac{1}{N}$. Let the time step be $\Delta t$ and let $\boldsymbol{P_{i}^{n}}$ be the set of primitive variables in cell $i$ at time $t_{n}=n\Delta t$. Let the level set function $f$ have root $x_{n}$ where x_{n}\in\left[\text{\left(i+\frac{1}{2}\right)}\Delta x,\text{\left(i+\frac{3}{2}\right)}\Delta x\right]. Thus, at time $t_{n}$ the interface lies between the cells with primitive variables $\boldsymbol{P_{i}^{n},P_{i+1}^{n}}$. Define two sets of primitive variables:

[TABLE]

[TABLE]

where:

[TABLE]

All cells in $\boldsymbol{P^{\left(1\right)}}$ to the left of the interface have the same state variables as those of $\boldsymbol{P^{n}}$. All cells to the right have the same pressure and velocity as their counterparts in $\boldsymbol{P^{n}}$, but the same entropy as $\boldsymbol{P_{i}^{n}}$. This affects their density. The situation is analogous for $\boldsymbol{P^{\left(2\right)}}$. This is demonstrated in Figure 1 on page 1.

$\boldsymbol{P^{\left(1\right)}},\boldsymbol{P^{\left(2\right)}}$ are stepped forward by time step $\Delta t$ using a standard Eulerian method. $f$ is advected using (20), taking the velocity in each cell to be that of $\boldsymbol{P^{n}}$. Now let $f\left(x_{n+1}\right)=0$ where $x_{n+1}\in\left[\left(k+\frac{1}{2}\right)\Delta x,\left(k+\frac{3}{2}\right)\Delta x\right]$ for some $k$. Define:

[TABLE]

The rationale behind the original GFM is that in most applications, pressure and velocity are continuous across the interface, and thus the ghost cells may take the real pressure and velocity values. Entropy is generally discontinuous at a contact discontinuity, resulting in large truncation errors if a standard finite difference scheme is used to solve the system. Thus, entropy is extrapolated as a constant from the interface boundary cell into the ghost region.

Fedkiw et al. advised to use the isobaric fix technique. This involves setting the entropy of cell $i$, and all cells in the right ghost region, to that of cell $i-1$, and setting the entropy of cell $i+1$, and all cells in the left ghost region, to that of cell $i+2$. This is demonstrated in Figure 2 on page 2.

Effectively, the ghost regions behave like they are composed of the same fluid as the regions they extend (as they have the same entropy), facilitating calculation of the next time step, but they have the same pressure and velocity profiles as the real fluids they replace, meaning the boundary conditions at the interface are upheld.

1.3.3 The Riemann Ghost Fluid Method

The Riemann Problem in its general form is the solution of the following initial value problem. Given a set of variables $\boldsymbol{P}$ dependent on space and time, and a hyperbolic set of equations which govern their spatio-temporal evolution, $\boldsymbol{P}\left(x,t\right)$ is sought for $t>0$, given the initial condition:

[TABLE]

This problem is denoted by $RP\left(\boldsymbol{P_{L}},\boldsymbol{P_{R}}\right)$. Exact solvers exist for the Riemann Problem for various sets of governing equations, such as the Euler equations [65], the equations of non-linear elasticity [5], or the shallow water equations [1], among others. There also exist approximate solvers for general conservative [43, 38] or non-conservative [15] hyperbolic systems of PDEs. The references given here form a very small sample of the work that has been done in this area.

The solution of the Riemann Problem usually takes the form of a set of waves, between which $\boldsymbol{P}$ is constant. The waves can either be a contact discontinuity (across which pressure and velocity are continuous), a shock (across which all variables may be discontinuous), or a rarefaction (along which the variables vary continuously between their values on either side of the wave). The number and form of the waves are determined by the governing equations and the initial conditions. The states of the variables either side of the contact discontinuity in the middle are known as the star states. This qualitative description is depicted in Figure 3 on page 3.

Liu et al. [37] demonstrated that the original GFM fails to resolve strong shocks at material interfaces. This is because the method effectively solves two separate single-fluid Riemann problems. The waves present in these Riemann problems do not necessarily correspond to those in the real Riemann problem across the interface. The Riemann Ghost Fluid Method of Sambasivan et al. [59] aims to rectify this.

Given $\boldsymbol{P^{n}}$ and $x_{n}\in\left[\left(i+\frac{1}{2}\right)\Delta x,\left(i+\frac{3}{2}\right)\Delta x\right]$, the ghost cells for fluid 1 are populated with the left star state of $RP\left(\boldsymbol{P_{i-1}^{n}},\boldsymbol{P_{i+2}^{n}}\right)$, and the ghost cells for fluid 2 are populated with the right star state. $RP\left(\boldsymbol{P_{i-1}^{n}},\boldsymbol{P_{i+2}^{n}}\right)$ is solved rather than $RP\left(\boldsymbol{P_{i}^{n}},\boldsymbol{P_{i+1}^{n}}\right)$, as $\boldsymbol{P_{i}^{n}},\boldsymbol{P_{i+1}^{n}}$ often contain errors generated by the fact that they lie on the material interface. $\boldsymbol{P^{n+1}}$ is then generated as before from the newly formed $\boldsymbol{P^{\left(1\right)}},\boldsymbol{P^{\left(2\right)}}$. This process is demonstrated in Figure 4 on page 4.

1.4 Finite Volume Scheme

The presented system is solved with a finite volume scheme in this study, using a regularized square meshing. Once the values of the ghost cells have been calculated, we require a numerical method to calculate the value of the grid at the next time step. In this study, the finite volume method presented in [32, 33] was used. It is outlined here for completeness.

Note that (1a), (1b), (1c), (1d), (1e) can be written in the following form:

[TABLE]

As described in [64], a method to solve a system such as this is to solve the following subsystems:

[TABLE]

[TABLE]

Let $H^{\delta t},S^{\delta t}$ be the operators that take data $\boldsymbol{Q}\left(x,t\right)$ to $\boldsymbol{Q}\left(x,t+\delta t\right)$ under systems (27a) and (27b) respectively. A second-order scheme (in time) for solving the full set of PDEs over time step $\left[0,\Delta t\right]$ is obtained by calculating $\boldsymbol{Q_{\Delta t}}$ using a Strang splitting:

[TABLE]

In this study, the homogeneous subsystem is be solved using a WENO reconstruction of the data, followed by a finite volume update, and the temporal ODEs will be solved with appropriate ODE solvers. The WENO method was chosen due to the arbitrarily high-order spatial reconstructions it is able to produce.

Noting that $\frac{d\rho}{dt}=0$ over the ODE time step, the operator $S$ entails solving the following systems:

[TABLE]

[TABLE]

These systems are solved separately, and thus he second-order Strang splitting becomes:

[TABLE]

where $D^{\delta t},T^{\delta t}$ are the operators solving the distortion and thermal impulse ODEs respectively, over time step $\delta t$.

The constraint (18) is enforced by rescaling the singular values of the distortion all by the same factor at each timestep, to ensure that $\det A=\frac{\rho}{\rho_{0}}$.

1.4.1 The Homogeneous System

A WENO reconstruction of the cell-averaged data is performed at the start of the time step (as described in [19]). Focusing on a single cell $C_{i}$ at time $t_{n}$, we have $\boldsymbol{w^{n}}\left(\boldsymbol{x}\right)=\boldsymbol{w^{n}}_{p}\Psi_{p}\left(\boldsymbol{\chi}\left(\boldsymbol{x}\right)\right)$ in $C_{i}$ where $\Psi_{p}$ is a tensor product of basis functions in each of the spatial dimensions. The flux in $C$ is approximated by $\boldsymbol{F}\left(\boldsymbol{x}\right)\approx\boldsymbol{F}\left(\boldsymbol{w}_{p}\right)\Psi_{p}\left(\boldsymbol{\chi}\left(\boldsymbol{x}\right)\right)$. $\boldsymbol{w}_{p}$ are stepped forwards half a time step using the update formula:

[TABLE]

i.e.

[TABLE]

where $\boldsymbol{\chi_{p}}$ is the node corresponding to $\Psi_{p}$.

Integrating (27a) over $C$ gives:

[TABLE]

where

[TABLE]

[TABLE]

where $V$ is the volume of $C$ and $\boldsymbol{Q^{-},Q^{+}}$ are the interior and exterior extrapolated states at the boundary of $C$, respectively.

Note that (27a) can be rewritten as:

[TABLE]

where $\boldsymbol{M}=\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{Q}}+\boldsymbol{B}$. Let $\boldsymbol{n}$ be the normal to the boundary at point $\boldsymbol{s}\in\partial C$. For the GPR model, $\hat{M}=\boldsymbol{M}\left(\boldsymbol{Q}\left(\boldsymbol{s}\right)\right)\cdot\boldsymbol{n}$ is a diagonalizable matrix with decomposition $\hat{M}=\hat{R}\hat{\Lambda}\hat{R}^{-1}$ where the columns of $\hat{R}$ are the right eigenvectors and $\hat{\Lambda}$ is the diagonal matrix of eigenvalues. Define also $\boldsymbol{\hat{F}}=\boldsymbol{F}\cdot\boldsymbol{n}$ and $\hat{B}=\boldsymbol{B}\cdot\boldsymbol{n}$. Using these definitions, the interface terms arising in the FV formula have the following form:

[TABLE]

$\tilde{M}$ is chosen to either correspond to a Rusanov/Lax-Friedrichs flux (see [64]):

[TABLE]

or a Roe flux (see [18]):

[TABLE]

or a simplified Osher–Solomon flux (see [18, 17]):

[TABLE]

where

[TABLE]

$\tilde{B}$ takes the following form:

[TABLE]

$\boldsymbol{P_{i}^{n+\frac{1}{2}}},\boldsymbol{D_{i}^{n+\frac{1}{2}}}$ are calculated using an $N+1$-point Gauss-Legendre quadrature, replacing $\boldsymbol{Q}\left(\boldsymbol{x},t_{n+\frac{1}{2}}\right)$ with $\boldsymbol{w^{n+\frac{1}{2}}}\left(\boldsymbol{x}\right)$.

1.4.2 The Thermal Impulse ODEs

In [32] it was shown that (29b) has the following analytical solution:

[TABLE]

where

[TABLE]

[TABLE]

1.4.3 The Distortion ODEs

Take the following singular value decomposition:

[TABLE]

Denote the singular values of $A$ by $a_{1},a_{2},a_{3}$ and define:

[TABLE]

Define also:

[TABLE]

[TABLE]

In [32] it was shown that for Newtonian fluids, after timestep $\Delta t$, $x_{i}$ become:

[TABLE]

where

[TABLE]

[TABLE]

These new values for $x_{i}$ are used to calculate the value of $A$ after timestep $\Delta t$.

In [33] it was shown that for elastoplastic solids, the same procedure can be undertaken, where now:

[TABLE]

where

[TABLE]

[TABLE]

1.4.4 Time Step

Let $\Lambda_{i}^{n}$ be the set of eigenvalues of the GPR system evaluated at $\boldsymbol{Q_{i}^{n}}$ (given explicitly in Section 2.1). $C_{cfl}<1$ is a constant (usually taken to be $0.9$, unless the problem being simulated is particularly demanding, requiring a lower value). The eigenvalues determine the speed of propagation of information in the solution to the Riemann Problem at the cell interfaces, and the time step is chosen to ensure that the characteristics do not enter into other cells between $t_{n}$ and $t_{n+1}$:

[TABLE]

2 A Riemann Ghost Fluid Method for the GPR Model

2.1 Eigenstructure of the GPR Model

Take a hyperbolic system of the following form, noting that the GPR model takes this form:

[TABLE]

Let $\left\{\boldsymbol{l_{i}}\right\}$ be the set of left eigenvectors of the system matrix $M$, in other words:

[TABLE]

Along characteristics corresponding to $\lambda_{i}$, we have:

[TABLE]

The second line of this equation is crucial to Riemann Ghost Fluid Method presented in this study. As such, we must now investigate the eigenvalues and eigenvectors of $M$.

2.1.1 Eigenvalues

Considering the primitive system matrix (142), it is clear that the eigenvalues of the GPR system in the first spatial axis consist of $v_{1}$ repeated 8 times, along with the roots of:

[TABLE]

where

[TABLE]

[TABLE]

By the properties of block matrices444If $A$ is invertible, \det\left(\begin{array}[]{cc}A&B\\ C&D\end{array}\right)=\det\left(A\right)\det\left(D-CA^{-1}B\right), the remaining eigenvalues are $v_{1}$ and the roots of $\left|\left(v_{1}-\lambda\right)^{2}I-\Xi_{1}\Xi_{2}\right|=0$. Thus, $\lambda_{i}=v_{1}\pm\sqrt{\tilde{\lambda_{i}}}$ where the $\tilde{\lambda_{i}}$ are the eigenvalues of the following matrix:

[TABLE]

where $\Omega$ is given shortly. Similar results hold for the other two spatial directions. In general it is not possible to express the eigenvalues of $\Xi$ in terms of the eigenvalues of its submatrices. Note, however, that if $c_{t}=0$ then one of the eigenvalues is [math] and the remaining eigenvalues can be found analytically, using the form given in the appendix of [16].

It is straightforward to verify the following:

[TABLE]

The quantity $\Omega$ is named here the acoustic tensor, due to its similarity to the acoustic tensor described in [4]:

[TABLE]

where $E_{ij}^{d}=\delta_{id}\delta_{jd}$.

2.1.2 Eigenvectors (with Heat Conduction)

By hyperbolicity of the system, $\Xi$ can be expressed as:

[TABLE]

where $D$ is a diagonal matrix with positive diagonal entries. The eigenvectors corresponding to $\lambda_{i}=v_{1}\pm\sqrt{\tilde{\lambda_{i}}}$ take the form \left(\begin{array}[]{cccc}\hat{u}&0_{6}&\tilde{u}&0_{2}\end{array}\right)^{T} where $\boldsymbol{\hat{u}}\in\mathbb{R}^{5},\boldsymbol{\tilde{u}}\in\mathbb{R}^{4}$ satisfy:

[TABLE]

Thus, $\Xi_{2}\boldsymbol{\tilde{u}}=\pm\sqrt{\tilde{\lambda_{i}}}\boldsymbol{\hat{u}}$ and $\Xi_{1}\boldsymbol{\hat{u}}=\pm\sqrt{\tilde{\lambda_{i}}}\boldsymbol{\tilde{u}}$. Combining these results, $\Xi\boldsymbol{\tilde{u}}=\tilde{\lambda_{i}}\boldsymbol{\tilde{u}}$. Thus, $\boldsymbol{\tilde{u}}$ is a right eigenvector of $\Xi$ and, taking the form $Q^{-1}\boldsymbol{e_{i}}$ for some $i=1\ldots 4$.

The four eigenvectors corresponding to eigenvalues of the form $v_{1}+\sqrt{\tilde{\lambda_{i}}}$ are columns 1-4 of matrix $R$ in (66). Those corresponding to eigenvalues of the form $v_{1}-\sqrt{\tilde{\lambda_{i}}}$ are columns 5-8. By inspection (using the system matrix (142)), it can be verified that the remaining 9 eigenvectors (corresponding to eigenvalue $v_{1}$) are the remaining columns.

Note that the index $d$ appearing in these representations should be taken as $1,2,3$ for eigenvectors in directions $x,y,z$, respectively. $0_{m,n}$ is defined to be the 0-matrix of shape $\left(m,n\right)$ and $I_{n}$ the identity matrix of size $n$.

[TABLE]

where

[TABLE]

[TABLE]

A similar analysis yields the left eigenvectors as the rows of (68).

[TABLE]

2.1.3 Eigenvectors (without Heat Conduction)

If the system does not include the heat conduction terms, the eigenstructure of the system matrix changes. $\Xi_{1},\Xi_{2},\Xi$ now take the following values:

[TABLE]

[TABLE]

[TABLE]

Using the eigendecomposition $\Xi=Q^{-1}D^{2}Q$ as before, we have:

[TABLE]

By considering their products with the first 8 columns of $R$, two of the left eigenvectors corresponding the the 7th and 8th right eigenvectors must come in the form of the rows of the following matrix:

[TABLE]

where $W\in\mathbb{R}^{2,5}$ and $X,Y,Z\in\mathbb{R}^{2,3}$, and:

[TABLE]

[TABLE]

Solving this system for $X,Y,Z$:

[TABLE]

[TABLE]

Define:

[TABLE]

where

[TABLE]

[TABLE]

By the properties of block matrices:

[TABLE]

By the orthonormality of eigenvectors, we must have:

[TABLE]

Thus, it is straightforward to confirm that:

[TABLE]

[TABLE]

Thus, we have:

[TABLE]

Finally, combining the preceding results with (73), we have:

[TABLE]

2.2 Solving the Riemann Problem

Barton et al. have presented an RGFM for the equations of non-linear elasticity [7, 3]. Owing to the similarity of the structure of the non-linear elasticity equations to those of the GPR model (differing only in the presence of source terms, the form of the shear stress tensor, and possibly the EOS), their method is built upon here. The resulting method is named the GPR-RGFM.

The Riemann Problem of the GPR model takes the form demonstrated in Figure 5 on page 5. Assuming all waves are distinct, there are four waves on either side of the contact discontinuity. On each side, one wave corresponds to the thermal impulse (manifesting as a heat wave) and the other three correspond to the distortion components in the axis in which the Riemann Problem is considered (manifesting as two shear waves and one longitudinal pressure wave). It is important to note that - owing to the source terms - the star states are not constant in the spacetime region in which they reside, so the method presented here produces only an approximation to them.

The method is presented here along the first spatial axis. It can easily be adapted along any axis by taking the components of all relevant vector quantities (velocity, distortion, and thermal impulse) in the direction normal to the interface.

Denote the vector of primitive variables by $\boldsymbol{P}$. Take the set of left eigenvectors $L$ (68) with eigenvalues $\left\{\lambda_{i}\right\}$ . From (54), we have the standard set of relations along characteristics (curves along which $\frac{dx}{dt}=\lambda_{i}$):

[TABLE]

In what follows, we enact an operator splitting of the two processes present in the system (94):

[TABLE]

[TABLE]

$\boldsymbol{P^{*K}}$ is now sought, where $K=L$ or $K=R$, denoting the left or right sides of the interface, respectively. Take the following linearization:

[TABLE]

13 relations from (95a) are taken: 4 regarding the 4 sets of characteristics traveling into the contact discontinuity from side $K$ (with speeds greater or less than $v$ for $K=L$ or $K=R$, respectively), and 9 relating to the contact discontinuity itself. This is demonstrated in Figure 6 on page 6. 4 more relations must be derived to solve the system for $\boldsymbol{P^{*K}}$.

Define the total stress tensor as:

[TABLE]

The values of $\Sigma,T$ under variables $\boldsymbol{P^{*K}}$ are obtained by expanding the following Taylor series:

[TABLE]

[TABLE]

Thus, we have:

[TABLE]

[TABLE]

These are the extra required relations. Thus we have:

[TABLE]

where $\hat{L}^{K}$ takes the form found in (115), with $\xi=-1$ for $K=R$ and $\xi=1$ for $K=L$, and:

[TABLE]

The inverse of $\hat{L}^{K}$ takes the form found in (116).

$\hat{L}^{K},\hat{L}^{K}{}^{-1}$ are evaluated at $\boldsymbol{P^{K}}$. It remains to find expressions for $\Sigma^{*}$ and $T^{*}$ in terms of $\boldsymbol{P^{L}},\boldsymbol{P^{R}}$ to close the system. The obtained values depend on the boundary conditions chosen, as explained below.

2.2.1 Boundary Conditions

Stick Boundary Conditions

The following boundary conditions are used:

[TABLE]

[TABLE]

Taking the relevant rows of $\boldsymbol{P^{*K}}=\boldsymbol{P^{K}}+\hat{L}^{K}{}^{-1}\boldsymbol{c^{K}}$:

[TABLE]

Thus:

[TABLE]

Slip Boundary Conditions

The following boundary conditions are used:

[TABLE]

[TABLE]

Taking the relevant rows of $\boldsymbol{P^{*K}}=\boldsymbol{P^{K}}+\hat{L}^{K}{}^{-1}\boldsymbol{c^{K}}$:

[TABLE]

where

[TABLE]

Thus:

[TABLE]

where

[TABLE]

Vacuum Boundary Conditions

The following boundary conditions are used:

[TABLE]

[TABLE]

Taking the relevant row of $\boldsymbol{P^{*K}}=\boldsymbol{P^{K}}+\hat{L}^{K}{}^{-1}\boldsymbol{c^{K}}$:

[TABLE]

As $q_{1}^{*}=0$ implies that $J_{1}^{*}=0$, we have:

[TABLE]

Iteration

(100) is solved for $\boldsymbol{P^{*K}}$, which is taken to be the star state if the following conditions are satisfied:

[TABLE]

[TABLE]

where

[TABLE]

These convergence criteria are chosen so that the variables required to be less than $TOL$ are dimensionless. At every iteration, (95b) is solved using the ODE solvers presented in [32, 33].

2.2.2 Linear Conditions with Heat Conduction

We now obtain $\hat{L}^{K}$ and its inverse in order to solve (100). Replacing the first four lines of (68) with the conditions (99a), (99b), we have:

[TABLE]

Thus, the inverse of the left-eigenvector matrix becomes:

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

and also:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By inversion of block matrices555\left(\begin{array}[]{cc}A&B\\ C&D\end{array}\right)^{-1}=\left(\begin{array}[]{cc}\left(A-BD^{-1}C\right)^{-1}&-\left(A-BD^{-1}C\right)^{-1}BD^{-1}\\ -D^{-1}C\left(A-BD^{-1}C\right)^{-1}&D^{-1}\left(I+C\left(A-BD^{-1}C\right)^{-1}BD^{-1}\right)\end{array}\right):

[TABLE]

where

[TABLE]

2.2.3 Linear Conditions without Heat Conduction

If the heat conduction terms are dropped from the GPR model, the eigenstructure of the system changes, along with the solution of the linear conditions. $\Xi$ retains the same definition, but is now a $3\times 3$ matrix (comprising the top-left corner of $\Xi$ under heat conduction). Thus, $Q,D$ are also $3\times 3$ matrices. Taking the eigenvectors (93), the linear conditions become:

[TABLE]

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

By inversion of block matrices:

[TABLE]

where

[TABLE]

3 Results

The GPR-RGFM is now assessed. The first fives tests in this chapter are Riemann problems that have appeared elsewhere in the literature. Reference solutions to these problems have been calculated by various methods, as described for each test individually. The sixth test is new; it assess the ability of the GPR-RGFM to correctly model heat conduction across interfaces. The last two tests are well-known two-dimensional problems, to demonstrate the applicability of the method to multiple dimensions. The stiffened gas EOS parameters for three commonly-used fluids are given in Table 1 on page 1.

3.1 Helium Bubble

The interface between two different gases is now modeled. As in Test B of Wang et al. [66], a bubble of helium - surrounded by air - initially occupies the region $x\in\left[0.4,0.6\right]$. A shock front in the air, initially at $x=0.05$, travels towards the helium bubble. The initial conditions are given in Table 2 on page 2 and the EOS parameters for each material are given in Table 1 on page 1.

200 cells are used. Reference solutions are computed using the exact solver for mixed ideal gas Riemann problems under the Euler equations (presented in [64]). The results for times $t=7\times 10^{-4}$ and $t=14\times 10^{-4}$ are displayed in Figure 7 on page 7. In the former, the shock is about to hit the helium bubble (corresponding to the region of low density). In the latter, the shock has traveled through the helium bubble, compressing it slightly, and the bubble itself has moved almost 0.1 spatial units to the right. There is good correspondence with the results in [66] and the sharp discontinuity in density is maintained.

3.2 Water-Air Shock Tube

This test comprises an interface between water and air, with initial data taken from Chinnayya et al. [12] (see Table 3 on page 3). The aim of this test is to evaluate the ability of the GPR-RGFM at capturing interfaces between qualitatively different fluids. The water is initially at high pressure, and the air at atmospheric pressure. Due to the large difference in state variables and qualitative characteristics of the two fluids, this is an example of a test with which the original GFM for the Euler equations does not perform well.

The results using the GPR-RGFM with 200 cells are shown in Figure 8 on page 8, along with the exact solution to the Euler equations (computed using the extension to the stiffened gas equations of the exact Riemann solver presented in [64]). As can be seen, the material interface is captured well, with the correct intermediate density found by the numerical method.

3.3 PBX9404-Copper Shock Tube

This test is taken from [7], with the aim of testing the ability of the GPR-RGFM to model interfaces between fluids and solids. High pressure, reacted PBX9404 is in contact with copper at position $x=0.5$ on domain $x\in\left[0,1\right]$, with both materials initially at rest. The pressure of the PBX is initially $18.9GPa$, and the entropy of the copper is initially [math]. The PBX follows an ideal gas EOS, with parameters $\rho_{0}=1840$, $\gamma=2.85$, $c_{s}=1$, $\mu=10^{-2}$. The copper follows the Godunov-Romenski EOS, with parameters $\rho_{0}=8930$, $c_{v}=390$, $T_{0}=300$, $c_{0}=3939$, $\alpha=1$, $\beta=3$, $\gamma=2$, $b_{0}=2141$. The test is run until time $t=0.5\times 10^{-6}$, using 500 cells.

The exact solution to this test is calculated using the iterative solver described in [4]. The error in the wavespeeds is calculated from the residual error in the traction and velocities across the central contact, as required by the Rankine–Hugoniot conditions and boundary conditions. The wavespeeds are found by iteratively reducing the residual using the Newton–Raphson method.

Plots for density, velocity, and total stress are given in Figure 9 on page 9. As can be seen, the GPR-RGFM is able to reproduce the solution to high fidelity, with a perfectly sharp discontinuity in the density, and a very well resolved discontinuity in the total stress.

3.4 Aluminium in Vacuum

This test is taken from [3]. The initial conditions of the test consist of a slab of aluminium, initially with velocity \left(\begin{array}[]{ccc}2&0&0.1\end{array}\right), meeting a vacuum at point $x=0.5$, on the domain $x\in\left[0,1\right]$. The distortion of the aluminium is initially given by:

[TABLE]

The initial density of the aluminium is thus given as $\rho=\rho_{0}\det\left(A\right)$. The aluminium is modeled using the Godunov-Romenski EOS, with parameters $\rho_{0}=2.71$, $c_{v}=9\times 10^{-4}$, $T_{0}=300$, $c_{0}=5.037$, $\alpha=1$, $\beta=3.577$, $\gamma=2.088$, $b_{0}=3.16$, $c_{t}=2$, $\kappa=204$.

The test was run until time $t=0.06$, using 500 cells. The results of solving this problem with the GPR-RGFM, not including thermal conduction (as in [3]), are given in Figure 10 on page 10. The results of solving the problem, including thermal conduction, are given in Figure 11 on page 11. The exact solutions are calculated using the iterative method presented in [4], as described in the previous test.

As can be seen, in both cases, the GPR-RGFM is able to accurately capture the longitudinal wave and the two transverse shock waves that propagate to the left side of the initial point of contact. Note that at $t=0.06$, the vacuum occupies the region $\left[\sim 0.65,1\right]$. As this region is empty, the plots in the aforementioned figures are shown over the interval $\left[0,0.7\right]$, to give greater resolution to the region of interest.

Without thermal conduction, the interface suffers from a “heating error” of the same kind discussed in [3], manifesting itself as a slight undershoot in the density of the metal at the interface. Note that by incorporating thermal conduction into the numerical method, this heating error completely disappears, without the use of an entropy fix (as in [3]). It must be noted that, in this case, the waves in the state variables now appear to be slightly more diffused than the reference solution. This is the expected effect of incorporating the phenomenon of thermal conduction into this physical problem.

3.5 Heat Conduction in a Gas

This test is based on the Heat Conduction in a Gas Test of Dumbser et al. [16]. Two ideal gases at different temperatures are initially in contact at position $x=0$. The initial conditions for this problem are given in Table 4 on page 4.

The material parameters are taken to be: $\gamma=1.4$, $c_{v}=2.5$, $\rho_{0}=1$, $p_{0}=1$, $c_{s}=1$, $c_{t}=1$, $\mu=10^{-2}$, $\kappa=10^{-2}$. An interface is initially placed between the two volumes of air at $x=0.5$. The final time is taken to be $t=1$, and 200 cells are used. Results are displayed in Figure 12 on page 12, using the results from [16] as a reference. The material interface is denoted by a dashed vertical line.

The temperature curve generated using the GPR-RGFM matches very well the reference solution. The interface has moved to $x=0.53756$, as is to be expected, as the cooler gas on the left expands as it heats up, and the hotter gas on the right contracts as it cools. Initially, the mass of the left volume is 1 and the right volume is 0.25. At $t=1$, these masses are 0.9997 and 0.2503, respectively. Thus, mass on either side is conserved to a good approximation. Although the GPR-RGFM results for the heat flux match the reference solution well over most of the domain, there are aberrations in a small region around the interface. Although this doesn’t affect the temperature curve, these discrepancies are undesirable, and possible methods to rectify them are discussed in Section 4.

3.6 Taylor Bar

This follows a similar form to that found in [9, 39]. A bar of aluminium of dimensions $100\times 500$ travels towards a solid wall at speed $0.015$. The surrounding environment is a vacuum. The aluminium bar is modelled by the shock Mie-Gruneisen equation of state, with parameters $\rho_{0}=2.785,$ $c_{v}=9\times 10^{-4}$, $c_{0}=0.533$, $\Gamma_{0}=2$, $s=1.338$. The aluminium also follows a plasticity law with parameters $b_{0}=0.305$, $\sigma_{Y}=0.003$, $\tau_{0}=1$, $n=20$. The domain has dimensions $300\times 510$, with $\Delta x,\Delta y=1$.

The density and plastic deformation of the bar at times $t=0.0025$ and $t=0.005$ are shown in Figure 13 on page 13. Unfortunately there are no experimental results for this test, but the reader is asked to note the good agreement here with the results found in [39]. In that study, the boundary between the bar and the vacuum is captured using a Lagrangian scheme, and it is reassuring that the same behaviour is captured here with a characteristically different numerical method.

3.7 Aluminum Plate Impact

This test follows the form found in Michael & Nikiforakis [42] (based on the original formulation found in [31]). An aluminum projectile impacts upon an aluminum plate at speed $400$. The domain is $\left[0,0.03\right]\times\left[0,0.04\right]$, with the projectile initially occupying $\left[0.001,0.006\right]\times\left[0.014,0.026\right]$, and the plate occupying $\left[0.006,0.028\right]\times\left[0.003,0.037\right]$. We have $\Delta x,\Delta y=10^{-4}$. The surroundings are taken to be a vacuum. The aluminium follows a Godunov-Romenski EOS with parameters $\rho_{0}=2710$, $c_{v}=900$, $T_{0}=300$, $c_{0}=5037$, $\alpha=1$, $\beta=3.577$, $\gamma=2.088$, $b_{0}=3160$, $\sigma_{Y}=4\times 10^{8}$, $\tau_{0}=1$, $n=100$. Gauges are placed initially at $x=0.0078125$, $0.0114375$, $0.0150625$, $0.0186875$, $0.0223125$ to measure the state variables over time, and these gauges are permitted to move with the local velocity of the material. The test is run until time $t=5\times 10^{-6}$.

The pressure contours throughout the aluminium at various times are shown in Figure 14 on page 14. Despite relying on a slightly different plasticity model to that found in [42], it can be seen that these plots are in very good agreement with those found in the aforementioned publication. Note that release waves can be seen on the sides of the plate, in agreement with Michael & Nikiforakis.

Plots over time of the $x$-velocity, pressure, density, and total stress - as measured by the gauges - are given in Figure 15 on page 15. Note the good agreement between these plots and those found in [42, 31], both in terms of their qualitative shape, and the arrival times of the waves that they represent. One can clearly see the separation between the elastic precursor wave and the trailing plastic wave in the impacted plate, and the subsequent return waves that are generated once these waves reach the end of the plate. This implies that the GPR-RGFM has correctly captured the aluminium-vacuum interface.

3.8 Confined C4 Detonation without Backplate

This test is a variation of that found in [42]. A steel bar of length 0.03 and width 0.018 impacts upon a steel plate of depth 0.003, which is covering a region of depth 0.009 composed of C4 . The bar is initially traveling with speed 700. The system is surrounded by air.

The steel is modeled using a shock Mie-Gruneisen EOS, with parameters $\rho_{0}=7870$, $c_{v}=134$, $c_{0}=4569$, $\Gamma_{0}=2.17$, $s=1.49$, $c_{s}=3235$, $\sigma_{Y}=0.53\times 10^{9}$, $\tau_{0}=1$, $n=10$. The C4 is modeled using a JWL EOS, with parameters $\rho_{0}=1601$, $c_{v}=2.487\times 10^{6}/1601$, $\Gamma_{0}=0.8938$, $A=7.781\times 10^{13}$, $B=-5.031\times 10^{9}$, $R_{1}=11.3$, $R_{2}=1.13$, $c_{s}=1487$. The air is modeled using an ideal gas EOS, with parameters $\rho_{0}=1.18$, $c_{v}=718$, $\gamma=1.4$, $c_{s}=50$, $\mu=1.85\times 10^{-5}$.

The reaction of the C4 is captured using the ignition and growth model [34], where total energy $E$ is modified to include the term:

[TABLE]

with $\lambda$ being the volume fraction of unreacted material, governed by the dynamical equation:

[TABLE]

[TABLE]

In this test, the parameters are taken to be $Q=9\times 10^{9}/1601$, $I=4\times 10^{6}$, $G_{1}=1.4\times 10^{-20}$, $a=0.0367$, $b=2/3$, $c=2/3$, $d=1/3$, $x=7$, $y=2$, $\phi_{I}=0.022$, $\phi_{G1}=1$.

Figure 16 on page 16 displays the resulting pressure and C4 concentration at times $t=2.4\times 10^{-6}$ and $t=4.9\times 10^{-6}$. As can be seen, the kinetic energy of the bar is correctly transmitted to the steel plate, with the plate deforming in a manner qualitatively identical to that found in [42]. This energy is in turn transmitted to the C4, leading to an exothermic reaction and a symmetrical wavefront that travels through the material. The C4 concentration is depleted to 0.93 at time $t=2.4\times 10^{-6}$ and to 0.915 at time $t=4.9\times 10^{-6}$.

3.9 Confined C4 Detonation

This test is identical to the previous test, except a steel plate of depth 0.003 is now placed behind the C4, so that the explosive is entirely confined. As can be seen from Figure 17 on page 17, the kinetic energy of the bar is once again correctly transmitted to the steel plate and C4, with the same deformation occurring in the first steel plate. At time $t=2.4\times 10^{-6}$ we see the wave in the C4 both partially rebounding off the backplate back into the reactant, and partially traveling on through the backplate. At the earlier time, the reactant concentration has been depleted to 0.88, and at the later time to 0.865. This corroborates the results of [42], in that the presence of the backplate accelerates the reactive processes.

3.10 Confined C4 Detonation with Air Gap

This problem is designed to test the ability of the framework presented in this paper to capture the interaction of widely varying media. It is identical to the previous problem, except an air gap is now placed between the first steel plate and the C4. The air has the same EOS parameters as the surrounding air.

As can be seen in Figure 18 on page 18, the rod displaces the air (with the numerical method coping with contact of the region representing the plate with the region representing the C4). The displacement of the air enables the kinetic energy of the rod to be transmitted through the plate and into the C4, as before. The earlier time of $t=2.4\times 10^{-6}$ corresponds with the instant after the plate makes contact with the C4. The concentration is depleted to 0.9998 at $t=2.4\times 10^{-6}$ and 0.894 at $t=4.9\times 10^{-6}$. The latter value is lower than the corresponding value in the previous test, as the reaction has been delayed by the presence of the air gap.

3.11 Convergence Study

In order to determine the order of convergence of the Riemann Ghost Fluid Method presented in this study, the tests from 3.2, 3.3, and 3.4 were run with cell counts of $50,100,150,200,250,300$, and the error as compared with the exact solutions were calculated at the interfaces. The reason for choosing these two tests was to incorporate a range of different interface types (gas-liquid, fluid-solid, and solid-vacuum). The results are displayed in Table 5 on page 5. As can be seen, in all tests, the convergence rate is roughly first-order.

4 Discussion

The Riemann Ghost Fluid Method presented in this study has been demonstrated to be an effective way of accurately simulating the interfaces between several different materials (in all three main phases of matter, plus vacuum), described by the GPR model. Unlike in many existing implementations, heat conduction across the interface was simulated accurately, leading to the representation of physical phenomena that are often overlooked, and to the redundancy of numerical techniques that are sometimes used to enforce more empirically-accurate results (such as entropy and temperature fixes).

The framework presented here greatly simplifies the conceptual framework required for multimaterial interactions. Implementation should be easier and quicker, and future work can be more focused on a single model, rather than several fundamentally different frameworks.

4.1 Limitations

Throughout this study, the various fluids have been assumed to be immiscible. Whilst this is a common assumption in situations where mixing is low or practically non-existent, there are many problems which may require it. An area of further research would be the implementation of a mixture model such as that proposed by Romenski et al. [55, 54], which uses the same thermal conduction system as the GPR model.

The truncation of the Taylor series expansions (98a) and (98b) used to find the star states of the heat flux and the viscous stress tensor implicitly assumes that there are only small differences between the side states and the star states of the variables upon which $q_{1}^{*},\boldsymbol{\sigma_{1}^{*}}$ depend ($\rho$, $p$, $J_{1}$, and $\boldsymbol{A_{1}}$). If this is the case, higher order terms of the expansion can be neglected. If it is not, however, the method may fail. The linearised nature of the GPR-RGFM solver also implicitly assumes that all waves of interest present in the Riemann Problem are shocks. Thus, strong rarefactions may cause the method to fail.

4.2 Potential Improvements

As noted in Section 3.5, the GPR-RGFM method does not necessarily ensure the continuity of the normal component of heat flux across interfaces that feature discontinuous temperature jumps. This is despite the method accurately modeling both the heat conduction across the interface over time, and the corresponding evolution to thermal equilibrium between the two materials. The reason for this is that the star states produced by the linearised solver presented in Section 2.2 represent the state of the system at a time slightly beyond the point in time at which they are applied to neighboring ghost cells. With more simple systems of PDEs - such as the Euler equations - this often doesn’t matter, as the star states are constant in time, or their time evolution is easily calculated. Owing to the source terms in the GPR equations, however, the star states evolve in a manner for which an analytical solution is not available. Thus, when the star states are applied to the ghost cells, they contain higher heat fluxes and slightly different temperatures to the actual values at the interface at that moment in time, leading to the slight aberration apparent in the heat flux in Figure 12 on page 12.

A possible solution to this is as follows: Take materials $L$ and $R$ either side of an interface. Given the states straddling the interface, $Q_{L},Q_{R}$, derive $Q_{L}^{*}$ using the procedure outlined in Section 2.2. Then, by inverting this procedure, derive a state $\tilde{Q_{R}}$ such that if $Q_{L},\tilde{Q_{R}}$ are both states for material $L$ (rather than for materials $L$ and $R$, respectively), following the procedure in Section 2.2 obtains the same star state, $Q_{L}^{*}$. The derivation of such an inverse mapping should be feasible, possibly with the addition of some physical constraints. In this way, $\tilde{Q_{R}}$ represents a state for the righthand cell at the current time, which - if cell $R$ were occupied by material $L$ - would result in the same state on the lefthand side of the interface at the end of the current timestep as if the righthand cell were occupied by $Q_{R}$ and material $R$.

Another clear improvement to the GPR-RGFM method presented would be to use a better Riemann solver than the iterative, linearised solver devised in Section 2.2. Let $L$ be the matrix of left eigenvectors of the primitive system. As noted previously, the solver relies upon the fact that each of the following relations holds along the characteristic to which it corresponds:

[TABLE]

Methods to accurately integrate (130) from the left and right interface boundary states to their respective star states warrant further research.

Alternatively, a completely different approximate Riemann solver could be employed, such as the universal HLLEM solver of Dumbser et al. [15]. This path-conservative formulation of the HLLEM solver works for general non-conservative systems (such as the GPR model) and is simple to implement. It’s based upon a new path-conservative HLL method (building on the original method of Harten, Lax, and van Leer [27]) but is claimed to be able to represent linearly degenerate intermediate waves “with a minimum of smearing” by evaluating the eigenvalues and eigenvectors of the intermediate characteristic fields (given in Section 2.1).

There are iterative exact Riemann solvers for the equations of non-linear elasticity (to which the GPR model reduces as $\tau_{1}\rightarrow\infty$). Thus, they will work for applications of the GPR model to solids problems (and perhaps to very viscous fluids problems too). Although these solvers are computationally expensive, they are only used once at each material interface point at each time step, and thus the added accuracy that they provide may be desirable. There are two ways to formulate the equations of non-linear elasticity: one in which the deformation tensor (the analogue of the inverse of the GPR model’s distortion tensor) is evolved in time, and one in which its inverse (the analogue of $A$) is evolved instead. Miller’s exact solver [43] uses the first formulation and the solver of Barton et al. [5] uses the second. The former can be used to evolve $A^{-1}$, from which $A$ can be calculated. Unfortunately, both solvers critically assume that the source terms of the system vanish, and so are unlikely to produce the correct boundary conditions for the GPR-RGFM when modeling relatively inviscid fluids. It should also be noted that they cannot be used for problems involving heat conduction across material interfaces, and they do not take the thermal conduction subsystem of the GPR model into account.

5 References

6 Acknowledgments

The authors acknowledge financial support from the EPSRC Centre for Doctoral Training in Computational Methods for Materials Science under grant EP/L015552/1.

7 Appendix

Taking the ordering $\boldsymbol{P}$ of primitive variables in (144), note that (1e), (1b), (1c), (1d) can be stated as:

[TABLE]

[TABLE]

where the total stress tensor $\Sigma=pI+\rho A^{T}\psi$. Note that:

[TABLE]

Thus, the energy equation becomes:

[TABLE]

Simplifying:

[TABLE]

We have666

$\displaystyle\frac{p-\rho^{2}E_{\rho}}{\rho E_{p}}$ $\displaystyle=\frac{\rho^{2}\left.E_{\rho}\right|_{s}-\rho^{2}\left.E_{\rho}\right|_{p}}{\rho\left.E_{p}\right|_{\rho}}=\rho\frac{\left.E_{\rho}\right|_{s}-\left(\left.E_{\rho}\right|_{s}+\left.E_{s}\right|_{\rho}\left.s_{\rho}\right|_{p}\right)}{\left.E_{s}\right|_{\rho}\left.s_{p}\right|_{\rho}}$

(135)

$\displaystyle=\rho\frac{-\left.s_{\rho}\right|_{p}}{\left.s_{p}\right|_{\rho}}=\rho\left.\frac{\partial p}{\partial\rho}\right|_{s}$

777

$\frac{c_{t}^{2}T}{\rho E_{p}}=\frac{c_{t}^{2}T}{\rho c_{v}T_{p}}=\frac{\rho c_{h}^{2}}{T_{p}}$

(136)

888

$\left.\frac{\partial E}{\partial A}\right|_{\rho,p}=\left(c_{s}^{2}-\frac{\rho}{\Gamma}\frac{\partial c_{s}^{2}}{\partial\rho}\right)\frac{\psi}{c_{s}^{2}}=\left(1-2\frac{\rho^{2}}{\rho\Gamma}\frac{\partial\log\left(c_{s}\right)}{\partial\rho}\right)\psi$

(137)

$\displaystyle\frac{\partial\sigma}{\partial\rho}$ $\displaystyle=\frac{\partial}{\partial\rho}\left(-\rho c_{s}^{2}A^{T}\frac{\psi}{c_{s}^{2}}\right)=-c_{s}^{2}A^{T}\frac{\psi}{c_{s}^{2}}-\rho\frac{\partial c_{s}^{2}}{\partial\rho}A^{T}\frac{\psi}{c_{s}^{2}}$

(138)

$\displaystyle=\frac{\sigma}{\rho}+2\frac{\partial\log\left(c_{s}\right)}{\partial\rho}\sigma$

:

[TABLE]

[TABLE]

The full system then becomes:

[TABLE]

[TABLE]

Thus, the GPR system can be written in the following form:

[TABLE]

where the first component of $M$ is given on (142) for illustrative purposes.

[TABLE]

where we have:

[TABLE]

[TABLE]

[TABLE]

[TABLE]