Closure Models for Lagrangian Gas Dynamics and Elastoplasticity Equations in Multimaterial Cells

Abstract

Mixed cells (multicomponent cells) emerging in the development of Lagrangian‐Eulerian (ALE) or Eulerian numerical techniques for solving the gas dynamics and elastoplasticity equations in multicomponent media contain either interfaces between materials or a mixture of materials. There is a problem of correctly approximation of the equations in such cells and the ALE code accuracy and performance depend on how the problem is resolved. Many approximation methods use the equation splitting into two stages, one of which consists in solving a given equation in Lagrangian variables. If mixed cells are simulated, the system of equations describing the gas dynamics and elastoplasticity is unclosed and there is a need to introduce additional closure relations that will allow determining the thermodynamic parameters of components using the available data for the mixture of components, as a whole. The chapter presents a review of the equation closure methods and results of the methods verification using several test problems having exact solutions.

1. Introduction

Mixed cells in arbitrary Lagrangian‐Eulerian (ALE) or Eulerian methods contain interfaces between different materials or a mixture of materials. In the next section, we will not distinguish these two methods, considering that both methods solve the advection equation, including the vicinity of mixed cells. Most of these methods use a two‐stage approximation of equations. The first stage considers gas dynamics or elastoplasticity equations without convective terms. The convective transfer comes into play at the second stage. Among many similar methods, we consider only the ALE methods that contain Lagrangian gas dynamics and elastoplasticity in the pure form, and the problem of mixed cells at that particular stage is the subject of research reported here. Note that mixed cells can be present even in purely Lagrangian techniques, and the problems related to their presence should also be addressed in this case.

Here, we will generally use the term “Lagrangian gas dynamics” (or simply “gas dynamics”), bearing in mind that, in the case of elastoplasticity, this will also involve equation terms related to the stress tensor deviator. Historically, several approaches to the problem of mixed cells in gas dynamics associated with materials distinction in such cells have been considered. In this chapter, we consider only the single‐velocity model of matter. The major approach that has become predominant these days uses complete thermodynamic distinction of materials.¹

Early in the development of Eulerian methods, a smaller number of parameters have been used to identify the materials; for example, in [1, 2], mass fractions of the materials and average energy of matter were employed. Accordingly, other closing relations were used, the required number of which in this case is plus one compared with the complete materials distinction. The models thus considered include the “isobar-isothermal” and the “isobar-isodQ” models, which, although successful in some respect, in the general case failed to deliver acceptable accuracy of results.

Next, we will use the term “material” meaning that, mathematically, an interface can also divide identical materials; moreover, one of the materials can be vacuum and/or a perfectly rigid body.

Thermodynamic parameters in gas dynamics include density, internal energy, and pressure. If other processes are modeled, the number of parameters increases; for example, for elastoplasticity, additional parameters will include components of the stress tensor deviator. In addition to thermodynamic parameters, volume fractions of constituent materials are introduced in each mixed cell that can be used to determine the geometric location of the interface inside a mixed cell, which is used in some models.²

The problem of identifying the contact location based on the material volume fractions is beyond the scope of this study; it is a separate problem discussed in dedicated studies (see, e.g. [1, 3–11]).

This approach to materials identification allows one to model mixed cells containing not only contacting but also intermingled materials. When mixed cells are used for gas dynamics equations, additional closing relations are needed, which in fact define the interaction of materials inside the cell (the subcell interaction). Most of the known models manage with information about volume (or mass) fractions of the materials and their thermodynamic states [12–26]. Such models can be divided into two classes according to the number of computational stages involved.³

Our classification and description of models is limited to the case of two materials in a cell, although many formulas mentioned in this chapter are also suitable for their larger number. For this reason, some models developed specifically for the case of several materials in a cell are left beyond the scope of our review.

The first class of models is based on introducing closure models at a single stage, while the second class includes two‐stage models, in which the second stage is in fact complementary to the first one and involves additional interaction between materials inside a mixed cell (so‐called subcell interaction).

Next, we often use the terms “model” and “method” without distinction. One should note here that a method is understood to be an algorithm implemented in the form of a program and based on some physical model.

Basic single‐stage closure methods include the following:

1. Method based on the model of equal pressures of constituent materials (the P method) [12],

2. Tipton's method [13],

3. Delov's method based on the acoustic Riemann solver [14] (this approach is also used in the DSS [22] and KSR [23] methods developed later),

4. The K&S method based on a local Riemann problem [15]. This method is not described in this work because of its impracticability, as noted in [15], but its test simulation results are given for comparison.

Note that these four models and methods developed on their basis are relaxation with respect to pressure. In a number of studies, nonrelaxation methods have been proposed, which use the following assumptions (models):

1. Equal velocity divergences ($\nabla ·\mathbf{u}$) of the materials [17].

2. Equal pressure changes (Δp) of the materials [18],

3. Equal velocities (Δu) of the materials behind a weak shock wave [19].

Two‐stage models include the stage of subcell interactions between the materials in the nonequilibrium state; so the first stage here can only use models 5–7. This approach for closure models has been proposed independently in [20, 21]. The subcell pressure relaxation method in [20] is versatile and it is used in combination with models 5–7, denoted as the $\nabla ·\mathbf{u}$‐PR, Δp‐PR and Δu‐PR (pressure relaxation) methods.

All the above‐mentioned methods do not employ the contact location inside a mixed cell. However, there are methods that make essential use of the information about the contact location. A method of this kind was first proposed in [21] and then developed in the “interface‐aware subscale dynamics” IA‐SSD method [24, 25] for the multimaterial case. It offers a two‐stage model, the first stage of which employs the $\nabla ·\text{u}$ model. At the second stage, driven by the materials’ individual pressures, the interface between the materials moves normally to it. The interface is reconstructed based on the volume fractions, and its motion is accomplished based on the solution of an acoustic Riemann problem (model 3).

Let us point out one common feature (associated with the assumptions made in the models) of all the above closure methods. Velocity in the methods is normal to the interface (definite or imaginary) irrespective of interface location relative to the vector of velocity. In fact, they are isotropic in the sense that compression (expansion) ratios of materials are assumed to be equal in all directions. One can mention a number of other closure methods that employ algorithms similar to those used in the above‐mentioned models [27–33]. This property of the methods is quite acceptable for most applications, but there are problems (see below), as applied to which it results in significant errors in simulations.

In [34], anisotropic closure methods, ACM‐1 and ACM‐2, are proposed, which are an extension of models 5–7. They possess all the advantages of methods 5–7, which are central to the EGAK code [35] when modeling flows, for which one can assume that they are isotropic, but have an important advantage when modeling more complex flows.

Apart from the basic closure method, mixed cells require additional relations to address the ways of pressure and artificial viscosity calculations for the whole cell and artificial viscosity calculations for the materials. Six approaches to calculate the artificial viscosity of materials are discussed in [36].

2. Finite difference approximation of elastoplasticity equations

2.1. Initial equations for multimaterial elastoplasticity

The initial set of equations solved at the Lagrangian stage for 2D elastoplastic flows is the following:

$\frac{d\mathbf{u}}{dt}=-\frac{1}{\rho }\nabla ·\mathbf{T},$E1

$\frac{d{\rho }_{\xi }}{dt}=-{\rho }_{\xi }\nabla ·{\mathbf{u}}_{\xi },$E2

$\frac{d{\beta }_{\xi }}{dt}={\beta }_{\xi }\left(\nabla ·{\mathbf{u}}_{\xi }-\nabla ·\mathbf{u}\right),$E3

$\frac{d{e}_{\xi }}{dt}=\frac{1}{{\rho }_{\xi }}Sp({\mathbf{T}}_{\xi }{\mathbf{D}}_{\xi }),$E4

$\frac{d\mathbf{r}}{dt}=\mathbf{u}.$E5

In this set of equations: u(u_x, u_y) is the velocity, ρ is the density, T is the stress tensor, D is the strain rate tensor, and е is the specific internal energy, βis the volume fraction of the material (${\beta }_{\xi }=\frac{M_{\xi }}{V}$), r(x,y) is the radius vector. The subscript ξ is the material index; also note that in the expression for the velocity divergence (or simply divergence for short) it relates to the divergence as a whole, rather than to the velocity. Bold type here and below is used to indicate the vector, tensor, and deviator.

Eq. (3) can be derived from the equation of continuity (2), in which the materials’ density is substituted with its expression in the form of ${\rho }_{\xi }=\frac{M_{\xi }}{V_{\xi }}$, where M_ξ and V_ξ are the mass and volume of the materials in the Lagrangian cell. We obtain $\frac{d{V}_{\xi }}{dt}=V_{\xi }\nabla ·{\mathbf{u}}_{\xi }$ and then introduce an expression for V_ξ in terms of volume fractions V_ξ = β_ξV into this equation. Thus, Eq. (3) is a consequence of Eq. (2), and we give it here solely for the purpose of empathizing that the volume fractions for multimaterial matter should also be updated to tⁿ⁺¹. In the single‐material case, Eqs. (1)–(5) come down to usual Lagrangian gas dynamics equations, because Eq. (3) is not present in this case, and the quantities in other equations are written without material indices.

Stress and strain rate tensors are expressed as follows:

$\mathbf{T}=\left[\begin{array}{ccc}{T}_{xx}& T_{xy}& 0\\ T_{yx}& T_{yy}& 0\\ 0& 0& {\mathbf{T}}_{\phi }\end{array}\right],\mathbf{D}=\left[\begin{array}{ccc}{d}_{xx}& d_{xy}& 0\\ d_{yx}& d_{yy}& 0\\ 0& 0& d_{\phi }\end{array}\right].$E6

The stress tensor is represented as a sum of the spherical part (pressure p) and the deviator S(S_xx, S_yy , S_xy , S_φ). Deviator components are defined by the relation S_ij = T_ij - δ_ijp.

For the materials, we define equations of state

$p_{\xi }=P_{\xi }\left({\rho }_{\xi },e_{\xi }\right),$E7

and equations to express the deviator S_ξ as a function of the strain rate tensor D_ξ

$f_{\xi }({\mathbf{S}}_{\xi },{\mathbf{D}}_{\xi })=0.$E8

The specific form of Eq. (8) is determined by the model of matter adopted.

EGAK uses decomposition in physical processes, in which the pressure‐related terms are approximated at the Lagrangian gas dynamics stage, and the terms related to the stress tensor deviator, at the other stage of the computation. In the present technique, materials can be both different substances with their equations of state and vacuum.

2.2. Finite difference approximation of elastoplasticity equations

EGAK uses a quadrangular mesh with node‐centered velocities and all the other quantities (${\rho }_{\xi }^{}{,}_{}{\beta }_{\xi }^{}{,}_{}{e}_{\xi }^{},p_{\xi }^{}$, S_ξ) defined at cell centers and for each material individually. Also note that for the purpose of program implementation, pressure in Eqs. (1)–(5) is replaced with a sum of pressure and artificial (computational) viscosity for matter as a whole and for the materials, p → p + q and p_ξ → p_ξ + q_ξ, respectively. Known quantities (basic variables) in Eqs. (1)–(5) in the 2D case include ${\mathrm{u}}_x{,}_{}{\mathrm{u}}_y{,}_{}{\rho }_{\xi }^{}{,}_{}{\beta }_{\xi }^{}{,}_{}{e}_{\xi }^{},p_{\xi }^{}$, S_ξ, and the quantities p, S, q, q_ξ, $\nabla ·{\mathbf{u}}_{\xi }$, D_ξ need to be determined. In the following formulas, the subscript means discretization in space and the material number ξ in the multimaterial case, and the superscript denotes discretization in time. Cell‐centered quantities are marked with a semi‐integer superscript (for example, i +1/2), and the node‐centered ones, with an integer superscript (i). If in a specific formula it is clear from the context that the superscripts are the same for all the quantities, they are omitted. Cell masses in Lagrangian gas dynamics remain constant in the course of calculations, so they have no temporal indexing.

