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 \cdot \mathbf{T},\tag{1}$$

$$\frac{d\rho\_{\xi}}{dt} = -\rho\_{\xi}\nabla \cdot \mathbf{u}\_{\xi},\tag{2}$$

$$\frac{d\boldsymbol{\beta}\_{\boldsymbol{\xi}}}{dt} = \boldsymbol{\beta}\_{\boldsymbol{\xi}} (\boldsymbol{\nabla} \cdot \mathbf{u}\_{\boldsymbol{\xi}} - \boldsymbol{\nabla} \cdot \mathbf{u}),\tag{3}$$

$$\rho \frac{de\_{\xi}}{dt} = \frac{1}{\rho\_{\xi}} \text{Sp}(\mathbf{T}\_{\xi} \mathbf{D}\_{\xi}),\tag{4}$$

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

In this set of equations: u(ux, uy) 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 (βξ <sup>¼</sup> <sup>M</sup><sup>ξ</sup> <sup>V</sup> ), 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 ρξ <sup>¼</sup> <sup>M</sup><sup>ξ</sup> V<sup>ξ</sup> , where M<sup>ξ</sup> and V<sup>ξ</sup> are the mass and volume of the materials in the Lagrangian cell. We obtain dV<sup>ξ</sup> dt ¼ Vξ∇ � u<sup>ξ</sup> and then introduce an expression for V<sup>ξ</sup> in terms of volume fractions V<sup>ξ</sup> = βξ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 <sup>n</sup>+1. 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:

EGAK code [35] when modeling flows, for which one can assume that they are isotropic, but

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

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

∇ � T, (1)

SpðTξDξÞ, (4)

dt <sup>¼</sup> <sup>u</sup>: (5)

, where M<sup>ξ</sup> and V<sup>ξ</sup> are the mass and

dt ¼ Vξ∇ � u<sup>ξ</sup> and then introduce an

<sup>n</sup>+1. In the single-material

dt <sup>¼</sup> <sup>−</sup>ρξ<sup>∇</sup> � <sup>u</sup>ξ, (2)

dt <sup>¼</sup> βξð<sup>∇</sup> � <sup>u</sup>ξ−<sup>∇</sup> � <sup>u</sup>Þ, (3)

have an important advantage when modeling more complex flows.

2.1. Initial equations for multimaterial elastoplasticity

2. Finite difference approximation of elastoplasticity equations

du dt <sup>¼</sup> <sup>−</sup> <sup>1</sup> ρ

dρξ

dβξ

de<sup>ξ</sup> dt <sup>¼</sup> <sup>1</sup> ρξ

dr

In this set of equations: u(ux, uy) 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

Eq. (3) can be derived from the equation of continuity (2), in which the materials' density is

expression for V<sup>ξ</sup> in terms of volume fractions V<sup>ξ</sup> = βξ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

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.

<sup>V</sup> ), 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

V<sup>ξ</sup>

discussed in [36].

72 Lagrangian Mechanics

following:

