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

αξ ¼ βξ �

where Δe<sup>ξ</sup> is the energy increment for the constituent material, and Δe, for the whole cell.

τg nþ1=2 ξ ρn ξ

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

<sup>∇</sup> � <sup>u</sup><sup>n</sup>þ1=<sup>2</sup>

ρ

Using Eq. (23) and the given ways of finding ∇ � u<sup>ξ</sup> from Eq. (31) one can obtain the values of

sents the change in the materials' internal energy as a result of their pressure relaxation. If no

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.

where the quantity λξ is determined by the chosen model of distributing the total divergence

<sup>ξ</sup> that represent additional changes in the materials' internal energy to meet the energy

where <sup>ρ</sup> <sup>¼</sup> <sup>X</sup>βξρξ:

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

Δe<sup>ξ</sup> ¼ −

We insert their expression (30) into Eq. (29) for Δe<sup>ξ</sup> and, using Eq. (28), obtain

<sup>ξ</sup> <sup>¼</sup> gnþ1=2, it follows from Eq. (31) that <sup>X</sup>αξΔ<sup>e</sup>

=

of the mixed cell to the constituent materials from the relation

τ <sup>X</sup>βξ g nþ1=2 <sup>ξ</sup> <sup>∇</sup> � <sup>u</sup><sup>n</sup>þ1=<sup>2</sup> ξ

has the following form:

76 Lagrangian Mechanics

balance requirement.

closing conditions.

nþ1=2

pressure relaxation is used, i.e., Δe

Let us write the expression for ψξ as

Given that g

Δe = −

<sup>τ</sup>gnþ1=<sup>2</sup><sup>∇</sup> � <sup>u</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>ρ</sup> <sup>¼</sup> <sup>−</sup>

ρξ ρ ,

<sup>Δ</sup><sup>e</sup> <sup>¼</sup> <sup>X</sup>αξΔeξ, (29)

<sup>þ</sup>Xαξ

=

<sup>ξ</sup> ¼ 0, then the definitions of ψξ follow directly from the

ψξ ¼ βξλξ, (32)

∇ � u<sup>ξ</sup> ¼ λξ∇ � u: (33)

<sup>ξ</sup> : (30)

Δe =

<sup>ξ</sup>: (31)

<sup>ξ</sup> ¼ 0. Thus, this term repre-

(28)

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} \tag{34}$$

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{split} \Delta E &= (e^{n+1/2} - e^n)M = -\frac{\tau (p^{n+1/2} + q^n)}{\rho^n} \nabla \cdot \mathbf{u}^{n+1/2} M \\ &= -\tau (p^{n+1/2} + q^n) \nabla \cdot \mathbf{u}^{n+1/2} V^n = (p^{n+1/2} + q^n) \mu \Delta V, \end{split} \tag{35}$$

where <sup>μ</sup> <sup>¼</sup> <sup>−</sup><sup>∇</sup> � <sup>u</sup><sup>n</sup>þ1=<sup>2</sup>V<sup>n</sup>τ=Δ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{aligned} \Delta E\_{\xi} &= (p^{n+1/2} + q^n) \mu \Delta V\_{\xi}, \\ \Delta V &= \sum\_{\xi} \Delta V\_{\xi}, \end{aligned} \tag{36}$$

$$P\left(\frac{M\_{\xi}}{V\_{\xi}^n + \Delta V\_{\xi}}, \frac{E\_{\xi} + \Delta E\_{\xi}}{M\_{\xi}}\right) = p^{n+1/2}.$$

The system (Eq. (36)) contains 2N + 1 equation in 2N + 1 unknown ΔVξ, ΔEξ, and pn+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}^{\imath+1/2} = e\_{\xi} + \frac{\Delta E\_{\xi}}{M\_{\xi}}, \quad \beta\_{\xi}^{\imath+1/2} = \frac{V\_{\xi} + \Delta V\_{\xi}}{V + \Delta V}. \tag{37}$$

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 <sup>r</sup><sup>n</sup>þ1=<sup>2</sup> <sup>¼</sup> <sup>r</sup><sup>n</sup> <sup>þ</sup> <sup>τ</sup>=<sup>2</sup> � <sup>u</sup><sup>n</sup>, and in Eq. (13), χ ¼ 1:0. Then,

$$p\_{\xi}^{n+1/2} = p\_{\xi}^{n} \cdot p\_{\xi}^{n} \cdot (c\_{\xi}^{n})^{2} \cdot \pi \cdot \nabla \cdot \mathbf{u}\_{\xi}^{n}. \tag{38}$$

This method assumes that the pressures of the constituent materials at the half-timestep should equal the average pressure pnþ1=2, which requires

$$p^{n+1/2} = p^{n+1/2}\_{\\\\\xi} + R^{n+1/2}\_{\\\xi}, \\
\tag{39}$$

where R<sup>n</sup>þ1=<sup>2</sup> <sup>ξ</sup> is different for different materials. The following equation is used to determine this quantity:

$$R^{n+1/2}\_{\xi} = -\rho^n\_{\xi} \cdot c^n\_{\xi} \cdot h^n \cdot \nabla \cdot \mathbf{u}^n\_{\xi},\tag{40}$$

where hn is the characteristic mesh spacing.

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

$$\mathbf{p}^{n+1/2} = \mathbf{p}\_{\boldsymbol{\xi}}^{n} - \tilde{\mathbf{B}}\_{\boldsymbol{\xi}}^{n} \nabla \cdot \mathbf{u}\_{\boldsymbol{\xi}}^{n},\tag{41}$$

