Numerical Simulation of a High-Speed Impact of Metal Plates Using a Three-Fluid Model

Abstract

The process of wave formation at the contact boundary of colliding metal plates is a fundamental basis of explosive welding technology. In this case, the metals are in a pseudo-liquid state at the initial stages of the process, and from a mathematical point of view, a wave formation process can be described by compressible multiphase models. The work is devoted to the development of a three-fluid mathematical model based on the Baer–Nunziato system of equations and a corresponding numerical algorithm based on the HLL and HLLC methods, stiff pressure, and velocity relaxation procedures for simulation of the high-speed impact of metal plates in a one-dimensional statement. The problem of collision of a lead plate at a speed of 500 m/s with a resting steel plate was simulated using the developed model and algorithm. The problem statement corresponded to full-scale experiments, with lead, steel, and ambient air as three phases. The arrival times of shock waves at the free boundaries of the plates and rarefaction waves at the contact boundary of the plates, as well as the acceleration of the contact boundary after the passage of rarefaction waves through it, were estimated. For the case of a 3-mm-thick steel plate and a 2-mm-thick lead plate, the simulated time of the rarefaction wave arrival at the contact boundary constituted 1.05 μs, and it was in good agreement with the experimental value 1.1 μs. The developed numerical approach can be extended to the multidimensional case for modeling the instability of the contact boundary and wave formation in the oblique collision of plates in the Eulerian formalism.

1. Introduction

A significant number of studies have been devoted to the problem of a high-speed impact of metal plates, both from the point of view of fundamental issues of wave formation and the development of instability of the contact boundary [1,2,3], and from the practical point of view of optimizing the explosion welding process [4,5,6,7]. A historical overview of fundamental studies on the phenomenon of wave formation that occurs during the oblique high-speed impact of metal plates can be found in [1,2].

At the initial stage of the impact process, metals behave as immiscible pseudo-liquids [8], so it is appropriate to consider this problem using a heterogeneous media mechanics approach or a diffuse interface approach. There are few such studies available in the literature [9,10], and this paper is intended to partially fill this gap. The problem of high-speed impact of metal plates, in general, is at least a three-phase or three-fluid problem (material of one plate, material of another plate, ambient medium) when each phase is compressible, and its volume fraction ranges from zero to one. So, the diffuse interface models for the impact problem are generally based on the Baer–Nunziato (BN) system of equations [11], although other multiphase models can also be used [12,13]. There are two main classes of the three-fluid models for the BN equations that take their origin from the papers [14,15].

The multiphase nature of the process also makes it hard to study the development of the instability of the contact boundary in comparison with common single-phase phenomena. For example, one of the first models of wave formation during high-speed impact was proposed in [16]. The model considered the oblique interaction of a liquid jet with a moving surface covered with soft silicon putty, that is, actually the two-fluid model in the current terminology. As another example, in [17], the development of the Richtmayer–Meshkov instability was studied when a shock wave (SW) interacted with the contact surface, not between pure gas components, but with the curved interphase surface of a dense particle cloud.

It is worth noting that the problems of the high-speed collision were historically solved in the Lagrangian formalism, or, at least at some stage of the numerical algorithm, there was a Lagrangian stage [18]. In the diffuse interface methods within the framework of the Eulerian formalism, the key element is the calculation of parameters in the so-called “mixture cells”, in which various interacting materials are simultaneously present. In an early method of diffuse boundaries from [18], when describing the procedure for calculating parameters in “mixture cells”, the concept of the volume fraction of the phase was not introduced at all. The method was described as a set of heuristic numerical rules that operated with the donor and acceptor cells. The procedures for the stable calculation of “mixture cells” theoretically correspond to the stage of pressure relaxation procedure in the methods for solving the BN equations [9]. The Eulerian model can describe, for example, the generation of new interfaces without having to re-mesh the domain or at least destroy and create cells as in Lagrangian models [10]. The need to manipulate the computational grid complicates the use of Lagrangian techniques in the multidimensional case. It is the Eulerian methods that are of interest to us in this work because the further development of the methods of this work will be associated with the study of the instability of the interphase boundary, its possible significant deformations, which can create problems for Lagrangian methods due to the complexities of grid rearrangement. However, as was noted in [10], with the thickness of “the mixture region” being increased in time, the diffuse interface methods can only be used for short times. Hence, high-velocity processes (impacts, explosive phase transitions, etc.) are a natural application of such methods.

So, this paper continues our previous studies [19,20,21] and has two main goals. First, the goal is to develop a three-fluid model based on the BN system of equations and the computational algorithms Harten–Lax–van Leer (HLL) and Harten–Lax–van Leer-Contact (HLLC) (the latter method is particularly preferred because of much more accurate description of the interfacial boundaries) for the Eulerian simulation of the plates collision process. Secondly, the aim of the work is to qualitatively reproduce the main features of the full-scale experiment [8] on the high-speed normal impact of two metal plates. Note more recent studies on the planar impact of metal plates. For example, in [22], the collision of 0.4-mm-thick aluminum flyer plates with velocities of 600–700 m/s and aluminum samples with a thickness of 2.85–3 mm was considered. Our attention to the results from [8] was connected to the interesting and important mechanical effect of the development of interfacial boundary instability, reported in [8]. This effect is in close connection with the wave formation process underlying explosive welding technology.

2. Physical Statement of the Problem

We shall consider the interaction of a lead plate with the thickness h_lead and initial density ρ_lead,0 = 11300 kg/m³ with a steel plate with the thickness h_steel and initial density ρ_steel,0 = 7900 kg/m³. The lead plate is thrown in the direction of the steel one normal to the contact surface with a velocity of 500 m/s. Metal plates are surrounded by a layer of air 2 mm thick (Figure 1). It is assumed that at the initial stage of the impact during roughly the first 10 µs, the metals behave as pseudo-fluids [8] so the material of each plate can be considered as a compressible fluid (as well as the surrounding air) in a three-fluid BN-like model. With such a multiphase flow model, the same equations are solved everywhere with the same numerical scheme. This is achieved by adding a negligible quantity of the other phase in pure phases. The initial pressure is 10⁵ Pa everywhere.