(βξ <sup>¼</sup> <sup>M</sup><sup>ξ</sup>

vector, tensor, and deviator.

substituted with its expression in the form of ρξ <sup>¼</sup> <sup>M</sup><sup>ξ</sup>

volume of the materials in the Lagrangian cell. We obtain dV<sup>ξ</sup>

volume fractions for multimaterial matter should also be updated to t

$$\mathbf{T} = \begin{bmatrix} T\_{xx} & T\_{xy} & 0 \\ T\_{yx} & T\_{yy} & 0 \\ \mathbf{0} & \mathbf{0} & \mathbf{T}\_{\phi} \end{bmatrix}, \quad \mathbf{D} = \begin{bmatrix} d\_{xx} & d\_{xy} & 0 \\ d\_{yx} & d\_{yy} & 0 \\ \mathbf{0} & \mathbf{0} & d\_{\phi} \end{bmatrix}. \tag{6}$$

The stress tensor is represented as a sum of the spherical part (pressure p) and the deviator S(Sxx, Syy , Sxy , Sφ). Deviator components are defined by the relation Sij = Ti<sup>j</sup> - δijp.

For the materials, we define equations of state

$$p\_{\xi} = P\_{\xi}(\rho\_{\xi}, e\_{\xi}),\tag{7}$$

and equations to express the deviator S<sup>ξ</sup> as a function of the strain rate tensor D<sup>ξ</sup>

$$f\_{\vec{\xi}}(\mathbf{S}\_{\vec{\xi}}, \mathbf{D}\_{\vec{\xi}}) = 0. \tag{8}$$

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 (ρξ, βξ,eξ, pξ, 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<sup>ξ</sup> ! p<sup>ξ</sup> + qξ, respectively. Known quantities (basic variables) in Eqs. (1)–(5) in the 2D case include ux, uy, ρξ, βξ,eξ, pξ, Sξ, and the quantities p, S, q, qξ, ∇�uξ, D<sup>ξ</sup> 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 <sup>n</sup> and that we seek to update their values at time t <sup>n</sup>+1 = t <sup>n</sup> + τ, 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 + \boldsymbol{\tau} \cdot \mathbf{u}^n,\tag{9}$$

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

$$p\_{\xi}^{n} = P(\rho\_{\xi}^{n}, e\_{\xi}^{n}),\tag{11}$$

$$\nabla \cdot \mathbf{u}^n = \frac{(V^{n+1/2} - V^n)}{(\tau \cdot V^n)},\tag{12}$$

$$p\_{\xi}^{n+1/2} = p\_{\xi}^{n} - \chi \cdot \rho\_{\xi}^{n} \cdot \left(\boldsymbol{c}\_{\xi}^{n}\right)^{2} \cdot \boldsymbol{\pi} \cdot \overline{\nabla \cdot \mathbf{u}}\_{\xi}^{n}.\tag{13}$$

In Eq. (13), <sup>χ</sup> <sup>¼</sup> <sup>0</sup>:6 (this value was chosen in [17]), cn <sup>ξ</sup> is the speed of sound.

Full timestep

$$M\_{i,j} = 0.25 \cdot \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} \cdot \beta\_{\xi} \cdot V,\tag{14}$$

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

$$\Rightarrow \mathbf{u}\_{i,j}^{n+1} = \mathbf{u}\_{i,j}^{n} - \left(\frac{\boldsymbol{\pi}}{M\_{i,j}}\right) \cdot \left(\nabla \overline{\mathbf{g}}^{n+1/2} + \overline{\nabla \cdot \boldsymbol{S}}^{n}\right)\_{i,j},\tag{16}$$

$$\mathbf{u}^{n+1/2} = (\mathbf{u}^n + \mathbf{u}^{n+1})/2,\tag{17}$$

$$\mathbf{r}^{n+1} = \mathbf{r}^n + \boldsymbol{\tau} \cdot \mathbf{u}^{n+1},\tag{18}$$

$$V\_{\xi}^{n+1} = V\_{\xi}^{n} + \tau \cdot V\_{\xi}^{n+1} \overline{\nabla \cdot \mathbf{u}\_{\xi}}^{n+1},\tag{19}$$

$$
\rho\_{\xi}^{n+1} = \frac{M\_{\xi}}{V\_{\xi}^{n+1}},
\tag{20}
$$

$$
\beta\_{\xi}^{n+1} = \frac{V\_{\xi}}{V^{n+1}},
\tag{21}
$$

$$\boldsymbol{\varepsilon}\_{\boldsymbol{\xi}}^{n+1} = \boldsymbol{\varepsilon}\_{\boldsymbol{\xi}}^{n} - \frac{\boldsymbol{\pi}}{\rho\_{\boldsymbol{\xi}}^{n}} \cdot \left[ \boldsymbol{g}\_{\boldsymbol{\xi}}^{n+1/2} \cdot \overline{\boldsymbol{\nabla} \cdot \mathbf{u}\_{\boldsymbol{\xi}}^{n+1/2}} - \boldsymbol{S}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{x}}^{n} \overline{\boldsymbol{d}}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{x}}^{n+1/2} - \boldsymbol{S}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{y}}^{n} \overline{\boldsymbol{d}}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{y}}^{n+1/2} + \left( \boldsymbol{S}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{x}}^{n} + \boldsymbol{S}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{y}}^{n} \right) \overline{\boldsymbol{d}}\_{\boldsymbol{\xi},\boldsymbol{\varepsilon}\mathbf{y}}^{n+1/2} \right]. \tag{22}$$

