4. Artificial viscosity

the mutual orientation of the velocity and the interface. The total divergence is written as the

obtained by velocity decomposition to two components: along the interface (uτ) and normal to

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

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

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

> ∇ � u ¼ ∇ � u<sup>n</sup> þ ∇ � u<sup>τ</sup> ∇ � u<sup>ξ</sup> ¼ ∇ � uξ<sup>n</sup> þ ∇ � uξτ

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

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

∇ � uξ<sup>n</sup> ¼ λξ∇ � u<sup>n</sup>

X k

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

!<sup>−</sup><sup>1</sup>

βk ck

λξ ¼ c<sup>ξ</sup>

<sup>ξ</sup><sup>n</sup>, to the divergence ∇ � uξ<sup>n</sup>:

∇ � u<sup>n</sup>

A ¼ A<sup>0</sup> �

∇ � u ¼ ∇ � u<sup>τ</sup> þ ∇ � un, (88)

<sup>∇</sup> � <sup>u</sup> , (89)

: (90)

: (92)

∇ � uξτ ¼ ∇ � uτ: (91)

sum of two components:

90 Lagrangian Mechanics

it (un). We also assume that

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

constant A in Section 3.1.7.1, i.e., A<sup>0</sup> = 1.

formula (87).

location):

relation

follows from it that

contribution,<sup>∇</sup> � <sup>u</sup>′

3.2.3. The ACM-2 model

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{cases} q^{\boldsymbol{\eta}} = \mathbb{C}\_{1} \cdot \rho^{\boldsymbol{\eta}} (h^{\boldsymbol{\eta}} \nabla \cdot \mathbf{u}^{\boldsymbol{\eta}})^{2} + \mathbb{C}\_{0} \cdot \rho^{\boldsymbol{\eta}} \cdot \boldsymbol{c} \cdot h^{\boldsymbol{\eta}} \nabla \cdot \mathbf{u}^{\boldsymbol{\eta}} & \text{if } \nabla \cdot \mathbf{u}^{\boldsymbol{\eta}} < 0 \\ q^{\boldsymbol{\eta}} = 0 & \text{if } \nabla \cdot \mathbf{u}^{\boldsymbol{\eta}} \succeq 0 \end{cases} \tag{95}$$

In Eq. (95), C<sup>1</sup> = 1 and C<sup>0</sup> = 0.2 are the fixed coefficients. In addition, expression (95) contains the characteristic dimension <sup>h</sup> and the divergence <sup>∇</sup> � <sup>u</sup><sup>n</sup> 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., <sup>q</sup> <sup>¼</sup> <sup>X</sup>ψξqξ, ψξ <sup>¼</sup> βξλξ, 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\_{\vec{\xi}} = -\frac{q\_{\vec{\xi}}^n}{\rho\_{\vec{\xi}}^n} \cdot \lambda\_{\vec{\xi}}^n \cdot \tau \cdot (\nabla \cdot \mathbf{u}^{n+1/2}).\tag{96}
$$

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 \cdot \Delta \rho \approx -\rho \cdot c^2 \cdot \tau \cdot \nabla \cdot \mathbf{u}.\tag{97}
$$

considering that Δρ≈−ρ � ∇ � u � τ.

Using the EOS p = P(r,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)\_\varepsilon \Delta \rho + \left(\frac{\partial p}{\partial \varepsilon}\right)\_\rho \Delta \varepsilon. \tag{98}
$$

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

$$
\Delta e = -\frac{p}{\rho} \cdot \tau \cdot \nabla \cdot \mathbf{u}.\tag{99}
$$

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

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

$$
\left(\frac{\partial p}{\partial \rho}\right)\_\epsilon + \left(\frac{\partial p}{\partial \varepsilon}\right)\_\rho \cdot \frac{p}{\rho^2} = \varepsilon^2. \tag{100}
$$

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

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

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., <sup>q</sup> <sup>¼</sup> <sup>X</sup>ψξqξ, ψξ <sup>¼</sup> βξλξ, 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

In accordance with the difference scheme of EGAK, the viscosity-related change in the mate-

The viscosity-related change in the materials' pressure at a timestep can be obtained as follows.

� Δρ ≈ −ρ � c

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

∂p ∂e � �

ρ

e Δρ þ

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

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

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

Þ: (96)

<sup>2</sup> � <sup>τ</sup> � <sup>∇</sup> � <sup>u</sup>: (97)

Δe: (98)

<sup>ρ</sup> � <sup>τ</sup> � <sup>∇</sup> � <sup>u</sup>: (99)

redistribution of dissipated energy between the materials.

rials' specific internal energy over a timestep is given by the formula

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

<sup>Δ</sup><sup>p</sup> <sup>≈</sup> <sup>∂</sup><sup>p</sup> ∂ρ � �

qn ξ ρn ξ � λn

S

<sup>Δ</sup><sup>p</sup> <sup>≈</sup> <sup>∂</sup><sup>p</sup> ∂ρ � �

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

EGAK on methods 5–7 and ACM-1).

For adiabatic flows, it holds true that

considering that Δρ≈−ρ � ∇ � u � τ.

represented as

92 Lagrangian Mechanics

we obtain

$$
\Delta e = -\frac{p+q}{\rho} \cdot \boldsymbol{\tau} \cdot \nabla \cdot \mathbf{u}.\tag{101}
$$

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

$$
\Delta p^{\varpi -} \left[ \rho \cdot c^2 + \left( \frac{\partial p}{\partial c} \right)\_\rho \cdot \frac{q}{\rho} \right] \tau \cdot \nabla \cdot \mathbf{u}. \tag{102}
$$

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 \mathbf{e}\_{\xi}}\right)\_{\rho} \cdot \frac{q\_{\xi}^{n}}{\rho\_{\xi}^{n}} \cdot \lambda\_{\xi}^{n} \cdot \pi \cdot (\nabla \cdot \mathbf{u}\_{\xi}^{n+1/2}).\tag{103}
$$

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).

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.


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