3. Mathematical Model

The mathematical model is based on the approach [14] and describes a flow of three immiscible compressible fluids:

$$U_t+F_x\left(u\right)=H\left(U,{\left({\alpha }_{air}\right)}_x,{\left({\alpha }_{steel}\right)}_x\right)+\mathbf{P},$$

(1)

$$U=\left[\begin{array}{c}{\alpha }_{air}\\ {\alpha }_{air}{\rho }_{air}\\ {\alpha }_{air}{\rho }_{air}{v}_{air}\\ {\alpha }_{air}{\rho }_{air}{E}_{air}\\ {\alpha }_{steel}\\ {\alpha }_{steel}{\rho }_{steel}\\ {\alpha }_{steel}{\rho }_{steel}{v}_{steel}\\ {\alpha }_{steel}{\rho }_{steel}{E}_{steel}\\ {\alpha }_{lead}{\rho }_{lead}\\ {\alpha }_{lead}{\rho }_{lead}{v}_{lead}\\ {\alpha }_{lead}{\rho }_{lead}{E}_{lead}\end{array}\right], F=\left[\begin{array}{c}0\\ {\alpha }_{air}{\rho }_{air}{v}_{air}\\ {\alpha }_{air}({\rho }_{air}{v}_{air}^2+p_{air})\\ {\alpha }_{air}{v}_{air}({\rho }_{air}{E}_{air}+p_{air})\\ 0\\ {\alpha }_{steel}{\rho }_{steel}{v}_{steel}\\ {\alpha }_{steel}({\rho }_{steel}{v}_{steel}^2+p_{steel})\\ {\alpha }_{steel}{v}_{steel}({\rho }_{steel}{E}_{steel}+p_{steel})\\ {\alpha }_{lead}{\rho }_{lead}{v}_{lead}\\ {\alpha }_{lead}({\rho }_{lead}{v}_{lead}^2+p_{lead})\\ {\alpha }_{lead}{v}_{lead}({\rho }_{lead}{E}_{lead}+p_{lead})\end{array}\right],$$

(2)

$$H=\left[\begin{array}{c}-v_i{({\alpha }_{air})}_x\\ 0\\ p_i{({\alpha }_{air})}_x\\ p_i{v}_i{({\alpha }_{air})}_x\\ -v_i{({\alpha }_{steel})}_x\\ 0\\ p_i{({\alpha }_{steel})}_x\\ p_i{v}_i{({\alpha }_{steel})}_x\\ 0\\ p_i{({\alpha }_{lead})}_x\\ p_i{v}_i{({\alpha }_{lead})}_x\end{array}\right], \mathbf{P}=\left[\begin{array}{c}\mu (p_{air}-p_{lead})\\ 0\\ \lambda (v_{lead}-v_{air})\\ \lambda v_i(v_{lead}-v_{air})-\mu p_i(p_{air}-p_{lead})\\ \mu (p_{steel}-p_{lead})\\ 0\\ \lambda (v_{lead}-v_{steel})\\ \lambda u_i(v_{lead}-v_{steel})-\mu p_i(p_{steel}-p_{lead})\\ 0\\ \lambda (v_{air}-v_{lead})\\ \lambda u_i(v_{air}-v_{lead})-\mu p_i(p_{lead}-p_{air})\end{array}\right],$$

(3)

$${\alpha }_{air}+{\alpha }_{steel}+{\alpha }_{lead}=1, E_k=\frac{v_k^2}{2}+e_k\left(p_k,{\rho }_k\right)=\frac{v_k^2}{2}+\frac{p_k+{\gamma }_k{P}_{0k}}{{\rho }_k\left({\gamma }_k-1\right)},$$

(4)

$$p_i={\displaystyle \sum _{k=1}^3{\alpha }_k{p}_k}, v_i={\displaystyle \sum _{k=1}^3{\alpha }_k{\rho }_k{v}_k}/{\displaystyle \sum _{k=1}^3{\alpha }_k{\rho }_k}.$$

(5)

Here U denotes the vector of “conservative” variables, F is the differential flux vector, H is the differential source term, and P is the algebraic source term, connected to the relaxation processes. The notation is conventional: t is the time, x is the spatial coordinate, α is the volume fraction, ρ is the true density, v is the velocity, p is the pressure, e is the specific internal energy, and E is the specific total energy. Parameters with indexes k = 1, 2, 3 correspond to the air, steel, and lead phases, respectively. Further, we use either these indexes or direct subscripts “air”, “steel”, and “lead” for clarity. Pressure and velocity with the “i” index correspond to the parameters at the fluid interfaces.

For each phase, the stiffened gas equation of state (EOS) with the parameters γ and P₀ is used. In general, the properties of the metals EOS are extremely important for quantitatively and even qualitatively correct simulation of the characteristics of the impact process especially if the parameters of the impact imply, for instance, phase transitions. In [23], the hypervelocity impact problem was simulated using three different EOSs for lead that provided the effects of melting in strong SWs, evaporation in rarefaction waves (RW), and spallation. The specific features in the simulation results with the use of different EOSs were discussed. For the BN-type models, Mie–Gruneisen EOS (more general than stiffened gas EOS) can be applied [9,24,25] and the influence of this type of EOS on the simulation results is the subject of further study.

Parameter μ is responsible for pressure relaxation, and λ is for velocity relaxation. In the present work, we assume the instantaneous pressure and velocity equilibrium at the interface boundaries:

$$p_{air}=p_{lead}=p_{steel},$$

(6)

$$v_{air}=v_{lead}=v_{steel}.$$

(7)