In Eq. (22), g nþ1=2 <sup>ξ</sup> ¼ p nþ1=2 <sup>ξ</sup> <sup>þ</sup> qn <sup>ξ</sup>. The methods to calculate the materials' artificial viscosity qn <sup>ξ</sup> 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{aligned} p^{n+1/2} &= \sum\_{\xi} \psi\_{\xi} p\_{\xi}^{n+1/2}, \\ q^{n} &= \sum\_{\xi} \psi\_{\xi} q\_{\xi}^{n}, \\ \mathbf{S}^{n} &= \sum\_{\xi} \psi\_{\xi} \mathbf{S}\_{\xi}^{n}. \end{aligned} \tag{23}$$

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 ∇ � uξ, ψξ, 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{aligned} V &= \sum\_{\xi} V\_{\xi} \\ \text{or} \\ \sum \beta\_{\xi} &= 1, \end{aligned} \tag{24}$$

the consequence of which is the relation

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

> nþ1=2 <sup>i</sup>þ1,jþ<sup>1</sup>,<sup>r</sup>

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

<sup>ξ</sup>−<sup>χ</sup> � <sup>ρ</sup><sup>n</sup>

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

Mξ,i−1=2,j−1=<sup>2</sup> þ Mξ,i−1=2,j−1=<sup>2</sup> þ Mξ,iþ1=2,j−1=<sup>2</sup> þ Mξ,iþ1=2,jþ1=<sup>2</sup>

pn <sup>ξ</sup> <sup>¼</sup> <sup>P</sup>ðρ<sup>n</sup> ξ,e n

nþ1=2 <sup>i</sup>,jþ<sup>1</sup> ,<sup>r</sup>

� �

<sup>−</sup>V<sup>n</sup><sup>Þ</sup>

� �, <sup>M</sup><sup>ξ</sup> <sup>¼</sup> ρξ � βξ � <sup>V</sup>,

þ ∇ � S <sup>n</sup> � �

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

nþ1=2

� ∇g

<sup>n</sup> <sup>þ</sup> <sup>τ</sup> � <sup>u</sup><sup>n</sup>þ<sup>1</sup>

<sup>ξ</sup>,yyd nþ1=2 <sup>ξ</sup>,yy <sup>−</sup> <sup>2</sup>Sn

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

h i

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

<sup>ξ</sup> <sup>þ</sup> <sup>τ</sup> � <sup>V</sup><sup>n</sup>þ<sup>1</sup>

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

ρ<sup>n</sup>þ<sup>1</sup> <sup>ξ</sup> <sup>¼</sup> <sup>M</sup><sup>ξ</sup> V<sup>n</sup>þ<sup>1</sup> ξ

β<sup>n</sup>þ<sup>1</sup> <sup>ξ</sup> <sup>¼</sup> <sup>V</sup><sup>ξ</sup> nþ1=2 <sup>i</sup>þ1,<sup>j</sup> ,<sup>r</sup>

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

<sup>ξ</sup> is the speed of sound.

i,j

i,j

Þ=2, (17)

, (18)

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

<sup>ξ</sup>,xx <sup>þ</sup> Sn ξ, yy � �

d nþ1=2 ξ,ϕ

:

(22)

<sup>ξ</sup> are

, (20)

<sup>V</sup><sup>n</sup>þ<sup>1</sup> , (21)

<sup>ξ</sup>,xyd nþ1=2 <sup>ξ</sup>,xy <sup>þ</sup> Sn

<sup>ξ</sup>. The methods to calculate the materials' artificial viscosity qn

nþ1=2 i,j

<sup>i</sup>þ1=2,jþ1=<sup>2</sup> <sup>¼</sup> <sup>V</sup> <sup>r</sup>

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

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

Mi,<sup>j</sup>

) <sup>u</sup><sup>n</sup>þ<sup>1</sup>

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

nþ1=2 <sup>ξ</sup> <sup>þ</sup> qn <sup>ð</sup>u<sup>n</sup>þ<sup>1</sup> <sup>i</sup>,<sup>j</sup> −u<sup>n</sup> i,j Þ <sup>τ</sup> <sup>¼</sup> <sup>∇</sup><sup>g</sup>

<sup>i</sup>,<sup>j</sup> <sup>¼</sup> <sup>u</sup><sup>n</sup> i,j <sup>−</sup> <sup>τ</sup> Mi,<sup>j</sup> � �

