1 Introduction

Solid-fluid interaction that exhibits large deformations occurs in various applications, such as boundary layer flows^[1], penetration^[2], and explosion^[3]. Compared with experiments^[4-6], the application of numerical simulations^[7-8] is a common approach to obtaining more details regarding the flows. In numerical simulations of such phenomena, the accurate capture of the interface between different materials is a big challenge. The Lagrangian framework^[9] is widely used to simulate multi-material flows because the interface can be computed naturally by Lagrangian methods with minimal grids. However, it fails when treating problems that involve large deformations due to the mesh distortion^[10], and a fully Eulerian method becomes more suitable for such simulations. The most commonly used Eulerian framework for solid mechanics can be found in reports by Godunov and Romenskii^[11] and Plohr and Sharp^[12], where the shear stress is controlled by the deformation gradient tensor or its inverse. In multi-material Eulerian computation where the grid is fixed, the material interface may be inside a cell, leading to a mixture cell. An interface model should be considered to treat the mixture cell.

Two typical interface models have been reported in the literature: the sharp interface model (SIM) and the diffuse interface model (DIM). Compared with SIMs, such as the volume of fluid method^[13] and the level-set method^[14], the DIM^[15] has been rapidly developed in recent years because it is simpler and more conservative. In the DIM, flows in the whole domain are viewed as a mixture, and the same equations are solved. Most DIMs were based on the seven-equation models of Baer and Nunziato^[16] and Saurel and Abgrall^[17]. A series of reduced models for the seven-equation models were investigated by Saurel et al.^[18], Kapila et al.^[19], and Murrone^[20], assuming pressure equilibrium or velocity equilibrium. Direct extensions of Murrone's model and Saurel's model to solid-fluid interaction have been investigated by Favrie et al.^[21], Favrie and Gavrilyuk^[22], and Ndanou et al.^[23]. However, it was found that some intrinsic deficiencies of the DIM, which can be neglected for fluid-fluid interaction, appeared to have a significant effect on solid-fluid interaction.

A well-known deficiency of most existing DIMs for solid-fluid interaction is that they cannot handle problems with pure phases. If the volume fraction of a certain material takes the value of zero or unity, one of the phase-density equations is degenerated to a "0 = 0" equation because volume fractions are involved in conservative variables of phase-density equations. An artificial quantity of another material must be added into the pure phase to avoid the degeneracy of the model. However, the physical properties and wave structures of the solid and fluid are quite different. Even a small quantity of solids in the gas phase may significantly affect the physical structures. First, this small quantity of solids can generate non-zero shear stresses^[22], and shear waves are introduced, which is unacceptable in pure gas. Second, the sound speeds of gas phases are dramatically altered by this small quantity of solids, which was reported by Kluth and Després^[24] in a comparison of the Lagrangian result and Eulerian result in an impact problem. Therefore, it is necessary to propose a model that is able to handle zero/unity volume fractions.

The difficulty in treating the zero/unity volume fraction stems from the phase-density equations. To overcome this, several attempts have been made where mixture-density equations were used instead of phase-density equations in fluid-fluid interaction. Abgrall and Karni^[25] adopted this formulation and added a nonconservative Γ equation to describe the interface (, and γ is the ratio of specific heats). Shyue^[26] showed that the Γ equation was identical to the nonconservative volume-fraction equation under pressure-equilibrium hypotheses. Extensions for different equations of state were reported by Shyue^[27-28]. The above models aimed to avoid the use of the information of phase densities. However, the phase densities are essential for the Mie-Grüeisen equation of state (EOS), which is widely used to describe metals. A new formulation of phase-density equations, which is able to treat zero/unity volume fractions, is therefore needed.

In this paper, it is proven that a new formulation of phase-density equations unrelated to volume fractions can be derived if the compaction term in the volume-fraction equation is neglected. A transport DIM based on the new formulation of phase-density equations is proposed, which is able to treat zero/unity volume fractions. Consequently, zero/unity volume fractions do not result in degeneracy of the model. Mixture models for variables in the mixture region are then derived based on the pressure equilibrium assumption.

This paper is arranged as follows. In Section 2, the transport DIM is derived for the solid-fluid interaction. In Section 3, different numerical tests are presented to evaluate the proposed approach. The paper is concluded in Section 4.

2 Transport DIM for solid-fluid interaction 2.1 Original pressure-equilibrium DIM for solid-fluid interaction

The original pressure-equilibrium DIM for solid-fluid interaction^[22] can be written in the following form:

(1)