Suppose that we know all the basic quantities at time tⁿ and that we seek to update their values at time tⁿ⁺¹ = tⁿ + τ, where τ is the timestep chosen based on the requirement that the difference scheme should be stable (these issues are beyond the scope of this work). Let us write the difference scheme of EGAK for the multimaterial case.

First half‐timestep (determination of predicted pressure)

${\mathbf{r}}^{n+1/2}={\mathbf{r}}^n+\tau ·{\mathbf{u}}^n,$E9

$V_{i+1/2,j+1/2}^{n+1/2}=V\left({\mathbf{r}}_{i+1,j+1}^{n+1/2},{\mathbf{r}}_{i,j+1}^{n+1/2},{\mathbf{r}}_{i+1,j}^{n+1/2},{\mathbf{r}}_{i,j}^{n+1/2}\right),$E10

$p_{\xi }^n=P\left({\rho }_{\xi }^n,e_{\xi }^n\right),$E11

$\nabla ·{\mathbf{u}}_{}^n=\frac{\left(V_{}^{n+1/2}-V_{}^n\right)}{\left(\tau ·V_{}^n\right)},$E12

$p_{\xi }^{n+1/2}=p_{\xi }^n-\chi ·{\rho }_{\xi }^n·{\left(c_{\xi }^n\right)}^2·\tau ·{\overline{\nabla ·\mathbf{u}}}_{\xi }^n.$E13

In Eq. (13), $\chi =0.6$ (this value was chosen in [17]), $c_{\xi }^n$ is the speed of sound.

Full timestep

$M_{i,j}=0.25·{\displaystyle \sum _{\xi }\left(M_{\xi ,i-1/2,j-1/2}+M_{\xi ,i-1/2,j-1/2}+M_{\xi ,i+1/2,j-1/2}+M_{\xi ,i+1/2,j+1/2}\right)},M_{\xi }={\rho }_{\xi }·{\beta }_{\xi }·V,$E14

$M_{i,j}\frac{\left({\mathbf{u}}_{i,j}^{n+1}-{\mathbf{u}}_{i,j}^n\right)}{\tau }={\left({\overline{\nabla g}}^{n+1/2}+{\overline{\nabla ·S}}^n\right)}_{i,j},$E15

$\Rightarrow {\mathbf{u}}_{i,j}^{n+1}={\mathbf{u}}_{i,j}^n-\left(\frac{\tau }{M_{i,j}}\right)·{\left({\overline{\nabla g}}^{n+1/2}+{\overline{\nabla ·S}}^n\right)}_{i,j},$E16

${\mathbf{u}}^{n+1/2}=\left({\mathbf{u}}^n+{\mathbf{u}}^{n+1}\right)/2,$E17

${\mathbf{r}}^{n+1}={\mathbf{r}}^n+\tau ·{\mathbf{u}}^{n+1},$E18

$V_{\xi }^{n+1}=V_{\xi }^n+\tau ·V_{\xi }^{n+1}{\overline{\nabla ·\mathbf{u}}}_{\xi }^{n+1},$E19

${\rho }_{\xi }^{n+1}=\frac{M_{\xi }^{}}{V_{\xi }^{n+1}},$E20

${\beta }_{\xi }^{n+1}=\frac{V_{\xi }^{}}{V_{}^{n+1}},$E21

$\begin{array}{ll}{e}_{\xi }^{n+1}\hfill & =e_{\xi }^n-\frac{\tau }{{\rho }_{\xi }^n}·[g_{\xi }^{n+1/2}·{\overline{\nabla ·u}}_{\xi }^{n+1/2}-S_{\xi ,xx}^n{\overline{d}}_{\xi ,xx}^{n+1/2}-S_{\xi ,yy}^n{\overline{d}}_{\xi ,yy}^{n+1/2}-2S_{\xi ,xy}^n{\overline{d}}_{\xi ,xy}^{n+1/2}+\hfill \\ \hfill & \left(S_{\xi ,xx}^n+S_{\xi ,yy}^n\right){\overline{d}}_{\xi ,\phi }^{n+1/2}].\hfill \end{array}$E22

In Eq. (22), $g_{\xi }^{n+1/2}=p_{\xi }^{n+1/2}+q_{\xi }^n$. The methods to calculate the materials’ artificial viscosity $q_{\xi }^n$ are discussed in Section 4. The bar denotes the difference counterpart of a corresponding operator (the formulas are generally known, so we skip them). In the following, we will not use the bars assuming that all the operators are difference operators. We have not described the way of updating the materials’ stresses, i.e., the approximation of Eq. (16), as it is beyond the scope of this work; we only note here that these equations include components of the tensor D_ξ.

The formulas for total pressure, viscosity, and stress deviator are the following:

$\begin{array}{l}{p}_{}^{n+1/2}={\displaystyle \sum _{\xi }{\psi }_{\xi }^{}{p}_{\xi }^{n+1/2}},\\ q_{}^n={\displaystyle \sum _{\xi }{\psi }_{\xi }^{}{q}_{\xi }^n},\\ {\text{S}}_{}^n={\displaystyle \sum _{\xi }{\psi }_{\xi }^{}{\text{S}}_{\xi }^n.}\end{array}$E23

where the factor ψ_ξ is determined by the chosen closure model (see below).

Thus, the quantities that are not yet determined in Eqs. (14)–(22) include $\nabla ·{\mathbf{u}}_{\xi }$, ψ_ξ, q_ξ, and D_ξ. To calculate these, one needs to use some closure relations, being the consequences of different assumptions (models) about thermodynamic states of the materials in mixed cells.

When introducing the closure relations, one should fulfill some requirements resulting from the laws of conservation.

Requirement 1 is additivity of volumes (the law of conservation of “volume”)

$\begin{array}{l}V={\displaystyle \sum _{\xi }{V}_{\xi }}\\ \text{or}\\ \\ {\displaystyle \sum {\beta }_{\xi }}=1,\end{array}$E24

the consequence of which is the relation

$\begin{array}{l}\Delta V={\displaystyle \sum _{\xi }\Delta V_{\xi }},\\ \text{or}\\ \\ {\displaystyle \sum {\beta }_{\xi }\nabla ·{\text{u}}_{\xi }^{}}=\nabla ·\text{u}.\end{array}$E25

The natural extension of relation (25) is $\mathbf{D}={\displaystyle \sum {\beta }_{\xi }{\mathbf{D}}_{\xi }}$, which is fulfilled at

${\mathbf{D}}_{\xi }=\mathbf{D}\left(\frac{\nabla ·{\mathbf{u}}_{\xi }}{\nabla ·\mathbf{u}}\right).$E26

Formula (26) is used to determine D_ξ in Eq. (22) and when approximating Eq. (16).

Requiment 2 is additivity of energies (the law of conservation of energy)

$e={\displaystyle \sum {\alpha }_{\xi }{e}_{\xi }},$E27

where the mass fraction ${\alpha }_{\xi }=\frac{M_{\xi }}{M}$ is given by the following expression:

$\begin{array}{l}{\alpha }_{\xi }={\beta }_{\xi }·\frac{{\rho }_{\xi }}{\rho },\\ \text{where}\\ \rho ={\displaystyle \sum {\beta }_{\xi }{\rho }_{\xi }}.\end{array}$E28

The requirement (27) can also be written for increments of specific energies

$\Delta e={\displaystyle \sum {\alpha }_{\xi }\Delta e_{\xi }},$E29

where $\text{Δ}{e}_{\xi }$ is the energy increment for the constituent material, and Δe, for the whole cell.

Let us consider closure methods for the case of approximation of gas dynamics equations. In EGAK, the difference approximation of the energy equation (22) for the constituent materials has the following form:

$\Delta e_{\xi }=-\frac{\tau g_{\xi }^{n+1/2}}{{\rho }_{\xi }^n}\nabla ·{\mathbf{u}}_{\xi }^{n+1/2}.$E30

We insert their expression (30) into Eq. (29) for $\text{Δ}{\mathbf{e}}_{\xi }$ and, using Eq. (28), obtain

$-\frac{\tau g_{}^{n+1/2}\nabla ·{\mathbf{u}}_{}^{n+1/2}}{\rho }=-\frac{\tau {\displaystyle \sum {\beta }_{{}_{\xi }}^{}{g}_{\xi }^{n+1/2}\nabla ·{\mathbf{u}}_{\xi }^{n+1/2}}}{\rho }+{\displaystyle \sum {\alpha }_{{}_{\xi }}^{}\Delta e_{\xi }^{/}}.$E31

Using Eq. (23) and the given ways of finding $\nabla ·{\mathbf{u}}_{\xi }^{}$ from Eq. (31) one can obtain the values of $\Delta {\mathbf{e}}_{\xi }^{/}$ that represent additional changes in the materials’ internal energy to meet the energy balance requirement.

Given that $g_{\xi }^{n+1/2}=g_{}^{n+1/2}$, it follows from Eq. (31) that ${\displaystyle \sum {\alpha }_{\xi }\Delta e_{\xi }^{/}}=0$. Thus, this term represents the change in the materials’ internal energy as a result of their pressure relaxation. If no pressure relaxation is used, i.e., $\Delta e_{\xi }^{/}=0$, then the definitions of ψ_ξ follow directly from the closing conditions.

In the general case of nonequal pressures, the requirement Eq. (31) can be fulfilled, when some conditions imposed on the function ψ_ξ are satisfied. These conditions are discussed in Section 3.

Let us write the expression for ψ_ξ as

${\psi }_{\xi }={\beta }_{\xi }{\lambda }_{\xi },$E32

where the quantity ${\lambda }_{\xi }$ is determined by the chosen model of distributing the total divergence of the mixed cell to the constituent materials from the relation

$\nabla ·{\mathbf{u}}_{\xi }={\lambda }_{\xi }\nabla ·\mathbf{u}.$E33

3. Closure methods for gas dynamics equations

There is no single closure method for gas dynamics equations in mixed cells that might be suitable for all types of flow. Closing relations are quite numerous, but many of them are not used nowadays and are of historical interest only. Next, we consider the most frequently used closure methods, most of which have been implemented in EGAK.

3.1. Isotropic single‐stage closure methods

It will be convenient to introduce some closure methods if we consider the 1D problem, as shown in Figure 1. The i – 1/2th cell is mixed; it contains two materials; the interface is denoted by point A. Depending on the way of velocity definition at the point A, one can derive one or another method for calculating divergences (and densities) of the materials in the mixed cell.

Figure 1.

Computational mesh. Point A is the interface.

3.1.1. Method based on the equal pressure model