Conditions (6) and (7) correspond to the stiff relaxation assumption μ→+∞, λ→+∞ [9,11] (see Section 4.3 and Section 4.4).

4. Numerical Algorithm

4.1. Splitting Scheme

The computational algorithm was based on the Strang splitting principle:

$${\mathbf{U}}_j^{n+1}=L_p{L}_v{L}_h{\mathbf{U}}_j^n.$$

(8)

We denote vector grid function by the same letter U as an unknown function in the defining system (1) – (5) but with the spatial index j and the time index n. At the hyperbolic stage of the algorithm denoted in (8) as operator L_h, the defining system of Equations (1)–(5) was solved with P = 0. After that, a velocity relaxation procedure L_v was implemented. Finally, a pressure relaxation procedure L_p was carried out. Consider each of the stages in more detail.

4.2. Hyperbolic Step

The computational domain was a one-dimensional segment, which was divided into N uniform cells. The size of the computational cell was Δx. Two approaches for the hyperbolic step were implemented, namely HLL and HLLC schemes.

4.2.1. HLL Method

Consider the HLL method for the three-fluid model first [9,14]. To begin with, the initial system (1)–(5) is formally divided into the first and fifth equations and the remaining sub-system:

$$u_t+f_x\left(u\right)=h\left(u,{\left({\alpha }_{air}\right)}_x,{\left({\alpha }_{air}\right)}_x\right),$$

(9)

$$u=\left[\begin{array}{c}{\alpha }_{air}{\rho }_{air}\\ {\alpha }_{air}{\rho }_{air}{v}_{air}\\ {\alpha }_{air}{\rho }_{air}{E}_{air}\\ {\alpha }_{steel}{\rho }_{steel}\\ {\alpha }_{steel}{\rho }_{steel}{v}_{steel}\\ {\alpha }_{steel}{\rho }_{steel}{E}_{steel}\\ {\alpha }_{lead}{\rho }_{lead}\\ {\alpha }_{lead}{\rho }_{lead}{v}_{lead}\\ {\alpha }_{lead}{\rho }_{lead}{E}_{lead}\end{array}\right], f=\left[\begin{array}{c}{\alpha }_{air}{\rho }_{air}{v}_{air}\\ {\alpha }_{air}({\rho }_{air}{v}_{air}^2+p_{air})\\ {\alpha }_{air}{v}_{air}({\rho }_{air}{E}_{air}+p_{air})\\ {\alpha }_{steel}{\rho }_{steel}{v}_{steel}\\ {\alpha }_{steel}({\rho }_{steel}{v}_{steel}^2+p_{steel})\\ {\alpha }_{steel}{v}_{steel}({\rho }_{steel}{E}_{steel}+p_{steel})\\ {\alpha }_{lead}{\rho }_{lead}{v}_{lead}\\ {\alpha }_{lead}({\rho }_{lead}{v}_{lead}^2+p_{lead})\\ {\alpha }_{lead}{v}_{lead}({\rho }_{lead}{E}_{lead}+p_{lead})\end{array}\right], h=\left[\begin{array}{c}0\\ p_i{({\alpha }_{air})}_x\\ p_i{v}_i{({\alpha }_{air})}_x\\ 0\\ p_i{({\alpha }_{steel})}_x\\ p_i{v}_i{({\alpha }_{steel})}_x\\ 0\\ p_i{({\alpha }_{lead})}_x\\ p_i{v}_i{({\alpha }_{lead})}_x\end{array}\right].$$

(10)

Finite-volume approximations of the detached equations and the system (9), (10) are the following:

$$\begin{array}{cc}{\left({\alpha }_k\right)}_{h,j}^{n+1}=& {\left({\alpha }_k\right)}_j^n-\frac{\Delta t^n}{\Delta x}\left[\frac{{\left(v_i\right)}_j^n\left(S_{j+1/2}^{+}{\left({\alpha }_k\right)}_j^n-S_{j+1/2}^{-}{\left({\alpha }_k\right)}_{j+1}^n\right)+S_{j+1/2}^{+}{S}_{j+1/2}^{-}\left({\left({\alpha }_k\right)}_{j+1}^n-{\left({\alpha }_k\right)}_j^n\right)}{S_{j+1/2}^{+}-S_{j+1/2}^{-}}-\right.\\ & \left.-\frac{{\left(v_i\right)}_j^n\left(S_{j-1/2}^{+}{\left({\alpha }_k\right)}_{j-1}^n-S_{j-1/2}^{-}{\left({\alpha }_k\right)}_j^n\right)+S_{j-1/2}^{+}{S}_{j-1/2}^{-}\left({\left({\alpha }_k\right)}_j^n-{\left({\alpha }_k\right)}_{j-1}^n\right)}{S_{j-1/2}^{+}-S_{j-1/2}^{-}}\right], k=1,2,\end{array}$$

(11)

$$u_{h,j}^{n+1}=u_j^n-\frac{\Delta t^n}{\Delta x}\left[{\mathbf{f}}_{j+1/2}\left(u_j^n,u_{j+1}^n\right)-{\mathbf{f}}_{j-1/2}\left(u_j^n,u_{j+1}^n\right)\right]+\mathbf{h}\left(u_j^n,\Delta {\left({\alpha }_{air}\right)}_x,\Delta \left({\alpha }_{steel}\right)\right).$$

(12)

The numerical flux f_j+1/2 in (12) is calculated using the HLL scheme:

$${\mathbf{f}}_{j+1/2}\left(u_j^n,u_{j+1}^n\right)=\frac{S_{j+1/2}^{+}\mathbf{f}\left({\mathbf{u}}_j^n\right)-S_{j+1/2}^{-}\mathbf{f}\left({\mathbf{u}}_{j+1}^n\right)+S_{j+1/2}^{+}{S}_{j+1/2}^{-}\left({\mathbf{u}}_{j+1}^n-{\mathbf{u}}_j^n\right)}{S_{j+1/2}^{+}-S_{j+1/2}^{-}},$$