ρ₁ and ρ₂ are the phase densities. α₁ and α₂ represent the volume fractions of each phase, where α₁+α₂=1. All phases are assumed to have equal velocities u_k. The mixture stress is modeled as σ_ij=α₁σ_{ij, 1}+α₂σ_{ij, 2}. E=ε+u₁²/2+u₂²/2 is the total mixture energy, where the mixture internal energy ε is modeled as ε=α₁ε₁+α₂ε₂. p is the mixture pressure, and g_{ij, l}^e (l=1, 2) is the effective inverse deformation gradient tensor of each phase. L_ik^p is considered as the plastic effect, a_l (l=1, 2) is the sound speed of each phase, and Q represents the compaction effect.

It is evident that one of the phase-density equations in Eq. (1) can be degenerated to a "0 = 0" equation if α₁=0 or α₁=1. A new formulation of the phase-density equations is needed, which is able to treat zero/unity volume fractions, and this is derived in the next part of this section.

2.2 Transport DIM for solid-fluid interaction

Expanding the phase-density equations, we can get

(2)

Rearranging Eq. (2) gives

(3)

Considering the volume-fraction equation in Eq. (1), Eq. (3) can be reduced as follows:

Only if the compaction effects are neglected, and we set Q=0, we can get a formulation that is unrelated to the volume fraction,

(4)

It is seen that the form of Eq. (4) is the same as the mixture-density equation. The volume fraction does not appear in Eq. (4), and the phase density is only affected by the divergence of velocity. The variation of the volume fraction does not affect the phase density. Therefore, Eq. (4) is valid for arbitrary volume fractions, and zero/unity volume fractions do not lead to a degenerated equation. In this study, Eq. (4) is used to compute the phase density. Then, the new model can be written as follows:

(5)

In order to obtain a new formulation of phase-density equations that are valid for arbitrary volume fractions, the nonconservative term related to the compaction effect in the volume-fraction equation is omitted. Therefore, the new DIM is not able to describe cavitation and compaction effects. Hence, the volume fraction is mainly used to track the interface between two immiscible materials. From the definition of Murrone^[20], we refer to the proposed model as a transport DIM for solid-fluid interaction.

In this study, the Mie-Grüneisen EOS^[29] is used for solids, where

in which η=ρ/ρ₀. The stiffness gas EOS is used for fluids, where

The shear stresses are described by hyperelastic-viscoplastic constitutive equations^[22].

In numerical simulations, there is an artificial mixture region near the interface owing to numerical dissipation. The mixture density, Cauchy stress, pressure, internal energy, elastic energy, and sound speed in the artificial mixture region are to be modeled. The mixture density is defined by

The mixture Cauchy stress σ_ij is defined by

where σ_{ij, 1} and σ_{ij, 2} are computed as a function of g_{ij, l}. The total energy ρE is expressed as ρE=ρ u_ku_k/2+ρε^h+ρε^e. The mixing rules of the hydrodynamic internal energy and elastic energy are

The ratio of the current and initial densities can be obtained as follows:

In the above expressions, subscripts 1 and 2 represent the values of materials 1 and 2, respectively. Subscript 0 represents the value in the undeformed state, where the Cauchy stress is zero. If pressure equilibrium (p₁=p₂=p) is assumed, the mixture hydrodynamic internal energy can be written in the following form:

where

Then, the mixture pressure can be computed as

Another important variable is the mixture sound speed, which is essential in flux-vector splitting methods. The mixing rule for the sound speed of this system is modeled as ρc²=α₁ρ₁c₁²+α₂ρ₂c₂².

3 Numerical tests

In this section, we evaluate the proposed approach by performing four numerical tests. The impacts of the copper-gas mixture, tangential velocity discontinuity in pure gas, copper-aluminum impact, and two-dimensional projectile are considered. The governing equations are solved using the finite-difference method^[30].

3.1 Impact of copper-gas mixture

We consider a copper tube with a certain amount of porosity. The EOS of the copper is given by the Mie-Grüneisen model, which is characterized by the following parameters: ρ₀=8 930, a₀=3 940 m/s, Γ₀=2, and s=1.49. The shear modulus and yield stress are 45 GPa and 90 MPa, respectively. The EOS of the gas is given by the real gas model characterized by the following parameters: ρ₀=1, γ₀=1.4, and p_∞=0. The pressure p=10⁵ Pa is uniform throughout the entire domain. The velocity is taken as u=100 m/s on the left and u=-100 m/s on the right, separated at x=0.5. A mesh with 200 computational cells is used. The volume fraction of the air is α_gas=0.1. The elastic case and elastic-plastic case are considered simultaneously.

Figure 1 shows the computed results at time t=50 μs. The symbols correspond to the numerical solution of the original pressure-equilibrium model, while the lines correspond to the current model. E represents the elastic case, and P represents the elastic-plastic case. If plasticity is considered, the density after impact is higher. The phase densities computed by the proposed model and the original model agree well in both cases, verifying that the new model and the original model are identical for solid-gas mixtures.

3.2 Tangential velocity discontinuity in pure gas