Method 1 uses the assumption that the materials have equal pressures (as proposed in [1]); in addition, artificial viscosities are assumed to be equal, too:

$\begin{array}{l}{p}_{\xi }=p,\\ q_{\xi }=q.\end{array}$E34

For EGAK, method 1 based on the model (Eq. (34)) has been developed in [37].

This method first calculates the energy increment for the cell as a whole:

$\begin{array}{c}\Delta E=\left(e_{}^{n+1/2}-e^n\right)M=-\frac{\tau \left(p^{n+1/2}+q^n\right)}{{\rho }^n}\nabla ·{\mathbf{u}}_{}^{n+1/2}M\\ =-\tau \left(p^{n+1/2}+q^n\right)\nabla ·{\mathbf{u}}_{}^{n+1/2}{V}^n=\left(p^{n+1/2}+q^n\right)\mu \Delta V,\end{array}$E35

where $\mu =-\nabla ·{\mathbf{u}}_{}^{n+1/2}{V}^n\tau /\Delta V$.

By analogy with Eq. (35), the energy increment equation for the materials can be rewritten by adding respective material indices in the quantities ΔV and ΔE. Then, using expressions (25) and (34), one can write the following closed system of equations:

$\begin{array}{c}\Delta E_{\xi }=\left(p^{n+1/2}+q^n\right)\mu \Delta V_{\xi },\\ \Delta V={\displaystyle \sum _{\xi }\Delta V_{\xi }},\\ P\left(\frac{M_{\xi }}{V_{\xi }^n+\Delta V_{\xi }},\frac{E_{\xi }+\Delta E_{\xi }}{M_{\xi }}\right)=p^{n+1/2}.\end{array}$E36

The system (Eq. (36)) contains 2N + 1 equation in 2N + 1 unknown ΔV_ξ, ΔE_ξ, and p^n+1/2 and can be solved iteratively. Solving the system of equations gives updated values of specific energies and volume fractions of the materials:

$e_{\xi }^{n+1/2}=e_{\xi }^{}+\frac{\Delta E_{\xi }^{}}{M_{\xi }},{\beta }_{\xi }^{n+1/2}=\frac{V_{\xi }^{}+\Delta V_{\xi }^{}}{V+\Delta V}.$E37

Among the drawbacks of this method one should note its expensiveness because of the iterative methods that are needed for solving the system (Eq. (36)) with complex equations of state. Also note that the assumption about equal pressures can turn out to be inconsistent, for example, in problems associated with energy release at every timestep.

3.1.2. Tipton's method

The model underlying Tipton's method is close to the model (34). Let us consider it as applied to the difference scheme implemented in [22] and being slightly different from Eqs. (9)–(22). At the first half‐timestep, instead of Eq. (9) the method uses the equation ${\mathbf{r}}_{}^{n+1/2}={\mathbf{r}}_{}^n+\tau /2·{\mathbf{u}}_{}^n,$ and in Eq. (13), $\chi =1.0$. Then,

$p_{\xi }^{n+1/2}=p_{\xi }^n-{\rho }_{\xi }^n·{\left(c_{\xi }^n\right)}^2·\tau ·\nabla ·{\mathbf{u}}_{\xi }^n.$E38

This method assumes that the pressures of the constituent materials at the half‐timestep should equal the average pressure $p_{}^{n+1/2}$, which requires

$p_{}^{n+1/2}=p_{\xi }^{n+1/2}+R_{\xi }^{n+1/2},$E39

where $R_{\xi }^{n+1/2}$ is different for different materials. The following equation is used to determine this quantity:

$R_{\xi }^{n+1/2}=-{\rho }_{\xi }^n·c_{\xi }^n·h^n·{\overline{\nabla ·\mathbf{u}}}_{\xi }^n,$E40

where hⁿ is the characteristic mesh spacing.

Eqs. (38) and (39) can be combined into one equation.

$p_{}^{n+1/2}=p_{\xi }^n-{\tilde{B}}_{\xi }{}^n\nabla ·{\mathbf{u}}_{\xi }^n,$E41

where ${\tilde{B}}_{\xi }{}^n\equiv {\left(c_{\xi }^n\right)}^2{\rho }_{\xi }^n·\tau \left[1+\frac{h^n}{\left(\tau c_{\xi }^n\right)}\right].$

Eq. (41) together with the requirement (Eq. (25)) constitutes a closed set of algebraic equations with the unknowns $\nabla ·{\mathbf{u}}_{\xi }^n$ and $p_{}^{n+1/2}$. As mentioned in Section 3.1.1, the system can be solved; for ξ = 1, 2, the solution is the following:

$\begin{array}{l}\nabla ·{\mathbf{u}}_{\xi }^n=\frac{\left(p_{\xi }^n-{\overline{p}}^n\right)+{\overline{B}}^n·\nabla ·{\mathbf{u}}^n}{{\tilde{B}}_{\xi }^n},\\ p^{n+1/2}={\overline{p}}^n-{\overline{B}}_n\nabla ·{\mathbf{u}}^n,\end{array}$E42

where the barred terms denote average values of the quantities

${\overline{p}}_{}^n\equiv \frac{{\displaystyle \sum \left({\beta }_{\xi }^n{p}_{\xi }^n/{\tilde{B}}_{\xi }^n\right)}}{{\displaystyle \sum \left({\beta }_{\xi }^n/{\tilde{B}}_{\xi }^n\right)}}\text{and}{\overline{B}}_{}^n\equiv {\left[{\displaystyle \sum \left(\frac{{\beta }_{\xi }^n}{{\tilde{B}}_{\xi }^{}}\right)}\right]}^{-1}.$E43

Eq. (42) leads to an expression for the change in the volume fractions

$$\Delta {\beta }_{\xi }^{n+1/2}={\beta }_{\xi }^n\left(\frac{p_{\xi }^n-{\overline{p}}_{}^n}{{\tilde{B}}_{\xi }{}^n}\right)+{\beta }_{\xi }^n\left[\tau (\frac{{\overline{B}}^n}{{\tilde{B}}_{\xi }^n}-1)\right]\nabla ·{\mathbf{u}}_{}^n.$$E44

To update the quantities at tⁿ⁺¹, another assumption is made that increments of the quantities after a full timestep are twice the half timestep increments:

$\begin{array}{l}\Delta {\beta }_{\xi }^{n+1}=2\Delta {\beta }_{\xi }^{n+1/2}\\ \text{and}\\ {\beta }_{\xi }^{n+1}={\beta }_{\xi }^n+2\Delta {\beta }_{\xi }^{n+1}.\end{array}$E45

The difference energy equation for the materials in the cell is the following:

$e_{\xi }^{n+1}=e_{\xi }^n-p^{n+1/2}\frac{\Delta V_{\xi }^{n+1}}{M_{\xi }},$E46

where $\Delta V_{\xi }^{n+1}=\Delta {\beta }_{\xi }^{n+1/2}·V_{}^{n+1}$.

3.1.3. Delov's method

Method 3 based on the acoustic Riemann solver is proposed in [14]. Let us consider this method as applied to the one‐dimensional problem (see Figure 1) with a single velocity component u. We consider the following acoustic Riemann problem in a mixed cell:

$\begin{array}{c}p=p_1^n,u=u_{i-1}^{n+1/2}\text{for}m<m_{А},\\ p=p_2^n,u=u_i^{n+1/2}\text{for}m>m_{А},\end{array}$E47

where m is the mass variable (m_A is the cell mass).

The set of equations for the quantities similar to the Riemann invariants along the advanced and retarded characteristics has the following form:

$\begin{array}{c}{u}_A^{n+1/2}+\frac{p_A^{n+1/2}}{{\left(\rho c\right)}_1^n}=u_{i-1}^{n+1/2}+\frac{p_1^n}{{\left(\rho c\right)}_1^n},\\ u_A^{n+1/2}-\frac{p_A^{n+1/2}}{{\left(\rho c\right)}_2^n}=u_i^{n+1/2}-\frac{p_2^n}{{\left(\rho c\right)}_2^n}.\end{array}$E48

The solution to this set will be the following expression for the velocity $u_A^{n+1/2}$:

$u_A^{n+1/2}=\frac{u_{i-1}^{n+1/2}·{\left(\rho c\right)}_1^n+u_i^{n+1/2}·{\left(\rho c\right)}_2^n+p_1^n-p_2^n}{{\left(\rho c\right)}_1^n+{\left(\rho c\right)}_2^n}.$E49

Now, let us write the equation of continuity for the problem at issue. By replacing the divergence with a corresponding expression, we obtain the following equation:

$\frac{1}{{\rho }_{}^{n+1/2}}-\frac{1}{{\rho }_{}^n}=\tau ·\frac{u_i^{n+1/2}-u_{i-1}^{n+1/2}}{M},$E50

where $M=\rho h$ is the linear cell mass.

A similar equation for the materials is obtained if one of the velocities is replaced with a velocity at the point A and a respective index is attached to density and mass. After substituting the expressions for velocity (49) into the equation of continuity (50), we obtain

$\frac{1}{{\rho }_{\left\{\begin{array}{l}1\\ 2\end{array}\right\}}^{n+1/2}}-\frac{1}{{\rho }_{\left\{\begin{array}{l}1\\ 2\end{array}\right\}}^n}=\frac{{\tau }^{}}{M_{\left\{\begin{array}{l}2\\ 1\end{array}\right\}}}·\left[\frac{{\left(\rho c\right)}_{\left\{\begin{array}{l}2\\ 1\end{array}\right\}}^n}{{\left(\rho c\right)}_1^n+{\left(\rho c\right)}_2^n}·\left(u_i^{n+1/2}-u_{i-1}^{n+1/2}\right)+\frac{p_1^n-p_2^n}{{\left(\rho c\right)}_1^n+{\left(\rho c\right)}_2^n}\right].$E51

Thus, as distinct from models 1–3, the present model does not employ any equilibration algorithm for pressure relaxation in this method, because the volume changes of the materials are also controlled by their pressures.

The extension of Eq. (51) in the multimaterial case without constraint of the number of materials can be written as follows:

$\frac{1}{{\rho }_{\xi }^{n+1/2}}=\frac{1}{{\rho }_{\xi }^n}+\frac{{\lambda }_{\xi }^n}{{\alpha }_{\xi }^n}·\left(\frac{1}{{\rho }_{}^{n+1/2}}-\frac{1}{{\rho }_{}^n}\right)+\frac{\omega ·\tau }{{\beta }_{\xi }^n·{\rho }_{\xi }^n·h_{}^n}\left(p_{\xi }^n-p_{\sum }^n\right)·\frac{1}{{\left(\rho с\right)}_{\sum }^n},$E52

where

$p_{\sum }^n=\frac{1}{N}{\displaystyle \sum _{\xi =1}^N{p}_{\xi }^n},$E53

${\left(\rho c\right)}_{\sum }^n=\frac{1}{N}{\displaystyle \sum _{\xi =1}^N{\left(\rho c\right)}_{\xi }^n},$E54

${\lambda }_{\xi }^n=\frac{1}{\left(N-1\right)}\left(1-\frac{{\left(\rho c\right)}_{\xi }^n}{{\displaystyle \sum _{\zeta =1}^N{\left(\rho c\right)}_{\zeta }^n}}\right),$E55