(13)

$$S_{j+1/2}^{+}=\underset{k=1,2,3}{\max }\left(0,{\left(v_k\right)}_j^n+{\left(c_k\right)}_j^n,{\left(v_k\right)}_{j+1}^n+{\left(c_k\right)}_{j+1}^n\right), S_{j+1/2}^{-}=\underset{k=1,2,3}{\min }\left(0,{\left(v_k\right)}_j^n-{\left(c_k\right)}_j^n,{\left(v_k\right)}_{j+1}^n-{\left(c_k\right)}_{j+1}^n\right),$$

(14)

where c is the speed of sound:

$$c_k=\sqrt{{\gamma }_k\frac{p_k+P_{0k}}{{\rho }_k}}, k=1,2,3.$$

(15)

Spatial derivatives of α_air and α_steel in the non-conservative term in (12) are approximated as:

$$\Delta {\left({\alpha }_k\right)}_x=\frac{1}{\Delta x}\left(\frac{S_{j+1/2}^{+}{\left({\alpha }_k\right)}_j^n-S_{j+1/2}^{-}{\left({\alpha }_k\right)}_{j+1}^n}{S_{j+1/2}^{+}-S_{j+1/2}^{-}}-\frac{S_{j-1/2}^{+}{\left({\alpha }_k\right)}_{j-1}^n-S_{j-1/2}^{-}{\left({\alpha }_k\right)}_j^n}{S_{j-1/2}^{+}-S_{j-1/2}^{-}}\right), k=1,2.$$

(16)

The latter approximation, consistent with the approximation of the divergent part of the system (9), (10) ensures the so-called “p-v” condition [26]. This condition claims that a multiphase system, uniform in velocity and pressure, should remain uniform through its evolution.

The time step is calculated dynamically from the condition:

$$\Delta t^n=\mathrm{CFL}\cdot {\min }_{k,j}\left(\frac{\Delta x}{\left|{\left(v_k\right)}_j^n\right|+{\left(c_k\right)}_j^n}\right),$$

(17)

where CFL is a coefficient between 0 and 1.

4.2.2. HLLC Method

To improve the quality of simulations, the HLLC method from [27,28] was extended for the three-fluid model. From the general point of view, the methods in [27,28] are less accurate than later developments [24,25,29] because they are actually a straightforward generalization of the common HLLC for the Euler equations (meanwhile, the approximation of non-conservative differential terms and equations for the volume fractions evolution are of great importance) and do not use the elements of the Riemann problem solution for the BN equations. However, the methods in [27,28] are much easier to implement and to extend to the three-fluid case.

Again as for the HLL method (11)–(17), we shall consider approximations of the first and the fifth equations of the system (1) and sub-system (9). Equations of volume fraction evolution are approximated as:

$${\left({\alpha }_k\right)}_{h,j}^{n+1}={\left({\alpha }_k\right)}_j^n-\frac{\Delta t^n}{\Delta x}\left[{\phi }_{k,j+1/2}-{\phi }_{k,j-1/2}\right], k=1,2,$$

(18)