In this case, a tangential velocity discontinuity is presented in pure gas. The initial discontinuity is located at x=0.5, and the tangential velocity is taken to v=100 m/s on the left and v=-100 m/s on the right. A small amount of copper (α=10^-6) is added in the gas phase for the original model, and pure gas is adopted for the current model. The EOS of the copper is given by the Mie-Grüneisen model, which is characterized by the following parameters: ρ₀=8 930, a₀=3 940 m/s, Γ₀=2, and s=1.49. The shear modulus is 45 GPa, and the solid is elastic. The tangential velocity and the shear stress at time t=1 ms with 200 cells are shown in Fig. 2. No shear stresses should be created in the gas phase. Therefore, a steady tangential velocity discontinuity should be retained. However, nonphysical shear waves are created by the artificial quantity of solids in the original model. However, the tangential velocity discontinuity is well maintained by the current model, and the shear stress of the mixture remains zero all the time.

3.3 Copper-aluminum impact

This test aims to show the advantage of the proposed model when there is the shock-interface interaction. The copper-aluminum impact is considered. The EOS of the copper is given by the Mie-Grüneisen model, which is characterized by the following parameters: ρ₀=8 930, a₀=3 940 m/s, Γ₀=2, and s=1.49. The shear modulus and the yield stress are 45 GPa and 90 MPa, respectively. The EOS of the aluminum is given by the Mie-Grüneisen model, which is characterized by the following parameters: ρ₀=2 785, a₀=5 328 m/s, Γ₀=2, and s=1.338. The shear modulus and the yield stress are 27.6 GPa and 300 MPa, respectively. The pressure p=10⁵ Pa is uniform throughout the entire domain. The initial interface is located at x=0.5, and the velocity is taken as u=100 m/s on the left and u=-100 m/s on the right. A small amount of one material (α=10^-6) is added to the other material. Three meshes with 200, 1 000, and 10 000 grids are considered to test the convergence of the models. Only the elastic-plastic case is compared in this case. The phase densities of each phase at time t=50 μs are shown in Fig. 3. The left figures are the results of the current model, while the right figures are the results of the original model. The phase densities computed by the current model converge rapidly with the increasing mesh number, while those of the original model hardly converge. At coarse grids, large oscillations appear for the phase densities computed by the original model near the interface. This test shows that the phase densities computed by the new model are more accurate when the shock-interface interaction occurs.

3.4 Two-dimensional projectile impact

This problem determines the impact of a two-dimensional aluminum projectile against a rigid wall. The initial projectile is a rectangular plate of length L = 5 m and thickness H = 1 m. The EOS for the solid phases is given by the Mie-Grüneisen model, and is characterized by the following parameters: ρ₀=2 785, a₀=5 328 m/s, Γ₀=2, and s=1.338. The shear modulus is 27.6 GPa, and the yield strength is 330 MPa. The EOS of the air is the real gas with ρ₀=1.0, γ₀=1.4, and p_∞=0.

Only the right half of the rod is computed. The initial computational domain is [x, y]∈[0, 2]×[0, 6]. The initial velocity is given by V₀=-150 m/s. The boundary conditions consist of a free boundary for the top and right boundaries. The left boundary is a symmetric boundary, while the bottom is a wall. The stopping time is t=8 ms, at which point the vertical velocity decreases to nearly zero. During the impact, the kinetic energy is entirely converted into the internal energy through plastic dissipation. Maire and Rebourcet^[29] and Kluth and Després^[24] simulated the same problem, and their final lengths were L=4.62 m and L= 4.34 m, respectively. In our simulation, the contour line of the solid volume fraction with α=0.5 is treated as the interface between the solid and gas. The length tends towards a limit that is approximately equal to L_final=4.53 m, which is close to the former results. The velocities in the x-direction of the current model and the original model at t=2.5 ms are compared in Fig. 4. A small quantity of solids (α=10^-6) is added in the gas phase for the original model, and pure gas is adopted for the current model. The dashed lines represent the interface, where α₁=0.5 is chosen to indicate the interface. The figure shows that the wave speed of the original model is larger than that of the current model owing to the addition of artificial solids. Oscillations can also be observed for the original model, which are suppressed by the current model.

4 Conclusions

An artificial quantity of solids must be added to the air phase in traditional DIMs to avoid degeneracy of the model for solid-fluid interaction. However, nonphysical waves will be generated because of the large difference between the properties of the solid and fluid. A transport DIM, which can resolve this difficulty, is proposed in this study. The main novelty of this work is to present a new formulation for phase densities, which is unrelated to volume fractions. Consequently, the proposed model is able to treat zero/unity volume fractions, which can prevent shear waves from traveling in pure gas. Numerical tests show that the proposed model performs better than the original model.

Acknowledgements  The authors are grateful to Prof. Jiequan LI from the Institute of Applied Physics and Computational Mathematics for his valuable suggestions.