where hⁿ is the characteristic mesh spacing, and ω∼1 is the factor introduced to improve stability conditions of the difference scheme.

From Eq. (55), using Eq. (12), one can obtain an expression for divergences:

$\nabla ·{\mathbf{u}}_{\xi }^n=\frac{{\lambda }_{\xi }^n}{{\beta }_{\xi }^n}·\nabla ·{\mathbf{u}}^n+\omega ·\frac{\left(p_{\xi }^n-p_{\sum }^n\right)}{{\beta }_{\xi }^n·h^n·{\left(\rho с\right)}_{\sum }^n}.$E56

Now, let us consider the quantity ${\lambda }_{\xi }^n$. The change in the materials’ volume, as we see from Eq. (55), can be written as

$V_{\xi }^{n+1}-V_{\xi }^n=\frac{{\lambda }_{\xi }^n}{{\beta }_{\xi }^n}·\left(V_{}^{n+1}-V_{}^n\right)+\omega ·\tau ·\frac{\left(p_{\xi }^n-p_{\sum }^n\right)}{{\left(\rho с\right)}_{\sum }^n}.$E57

From this, using the requirement (Eq. (25)), we obtain that the equality ${\displaystyle \sum _{\xi =1}^N{\lambda }_{\xi }^n{\left({\beta }_{\xi }^n\right)}^{-1}}=1$ must be fulfilled. One can easily see that in the 1D case, for ω = 1 and N = 2, Eq. (52) includes Eq. (51).

Now, let us show that the quantity ${\lambda }_{\xi }^n{\left({\beta }_{\xi }^n\right)}^{-1}$ can be taken as a function ${\psi }_{\xi }$ (${\psi }_{\xi }={\lambda }_{\xi }^n{\left({\beta }_{\xi }^n\right)}^{-1}$) for determining the average pressure using Eq. (23). For this purpose, it will suffice to consider the case when the materials have equal pressures. Indeed, in this case, the second term in (56) is equal to zero, and to meet the energy additivity requirement (31), it will be sufficient to assume that in Eq. (23)

${\psi }_{\xi }={\lambda }_{\xi }^n/{\beta }_{\xi }^n.$E58

If the materials have different pressures, when we use Eq. (58), the right‐hand member of the energy equation should be corrected by $\Delta {e^{\prime }}_{\xi }$, for example, in the form of

$\Delta {e^{\prime }}_{\xi }=\frac{\omega ·\tau }{{\alpha }_{\xi }{\rho }_{}^n·h_{}^n}·\frac{{\psi }_{\xi }^{}}{{\left(\rho c\right)}_{\Sigma }^n}·{\displaystyle \sum _{\xi }{p}_{\xi }^n·\left(p_{\xi }^n-p_{\Sigma }^n\right)}.$E59

This method provides good results in Lagrangian calculations, when the materials’ volume fractions in the cell are invariable and close to each other, i.e., at $1\gg {\beta }_{\xi }^{}\gg 0$. However, the values of β in ALE calculations can range freely within 0 < β < 1 and, at β close to zero, the method gives a noticeable error associated with the presence of division by β in Eq. (58). This situation is physically attributed to the fact that, in Riemann problem calculations, waves travel some distance that must be equal to or less than the size of the region occupied by each of the materials. As the choice of the timestep is based on the Courant constraint, the necessary requirement can be violated, which leads to unphysical results (e.g., a negative updated value of material volume). To fix this, additional constraints on volume increments are required.

3.1.4. Method 5 based on the equal compressibility model

Method 5 rests upon the most frequently used model of equal compressibility of materials, or to put it differently, of their equal divergences. The model has been proposed in [17]; it has also been considered in [38, 39] and is formulated as follows:

$\nabla ·{\mathbf{u}}_{\xi }^n=\nabla ·{\mathbf{u}}_{}^n.$E60

Naturally, it is also assumed that $\nabla ·{\mathbf{u}}_{\xi }^{n+1}=\nabla ·{\mathbf{u}}_{}^{n+1}$.

This method, which appears at first glance to be strange, results from the assumption that the velocity at the point A (Figure 1) is determined by distance‐linear interpolation between nodal velocities. In the 2D case, this method is generalized on the assumption of volume‐linear interpolation.

In this method, λ_ξ = 1 and formula (25) gives

${\psi }_{\xi }={\beta }_{\xi },$E61

which ensures the fulfillment of the requirement (31) at $\Delta {e^{\prime }}_{\xi }=0$. As a result, formula (23) transforms to

$p={{\displaystyle \sum {\beta }_{\xi }p}}_{\xi }^{}.$E62

Thus, all the quantities we need to solve Eqs. (9)–(22) are defined.

Formula (62), which is natural for a homogeneous mixture of ideal gases, has a simple interpretation for a heterogeneous mixture. Let us consider a mixed cell containing two materials of volume V₀ (see Figure 2).

Figure 2.

Graphic illustration of Formula (62).

Let us assign p_ξ to the centers of their volumes. If the volume is used as a variable, then the linear interpolation of pressure over the cell volume can be written as

$p(V)=2·\frac{p_2-p_1}{V_o}\left(V-\frac{V_1}{2}\right)+p_1.$E63

If we insert the value of $V=V_0/2$ into Eq. (63), we will obtain the very same Formula (62). It is easy to demonstrate that for the case of an arbitrary number of materials, the value of pressure at the point $V_0/2$ will also be defined by Formula (62) (first, two materials are considered, and then other materials are added one by one). Thus, in the heterogeneous case, for linear interpolation of pressure between the materials’ pressures, Formula (62) defines the pressure at the volume center of the mixed cell. As the approximation of the momentum equation uses cell‐related pressures, Formula (62) is consistent enough for determining the average pressure in the mixed cell.

Strengths and weaknesses of method 5 are evident. It is easy to implement and inexpensive in operation, but it can lead to nonphysical states of the materials. The point is that the materials in the mixed cell in calculations by this method are compressed uniformly, which leads to different pressures of the materials, which do not relax with time (see the test calculations in Section 5). Nevertheless, the method delivers quite acceptable results when used for flows with distinct interfaces.

3.1.5. Method 6 based on the equal pressure increments model

Bondarenko and Yanilkin [18] proposed a closure method based on the equality of pressure increments of the materials in the mixed cell. This model is mathematically expressed as

${\rho }_{\xi }{c}_{\xi }^2\nabla ·{\mathbf{u}}_{\xi }={\rho }_{\zeta }{c}_{\zeta }^2\nabla ·{\mathbf{u}}_{\zeta }.$E64

The condition (64) is derived as follows. For adiabatic flows,

$\frac{\partial p}{\partial t}=\frac{\partial p}{\partial \rho }·\frac{\partial \rho }{\partial t}=c_{}^2\rho \nabla ·\mathbf{u}.$E65

Going over in Eq. (65) to the materials and claiming that $\frac{\partial p_{\xi }}{\partial t}=\frac{\partial p_{\zeta }}{\partial t}$, we obtain the condition (64).

The set of algebraic equations (64), Eq. (25) is closed, and its solution is given by

$\nabla ·{\mathbf{u}}_{\xi }^n={\lambda }_{\xi }\nabla ·{\mathbf{u}}^n,$E66

where

${\lambda }_{\xi }^{}={\left(\frac{{\beta }_k^{}}{{\rho }_k^{}{c}_k^2}\right)}^{-1}·{\left({\rho }_{\xi }^{}{c}_{\xi }^2\right)}^{-1}.$E67

When this model is used in calculations for adiabatic flows, pressures will stay approximately equal, if the materials’ initial pressures in the cell are equal. However, in some cases, pressures may turn out to be different: if energy release is specified for one of the materials; behind a strong shock in the mixed cell, because the condition of equal pressure increments is incorporated in the adiabatic approximation; after calculating convective fluxes at the Eulerian stage, etc. In this case, the use of this model in its original form is associated with some problems.

One of them is the following. Let $p_{\xi }^n\ne p_{\varsigma }^n$. For an ideal gas, the following estimate is true:

$\nabla ·{\text{u}}_{\xi }\sim {\left({\rho }_{\xi }{c}_{\xi }^2\right)}^{\text{‐}1}\sim {\left({\gamma }_{\xi }{p}_{\xi }\right)}^{\text{‐}1}.$E68

It follows from Eq. (68) that at close values of γ_ξ, the values of $\nabla ·{\mathbf{u}}_{\xi }$ are inversely proportional to p_ξ. As a result, the material with a lower pressure relaxes more actively on relief, with an opposite trend in compression. However, according to the physics of the process, pressure relaxation should take place in any case. The predicted pressure for the lower pressure material can even turn out to become negative. To fix this flaw, in the case of cell expansion, method 5 uses a requirement that relative rather than absolute pressure increments of the materials should be equal. Let us write the modified equation (Eq. (65)) as

$\nabla ·{\mathbf{u}}_{\xi }\approx -\frac{\Delta p_{\xi }}{p_{\xi }}·\frac{p_{\xi }}{\tau p_{\xi }\text{d}{p}_{\xi }/\text{d}{p}_{\xi }},$E69

require that the condition $\Delta p_{\xi }^{}/p_{\xi }^{}=\Delta p_{\zeta }^{}/p_{\zeta }^{}$ is fulfilled, and obtain formula (70):

${\lambda }_{\xi }=-\frac{p_{{}_{\xi }}}{p_{{}_{\xi }}{c}_{{}_{\xi }}^2}·\frac{1}{{\displaystyle \sum {\beta }_{\text{k}}{p}_{\text{k}}/({\rho }_{\text{k}}{c}_{\text{k}}^2)}}.$E70

This value of λ_ξ will also be used in Eq. (32) to determine ψ_ξ, which meets the requirement (31) at $\Delta {e^{\prime }}_{\xi }=0$.

3.1.6. Method 7 based on the equal velocity increments model

This model has been proposed in [19]. The assumption that the materials’ velocity increments are equal that is the consequence of the fact that the set of gas dynamics equations is solved in the one‐velocity approximation. Since the materials’ velocities in the mixed cell are equal at any time, it is natural that the changes in the materials’ velocities at every timestep will also be equal. One can treat changes in physical quantities over a timestep as those resulting from the propagation of some perturbations. If one assumes that these perturbations are plane acoustic waves, for which $\Delta \rho /\rho =\Delta u/c$, then for $\nabla ·{\mathbf{u}}_{\xi }$ one can write the following expression:

$\nabla ·{\mathbf{u}}_{\xi }\approx -\frac{\Delta {\rho }_{\xi }}{{\rho }_{\xi }\tau }=-\frac{\Delta u_{\xi }}{c_{\xi }\tau }.$E71

Considering that the materials’ velocity increments $\Delta u_{\xi }$ in the mixed cell are assumed to be equal, for $\nabla ·{\mathbf{u}}_{\xi }$ we obtain the following relation (72):

$c_{\xi }\nabla ·{\mathbf{u}}_{\xi }=c_{\zeta }\nabla ·{\mathbf{u}}_{\zeta }.$E72

The set of algebraic equations (72) and (25) is closed, and its solution can be given by

$\nabla ·{\mathbf{u}}_{\xi }^n={\lambda }_{\xi }\nabla ·{\mathbf{u}}^n,$E73

where