where B~ <sup>ξ</sup> n <sup>≡</sup>ðcn ξÞ 2 ρn <sup>ξ</sup> � <sup>τ</sup> <sup>1</sup> <sup>þ</sup> hn <sup>ð</sup>τc<sup>n</sup> ξÞ h i: Eq. (41) together with the requirement (Eq. (25)) constitutes a closed set of algebraic equations with the unknowns <sup>∇</sup> � <sup>u</sup><sup>n</sup> <sup>ξ</sup> and p<sup>n</sup>þ1=2. As mentioned in Section 3.1.1, the system can be solved; for ξ = 1, 2, the solution is the following:

$$\nabla \cdot \mathbf{u}\_{\boldsymbol{\xi}}^{n} = \frac{(p\_{\boldsymbol{\xi}}^{n} - \overline{p}^{n}) + \overline{B}^{n} \cdot \nabla \cdot \mathbf{u}^{n}}{\tilde{B}^{n} \boldsymbol{\xi}},$$

$$p^{n+1/2} = \overline{p}^{n} - \overline{B}^{n} \nabla \cdot \mathbf{u}^{n},\tag{42}$$

where the barred terms denote average values of the quantities

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

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

$$
\Delta \beta\_{\boldsymbol{\xi}}^{n+1/2} = \beta\_{\boldsymbol{\xi}}^{n} \left( \frac{p\_{\boldsymbol{\xi}}^{n} - \overline{p}^{n}}{\check{B}\_{\boldsymbol{\xi}}^{n}} \right) + \beta\_{\boldsymbol{\xi}}^{n} \left[ \tau \left( \frac{\overline{B}^{n}}{\check{B}\_{\boldsymbol{\xi}}^{n}} - 1 \right) \right] \nabla \cdot \mathbf{u}^{n}. \tag{44}$$

To update the quantities at t <sup>n</sup>+1, another assumption is made that increments of the quantities after a full timestep are twice the half timestep increments:

$$\begin{aligned} \Delta \theta\_{\xi}^{\mathfrak{n}+1} &= 2\Delta \theta\_{\xi}^{\mathfrak{n}+1/2} \\ \text{and} \\ \theta\_{\xi}^{\mathfrak{n}+1} &= \theta\_{\xi}^{\mathfrak{n}} + 2\Delta \theta\_{\xi}^{\mathfrak{n}+1} . \end{aligned} \tag{45}$$

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

$$
\sigma\_{\xi}^{n+1} = \sigma\_{\xi}^{n} - p^{n+1/2} \frac{\Delta V\_{\xi}^{n+1}}{M\_{\xi}},\tag{46}
$$

where ΔV<sup>n</sup>þ<sup>1</sup> <sup>ξ</sup> <sup>¼</sup> Δβ<sup>n</sup>þ1=<sup>2</sup> <sup>ξ</sup> � <sup>V</sup><sup>n</sup>þ<sup>1</sup> .

#### 3.1.3. Delov's method

<sup>Δ</sup>E<sup>ξ</sup> ¼ ðpnþ1=<sup>2</sup> <sup>þ</sup> <sup>q</sup><sup>n</sup>ÞμΔVξ, <sup>Δ</sup><sup>V</sup> <sup>¼</sup> <sup>X</sup> ξ

> E<sup>ξ</sup> þ ΔE<sup>ξ</sup> M<sup>ξ</sup>

The system (Eq. (36)) contains 2N + 1 equation in 2N + 1 unknown ΔVξ, ΔEξ, and pn+1/2 and can be solved iteratively. Solving the system of equations gives updated values of specific

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

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 <sup>r</sup><sup>n</sup>þ1=<sup>2</sup> <sup>¼</sup> <sup>r</sup><sup>n</sup> <sup>þ</sup> <sup>τ</sup>=<sup>2</sup> � <sup>u</sup><sup>n</sup>,

This method assumes that the pressures of the constituent materials at the half-timestep

nþ1=2 <sup>ξ</sup> <sup>þ</sup> Rnþ1=<sup>2</sup>

<sup>ξ</sup> � c n

> <sup>ξ</sup>−B<sup>~</sup> <sup>ξ</sup> n <sup>∇</sup> � <sup>u</sup><sup>n</sup>

pnþ1=<sup>2</sup> <sup>¼</sup> pn

<sup>ξ</sup> is different for different materials. The following equation is used to determine

<sup>ξ</sup> � hn � <sup>∇</sup> � <sup>u</sup><sup>n</sup>

ξ−ρ<sup>n</sup> <sup>ξ</sup> � ðc n ξÞ

pnþ1=<sup>2</sup> <sup>¼</sup> <sup>p</sup>

R<sup>n</sup>þ1=<sup>2</sup> <sup>ξ</sup> <sup>¼</sup> <sup>−</sup>ρ<sup>n</sup>

,

ΔE<sup>ξ</sup> M<sup>ξ</sup> , β nþ1=2

!

<sup>P</sup> <sup>M</sup><sup>ξ</sup> Vn <sup>ξ</sup> þ ΔV<sup>ξ</sup>

example, in problems associated with energy release at every timestep.

p nþ1=2 <sup>ξ</sup> <sup>¼</sup> pn

should equal the average pressure pnþ1=2, which requires

where hn is the characteristic mesh spacing.