> r <sup>n</sup>þ<sup>1</sup> <sup>¼</sup> <sup>r</sup>

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

<sup>ξ</sup> <sup>−</sup> <sup>S</sup><sup>n</sup>

<sup>ξ</sup>, xxd nþ1=2 <sup>ξ</sup>, xx <sup>−</sup> Sn

In Eq. (13), <sup>χ</sup> <sup>¼</sup> <sup>0</sup>:6 (this value was chosen in [17]), cn

Full timestep

74 Lagrangian Mechanics

Mi,<sup>j</sup> ¼ 0:25 �

e nþ1 <sup>ξ</sup> ¼ e n <sup>ξ</sup> <sup>−</sup> <sup>τ</sup> ρn ξ � g nþ1=2

In Eq. (22), g

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

X ξ

<sup>n</sup> <sup>þ</sup> <sup>τ</sup> � <sup>u</sup><sup>n</sup>, (9)

<sup>ξ</sup>Þ, (11)

<sup>ð</sup><sup>τ</sup> � <sup>V</sup><sup>n</sup><sup>Þ</sup> , (12)

, (10)

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

, (15)

, (16)

(14)

$$\begin{aligned} \Delta V &= \sum\_{\xi} \Delta V\_{\xi}, \\ \text{or} \\\\ \sum \beta\_{\xi} \nabla \cdot \mathbf{u}\_{\xi} &= \nabla \cdot \mathbf{u}. \end{aligned} \tag{25}$$

The natural extension of relation (25) is <sup>D</sup> <sup>¼</sup> <sup>X</sup>βξDξ, which is fulfilled at

$$\mathbf{D}\_{\xi} = \mathbf{D} \left( \frac{\nabla \cdot \mathbf{u}\_{\xi}}{\nabla \cdot \mathbf{u}} \right) . \tag{26}$$

Formula (26) is used to determine D<sup>ξ</sup> in Eq. (22) and when approximating Eq. (16).

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

$$
\varepsilon = \sum \alpha\_{\xi} e\_{\xi},
\tag{27}
$$

where the mass fraction αξ <sup>¼</sup> <sup>M</sup><sup>ξ</sup> <sup>M</sup> is given by the following expression:

$$\begin{aligned} \alpha\_{\xi} &= \beta\_{\xi} \cdot \frac{\rho\_{\xi}}{\rho}, \\ \text{where} \\ \rho &= \sum \beta\_{\xi} \rho\_{\xi}. \end{aligned} \tag{28}$$

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

$$
\Delta \mathcal{e} = \sum \alpha\_{\xi} \Delta \varepsilon\_{\xi},
\tag{29}
$$

where Δe<sup>ξ</sup> 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 \mathbf{e}\_{\xi} = -\frac{\tau \mathbf{g}\_{\xi}^{n+1/2}}{\rho\_{\xi}^{n}} \nabla \cdot \mathbf{u}\_{\xi}^{n+1/2}. \tag{30}
$$

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

$$-\frac{\tau g^{n+1/2}\nabla \cdot \mathbf{u}^{n+1/2}}{\rho} = -\frac{\tau \sum \beta\_{\varepsilon} g\_{\varepsilon}^{n+1/2} \nabla \cdot \mathbf{u}\_{\xi}^{n+1/2}}{\rho} + \sum \alpha\_{\varepsilon} \Delta e\_{\xi}^{\prime}. \tag{31}$$

Using Eq. (23) and the given ways of finding ∇ � u<sup>ξ</sup> from Eq. (31) one can obtain the values of Δe = <sup>ξ</sup> that represent additional changes in the materials' internal energy to meet the energy balance requirement.

Given that g nþ1=2 <sup>ξ</sup> <sup>¼</sup> gnþ1=2, it follows from Eq. (31) that <sup>X</sup>αξΔ<sup>e</sup> = <sup>ξ</sup> ¼ 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., Δe = <sup>ξ</sup> ¼ 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},
\tag{52}
$$

where the quantity λξ is determined by the chosen model of distributing the total divergence of the mixed cell to the constituent materials from the relation

$$
\nabla \cdot \mathbf{u}\_{\xi} = \lambda\_{\xi} \nabla \cdot \mathbf{u}.\tag{33}
$$