$${\phi }_{k,j+1/2}=\left\{\begin{array}{l}{\left({\alpha }_k\right)}_j^n, \mathrm{if} \hspace{0.17em}{\left(v_i\right)}_j^n\ge 0,\\ {\left({\alpha }_k\right)}_{j+1}^n, \mathrm{else}.\end{array}\right.$$

(19)

The general finite-volume scheme (12) for the sub-system (9), (10) remains the same. The numerical flux f_j+1/2 is calculated using the HLLC scheme:

$${\mathbf{f}}_{j+1/2}\left(u_j^n,u_{j+1}^n\right)=\left\{\begin{array}{l}\mathbf{f}\left({\mathbf{u}}_j^n\right),\hspace{0.17em} \mathrm{if} \hspace{0.17em}{S}_{j+1/2}^{-}\ge 0,\\ \mathbf{f}\left({\mathbf{u}}_j^n\right)+S_{j+1/2}^{-}\left({\mathbf{q}}_j^n-{\mathbf{u}}_j^n\right),\hspace{0.17em} \mathrm{if} \hspace{0.17em}{S}_{j+1/2}^{-}<0\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}{S}_{j+1/2}^{*}\ge 0,\\ \mathbf{f}\left({\mathbf{u}}_{j+1}^n\right)+S_{j+1/2}^{+}\left({\mathbf{q}}_{j+1}^n-{\mathbf{u}}_{j+1}^n\right),\hspace{0.17em} \mathrm{if} \hspace{0.17em}{S}_{j+1/2}^{*}<0\hspace{0.17em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}{S}_{j+1/2}^{+}\ge 0,\\ \mathbf{f}\left({\mathbf{u}}_{j+1}^n\right),\hspace{0.17em} \mathrm{if} \hspace{0.17em}{S}_{j+1/2}^{+}\le 0.\end{array}\right.$$

(20)

$$S_{j+1/2}^{+}=\underset{k=1,2,3}{\max }\left({\left(v_k\right)}_j^n+{\left(c_k\right)}_j^n,{\left(v_k\right)}_{j+1}^n+{\left(c_k\right)}_{j+1}^n\right), S_{j+1/2}^{-}=\underset{k=1,2,3}{\min }\left({\left(v_k\right)}_j^n-{\left(c_k\right)}_j^n,{\left(v_k\right)}_{j+1}^n-{\left(c_k\right)}_{j+1}^n\right),$$

(21)

$$S_{j+1/2}^{*}=\frac{{\left(p_i\right)}_{j+1}^n-{\left(p_i\right)}_j^n+{\left({\rho }_i\right)}_j^n{\left(v_i\right)}_j^n\left(S_{j+1/2}^{-}-{\left(v_i\right)}_j^n\right)-{\left({\rho }_i\right)}_{j+1}^n{\left(v_i\right)}_{j+1}^n\left(S_{j+1/2}^{+}-{\left(v_i\right)}_{j+1}^n\right)}{{\left({\rho }_i\right)}_j^n\left(S_{j+1/2}^{-}-{\left(v_i\right)}_j^n\right)-{\left({\rho }_i\right)}_{j+1}^n\left(S_{j+1/2}^{+}-{\left(v_i\right)}_{j+1}^n\right)},$$

(22)

$${\rho }_i={\displaystyle \sum {\alpha }_k{\rho }_k},$$

(23)

$${\mathbf{q}}_j^n=\left[\begin{array}{c}{\mathbf{q}}_{j,air}^n\\ {\mathbf{q}}_{j,steel}^n\\ {\mathbf{q}}_{j,lead}^n\end{array}\right], {\mathbf{q}}_{j+1}^n=\left[\begin{array}{c}{\mathbf{q}}_{j+1,air}^n\\ {\mathbf{q}}_{j+1,steel}^n\\ {\mathbf{q}}_{j+1,lead}^n\end{array}\right],$$

(24)

$${\mathbf{q}}_{j,k}^n=\left[\begin{array}{c}{C}_{j,k}\\ C_{j,k}{S}_{j+1/2}^{*}\\ C_{j,k}\left[{\left(p_k\right)}_j^n/{\left({\rho }_k\right)}_j^n+\left(S_{j+1/2}^{*}-{\left(v_k\right)}_j^n\right)\left(S_{j+1/2}^{*}+{\left(p_k\right)}_j^n/{\left({\rho }_k\right)}_j^n/\left(S_{j+1/2}^{-}-{\left(v_k\right)}_j^n\right)\right)\right]\end{array}\right], k=1,2,3,$$

(25)

$${\mathbf{q}}_{j+1,k}^n=\left[\begin{array}{c}{C}_{j+1,k}\\ C_{j+1,k}{S}_{j+1/2}^{*}\\ C_{j+1,k}\left[{\left(p_k\right)}_{j+1}^n/{\left({\rho }_k\right)}_{j+1}^n+\left(S_{j+1/2}^{*}-{\left(v_k\right)}_{j+1}^n\right)\left(S_{j+1/2}^{*}+{\left(p_k\right)}_{j+1}^n/{\left({\rho }_k\right)}_{j+1}^n/\left(S_{j+1/2}^{+}-{\left(v_k\right)}_{j+1}^n\right)\right)\right]\end{array}\right], k=1,2,3,$$

(26)

$$C_{j,k}=\frac{{\left({\alpha }_k\right)}_j^n\left[S_{j+1/2}^{-}-{\left(v_{air}\right)}_j^n\right]}{S_{j+1/2}^{-}-S_{j+1/2}^{*}}, C_{j+1,k}=\frac{{\left({\alpha }_k\right)}_{j+1}^n\left[S_{j+1/2}^{+}-{\left(v_{air}\right)}_{j+1}^n\right]}{S_{j+1/2}^{+}-S_{j+1/2}^{*}}, k=1,2,3.$$

(27)

Spatial derivatives of α_air and α_steel in the non-conservative term in (12) are approximated as:

$$\Delta {\left({\alpha }_k\right)}_x=\frac{1}{\Delta x}\left({\psi }_{k,j+1/2}-{\psi }_{k,j-1/2}\right), k=1,2,$$

(28)

$${\psi }_{k,j+1/2}=\left\{\begin{array}{l}{\left({\alpha }_k\right)}_j^n, \mathrm{if} \hspace{0.17em}{S}_{j+1/2}^{*}\ge 0,\\ {\left({\alpha }_k\right)}_{j+1}^n, \mathrm{else}.\end{array}\right.$$

(29)

Simulations in Section 5 below show that the HLLC method (18)–(29) was less robust than the HLL one. The initial volume fraction of air α_air, i.e., the small parameter, in the regions that corresponded to the steel and lead plates, and in the HLLC simulations, should be tuned more carefully than in the HLL simulations. For example, the HLLC simulation failed if α_air was equal to 10^–6.

4.3. Velocity Relaxation

For the velocity relaxation, the following system of ordinary differential equations (ODE) for each phase k = 1, 2, 3 is solved in each computational cell [14]:

$$\left\{\begin{array}{l}\frac{\partial {\alpha }_k}{\partial t}=0,\\ \frac{\partial ({\alpha }_k{\rho }_k)}{\partial t}=0,\\ \frac{\partial ({\alpha }_k{\rho }_k{v}_k)}{\partial t}=\lambda \left(v_m-v_k\right),\\ \frac{\partial ({\alpha }_k{\rho }_k{E}_k)}{\partial t}=\lambda u_i\left(v_m-v_k\right),\end{array}\right.$$

(30)

where m is any index not equal to k. To achieve velocity equilibrium (7), the relaxation coefficient λ in (30) is assumed to tend to infinity. Such a proposal leads to the equation for equilibrium velocity:

$$\overline{v}=\frac{{\displaystyle \sum _k{({\alpha }_k{\rho }_k{v}_k)}_0}}{{\displaystyle \sum _k{({\alpha }_k{\rho }_k)}_0}},$$

(31)

where index “0” denotes values obtained after the hyperbolic step. After that, the specific internal energy correction is required for each phase:

$$e_k=e_{k,0}+\frac{1}{2}{\left(\overline{v}-v_{k,0}\right)}^2.$$

(32)

As the true densities of the phases do not change during velocity relaxation (as well as volume fractions), the correction of the internal energies (32) leads to the correction of phase pressures only.

4.4. Pressure Relaxation

At the pressure relaxation step, it is required to solve the following ODE system for each phase in each computational cell [14]:

$$\left\{\begin{array}{l}\frac{\partial {\alpha }_k}{\partial t}=\mu \left(p_k-p_m\right),\\ \frac{\partial \left({\alpha }_k{\rho }_k\right)}{\partial t}=0,\\ \frac{\partial \left({\alpha }_k{\rho }_k{v}_k\right)}{\partial t}=0,\\ \frac{\partial \left({\alpha }_k{\rho }_k{E}_k\right)}{\partial t}=\mu p_i\left(p_k-p_m\right),\end{array}\right.$$

(33)

where m is again any index not equal to k. The pressure relaxation procedure follows the same principle as velocity relaxation. When the relaxation parameter μ in (33) tends to infinity, all pressures after relaxation must be equal. Combining the volume fraction, mass, momentum, and energy equation, we get [14]

$$\frac{\partial e_k}{\partial t}=-p_i\frac{\partial }{\partial t}\left(\frac{1}{{\rho }_k}\right).$$

(34)

Trapezoidal approximation of (34) gives:

$${\overline{e}}_k-e_k^0=-\frac{1}{2}\left(\overline{p}+p_i^0\right)\left(\frac{1}{{\overline{\rho }}_k}-\frac{1}{{\rho }_k^0}\right),$$

(35)

where index “0” denotes values obtained after the velocity relaxation procedure, and values with bars are the relaxed ones. So, using the EOS for each phase, saturation constraint ${\alpha }_1+{\alpha }_2+{\alpha }_3=1$ and Formula (35) it is possible to obtain a system of non-linear equation for variables ${\overline{\rho }}_{air},{\overline{\rho }}_{steel},{\overline{\rho }}_{lead},\hspace{0.17em}\overline{p}$:

$$\{\begin{array}{c}2{\overline{\rho }}_{air}{\rho }_{air}^0\left(\frac{\overline{p}+{\gamma }_{air}{P}_{0air}}{{\overline{\rho }}_{air}\left({\gamma }_{air}-1\right)}-\frac{p_{air}^0+{\gamma }_{air}{P}_{0air}}{{\rho }_{air}^0\left({\gamma }_{air}-1\right)}\right)-\left(\overline{p}+p_i^0\right)\left({\overline{\rho }}_{air}-{\rho }_{air}^0\right)=0,\\ 2{\overline{\rho }}_{steel}{\rho }_{steel}^0\left(\frac{\overline{p}+{\gamma }_{steel}{P}_{0steel}}{{\overline{\rho }}_{steel}\left({\gamma }_{steel}-1\right)}-\frac{p_{steel}^0+{\gamma }_{steel}{P}_{0steel}}{{\rho }_{steel}^0\left({\gamma }_{steel}-1\right)}\right)-\left(\overline{p}+p_i^0\right)\left({\overline{\rho }}_{steel}-{\rho }_{steel}^0\right)=0,\\ 2{\overline{\rho }}_{lead}{\rho }_{lead}^0\left(\frac{\overline{p}+{\gamma }_{lead}{P}_{0lead}}{{\overline{\rho }}_{lead}\left({\gamma }_{lead}-1\right)}-\frac{p_{lead}^0+{\gamma }_{lead}{P}_{0lead}}{{\rho }_{lead}^0\left({\gamma }_{lead}-1\right)}\right)-(\overline{p}+p_i^0)({\overline{\rho }}_{lead}-{\rho }_{lead}^0)=0,\\ {\displaystyle \sum m_k/{\overline{\rho }}_k}-1=0,\end{array}$$

(36)

$$m_k={\alpha }_k^0{\rho }_k^0.$$

(37)

In the two-phase case, the analogous system can be reduced to one quadratic equation that can be solved analytically [30,31]. In the current case, as in the case of more complicated equilibrium conditions (for example, taking into account intergranular stresses in the gas—particles simulations [32]), the non-linear system (36)–(37) should be solved:

$$\mathbf{F}\left(\mathbf{X}\right)=0,$$

(38)

$$\mathbf{X}=\left[\begin{array}{c}{\overline{\rho }}_{air}\\ {\overline{\rho }}_{steel}\\ {\overline{\rho }}_{lead}\\ \overline{p}\end{array}\right], \mathbf{F}\left(\mathbf{X}\right)=\left[\begin{array}{c}2{\overline{\rho }}_{air}{\rho }_{air}^0\left(\frac{\overline{p}+{\gamma }_{air}{P}_{0air}}{{\rho }_{air}\left({\gamma }_{air}-1\right)}-\frac{p_{air}^0+{\gamma }_{air}{P}_{0air}}{{\rho }_{air}^0\left({\gamma }_{air}-1\right)}\right)-\left(\overline{p}+p_i^0\right)\left({\overline{\rho }}_{air}-{\rho }_{air}^0\right)\\ 2{\overline{\rho }}_{steel}{\rho }_{steel}^0\left(\frac{\overline{p}+{\gamma }_{steel}{P}_{0steel}}{{\rho }_{steel}\left({\gamma }_{steel}-1\right)}-\frac{p_{steel}^0+{\gamma }_{steel}{P}_{0steel}}{{\rho }_{steel}^0\left({\gamma }_{steel}-1\right)}\right)-\left(\overline{p}+p_i^0\right)\left({\overline{\rho }}_{steel}-{\rho }_{steel}^0\right)\\ 2{\overline{\rho }}_{lead}{\rho }_{lead}^0\left(\frac{\overline{p}+{\gamma }_{lead}{P}_{0lead}}{{\rho }_{lead}\left({\gamma }_{lead}-1\right)}-\frac{p_{lead}^0+{\gamma }_{lead}{P}_{0lead}}{{\rho }_{lead}^0\left({\gamma }_{lead}-1\right)}\right)-\left(\overline{p}+p_i^0\right)\left({\overline{\rho }}_{lead}-{\rho }_{lead}^0\right)\\ {\displaystyle \sum m_k/{\overline{\rho }}_k}-1\end{array}\right]$$

(39)

The solution of the system (38) – (39) is obtained with the Newton–Raphson method:

$${\mathbf{X}}^{\left(s+1\right)}={\mathbf{X}}^{\left(s\right)}-{\mathrm{J}}^{-1}\left({\mathbf{X}}^{\left(s\right)}\right)F\left({\mathbf{X}}^{\left(s\right)}\right),$$

(40)

$$\mathbf{J}\left(\mathbf{X}\right)=\frac{\partial \mathbf{F}}{\partial \mathbf{X}}=\left(\begin{array}{cccc}{A}_{air}& 0& 0& B_{air}\\ 0& A_{steel}& 0& B_{steel}\\ 0& 0& A_{lead}& B_{lead}\\ -m_{air}/{\overline{\rho }}_{air}^2& -m_{steel}/{\overline{\rho }}_{steel}^2& -m_{lead}/{\overline{\rho }}_{lead}^2& 0\end{array}\right),$$

(41)

$$A_k=-2\cdot \frac{p_k^0+{\gamma }_k{P}_{0k}}{{\gamma }_k-1}-\left(\overline{p}+p_i^0\right), B_k=\frac{2{\rho }_k^0}{{\gamma }_k-1}-\left({\overline{\rho }}_k+{\rho }_k^0\right), k=1,2,3,$$

(42)

where s is a current iteration number in the repetitive process (40)–(42).

5. Simulation Results

To obtain quantitatively valid characteristics of the collision process, the parameters of the stiffened-gas equations of state for the fluids in consideration must be calibrated using either experimental data or the data computed using the real-life wide range metal EOS. Parameters of the EOS for steel and lead were taken close to those found in our previous works [19,20]. Parameters were obtained by comparing the computed characteristics of the SWs formed at collision with the computed data obtained by the procedure described in [33], which used the wide range EOS for metals [34] (this procedure was implemented in [35]). We took γ = 3.0 and P₀ = 65.0 GPa for steel; and γ = 2.7 and P₀ = 15.5 GPa for lead. The maximal relative differences in SWs speed, pressure, and velocity at the contact surface and the densities behind the SW fronts in our simulations do not exceed 7% in comparison with the reference values. For air, the values γ = 1.4 and P₀ = 0 were used.

Simulations were carried out with a spatial resolution of 2.5 μm. For example, a plate with a thickness of 2 mm was resolved with the use of 800 computational cells. The CFL number in (17) was equal to 0.5. We set the non-penetrating condition on both boundaries of the computational area.

Consider the main stages of the process for the case of h_steel = 3 mm, h_lead = 2 mm. Figure 2 illustrates spatial distributions of metals’ volume fractions and pressure. The initial impact causes the formation of two SWs propagating in the opposite directions towards the free surfaces of the plates (see Figure 2a). At a time instant of about τ_steel = 0.6 μs after the beginning of plates interaction, SW in steel reaches a free boundary (see Figure 2b). As a result of an SW interaction with a free boundary of the plate, RW forms. It moves towards the interface between colliding plates (see Figure 2c). At a time instant of about T_steel = 1.1 μs, RW from steel reaches the interface and moves through it (see Figure 2d). Analogous processes occur in the lead plate. The times τ_steel, τ_lead, T_steel, and T_lead depend on the thicknesses of the plates. In the considered case, RW from the free boundary of the steel plate reaches the interface boundary between the plates faster than the RW from the free boundary of the lead plate. As is seen from Figure 2d, an RW from the steel plate interacts with a high-pressure region in the lead plate and the time instant T_lead has no sense in this case. As a result of an RW passing, the interface accelerates rapidly in the direction of a free boundary of the steel plate. This acceleration together with the density gradient in the opposite direction can cause a development of the Rayleigh–Taylor instability of the interface. This question was studied in [21] using three-dimensional Euler equations.

Figure 3 illustrates the dynamics of the velocity and the acceleration of the interface boundary between the plates. Both plots were obtained by numerical differentiation of the dependency x_1/2-t, where x_1/2 is a “center” of the diffuse interface between the plates (the coordinate of the computational cell center, where α_steel = α_lead ≈ 1/2 are accurate within the existence of a small amount of air, α_air ≈ 10^–5). The smooth region on the black curve between time instants of about 1.1 μs and 1.3 μs corresponds to the acceleration of the interfacial boundary between the plates due to the RW arrival. The interfacial boundary speed increased from about 195 m/s up to 425 m/s. The orange curve at the same time interval demonstrates a pronounced peak. Trembling of the acceleration curve up to about 1.1 μs has a numerical nature. It vanishes at the later time instants with the increase of the interface boundary speed after the RW passing and thus with the decrease of the effect of uncertainty of the x_1/2 definition due to the spatial discretization of the computational domain. Figure 3 gives quantitative estimations of the realized acceleration and the time of its action. These data can be used for the estimations of the possibility of the Rayleigh–Taylor instability development of the interfacial boundary. After the meeting of two RWs from the free boundaries of the plates after about 1.5 μs, the simulation crashed because of a negative pressure occurrence. The problem was not due to the numerical instability but rather due to going beyond the limits of applicability of the gas-dynamics model of the process in use. Large tensile forces that can lead to the formation of internal cracks in the physical experiment demand the usage of models like [10] at the following stages of the process of plates’ collision.

Following experiments [8], parametric simulations were performed for various thicknesses of a steel plate. The results in terms of time instants τ and T are presented in

Table 1. Times of arrival of an SW to the free boundary of a plate (τ) and arrival of an RW to the interfacial boundary between the plates (T) for the steel and lead plates; h_lead = 2 mm; simulations using the HLL method.

No.	hsteel, mm	τsteel, μs	τlead, μs	Tsteel, μs	Tlead, μs
1	2	0.38	0.8	0.70	–
2	3	0.54		1.05	–
3	4	0.72		1.38	1.38
4	4.5	0.82		–	1.38
5	5	0.91		–	1.38
6	10	1.85		–	1.38
7	20	×		–	1.38

. In this series of simulations, the HLL scheme was used as the results in terms of time instants τ and T were close for both HLL and HLLC methods. The second row in the table corresponds to the problem statement considered above (see Figure 2). A label “–” instead of the T value means that an RW from another plate reached the interfacial boundary first. A label “×” means that a time instant could not be measured because of the simulation crash after two RWs meeting. For a small thickness of the steel plate (up to 4 mm), the dynamics of the interfacial boundary motion is determined by an RW from the steel plate side. For h_steel = 4.5 mm, time instants τ_steel and τ_lead become almost equal. With the following increase of the h_steel, the process is determined by an RW from the lead plate. The experimental value of T_steel for simulation No. 2 is available and constitutes 1.1 μs [8]. The calculated value is 1.05 μs and it is very close to the experimental one.

6. Discussion

Let us compare the results of simulations using the HLL scheme (11)–(17) and the HLLC scheme (18)–(29) on the hyperbolic step. Figure 2a shows that SWs in both simulations are described almost identically. As expected, maximal differences were obtained in the description of the interfacial boundaries—between two plates and between free surfaces of the plates and ambient air. The non-moving free surface of the steel plate up to the moment of the SW arrival was not smeared at all in the HLLC simulation (see Figure 2a). This important property, valid for the HLLC method for the Euler equations solution, was inherited by the HLLC method for the BN equations [27,28]. Moving interfacial boundaries are smeared by both HLL and HLLC methods, but the HLL diffuse interface is an order of magnitude greater than the HLLC diffusive interface. Apparently, it is one of the reasons that the RW that occurs after SW arrives at the free boundary of the plate is described in the simulation using the HLL scheme with very large errors, compared with the HLLC approach. The pressure to the left of the tail of the RW in steel in Figure 2c,d does not fall to 1 atm in the HLL simulation. The RW is not localized with a long non-physical tail. On the contrary, the RW in the HLLC simulation has a profile close to the SW, which is typical for wave processes in condensed media.

Our further goal was to develop a two- and three-dimensional algorithm for studying the development of the instability of the interfacial boundary of the metal plates during oblique collision using a three-fluid model. The developed one-dimensional model of the process and the computational algorithm allows such extensions, for example, using a multidimensional algorithm from [14]. A two-dimensional three-fluid HLL approach was realized. Figure 4 demonstrates preliminary statement and simulation results for the problem of the oblique impact of the plates. The initial velocity of the lead plate was equal to 500 m/s and was directed normally to its surface (see Figure 4a). The initial angle between the surfaces of the plates was equal to 5°. All boundaries were free. The computational grid was uniform with a cell size equal to 10^–2 mm. Figure 4b demonstrates non-planar SWs, originating after an oblique impact, and the very beginning of the wave formation process at the interfacial boundary between the plates. The major problem we have faced in the two-dimensional simulation was the motion of the air in the gap between the plates being too fast due to the velocity relaxation procedure (30)–(32). Apparently, this part of the model was not very relevant to the real process in the multi-dimensional case. Note the up-to-date studies of such types of flow due to the body impact [36,37]. It was reported in [36] that the flow between a base and a cladding plate can affect the process of explosive welding. Further multi-dimensional simulations require the use of the HLLC scheme, at least an order of magnitude more detailed computational grids, and, probably, another problem formulation with the resolution of the angle of the lead plate for a more correct description of the dynamics of the contact point between the plates.

7. Conclusions

Thus, we adapted a multiphase model [14] from the class of diffuse interface approaches for the consideration of the problem of a high-speed collision of metal plates, which is a theoretical foundation for the problem of explosive welding. In addition to the common HLL method for solving BN equations [9], we adapted the HLLC method [27,28] for the three-phase case, which showed qualitatively better results than the HLL method. Nowadays, the most developed numerical approaches for the high-speed impact simulations include Arbitrary Lagrangian–Eulerian (ALE) methods [38], molecular dynamics methods [39,40], and smoothed particle hydrodynamics (SPH) methods [3,38]. As noticed in [38], a traditional pure Lagrangian analysis, such as in [1], is not able to accurately model the impact process due to excessive computational cells or elements in the finite element analysis (FEA) distortion near the contact zone. However, generally, the SPH method is less accurate than the pure Lagrangian FEA method especially when the deformation is not severe. However, only using the SPH method, in contrast to ALE, the authors of [38] focused on the jetting phenomenon (regarding a jet moving ahead of the collision point during the oblique impact of the plates) and the composition of the jetted material able to be simulated. These conclusions were confirmed in [3] where SPH simulations in Ansys Autodyn allowed authors to accurately reproduce the formation of the wave boundary, vortex zones, as well as the formation of a jet moving ahead of the collision point in the oblique impact of metal plates. In [39,40], the wave formation process during explosive welding was carried out using a LAMMPS molecular dynamics simulator. Apparently, the molecular dynamics method is the most physically relevant one, but the size of the bodies simulated by the molecular dynamics method is 4–5 orders of magnitude less than the size of the plates that are used in explosive welding experiments [40]. Additional scaling procedures should be used. So, in this paper, we studied the principal possibility of the simulation of the impact of the plates using the shock-capturing, algorithmically uniform, and thus attractive Eulerian three-phase diffuse interface method.

The proposed approach was applied to the problem of a planar impact of steel and lead plates of different thicknesses. The problem follows previous experimental studies [8]. The dynamics of shock and rarefaction waves in the samples were analyzed. Quantitative estimation of the acceleration of the interfacial boundary due to the passage of the rarefaction wave from a free boundary of one of the plates was obtained. For the case of the 3-mm-thick steel plate and the 2-mm-thick lead plate, the simulated time of the rarefaction wave arrival at the contact boundary constituted 1.05 μs, and it was in good agreement with the experimental value of 1.1 μs.

The preliminary results of the two-dimensional simulation of the oblique impact of the metal plates using the proposed approach are presented. The results are encouraging from the point of view of subsequent simulation of the wave formation process during the oblique impact.