<sup>ξ</sup> � <sup>τ</sup> <sup>1</sup> <sup>þ</sup> hn <sup>ð</sup>τc<sup>n</sup> ξÞ h i

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

:

energies and volume fractions of the materials:

3.1.2. Tipton's method

78 Lagrangian Mechanics

where R<sup>n</sup>þ1=<sup>2</sup>

this quantity:

where B~ <sup>ξ</sup> n <sup>≡</sup>ðcn ξÞ 2 ρn

and in Eq. (13), χ ¼ 1:0. Then,

e nþ1=2 <sup>ξ</sup> ¼ e<sup>ξ</sup> þ ΔVξ,

<sup>¼</sup> pnþ1=<sup>2</sup>:

<sup>ξ</sup> <sup>¼</sup> <sup>V</sup><sup>ξ</sup> <sup>þ</sup> <sup>Δ</sup>V<sup>ξ</sup>

<sup>2</sup> � <sup>τ</sup> � <sup>∇</sup> � <sup>u</sup><sup>n</sup>

<sup>V</sup> <sup>þ</sup> <sup>Δ</sup><sup>V</sup> : (37)

<sup>ξ</sup>: (38)

<sup>ξ</sup> , (39)

<sup>ξ</sup>, (40)

<sup>ξ</sup>, (41)

(36)

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{aligned} p &= p\_1^n, \quad \iota = \iota\_{i-1}^{n+1/2} \quad \text{for} \quad m < m\_A, \\ p &= p\_2^n, \quad \iota = \iota\_i^{n+1/2} \quad \text{for} \quad m > m\_A, \end{aligned} \tag{47}$$

where m is the mass variable (mA 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{aligned} \boldsymbol{u}\_A^{n+1/2} + \frac{\boldsymbol{p}\_A^{n+1/2}}{\left(\rho c\right)\_1^n} &= \boldsymbol{u}\_{i-1}^{n+1/2} + \frac{\boldsymbol{p}\_1^n}{\left(\rho c\right)\_1^n}, \\ \boldsymbol{u}\_A^{n+1/2} - \frac{\boldsymbol{p}\_A^{n+1/2}}{\left(\rho c\right)\_2^n} &= \boldsymbol{u}\_i^{n+1/2} - \frac{\boldsymbol{p}\_2^n}{\left(\rho c\right)\_2^n}. \end{aligned} \tag{48}$$

The solution to this set will be the following expression for the velocity unþ1=<sup>2</sup> <sup>A</sup> :

$$
\mu\_A^{n+1/2} = \frac{u\_{i-1}^{n+1/2} \cdot (\rho c)\_1^n + u\_i^{n+1/2} \cdot (\rho c)\_2^n + p\_1^n \neg p\_2^n}{\left(\rho c\right)\_1^n + \left(\rho c\right)\_2^n}.
\tag{49}
$$

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 \cdot \frac{u\_i^{n+1/2} - u\_{i-1}^{n+1/2}}{M},\tag{50}$$

where M ¼ ρ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}{\begin{pmatrix} 1\\ \rho \end{pmatrix}^{n+1/2}} - \frac{1}{\rho^n \begin{pmatrix} 1\\ 1\\ 2 \end{pmatrix}} = \frac{\tau}{M \begin{Bmatrix} 2\\ 1 \end{Bmatrix}} \cdot \begin{Bmatrix} \begin{pmatrix} \rho c \end{Bmatrix} \begin{pmatrix} 2\\ 1\\ \rho c \end{pmatrix} \\ \begin{Bmatrix} \rho c \end{Bmatrix} + \begin{Bmatrix} \begin{pmatrix} \rho c \end{Bmatrix} \begin{Bmatrix} 2\\ 1 \end{Bmatrix} \end{Bmatrix} + \begin{Bmatrix} \begin{Bmatrix} \mathbf{u}^{n+1/2}\_{1} - \mathbf{u}^{n+1/2}\_{i-1} \end{Bmatrix} + \begin{Bmatrix} p^n\_1 - p^n\_2\\ \rho c \end{Bmatrix} \end{Bmatrix} . \tag{51}$$

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}}{a\_{\xi}^{n}} \cdot \left(\frac{1}{\rho^{n+1/2}} - \frac{1}{\rho^{n}}\right) + \frac{\omega \cdot \tau}{\rho\_{\xi}^{n} \cdot \rho\_{\xi}^{n} \cdot h^{n}} (p\_{\xi}^{n} - p\_{\geq \cdot}^{n}) \cdot \frac{1}{(\rho c)\_{\Sigma}^{n}},\tag{52}$$

where

$$p\_{\Sigma}^{n} = \frac{1}{N} \sum\_{\xi=1}^{N} p\_{\xi}^{n},\tag{53}$$

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

Closure Models for Lagrangian Gas Dynamics and Elastoplasticity Equations in Multimaterial Cells http://dx.doi.org/10.5772/66858 81

$$
\lambda\_{\boldsymbol{\xi}}^{n} = \frac{1}{(N-1)} \begin{pmatrix} \boldsymbol{1} - \frac{(\rho c)\_{\boldsymbol{\xi}}^{n}}{N} \\ \boldsymbol{1} - \frac{\boldsymbol{\zeta}}{N} (\rho c)\_{\boldsymbol{\zeta}}^{n} \end{pmatrix}, \tag{55}
$$

where h<sup>n</sup> 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:

u<sup>n</sup>þ1=<sup>2</sup> <sup>A</sup> <sup>þ</sup> <sup>p</sup>