${\lambda }_{\xi }={\left(c_{\xi }{\displaystyle \sum _k\frac{{\beta }_k}{c_k}}\right)}^{-1}=\frac{c}{c_{\xi }},$E74

with

$$c={\left({\displaystyle \sum _k\frac{{\beta }_k}{c_k}}\right)}^{-1}.$$

This value of λ_ξ will also be used in Eq. (32) to determine ψ_ξ, which satisfies the requirement (31) at $\Delta {e^{\prime }}_{\xi }=0$.

3.1.7. Pressure relaxation methods

The use of models 5–7 as single‐stage methods in real simulations is associated with a number of problems that sometimes lead to inconsistent results. All these cases are related to the absence of the pressure relaxation mechanism for materials in mixed cells. The analysis shows that, even in model 6, despite the equal pressure increments of the materials at a timestep, equality of pressures at tⁿ⁺¹ is not always the case. The situation in the other two models is even worse.

For this reason, methods 5–7 that can be used as single‐stage methods are combined with subcell pressure relaxation. Next, we consider two known relaxation methods. In all the two‐stage isotropic closure models, cell divergence is redistributed among the materials only if it is different from zero. As for the second (subcell) stage of the models, it involves interactions between the materials if these are in the nonequilibrium state with no mandatory requirement that the divergence should be nonzero.

3.1.7.1. The PR method

Let us consider the PR pressure relaxation method proposed in [20]. As the materials occupy finite volumes in the mixed cell, equilibration of the materials’ pressures occurs not instantaneously (instantaneous pressure equilibration occurs only at the points of the surface, along which the materials contact each other), but over some time in several timesteps.

This method calculates $\nabla ·{\mathbf{u}}_{\xi }$ at the timestep in two stages:

$\nabla ·{\mathbf{u}}_{\xi }=\nabla ·{\mathbf{u}}_{\xi 1}+\nabla ·{\mathbf{u}}_{\xi 2}.$E75

In Eq. (75), $\nabla ·{\mathbf{u}}_{\xi \text{1}}$ is the material's divergence at the first stage obtained by one of the above methods. The second stage includes pressure relaxation of the materials. The second stage imposes the requirement that both $\nabla ·\mathbf{u}$ and the total internal energy, i.e., ΔE₂ = 0, should remain unchanged at that stage. Pressure relaxation is implemented by calculating additional divergences of the materials $\nabla ·{\mathbf{u}}_{\xi 2}$ by the formulas

$\nabla ·{\mathbf{u}}_{\xi 2}=-\frac{\Delta p_{\xi }}{\tau {\rho }_{\xi }{c}_{\xi }^2},$E76

where 

$\Delta p_{\xi }=A\frac{c\tau }{h}\left(p-p_{\xi }\right),$E77

where p is the average pressure. Expression (76) was derived using a known relation in the adiabatic approximation, $\Delta p=-\rho c^2\tau \nabla ·\mathbf{u}$. The factor $c\tau /h$, equal to the ratio of the timestep to the characteristic pressure equilibration time of the mixed cell $h/c$, determines the fraction of the materials’ pressure difference, by which the materials’ pressure equalizes. According to the meaning of expression (77), A ∼ 1; in this case, the materials’ pressures will not relax over a single timestep.

As p in the equilibration algorithm (not to confuse with the average pressure at the basic stage done by Eq. (23)) we take Eq. (62). Then, the requirement that the cell volume should be constant at this stage, ${\displaystyle \sum {\beta }_{\xi }\nabla ·{\mathbf{u}}_{\xi 2}}=0$, will be satisfied automatically. The choice of the formula for p is ambiguous. For example, the choice based on Eq. (23) may prove to be unbeneficial. Let us illustrate this with the following example. Suppose a mixed cell contains two ideal gases having the same EOS γ₁ = γ₂ = γ, but disparate pressures and volume fractions. Let p₁ = 1000, β₁ = 0.9 and p₂ = 1, β₂ = 0.1. Using Eqs. (23) and (25) combined with the relation ${\rho }_{\xi }{c}_{\xi }^2={\gamma }_{\xi }{p}_{\xi }^{}$, one can easily calculate that p = 10. Thus, the cell‐average pressure calculated by Eq. (23) is a hundred times lower than the pressure in the material occupying nearly the whole cell. Formula (62) is free of this flaw and, as shown above, has a certain mathematical basis.

This method results in the internal energy exchange between the materials. Indeed, let us represent the total change in the materials’ internal energies as

$\text{Δ}E=-P_{+}\text{Δ}{V}_{+}-P_{-}\text{Δ}{V}_{-},$E78

where $\begin{array}{l}{P}_{+}={\displaystyle \sum _{}{\beta }_{\xi }{p}_{\xi }}\text{ if }{p}_{\xi }>p;\\ P_{-}={\displaystyle \sum _{}{\beta }_{\xi }{p}_{\xi }}\text{ if }{p}_{\xi }\le p\end{array}$ and $\text{Δ}{V}_{+}$ and $\text{Δ}{V}_{-}$ are the volume changes of these materials. With pressure equilibration, the materials with $P_{+}$ expand, so $\text{Δ}{V}_{+}$ > 0 and $\text{Δ}{V}_{-}$ < 0. As far as it follows from the volume conservation condition that $\text{Δ}{V}_{+}=\left|\text{Δ}{V}_{-}\right|$, and by definition $P_{+}$>$P_{-}$, $\text{Δ}E$ in the pressure equilibration procedure by Eq. (78) will always be negative. This situation is associated with the fact that motion of interfaces causes internal (subcell) motion in the cell, and part of the cell's internal energy converts into the subcell kinetic energy. Since the subcell kinetic energy is not taken into account in the calculations, it is returned to the materials in the form of internal energy increments $\text{Δ}{e^{\prime }}_{\xi }$ in accordance with the equation

$-\frac{\tau }{{\rho }^n}{\displaystyle \sum _{\xi }{\beta }_{\xi }(p_{\xi }^{n+1/2}+q_{\xi }^n)\nabla ·{\mathbf{u}}_{\xi 2}^{n+1/2}}+{\displaystyle \sum _{\xi }{\alpha }_{\xi }\Delta {e^{\prime }}_{\xi }=0},$E79

The quantity $\nabla ·{\mathbf{u}}_{\xi 2}^{n+1}$ present in this expression is calculated by the formula

$\nabla ·{\mathbf{u}}_{\xi 2}^{n+1}=\nabla ·{\mathbf{u}}_{\xi 2}^n\left(1-Ac\tau /h\right).$

It remains to decide how to distribute the dissipated kinetic energy among the materials (formula (79) defines only the total dissipated energy $\Delta E={\displaystyle \sum _{\xi }{\alpha }_{\xi }\Delta {e^{\prime }}_{\xi }}$). In [20], it is assumed that $\Delta {e^{\prime }}_{\xi }=\Delta e^{\prime }$. In this case, for all the materials in the mixed cell we obtain from Eq. (79) that

$\Delta {e^{\prime }}_{\xi }=\Delta e^{\prime }=\frac{\tau }{{\rho }^n}{\displaystyle \sum _{\xi }{\beta }_{\xi }\left(p_{\xi }^{n+1/2}+q_{\xi }^n\right)\nabla ·{\mathbf{u}}_{\xi 2}^{n+1/2}}.$E80

Note that this pressure relaxation approach is universal, i.e., it is independent of the way, the total velocity divergence in the mixed cell is distributed among the materials. In EGAK, it is employed for the three above‐mentioned methods. However, its application to physically inconsistent generation of pressure difference between the materials through the basic closure relation may lead to excessive energy exchange between them, so for each specific method its consistency needs to be validated by test simulations.

3.1.7.2. Method of Barlow, Hill and Shashkov (BHS)

This method is described in detail in [25]; here we briefly summarize and illustrate the basic concept of the method for the case of two materials, which is sufficient for understanding the whole method. In its complete form, the method has been developed for the multimaterial case; for details see [25].

This method assumes that the total change in the material volume over a timestep is the sum of two terms:

$\Delta V_{\xi }=\Delta V_{1\xi }+\Delta V_{2\xi },$E81

where the subscripts 1 and 2 denote the two stages of the closure model.

The first stage employs model 5, which assumes that the materials’ divergences are equal and do not require the information on the contact location in the cell. The equality of divergences means that

$\Delta V_{1\xi }={\beta }_{1\xi }\Delta V.$E82