> unþ1=<sup>2</sup> <sup>A</sup> − p nþ1=2 A ðρcÞ n 2

unþ1=<sup>2</sup>

where M ¼ ρh is the linear cell mass.

80 Lagrangian Mechanics

� � <sup>−</sup> <sup>1</sup> ρn 1 2

are also controlled by their pressures.

materials can be written as follows:

where

1 ρ<sup>n</sup>þ1=<sup>2</sup> ξ

¼ 1 ρn ξ þ λn ξ αn ξ

� � <sup>¼</sup> <sup>τ</sup> M <sup>2</sup> 1 � � �

<sup>A</sup> <sup>¼</sup> unþ1=<sup>2</sup>

The solution to this set will be the following expression for the velocity unþ1=<sup>2</sup>

<sup>i</sup>−<sup>1</sup> � ðρcÞ

gence with a corresponding expression, we obtain the following equation:

1 <sup>ρ</sup><sup>n</sup>þ1=<sup>2</sup> <sup>−</sup> <sup>1</sup> n <sup>1</sup> <sup>þ</sup> unþ1=<sup>2</sup>

ðρcÞ n <sup>1</sup> þ ðρcÞ n 2

Now, let us write the equation of continuity for the problem at issue. By replacing the diver-

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 substitut-

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

The extension of Eq. (51) in the multimaterial case without constraint of the number of

þ

ξ¼1 pn

ξ¼1 ðρcÞ n

ω � τ βn <sup>ξ</sup> � <sup>ρ</sup><sup>n</sup>

<sup>ξ</sup> � hn <sup>ð</sup>pn

<sup>ξ</sup> <sup>−</sup> pn

<sup>P</sup>Þ � <sup>1</sup> ðρсÞ n

<sup>ξ</sup>, (53)

<sup>ξ</sup>, (54)

unþ1=<sup>2</sup> <sup>i</sup> <sup>−</sup>unþ1=<sup>2</sup> i−1

� unþ1=<sup>2</sup>

<sup>i</sup> <sup>−</sup>unþ1=<sup>2</sup> i−1 � �

<sup>ρ</sup><sup>n</sup> <sup>¼</sup> <sup>τ</sup> �

ing the expressions for velocity (49) into the equation of continuity (50), we obtain

ðρcÞ n 2 1 � �

ðρcÞ n <sup>1</sup> þ ðρcÞ n 2

� <sup>1</sup> <sup>ρ</sup><sup>n</sup>þ1=<sup>2</sup> <sup>−</sup> <sup>1</sup> ρn

� �

pn <sup>P</sup> <sup>¼</sup> <sup>1</sup> N X N

ðρcÞ n <sup>P</sup> <sup>¼</sup> <sup>1</sup> N X N

nþ1=2 A ðρcÞ n 1

<sup>¼</sup> <sup>u</sup><sup>n</sup>þ1=<sup>2</sup>

<sup>¼</sup> unþ1=<sup>2</sup> <sup>i</sup> <sup>−</sup> pn 2 ðρcÞ n 2 :

<sup>i</sup>−<sup>1</sup> <sup>þ</sup> pn

<sup>i</sup> � ðρcÞ

1 ðρcÞ n 1 ,

> n <sup>2</sup> <sup>þ</sup> pn <sup>1</sup> −pn 2

<sup>A</sup> :

<sup>M</sup> , (50)

<sup>þ</sup> pn

ðρcÞ n <sup>1</sup> þ ðρcÞ n 2

<sup>1</sup> − p<sup>n</sup> 2 3 7 7

<sup>P</sup> , (52)

5: (51)

: (49)

(48)

$$\nabla \cdot \mathbf{u}\_{\boldsymbol{\xi}}^{\boldsymbol{n}} = \frac{\lambda\_{\boldsymbol{\xi}}^{\boldsymbol{n}}}{\beta\_{\boldsymbol{\xi}}^{\boldsymbol{n}}} \cdot \nabla \cdot \mathbf{u}^{\boldsymbol{n}} + \boldsymbol{\omega} \cdot \frac{(p\_{\boldsymbol{\xi}}^{\boldsymbol{n}} - p\_{\boldsymbol{\Sigma}}^{\boldsymbol{n}})}{\beta\_{\boldsymbol{\xi}}^{\boldsymbol{n}} \cdot h^{\boldsymbol{n}} \cdot (\rho c)\_{\boldsymbol{\Sigma}}^{\boldsymbol{n}}}.\tag{56}$$

Now, let us consider the quantity λ<sup>n</sup> <sup>ξ</sup>. 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}} \cdot (V^{n+1} - V^{n}) + \omega \cdot \tau \cdot \frac{(p\_{\xi}^{n} - p\_{\Sigma}^{n})}{\left(\rho c\right)\_{\Sigma}^{n}}.\tag{57}$$

From this, using the requirement (Eq. (25)), we obtain that the equality X N ξ¼1 λn ξðβ<sup>n</sup> ξÞ <sup>−</sup><sup>1</sup> <sup>¼</sup> 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 λ<sup>n</sup> ξðβ<sup>n</sup> ξÞ <sup>−</sup><sup>1</sup> can be taken as a function ψξ (ψξ <sup>¼</sup> <sup>λ</sup><sup>n</sup> ξðβ<sup>n</sup> ξÞ −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}. \tag{58}
$$

If the materials have different pressures, when we use Eq. (58), the right-hand member of the energy equation should be corrected by Δe′ <sup>ξ</sup>, for example, in the form of

$$
\Delta \dot{\boldsymbol{e}}\_{\xi}^{\prime} = \frac{\boldsymbol{\omega} \cdot \boldsymbol{\tau}}{\alpha\_{\xi} \rho^{n} \cdot \boldsymbol{h}^{n}} \cdot \frac{\boldsymbol{\psi}\_{\xi}}{\left(\rho \boldsymbol{c}\right)^{n}\_{\boldsymbol{\Sigma}}} \cdot \sum\_{\xi} p\_{\xi}^{n} \cdot \left(p\_{\xi}^{n} - p\_{\boldsymbol{\Sigma}}^{n}\right). \tag{59}
$$

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≫βξ≫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 \cdot \mathbf{u}\_{\xi}^{\boldsymbol{\mu}} = \nabla \cdot \mathbf{u}^{\boldsymbol{\mu}}.\tag{60}$$

Naturally, it is also assumed that <sup>∇</sup> � <sup>u</sup><sup>n</sup>þ<sup>1</sup> <sup>ξ</sup> <sup>¼</sup> <sup>∇</sup> � <sup>u</sup><sup>n</sup>þ1.

Thismethod, 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},
\tag{61}
$$

which ensures the fulfillment of the requirement (31) at Δe′ <sup>ξ</sup> ¼ 0. As a result, formula (23) transforms to

$$p = \sum \beta\_{\xi} p\_{\underline{\xi}}.\tag{62}$$

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<sup>0</sup> (see Figure 2).

Figure 2. Graphic illustration of Formula (62).

Let us assign p<sup>ξ</sup> 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 \cdot \frac{p\_2 - p\_1}{V\_o} \left( V - \frac{V\_1}{2} \right) + p\_1. \tag{63}$$

If we insert the value of V ¼ V0=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 V0=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 \cdot \mathbf{u}\_{\xi} = \rho\_{\zeta} c\_{\zeta}^{2} \nabla \cdot \mathbf{u}\_{\zeta}. \tag{64}
$$

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

$$\frac{\partial p}{\partial t} = \frac{\partial p}{\partial \rho} \cdot \frac{\partial \rho}{\partial t} = c^2 \rho \nabla \cdot \mathbf{u}.\tag{65}$$

Going over in Eq. (65) to the materials and claiming that <sup>∂</sup>p<sup>ξ</sup> <sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>∂</sup>p<sup>ζ</sup> <sup>∂</sup><sup>t</sup> , we obtain the condition (64). The set of algebraic equations (64), Eq. (25) is closed, and its solution is given by

$$
\nabla \cdot \mathbf{u}\_{\xi}^{n} = \lambda\_{\xi} \nabla \cdot \mathbf{u}^{n}, \tag{66}
$$

where

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.

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

<sup>ξ</sup> <sup>¼</sup> <sup>∇</sup> � <sup>u</sup><sup>n</sup>: (60)

ψξ ¼ βξ, (61)

<sup>ξ</sup> ¼ 0. As a result, formula (23)

: (62)

þ p1: (63)

<sup>∇</sup> � <sup>u</sup><sup>n</sup>

the 2D case, this method is generalized on the assumption of volume-linear interpolation.

<sup>ξ</sup> <sup>¼</sup> <sup>∇</sup> � <sup>u</sup><sup>n</sup>þ1.

Thismethod, 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

<sup>p</sup> <sup>¼</sup> <sup>X</sup>βξp<sup>ξ</sup>

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

Let us assign p<sup>ξ</sup> to the centers of their volumes. If the volume is used as a variable, then the

V− V1 2 � �

p2− p<sup>1</sup> Vo

3.1.4. Method 5 based on the equal compressibility model

been considered in [38, 39] and is formulated as follows:

which ensures the fulfillment of the requirement (31) at Δe′

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

linear interpolation of pressure over the cell volume can be written as

pðVÞ ¼ 2 �

Naturally, it is also assumed that <sup>∇</sup> � <sup>u</sup><sup>n</sup>þ<sup>1</sup>

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

transforms to

82 Lagrangian Mechanics

volume V<sup>0</sup> (see Figure 2).

Figure 2. Graphic illustration of Formula (62).

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

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 pn <sup>ξ</sup> ≠ pn <sup>ς</sup> . For an ideal gas, the following estimate is true:

$$\nabla \cdot \mathbf{u}\_{\xi} \sim \left(\rho\_{\xi} c\_{\xi}^{2}\right)^{-1} \sim \left(\boldsymbol{\gamma}\_{\xi} \boldsymbol{p}\_{\xi}\right)^{-1}.\tag{68}$$

It follows from Eq. (68) that at close values of γξ, the values of ∇ � u<sup>ξ</sup> 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 \cdot \mathbf{u}\_{\xi} \approx -\frac{\Delta p\_{\xi}}{p\_{\xi}} \cdot \frac{p\_{\xi}}{\pi p\_{\xi} \mathbf{d} p\_{\xi} / \mathbf{d} p\_{\xi}},\tag{69}$$

require that the condition Δpξ=p<sup>ξ</sup> ¼ Δpζ=p<sup>ζ</sup> is fulfilled, and obtain formula (70):

$$
\lambda\_{\xi} = -\frac{p\_{\varepsilon}}{p\_{\varepsilon}c\_{\varepsilon}^{2}} \cdot \frac{1}{\sum \beta\_{\mathbf{k}} p\_{\mathbf{k}} / (\rho\_{\mathbf{k}} c\_{\mathbf{k}}^{2})}.\tag{70}
$$