The second (subcell) stage uses the model based on the acoustic Riemann problem (Delov's model), which calls for the reconstruction of the contact location in the mixed cell. In cell 1234, as shown in Figure 3, it is the segment AB.

Figure 3.

Contact location reconstructed after the first stage.

After the first stage, P₁ > P₂. Then, after the subcell stage, the contact will move to the location CD. The quantity of volume increment is represented by the quadrilateral ABCD determined based on the solution of the acoustic Riemann problem

$\Delta V_2=\frac{p_1-p_2}{{\rho }_1{c}_1+{\rho }_2{c}_2}{S}_{AB}{\tau }^n,$E83

where all the unindexed quantities are taken after the first stage, and S_AB is the area of the boundary between the materials.

Thus,

$\begin{array}{l}\text{Δ}{V}_{21}=\text{Δ}{V}_2,\\ \text{Δ}{V}_{22}=-\text{Δ}{V}_2.\end{array}$E84

The updated volumes of the materials calculated by formula (81) accounting for the volume increments must satisfy the following inequalities:

$V^{n+1}>V_{\xi }^{n+1}>0,$E85

which, however, is not always the case for the same reason as in Delov's method (see the remark at the end of Section 3.1.3). Therefore in [25], the authors introduce constraints on volume increments, which complicate the method significantly, especially in the multimaterial case.

3.2. Anisotropic closure models

3.2.1. Model fundamentals

Let us consider two limiting cases of contact location relative to the wave motion (shock, acoustic, or elastic wave) presented in Figure 4, in which cells contract or expand, i.e., the divergence is nonzero. In the first case (Figure 4a), most of motion is normal to the interface, so all the above models 1–7, each having its own accuracy, are suitable. In the second case (Figure 4b) most of motion occurs along the interface, while the lateral motion is insignificant and is therefore auxiliary. It means that the materials contract or expand tangentially to the interface; thus, equality of compressibilities, i.e., model 5, may be more consistent in this case. Indeed, calculations show that, for example, using models 6 and 7 for such flows in the elastoplastic case one can obtain a considerable error, while model 5 provides good accuracy.

Figure 4.

Two cases of contact location relative to wave motion.

The above fact implies that closure models 5–7 are inappropriate for modeling such flows. Thus, to ensure acceptable modeling accuracy for the two different types of flow (in different directions relative to the interface), different closure relations need to be used. For this purpose, two‐stage models are proposed.

3.2.2. The ACM‐1 model

At the first stage in the anisotropic model ACM‐1 [34], matter in the mixed cell moves as a whole, and all nonuniformities (including the interface) are assumed to be frozen. The freezing condition in terms of closing in the first approximation means that the materials’ divergences are equal.

The second stage includes relaxation of pressure (and stresses) concurrent with such motion. The work [34] suggests using the above pressure relaxation procedure at this stage with the degree of relaxation made dependent upon the mutual orientation of the velocity direction and the interface. If the velocity is normal to the interface, the pressure relaxation is the highest, and if the velocity is directed along the interface, it is the smallest.

In the implementation of the ACM‐1 model, two stages are used to calculate $\nabla ·{\mathbf{u}}_{\xi }$ at a timestep

$\nabla ·{\mathbf{u}}_{\xi }=\nabla ·{\mathbf{u}}_{\xi 1}+\nabla ·{\mathbf{u}}_{\xi 2}.$E86

In Eq. (86), $\nabla ·{\mathbf{u}}_{\xi 1}$ is calculated at the first stage in accordance with closure model 1:

$\nabla ·{\mathbf{u}}_{\xi 1}=\nabla ·\mathbf{u}.$E87

The second stage includes pressure relaxation of the materials in mixed cells according to the algorithm, basically similar to the algorithm described in Section 3.1.7.1. The only distinguishing feature is that for the ACM‐1 model, the coefficient A in Eq. (76) depends on the mutual orientation of the velocity and the interface. The total divergence is written as the sum of two components:

$\nabla ·\mathbf{u}=\nabla ·{\mathbf{u}}_{\tau }+\nabla ·{\mathbf{u}}_n,$E88

obtained by velocity decomposition to two components: along the interface (u_τ) and normal to it (u_n). We also assume that

$A=A_0·\frac{\nabla ·{\mathbf{u}}_n}{\nabla ·\mathbf{u}},$E89

where A₀ is some constant.

Thus, the coefficient A is a variable in this case. If the velocity is normal to the interface (Figure 4a), $\nabla ·{\mathbf{u}}_n=\nabla ·\mathbf{u}$ and A = A₀; this is the case with the highest pressure relaxation. If the velocity is directed along the interface (Figure 4b), $\nabla ·{\mathbf{u}}_n=0$ and A = 0, i.e., there is no pressure relaxation at all, and only the first stage of the closure method is in effect, i.e., formula (87).

The constant A₀ in Eq. (89) must be determined based on test simulation results. In [34], its value was determined in several simulations, and it proved to coincide with the value of the constant A in Section 3.1.7.1, i.e., A₀ = 1.

3.2.3. The ACM‐2 model

The anisotropic model ACM‐2 is formulated as follows. We divide the divergence of the entire cell and its materials into two components: normal and tangential (relative to the contact location):

$\begin{array}{l}\nabla ·\mathbf{u}=\nabla ·{\mathbf{u}}_n+\nabla ·{\mathbf{u}}_{\tau }\\ \nabla ·{\mathbf{u}}_{\xi }=\nabla ·{\mathbf{u}}_{\xi n}+\nabla ·{\mathbf{u}}_{\xi \tau }\end{array}.$E90

Along the interface, the materials are assumed to have equal compressibilities, i.e., the distribution of the corresponding divergence component among the materials is described by the relation

$\nabla ·{\mathbf{u}}_{\xi \tau }=\nabla ·{\mathbf{u}}_{\tau }.$E91

As for the divergence in the direction normal to the interface, it can be distributed to the materials using any of closure models 5–7. In this work, we use model 7, which, as shown in [40, 41], is the most widely applicable with respect to modeling different kinds of flow. It follows from it that

$\begin{array}{l}\nabla ·{\mathbf{u}}_{\xi n}={\lambda }_{\xi }\nabla ·{\mathbf{u}}_n\\ {\lambda }_{\xi }={\left(c_{\xi }{\displaystyle \sum _k\frac{{\beta }_k}{c_k}}\right)}^{-1}\end{array}.$E92

Once this part of divergence is distributed to the materials, relaxation of their pressures is carried out by the algorithm described in Section 3.1.7.1, which makes an additional contribution,$\nabla ·{{\mathbf{u}}^{\prime }}_{\xi n}$, to the divergence $\nabla ·{\mathbf{u}}_{\xi n}$:

${\widetilde{\nabla \cdot \mathrm{u}}}_{\xi n}=\nabla ·{\mathbf{u}}_{\xi n}+\nabla ·{{\mathbf{u}}^{\prime }}_{\xi n}.$E93

The ultimate formula for the distribution of $\nabla ·\mathbf{u}$ among the materials will be

$\nabla ·{\mathbf{u}}_{\xi }={\widetilde{\nabla \cdot \mathrm{u}}}_{\xi n}+\nabla ·{\mathbf{u}}_{\xi \tau }.$E94

4. Artificial viscosity

For mixed cells, closure relations for calculating divergences of constituent materials are insufficient. Such cells require an additional relation to determine the artificial viscosity for the cell as a whole and for each individual material. Two approaches can be used to calculate the viscosity of the materials.

The first approach is based on viscosity calculations directly for the materials based on their individual parameters. The viscosity for matter as a whole is then calculated based on the procedure used to calculate the average pressure from the materials’ pressures.

The second approach is based on the calculations of the viscosity for the cell as a whole based on the cell‐average parameters of matter and its distribution to the materials according to some assumptions on the way of its distribution.

The first approach is more expensive both in terms of research effort, and in terms of computations; therefore, the less complicated second approach has been developed more widely. This work mostly considers viscosity definition procedures based on the second approach.

4.1. Artificial viscosity for matter as a whole

The viscosity of matter as a whole in EGAK is calculated in the cells, and it is a combination of the von Neumann and Richtmyer type quadratic viscosity and linear viscosity

$\{\begin{array}{l}{q}^n=C_1·{\rho }^n{\left(h^n\nabla ·{\mathbf{u}}^n\right)}^2+C_0·{\rho }^n·с·h^n\nabla ·{\mathbf{u}}^n\text{ if }\nabla ·{\mathbf{u}}^n<0\\ q^n=0\text{ if }\nabla ·{\mathbf{u}}^n\ge 0\end{array}.$E95

In Eq. (95), C₁ = 1 and C₀ = 0.2 are the fixed coefficients. In addition, expression (95) contains the characteristic dimension h and the divergence $\nabla ·{\mathbf{u}}^n$ of the cell. Here, we will not address the issues of determining these quantities, as it does not matter for our purposes. Let us consider the approaches to viscosity distribution to materials.

4.2. Artificial viscosity of materials in mixed cells

The research summarized below has been carried out in [36]. Viscosity calculations for materials represent an ambiguous problem; for solving it correctly, it is usually insufficient to have data on the subcell behavior of the materials. The way of the materials’ viscosity representation governs the distribution of energy dissipated in the cell on shock propagation to the materials. The problem of energy distribution to the materials is a subcell problem that lacks information needed for obtaining the exact solution. In reality, the shock front width is generally much smaller than the mesh spacing, so the shock propagates through each material in the heterogeneous mixture practically independently. Shock parameters in each material are different, and they differ from average values in the mixed cells, so it is practically impossible to determine uniquely the fraction of dissipated energy for each material. Pressures and velocities of the materials are different behind the shock, and the processes of pressure and velocity relaxation provide additional redistribution of dissipated energy between the materials.

When considering the approaches to viscosity definition for the materials, let us characterize these approaches in terms of dissipated energy distribution among the materials and resulting changes in the materials’ pressures at a timestep. Before we consider different kinds of material viscosities, let us note that the cell‐average value of viscosity is determined using Eq. (23), i.e., $q={\displaystyle \sum {\psi }_{\xi }{q}_{\xi }},$ ${\psi }_{\xi }={\beta }_{\xi }{\lambda }_{\xi }$, where λ_ξ is a normalizing factor (see Section 3), which imposes a constraint on the closure methods (all the procedures described below have been studied using EGAK on methods 5–7 and ACM‐1).

In accordance with the difference scheme of EGAK, the viscosity‐related change in the materials’ specific internal energy over a timestep is given by the formula

$\Delta e_{\xi }=-\frac{q_{\xi }^n}{{\rho }_{\xi }^n}·{\lambda }_{\xi }^n·\tau ·\left(\nabla ·{\mathbf{u}}^{n+1/2}\right).$E96

The viscosity‐related change in the materials’ pressure at a timestep can be obtained as follows. For adiabatic flows, it holds true that

$\Delta p\approx {\left(\frac{\partial p}{\partial \rho }\right)}_S·\Delta \rho \approx -\rho ·c^2·\tau ·\nabla ·\mathbf{u}.$E97

considering that $\Delta \rho \approx -\rho ·\nabla ·\mathbf{u}·\tau$.

Using the EOS p = P(ρ,e), the total pressure change at a timestep in the general case can be represented as

$\Delta p\approx {\left(\frac{\partial p}{\partial \rho }\right)}_e\Delta \rho +{\left(\frac{\partial p}{\partial e}\right)}_{\rho }\Delta e.$E98

For adiabatic flows, the energy increment is calculated by the formula

$\Delta e=-\frac{p}{\rho }·\tau ·\nabla ·\mathbf{u}.$E99

Substituting this expression into Eq. (98) and comparing it with Eq. (96), considering Eq. (97), we obtain

${\left(\frac{\partial p}{\partial \rho }\right)}_e+{\left(\frac{\partial p}{\partial e}\right)}_{\rho }·\frac{p}{{\rho }^2}=c^2.$E100

On shock propagation, the energy increment is calculated by the formula

$\Delta e=-\frac{p+q}{\rho }·\tau ·\nabla ·\mathbf{u}.$E101

Using Eq. (100) from Eq. (98), we obtain the total pressure increment in the form of

$\Delta p\approx -\left[\rho ·c^2+{\left(\frac{\partial p}{\partial e}\right)}_{\rho }·\frac{q}{\rho }\right]\tau ·\nabla ·\mathbf{u}.$E102

It follows from Eq. (102) that the materials’ viscosity‐related pressure increment at a timestep equals

$\Delta p_{q\xi }\approx -{\left(\frac{\partial p_{\xi }}{\partial e_{\xi }}\right)}_{\rho }·\frac{q_{\xi }^n}{{\rho }_{\xi }^n}·{\lambda }_{\xi }^n·\tau ·\left(\nabla ·{\mathbf{u}}_{\xi }^{n+1/2}\right).$E103

Now, let us consider models to material viscosity definition. Six models have been explored. Table 1 provides their descriptions and formulas and changes in specific energy and pressure in accordance with Eqs. (96) and (103).

	Description of the approach to calculating ${\mathbf{q}}_{\mathbf{\xi }}$	Formula	$\mathbf{\Delta }{\mathbf{e}}_{\mathbf{\xi }}$	$\mathbf{\Delta }{\mathbf{p}}_{\mathbf{\xi }}$
1	Equal to the cell‐average viscosity	$q_{\xi }=q$	${\lambda }_{\xi }/{\rho }_{\xi }$	${\left(\frac{\partial p_{\xi }}{\partial e_{\xi }}\right)}_{\rho }·\frac{{\lambda }_{\xi }}{{\rho }_{\xi }}$
2	Viscosity with its quantities ${\rho }_{\xi },h_{\xi }={\beta }_{\xi }h,\nabla ·{\mathbf{u}}_{\xi }$	$q_{\xi }=q\frac{{\rho }_{\xi }{\beta }_{\xi }^2{\lambda }_{\xi }^2}{{\displaystyle \sum {\rho }_k{\beta }_k^3{\lambda }_k^3}}$	${\beta }_{\xi }^2{\lambda }_{\xi }^3$	${\left(\frac{\partial p_{\xi }}{\partial e_{\xi }}\right)}_{\rho }{\beta }_{\xi }^2{\xi }_{\xi }^3$
3	Proportional to material densities $q_{\xi }\sim {\rho }_{\xi }$	$q_{\xi }=q\frac{{\rho }_{\xi }}{{\displaystyle \sum {\rho }_k{\beta }_k{\xi }_k}}$	${\lambda }_{\xi }$	${\left(\frac{\partial p_{\xi }}{\partial e_{\xi }}\right)}_{\rho }{\lambda }_{\xi }$
4	Same energy increment $\Delta e_{\xi }=\Delta e_{\zeta }$	$q_{\xi }=q\frac{{\rho }_{\xi }}{\rho {\lambda }_{\xi }}$	$\Delta e_{\xi }=\Delta e_{\varsigma }$	${\left(\frac{\partial p_{\xi }}{\partial e_{\xi }}\right)}_{\rho }$
5	Same increment $\Delta p_{\xi }$$\Delta p_{\xi }=\Delta p_{\zeta }$	$q_{\xi }=q\frac{{\rho }_{\xi }}{{\lambda }_{\xi }{\left(\frac{\partial p_{\xi }}{\partial e_{\xi }}\right)}_{\rho }{\displaystyle \sum \frac{{\beta }_k{\rho }_k}{{\left(\partial p_k/\partial e_k\right)}_{\rho }}}}$	$\frac{1}{{\left(\partial p_{\xi }/\partial e_{\xi }\right)}_{\rho }}$	$\Delta p_{\xi }=\Delta p_{\varsigma }$
6	Proportional to $\Delta p_{\xi }$ in adiabatic approximation $q_{\xi }\sim △p_{\xi }$	$A\times {\rho }_{\xi }{c}_{\xi }^2={\left(\partial p_{\xi }/\partial e_{\xi }\right)}_{\rho }{q}_{\xi }/{\rho }_{\xi }$, where A is the proportionality factor	$\frac{{\lambda }_{\xi }{\rho }_{\xi }{c}_{\xi }^2}{{\left(\partial p_{\xi }/\partial e_{\xi }\right)}_{\rho }}$	${\lambda }_{\xi }{\rho }_{\xi }{c}_{\xi }^2$

Table 1.

Models to calculate the viscosity and specific energy and pressure increments of the materials.

Thus, in these six models to the materials’ viscosity definition, distribution of dissipated energy to the materials is differently dependent on their density, speed of sound, and divergence. The choice of one model or another depends both on the closure method, and on the modeled problem. Based on the test calculations in [36], the best performance was demonstrated by model 3.

5. Method for calculating mixed cells with vacuum

One of the materials available in EGAK is a zero‐pressure “vacuum.” For the case of vacuum, a special algorithm has been developed, which is the same for closure methods 1 and 5–7. The development of this algorithm was motivated by the fact that the pressure for vacuum is specified rather than controlled by the closure method.

In a mixed cell with vacuum, two cases are possible: $\nabla ·\mathbf{u}>0\text{ and }\nabla ·\mathbf{u}\le 0$.

The case of $\nabla ·\mathbf{u}>0$. In this case, it is assumed that

$\nabla ·{\mathbf{u}}_{\xi }=\nabla ·\mathbf{u}.$E104

The case of $\nabla ·\mathbf{u}\le 0$. In this case, the cell volume is assumed to decrease only due to a decrease in the vacuum volume, while the volumes of the other materials change only after the vacuum gets closed. This can be represented by the following formula:

$\begin{array}{l}\nabla ·{\mathbf{u}}_{\xi }=-\frac{2\left({\rho }_{\xi }^{n+1}-{\rho }_{\xi }^n\right)}{\left({\rho }_{\xi }^{n+1}+{\rho }_{\xi }^n\right)\tau },\\ {\rho }_{\xi }^{n+1}=\frac{V^n{\beta }_{\xi }^n{\rho }_{\xi }^n}{V^{n+1}{\beta }_{\xi }^{n+1}},\\ {\beta }_{\xi }^{n+1}=\frac{{\beta }_{\xi }^n\left(1-{\beta }_{vac}^{n+1}\right)}{\left(1-{\beta }_{vac}^n\right)},\\ {\beta }_{vac}^{n+1}=\frac{{\beta }_{vac}^n{V}^n-\min \left\{\Delta V;{\beta }_{vac}^n{V}^n\right\}}{V^{n+1}}.\end{array}$E105

The following is proposed for the anisotropic methods ACM‐1 and ACM‐2: we represent the total divergence, like in Eq. (88), as a sum of two items, $\nabla ·{\mathbf{u}}_{\tau }$ and $\nabla ·{\mathbf{u}}_n$. If the cell expands, i.e., if $\nabla ·\text{u}\ge 0$, then $\nabla ·{\mathbf{u}}_{\xi }=\nabla ·\mathbf{u}$.

When the cell contracts, i.e., when $\nabla ·\mathbf{u}<0$:

‐ if $\nabla ·{\mathbf{u}}_n<0$, then $\nabla ·{\mathbf{u}}_{\xi }=\nabla ·{\mathbf{u}}_{\xi \tau }^{}$, $\nabla ·{\mathbf{u}}_{vac}=\nabla ·\mathbf{u}-{\beta }_{\xi }·\nabla ·{\text{u}}_{\xi }/{\beta }_{vac}$;

‐ if $\nabla ·{\mathbf{u}}_n\ge 0$, then $\nabla ·{\mathbf{u}}_{\xi }=\nabla ·\mathbf{u}$.

Let us consider the cases illustrated in Figure 4. Suppose the cell is contracting. Then, if the velocity is normal to the interface (Figure 4a), we have $\nabla ·{\mathbf{u}}_n=\nabla ·\mathbf{u}$, $\nabla ·{\mathbf{u}}_{\tau }=0$, and as $\nabla ·{\mathbf{u}}_n<0$, we obtain $\nabla ·{\mathbf{u}}_{\text{vacuum}}=\nabla ·\mathbf{u}$, i.e., only the vacuum is contracting. If the velocity is directed along the interface (Figure 4b), $\nabla ·{\text{u}}_n=0$ and $\nabla ·{\mathbf{u}}_{\xi }=\nabla ·\mathbf{u}$.

6. Test problems and results

The author does not have results of testing all of the above closure methods on the problems modeled in the section, so below he basically presents the results for methods 1(P), 5($\nabla ·\mathbf{u}$), 6(Δp), 7(Δu), and, correspondingly, for the methods $\nabla ·\text{u}$‐PR, Δp‐PR, Δu‐PR, and ACM‐1 and ACM‐2, which have been developed with his direct participation. These methods have been tested on numerous problems, including those not included in this work (see [40, 41]). It does not seem possible to present results of all such calculations, so the author confines himself to three one‐dimensional and one two‐dimensional problem having analytical solutions. All the 1D calculations were done in Lagrangian variables, and the 2D one in Eulerian.

The following unified types of data processing are provided for all the calculations.

• Tables with quantities characterizing the order of convergence of the integral error of basic quantities in the calculations in the L₁ norm at a reference time.

• Tabulated values of basic quantities in mixed cells at a reference time.

The error is calculated by formula (106):

$\delta y={\Vert y_{\text{comp}}-y_{\text{exact}}\Vert }_1=A{h}^{\sigma },$E106

where h is the initial mesh spacing and y_comp and y_exact are the calculated and the exact value of the quantity at the cell center, respectively.

In EGAK calculations, mixed cells were of the same size as pure cells in pure‐cell calculations. In the mixed cell calculations, domain coordinates were shifted to the right at a distance of δx = h/2, where h is the mesh spacing in the corresponding calculation. In other studies, the size of mixed cells was doubled, but their number was less than one cell.

In addition, a two‐dimensional problem of elastic wave propagation in a thin plate is discussed, for which only EGAK results are presented and the analytical solution is given in [42]. For this problem, we have calculated the velocity of a longitudinal elastic wave in a plate. In [42], we have also solved this problem numerically using EGAK in the absence of mixed cells, which demonstrated good accuracy of the calculations in this setup. In [34], this problem has been solved numerically in the setup, in which interfaces are misaligned with grid lines thus producing mixed cells.

6.1. The water‐air shock tube problem

Setup. The problem has the following initial conditions [43]:

$\left(\gamma ,p_{\infty },\rho ,e,p,\mathbf{u}\right)=\{\begin{array}{l}(\text{4}\text{.4, 6}·{\text{10}}^{\text{8}}{\text{, 10}}^{\text{3}}\text{, 1}\text{.07}·{10}^6\text{, }{10}^9\text{, 0), if 0 }\le \text{ x }\le \text{ 0}\text{.7,}\\ (\text{1}\text{.4, 0, 50, 5}·{10}^4\text{, }{10}^6\text{, 0), if 0}\text{.7<x}\le \text{ 1}\text{.}\end{array}$E107

The EOS of water is $p=(\gamma -1)\rho e-\gamma p_{\infty }$, for which the squared speed of sound is calculated by the formula $c^2=\gamma (\gamma -1)\left(e-p_{\infty }/\rho \right)=\gamma \left(p+p_{\infty }\right)/\rho$.

The final time of the calculation is t = 2.2 · 10^-4. The exact solution to the problem has been obtained in [44]. The calculations were carried out on meshes having 250, 500, and 1000 cells.

Results. Table 2 presents the values of the factor A and the order of convergence of the integral error in the basic quantities at t = 2.2 · 10^-4, and Table 3 shows the exact and calculated values of the basic quantities for the closure methods, for which data are available in publications. Figure 5 shows the L₁ norm of the absolute error as a function of h.

Figure 5.

L₁ norm of the absolute error as a function of h.

Method	p ($\times {\mathbf{10}}^{-\mathbf{6}}$)		ρ		e ($\times {\mathbf{10}}^{-\mathbf{3}}$)		u
	A ($\times {\mathbf{10}}^{-\mathbf{2}}$)	σ	A ($\times {\mathbf{10}}^{-\mathbf{2}}$)	σ	A ($\times {\mathbf{10}}^{-\mathbf{2}}$)	σ	A ($\times {\mathbf{10}}^{-\mathbf{2}}$)	σ
Pure (EGAK)	5.13	0.84	1.92	0.86	1.3	0.84	6.0	0.91
$\nabla ·\mathbf{u}$	4	0.79	12.4	1.10	14.5	1.16	3.56	0.79
p	4.43	0.80	2.5	0.89	1.14	0.81	1.8	0.71
Δp	8.2	0.83	5.08	0.93	5.18	0.96	3.62	0.79
Δu	8.29	0.84	16.2	0.96	9.41	0.90	15.1	0.92
$\nabla ·\mathbf{u}$‐PR	4.1	0.79	2.28	0.89	1.03	0.80	1.57	0.70
Δp‐PR	7.11	0.87	4.20	0.88	5.78	0.92	6.97	0.88
Δu‐PR	6.26	0.85	2.49	0.83	2.04	0.83	4.27	0.84
Pure (Delov)	4.96	0.97	2.15	0.94	1.06	0.87	3.73	0.91
Delov	7.42	0.94	8.45	0.96	18.4	0.98	22.1	1.03
Pure (K&S, Tipton)	7.65	0.99	2.18	0.97	0.7	0.90	4.12	0.92
K&S	6.95	0.99	3.63	1.01	3.09	1.01	7.98	1.03
Tipton	4	0.79	12.4	1.10	14.5	1.16	3.56	0.79

Table 2.

Factor A and the rate of convergence σ with mesh refinement.

Method	p1 (× 10^-7)	p2 (× 10^-7)	ρ1 (× 10^-2)	ρ2 (× 10^-2)	e1 (× 10^-5)	e2 (× 10^-5)
Exact	1.599	1.599	8.050	2.204	9.704	1.813
Pure (EGAK)	1.599	1.599	7.993	1.535	9.773	2.605
p	1.595	1.595	7.764	1.377	10.062	2.896
$\nabla ·\mathbf{u}$	3.111	0.078	8.030	0.401	9.784	0.484
Δp	99.034	1.351	9.972	4.606	10.707	7.332
Δu	39.664	0.065	8.931	0.174	10.001	0.931
$\nabla ·\mathbf{u}$‐PR	1.594	1.594	7.579	1.504	10.306	2.650
Δp‐PR	1.599	1.599	7.455	0.471	10.477	8.479
Δu‐PR	1.599	1.599	7.395	1.104	10.562	3.620
Pure (Delov)	1.599	1.599	8.113	1.409	9.629	2.835
Delov	1.595	1.594	8.044	0.249	9.711	16.015
Pure (K&S, Tipton)	1.599	1.599	7.983	1.669	9.785	2.394
K&S	1.599	1.599	7.352	1.364	10.625	2.930
Tipton	1.599	1.599	7.315	0.591	10.680	6.765

Table 3.

Exact and calculated values of the basic quantities in mixed cells on a mesh having 1000 cells at t = 2.2×10^-4.

6.2. The mixed‐material shock transition problem

Setup: The problem has been proposed in [18]. The domain $-\text{2}<x<\text{1}$ contains a mixture of two ideal gases having the following parameters: ${\rho }_{\text{1}}^{\text{0}}=\text{1, }{e}_{\text{1}}^{\text{0}}=\text{0, }{\gamma }_{\text{1}}=\text{3, }{\beta }_{\text{1}}^{\text{0}}=\text{0}\text{.5}$ (material 1) and ${\rho }_{\text{2}}^{\text{0}}=\text{1, }{e}_{\text{2}}^{\text{0}}=\text{0, }{\gamma }_{\text{2}}=\text{1}\text{.2, }{\beta }_{\text{2}}^{\text{0}}=\text{0}\text{.5}$ (material 2). A constant velocity of u = 2 is given at the left boundary. Due to the specified boundary condition, a strong shock starts moving across the mixture. The problem has an analytical solution obtained in [20] assuming that the materials’ pressures are equal.

The values of densities are determined based on the condition that the shock is strong for each of the materials: ${\rho }_{\xi }=\left({\gamma }_{\xi }+1\right)/\left({\gamma }_{\xi }-1\right){\rho }_{\xi }^0$. It is implicitly supposed here that only one shock travels across the gases (additional waves reverberating between the interfaces are ignored). The volumes occupied by the materials behind the shock equal $V_{\xi }=V_{\xi }^0{\rho }_{\xi }^0/{\rho }_{\xi }$. The average density behind the shock front then equals

$\rho =\frac{{\displaystyle \sum M_{\xi }}}{{\displaystyle \sum V_{\xi }}}=\frac{{\displaystyle \sum V_{\xi }^0{\rho }_{\xi }^0}}{{\displaystyle \sum V_{\xi }}}=\frac{{\displaystyle \sum V_{\xi }^0}}{{\displaystyle \sum V_{\xi }\left({\gamma }_{\xi }-1\right)/\left({\gamma }_{\xi }+1\right)}}{\rho }^0.$E108

The laws of conservation of mass, momentum, and total energy for the shock (the Rankine‐Hugoniot relations) traveling across each of the materials make it possible to find the parameters of the gases behind the shock front:

$D=\rho /\left(\rho -{\rho }^0\right)\mathbf{u},$E109

${\mathbf{u}}^2=p\left(1/{\rho }^0-1/\rho \right),$E110

$e=0.5p\left(1/{\rho }^0-1/\rho \right).$E111

Using Eq. (109), we obtain the shock velocity, using Eq. (110), the average pressure, and using Eq. (111), the average energy of the mixture. The pressures of the materials are equal to the average pressure due to the assumption we made, and the energies of the materials can be obtained from the corresponding EOS.

This problem was calculated on a mesh having 600 cells. This problem is also of interest in terms of examining the effect of the approach to calculate the artificial viscosity for the materials.

Results: This problem differs from the two problems discussed above. First, there are no pure cells, so pure‐cell calculations are inapplicable in this case. Second, only some of the above dependences can be obtained for this problem, and these are presented below. In particular, it has almost no sense to perform convergence calculations for this problem, because the steady‐state solution in the mixed cells does not depend on the mesh spacing. Table 4 presents the values of the parameters in the mixed cell with x = 0.2 at t = 2 for the materials behind the shock obtained using Eqs. (109)–(111) and in the calculations (for the materials’ viscosities by approach 3). Table 4 shows the results obtained with different viscosities for the method $\nabla ·\text{u}$‐PR.

Method	D	p1	p2	ρ1	ρ2	e1	e2
Exact	2.839	5.677	5.677	2.0	11	1.419	2.581
p	2.997	5.992	5.992	1.694	13.422	1.769	2.232
$\nabla ·\mathbf{u}$	3.456	13.219	0.581	2.379	2.379	2.778	1.222
Δp	2.830	5.753	5.477	2.032	10.596	1.416	2.584
Δu	2.827	5.886	5.324	2.053	10.372	1.434	2.567
$\nabla ·\mathbf{u}$‐PR	2.859	5.715	5.715	1.956	11.253	1.461	2.539
Δp‐PR	2.837	5.668	5.668	2.011	10.939	1.409	2.591
Model 1	2.916	5.820	5.820	1.844	12.013	1.578	2.422
Model 2	2.965	5.909	5.909	1.762	12.720	1.677	2.323
Model 3	2.817	5.640	5.640	2.047	10.754	1.378	2.622
Model 4	2.859	5.715	5.715	1.956	11.253	1.461	2.539
Model 5	2.535	5.049	5.049	3.503	7.697	0.720	3.279
Model 6	3.11	6.207	6.207	1.545	15.587	1.991	3.110

Table 4.

Exact and calculated values of the basic quantities in the cell with × = 0.2 at t = 2 on a mesh having 600 cells (model 1 is the way of cell viscosity distribution to the materials for the closure method Δu‐PR).

6.3. 2D problem of elastic wave propagation in a plate

Next, we consider a 2D problem, in which a longitudinal elastic wave propagates in a thin plate and the wave velocity for which has been obtained analytically in [42]. The problem has been calculated using EGAK in the absence of mixed cells in [42] and with mixed cells in [34]. As a surrounding medium in the problem, we used air or vacuum.

Setup: In the calculations, a titanium projectile of length L = 10 cm flying at a velocity of v₀ = 0.01 km/s surrounded by air or vacuum hits a “rigid” wall. This produces an elastic wave in the projectile traveling toward its rear surface. Figure 6 shows a schematic drawing of the initial problem geometry; H = 1 is the thickness of the wall. The calculations were carried out on a fixed mesh having a mesh spacing of h = 0.2 cm. The parameters of the EOS and the model of matter for the materials are shown in Tables 5 and 6.

Figure 6.

Geometry of the problem of elastic wave propagation in a plate.

ρ0 (g/cm^3)	c0 (km/s)	n	Γ
4.5	4.842	3.4243	1.18

Table 5.

Mie‐Grüneisen EOS parameters.

	b (GPa)	k	c	m	Cv (kJ/(g×K))	Tm (K)	G (GPa)	ν
1.098	1.092	0.93	0.014	1.1	580 × 10^-6	1878	43	0.32

Table 6.

Johnson‐Cook model parameters.

The mesh was constructed in such a way that the projectile was initially surrounded by mixed cells containing titanium and vacuum (air) with a ratio β = 0.5. We also carried out calculations on an oblique mesh with a varied volume fraction. The field of volume fractions of titanium and a mesh fragment for this simulation are shown in Figure 7.

Figure 7.

Distribution of volume fractions in the simulation on a fragment of the oblique mesh.

Results: In this problem, there is certain difficulty determining the longitudinal wave velocity. To address this difficulty, the following approach has been proposed in [42]. Suppose the elastic wave front does not “smear” as it propagates in the material. As the projectile hits the “rigid” wall, the velocity of the projectile material behind the wave should be zero. Therefore, the rate of deceleration of the projectile's center of mass can be related to the elastic wave velocity. Figure 8 shows the time‐history plots for the velocity of the projectile's center of mass for three simulations. The plots demonstrate that these dependences are nicely approximated by a linear function (v = v₀ - At). The time it takes the step‐like elastic wave to travel all the way along the projectile is T = v₀/A. Here, T is the time, at which v = 0, which corresponds to the time of wave traveling all the way along the projectile. The longitudinal wave velocity is then defined as c_w = L/A. The error associated with the displacement of the projectile's rear end as the wave travels all the way along the projectile can be neglected, because the material's velocity is small compared to the wave velocity. Table 7 shows the values of longitudinal wave velocities for all simulations.

Figure 8.

Velocity of the plate's center of mass as a function of time in calculations with different closure conditions for mixed cells (simulations with air).

Closure method	Surrounding medium	cw (km/s)
Exact		5.3
Pure cells	Air	5.2
Δu‐PR	Vacuum	4.8
ACM‐1	Vacuum	5.1
ACM‐2	Vacuum	5.1
ACM‐1 (oblique mesh)	Vacuum	5.1
Δu‐PR	Air	4.7
ACM‐1	Air	5.1

Table 7.

Theoretical and calculated values of longitudinal elastic wave velocity in a plate.

The calculations of this 2D problem demonstrated that, for both of the anisotropic closure methods, the difference between the calculated elastic wave velocity and the exact solution is ∼4%, whereas for the method Δu‐PR it is ∼10%. No comparison with other methods was made, because the Δu‐PR method proved to be the most versatile among all the methods in EGAK as applied to a wide range of different problems.

6.4 Discussion of results and conclusions

The calculated data presented here and not included in this work demonstrate that all the methods under consideration have good convergence (the order of convergence is ∼1) to the exact solution with mesh refinement as applied to all 1D problems with interfaces. When comparing the methods, one should note that the order of convergence of calculations with closure methods is mostly governed by the order of convergence of the basic difference scheme. As for the error of the closure methods themselves, it is basically controlled by the value of the factor A in formula (106). The reader himself can choose the method he likes. However, two circumstances need to be mentioned, which are important when choosing the method. First, the methods differ in the amount of calculations. Second, the methods differ in the complications in program implementation associated with limitations in their applicability.

As for the 2D problem, the anisotropic methods have no alternative. They possess the same accuracy as the basic methods on the 1D problems, because they rest upon the same closure models, and are more accurate as applied to the 2D problem. Of the two anisotropic methods, it is worth giving preference to ACM‐1, because it is easier to implement.

Acknowledgments

The author sincerely thanks Goncharov E, Kolobyanin V, and Toporova O for the joint work on the closure models, Krayukhin A for problem 3 he proposed, and Kamm J and Shashkov M for consultations in the course of the work on the closure models and the setup of the 1D test problems.