This value of λξ will also be used in Eq. (32) to determine ψξ, which meets the requirement (31) at Δe′ <sup>ξ</sup> ¼ 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 Δρ=ρ ¼ Δu=c, then for ∇ � u<sup>ξ</sup> one can write the following expression:

$$\nabla \cdot \mathbf{u}\_{\xi} \approx -\frac{\Delta \rho\_{\xi}}{\rho\_{\xi} \pi} = -\frac{\Delta u\_{\xi}}{c\_{\xi} \pi} \,. \tag{71}$$

Considering that the materials' velocity increments Δu<sup>ξ</sup> in the mixed cell are assumed to be equal, for ∇ � u<sup>ξ</sup> we obtain the following relation (72):

$$
\mathbf{c}\_{\xi} \nabla \cdot \mathbf{u}\_{\xi} = \mathbf{c}\_{\zeta} \nabla \cdot \mathbf{u}\_{\zeta}. \tag{72}
$$

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

Closure Models for Lagrangian Gas Dynamics and Elastoplasticity Equations in Multimaterial Cells http://dx.doi.org/10.5772/66858 85

$$\nabla \cdot \mathbf{u}\_{\xi}^{n} = \lambda\_{\xi} \nabla \cdot \mathbf{u}^{n},\tag{73}$$

where

One of them is the following. Let pn

84 Lagrangian Mechanics

at Δe′

<sup>ξ</sup> ¼ 0.

<sup>ξ</sup> ≠ pn

should be equal. Let us write the modified equation (Eq. (65)) as

∇ � u<sup>ξ</sup> ∼ ðρξc

<sup>∇</sup> � <sup>u</sup>ξ≈ − <sup>Δ</sup>p<sup>ξ</sup>

require that the condition Δpξ=p<sup>ξ</sup> ¼ Δpζ=p<sup>ζ</sup> is fulfilled, and obtain formula (70):

λξ <sup>¼</sup> <sup>−</sup> <sup>p</sup><sup>ξ</sup> pξ c2 ξ �X

3.1.6. Method 7 based on the equal velocity increments model

equal, for ∇ � u<sup>ξ</sup> we obtain the following relation (72):

pξ

� <sup>p</sup><sup>ξ</sup> τpξdpξ=dp<sup>ξ</sup>

> 1 βkpk=ðρkc

This value of λξ will also be used in Eq. (32) to determine ψξ, which meets the requirement (31)

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 Δρ=ρ ¼ Δu=c, then for ∇ � u<sup>ξ</sup> one can write the following expression:

<sup>∇</sup> � <sup>u</sup>ξ≈ − Δρξ

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

ρξτ <sup>¼</sup> <sup>−</sup>

Considering that the materials' velocity increments Δu<sup>ξ</sup> in the mixed cell are assumed to be

Δu<sup>ξ</sup>

2 kÞ

2 ξÞ

It follows from Eq. (68) that at close values of γξ, the values of ∇ � u<sup>ξ</sup> 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

�<sup>1</sup> <sup>∼</sup> <sup>ð</sup>γξpξ<sup>Þ</sup>

<sup>ς</sup> . For an ideal gas, the following estimate is true:

: (68)

, (69)

: (70)

<sup>c</sup>ξτ : (71)

cξ∇ � u<sup>ξ</sup> ¼ cζ∇ � uζ: (72)

�1

$$
\lambda\_{\xi} = \left( c\_{\xi} \sum\_{k} \frac{\beta\_{k}}{c\_{k}} \right)^{-1} = \frac{c}{c\_{\xi}},
\tag{74}
$$

with

$$c = \left(\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 Δe′ <sup>ξ</sup> ¼ 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 <sup>n</sup>+1 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 twostage 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 ∇ � u<sup>ξ</sup> at the timestep in two stages:

$$
\nabla \cdot \mathbf{u}\_{\xi} = \nabla \cdot \mathbf{u}\_{\xi 1} + \nabla \cdot \mathbf{u}\_{\xi 2}. \tag{75}
$$

In Eq. (75), ∇ � uξ<sup>1</sup> 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 ∇ � u and the total internal energy, i.e., Δ E<sup>2</sup> = 0, should remain unchanged at that stage. Pressure relaxation is implemented by calculating additional divergences of the materials ∇ � uξ<sup>2</sup> by the formulas

$$\nabla \cdot \mathbf{u}\_{\xi2} = -\frac{\Delta p\_{\xi}}{\pi \rho\_{\xi} c\_{\xi}^{2}},\tag{76}$$

where

$$
\Delta p\_{\xi} = A \frac{c\tau}{h} (p - p\_{\xi}),
\tag{77}
$$

where p is the average pressure. Expression (76) was derived using a known relation in the adiabatic approximation, <sup>Δ</sup><sup>p</sup> <sup>¼</sup> <sup>−</sup>ρc<sup>2</sup>τ<sup>∇</sup> � <sup>u</sup>. The factor <sup>c</sup>τ=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, <sup>X</sup>βξ<sup>∇</sup> � <sup>u</sup>ξ<sup>2</sup> <sup>¼</sup> 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 γ<sup>1</sup> = γ<sup>2</sup> = γ, but disparate pressures and volume fractions. Let p<sup>1</sup> = 1000, β<sup>1</sup> = 0.9 and p<sup>2</sup> = 1, β<sup>2</sup> = 0.1. Using Eqs. (23) and (25) combined with the relation ρξc<sup>2</sup> <sup>ξ</sup> ¼ γξpξ, 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

$$
\Delta E = -P\_+ \Delta V\_+ - P\_- \Delta V\_-,\tag{78}
$$

where <sup>P</sup><sup>þ</sup> <sup>¼</sup> <sup>X</sup>βξp<sup>ξ</sup> if <sup>p</sup><sup>ξ</sup> <sup>&</sup>gt; <sup>p</sup>; <sup>P</sup><sup>−</sup> <sup>¼</sup> <sup>X</sup>βξp<sup>ξ</sup> if <sup>p</sup><sup>ξ</sup> <sup>≤</sup> <sup>p</sup> and ΔV<sup>þ</sup> and ΔV<sup>−</sup> are the volume changes of these materials.

With pressure equilibration, the materials with P<sup>þ</sup> expand, so ΔV<sup>þ</sup> > 0 and ΔV<sup>−</sup> < 0. As far as it follows from the volume conservation condition that ΔV<sup>þ</sup> ¼ jΔV−j, and by definition Pþ>P−, Δ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 Δe′ <sup>ξ</sup> in accordance with the equation

Closure Models for Lagrangian Gas Dynamics and Elastoplasticity Equations in Multimaterial Cells http://dx.doi.org/10.5772/66858 87

$$-\frac{\tau}{\rho^{\mathfrak{n}}}\sum\_{\xi} \beta\_{\xi} (p\_{\xi}^{\mathfrak{n}+1/2} + q\_{\xi}^{\mathfrak{n}}) \nabla \cdot \mathbf{u}\_{\xi \mathbf{2}}^{\mathfrak{n}+1/2} + \sum\_{\xi} a\_{\xi} \Delta e\_{\xi}^{\prime} = 0,\tag{79}$$

The quantity <sup>∇</sup> � unþ<sup>1</sup> <sup>ξ</sup><sup>2</sup> present in this expression is calculated by the formula

$$\nabla \cdot \mathbf{u}\_{\xi2}^{n+1} = \nabla \cdot \mathbf{u}\_{\xi2}^n (1 - Ac\tau/h).$$

It remains to decide how to distribute the dissipated kinetic energy among the materials (formula (79) defines only the total dissipated energy <sup>Δ</sup><sup>E</sup> <sup>¼</sup> <sup>X</sup> ξ αξΔe ′ <sup>ξ</sup>). In [20], it is assumed

that Δe′ <sup>ξ</sup> <sup>¼</sup> <sup>Δ</sup>e′ . In this case, for all the materials in the mixed cell we obtain from Eq. (79) that

$$
\Delta \dot{\boldsymbol{e}}\_{\cdot \xi}^{\cdot} = \Delta \dot{\boldsymbol{e}}^{\cdot} = \frac{\pi}{\rho^{\text{v}}} \sum\_{\xi} \beta\_{\xi} (\boldsymbol{p}\_{\xi}^{\text{n}+1/2} + \boldsymbol{q}\_{\xi}^{\text{v}}) \nabla \cdot \mathbf{u}\_{\xi \mathbf{2}}^{\text{n}+1/2}. \tag{80}
$$

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)

<sup>∇</sup> � <sup>u</sup>ξ<sup>2</sup> <sup>¼</sup> <sup>−</sup> <sup>Δ</sup>p<sup>ξ</sup>

<sup>Δ</sup>p<sup>ξ</sup> <sup>¼</sup> <sup>A</sup> <sup>c</sup><sup>τ</sup>

where p is the average pressure. Expression (76) was derived using a known relation in the adiabatic approximation, <sup>Δ</sup><sup>p</sup> <sup>¼</sup> <sup>−</sup>ρc<sup>2</sup>τ<sup>∇</sup> � <sup>u</sup>. The factor <sup>c</sup>τ=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

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, <sup>X</sup>βξ<sup>∇</sup> � <sup>u</sup>ξ<sup>2</sup> <sup>¼</sup> 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 γ<sup>1</sup> = γ<sup>2</sup> = γ, but disparate pressures and volume fractions. Let p<sup>1</sup> = 1000, β<sup>1</sup> = 0.9 and p<sup>2</sup> = 1, β<sup>2</sup> = 0.1. Using Eqs. (23) and (25) combined with

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

This method results in the internal energy exchange between the materials. Indeed, let us

With pressure equilibration, the materials with P<sup>þ</sup> expand, so ΔV<sup>þ</sup> > 0 and ΔV<sup>−</sup> < 0. As far as it follows from the volume conservation condition that ΔV<sup>þ</sup> ¼ jΔV−j, and by definition Pþ>P−, Δ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

<sup>ξ</sup> in accordance with the equation

represent the total change in the materials' internal energies as

<sup>P</sup><sup>þ</sup> <sup>¼</sup> <sup>X</sup>βξp<sup>ξ</sup> if <sup>p</sup><sup>ξ</sup> <sup>&</sup>gt; <sup>p</sup>;

<sup>P</sup><sup>−</sup> <sup>¼</sup> <sup>X</sup>βξp<sup>ξ</sup> if <sup>p</sup><sup>ξ</sup> <sup>≤</sup> <sup>p</sup>

internal energy increments Δe′

<sup>ξ</sup> ¼ γξpξ, one can easily calculate that p = 10. Thus, the cell-average pressure

ΔE ¼ −PþΔVþ−P−ΔV−, (78)

and ΔV<sup>þ</sup> and ΔV<sup>−</sup> are the volume changes of these materials.

where

86 Lagrangian Mechanics

single timestep.

the relation ρξc<sup>2</sup>

mathematical basis.

where

τρξc<sup>2</sup> ξ , (76)

<sup>h</sup> <sup>ð</sup>p−pξÞ, (77)

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},\tag{81}
$$

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. \tag{82}
$$

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<sup>1</sup> > P2. 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,\tag{83}
$$

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

Thus,

$$
\begin{array}{l}
\Delta V\_{21} = \Delta V\_{2}, \\
\Delta V\_{22} = -\Delta V\_{2}.
\end{array}
\tag{84}
$$

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

$$V^{n+1} > V^{n+1}\_{\\\\\xi} > 0,\tag{85}$$

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

After the first stage, P<sup>1</sup> > P2. 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

ρ1c<sup>1</sup> þ ρ2c<sup>2</sup>

where all the unindexed quantities are taken after the first stage, and SAB is the area of the

ΔV<sup>21</sup> ¼ ΔV2,

V<sup>n</sup>þ<sup>1</sup> > V<sup>n</sup>þ<sup>1</sup>

The updated volumes of the materials calculated by formula (81) accounting for the volume

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

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

SABτ<sup>n</sup>, (83)

<sup>Δ</sup>V<sup>22</sup> <sup>¼</sup> <sup>−</sup>ΔV2: (84)

<sup>ξ</sup> > 0, (85)

<sup>Δ</sup>V<sup>2</sup> <sup>¼</sup> <sup>p</sup>1<sup>−</sup> <sup>p</sup><sup>2</sup>

based on the solution of the acoustic Riemann problem

Figure 3. Contact location reconstructed after the first stage.

increments must satisfy the following inequalities:

boundary between the materials.

3.2. Anisotropic closure models

3.2.1. Model fundamentals

Thus,

88 Lagrangian Mechanics

case.

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 ∇ � u<sup>ξ</sup> at a timestep

$$
\nabla \cdot \mathbf{u}\_{\xi} = \nabla \cdot \mathbf{u}\_{\xi 1} + \nabla \cdot \mathbf{u}\_{\xi 2}. \tag{86}
$$

In Eq. (86), ∇ � uξ<sup>1</sup> is calculated at the first stage in accordance with closure model 1:

$$
\nabla \cdot \mathbf{u}\_{\xi1} = \nabla \cdot \mathbf{u}.\tag{87}
$$

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 \cdot \mathbf{u} = \nabla \cdot \mathbf{u}\_{\tau} + \nabla \cdot \mathbf{u}\_{n}, \tag{88}
$$

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

$$A = A\_0 \cdot \frac{\nabla \cdot \mathbf{u}\_n}{\nabla \cdot \mathbf{u}},\tag{89}$$

where A<sup>0</sup> is some constant.

Thus, the coefficient A is a variable in this case. If the velocity is normal to the interface (Figure 4a), ∇ � u<sup>n</sup> ¼ ∇ � u and A = A0; this is the case with the highest pressure relaxation. If the velocity is directed along the interface (Figure 4b), ∇ � u<sup>n</sup> ¼ 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<sup>0</sup> 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<sup>0</sup> = 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{aligned} \nabla \cdot \mathbf{u} &= \nabla \cdot \mathbf{u}\_n + \nabla \cdot \mathbf{u}\_\tau \\ \nabla \cdot \mathbf{u}\_\xi &= \nabla \cdot \mathbf{u}\_{\xi n} + \nabla \cdot \mathbf{u}\_{\xi \tau} \end{aligned} \tag{90}$$

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 \cdot \mathbf{u}\_{\xi\tau} = \nabla \cdot \mathbf{u}\_{\tau}.\tag{91}
$$

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{aligned} \nabla \cdot \mathbf{u}\_{\xi,\mathbf{u}} &= \lambda\_{\xi} \nabla \cdot \mathbf{u}\_{\eta} \\ \lambda\_{\xi} &= \left( c\_{\xi} \sum\_{k} \frac{\beta\_{k}}{c\_{k}} \right)^{-1} . \end{aligned} \tag{92}$$

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,<sup>∇</sup> � <sup>u</sup>′ <sup>ξ</sup><sup>n</sup>, to the divergence ∇ � uξ<sup>n</sup>:

Closure Models for Lagrangian Gas Dynamics and Elastoplasticity Equations in Multimaterial Cells http://dx.doi.org/10.5772/66858 91

$$
\widetilde{\nabla \cdot \mathbf{u}}\_{\xi n} = \nabla \cdot \mathbf{u}\_{\xi n} + \nabla \cdot \mathbf{u}\_{\xi n}^{'}.\tag{93}
$$

$$
\circ \widetilde{\nabla \cdot \mathbf{u}}\_{\xi n} = \nabla \cdot \mathbf{u}\_{\xi n} + \nabla \cdot \mathbf{u}\_{\xi n}^{'}.\tag{93}
$$

The ultimate formula for the distribution of ∇ � u among the materials will be

$$
\mathbf{\tilde{v}} \cdot \mathbf{u} = \mathbf{\tilde{u}} \cdot \mathbf{\tilde{u}}
$$

$$
\nabla \cdot \mathbf{u}\_{\xi} = \widehat{\nabla \cdot \mathbf{u}}\_{\xi \eta} + \nabla \cdot \mathbf{u}\_{\xi \tau}. \tag{94}
$$
